PLASMA SIMULATIONS OF EMISSION LINE REGIONS IN
HIGH ENERGY ENVIRONMENTS
By
Chris T. Richardson

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Astrophysics and Astronomy - Doctor of Philosophy
2013

ABSTRACT
PLASMA SIMULATIONS OF EMISSION LINE REGIONS IN
HIGH ENERGY ENVIRONMENTS
By
Chris T. Richardson

This dissertation focuses on understanding two di↵erent, but in each case extreme,
astrophysical environments: the Crab Nebula and emission line galaxies. These relatively local objects are well constrained by observations and are test cases of phenomena seen at high-z where detailed observations are rare. The tool used to study these
objects is the plasma simulation code known as Cloudy. The introduction provides
a brief summary of relevant physical concepts in nebular astrophysics and presents
the basic features and assumptions of Cloudy. The first object investigated with
Cloudy, the Crab Nebula, is a nearby supernova remnant that previously has been
subject to photoionization modeling to reproduce the ionized emission seen in the
nebula’s filamentary structure. However, there are still several unanswered questions:
(1) What excites the H2 emitting gas? (2) How much mass is in the molecular component? (3) How did the H2 form? (4) What is nature of the dust grains? A large
suite of observations including long slit optical and NIR spectra over ionized, neutral
and molecular gas in addition to HST and NIR ground based images constrain a
particularly bright region of H2 emission, Knot 51, which exhibits a high excitation
temperature of ⇠3000 K. Simulations of K51 revealed that only a trace amount of
H2 is needed to reproduce the observed emission and that H2 forms through an uncommon nebular process known as associative detachment. The final chapters of this
dissertation focus on interpreting the narrow line region (NLR) in low-z emission line

galaxies selected by a novel technique known as mean field independent component
analysis (MFICA). A mixture of starlight and radiation from an AGN excites the gas
present in galaxies. MFICA separates galaxies over a wide range of ionization into
subsets of pure AGN and pure star forming galaxies allowing simulations to reveal
the properties responsible for their observed variation in ionization. Emission line ratios can constrain the spectral energy distribution, excitation mechanism, abundances
and physical conditions present in these galaxies, while the large data set allows many
weaker emission lines to be used as consistency checks. By integrating over a wide
range of densities and radii from the excitation source, the variation in ionization
for AGN can be represented as change in the central concentration of clouds in the
NLR. Preliminary analysis from modeling star forming galaxies indicates that the
same interpretation might apply to galaxies without an AGN in which gas is excited
by starlight.

Copyright by
CHRIS T. RICHARDSON
2013

To my Bob.

v

ACKNOWLEDGMENTS

I would especially like to thank my adviser Jack Baldwin for his patience, his generosity, the clarity of his advice and his support of my career goals.

To Mom and Dad, Nan and Bob, you have all played a larger part in this dissertation
than you will ever know.

To Carolyn, thank you for your never ending support and reassurance. Traveling
through this process together has been unforgettable. By the way, we have cats.

To Brian O’Shea, thank you for your guidance, mentorship and candid advice throughout my graduate career.

To Gary Ferland, Ed Loh, James Allen, Paul Hewett, Charles Kuehn, Andy Fabian
and Philippe Salom´, thank you for your hard work and dedication to these projects.
e

To my friends, thank you for your support and years of laughter.

Thank you to the Michigan State University High Performance Computing Center
and the Institute for Cyber Enabled Research for providing the computational facilities and support necessary to run many of the models in this work.

A
Thanks to Christopher Waters for the L TEX class used to format this dissertation.

vi

TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

ix

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

x

1 A Brief Introduction to Nebular Astrophysics . . . . . . . . . . . . . . . . . . . .
1.1
Photoionization Equilibrium . . . . . . . . . . . . . . . . . . .
1.2
Thermal Equilibrium . . . . . . . . . . . . . . . . . . . . . . . .
1.3
Physical Conditions from Observations . . . . . . . . . . . . . .
1.4
Plasma Simulations with Cloudy . . . . . . . . . . . . . . . . .

1
1
4
6
8

2 The Nature of the H2 Emitting Gas in the Crab Nebula . . . . . . . . . . 10
2.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2
Observational Data . . . . . . . . . . . . . . . . . . . . . . . . 17
2.2.1
Narrow-Band Emission Line Images . . . . . . . . . . . . . 17
2.2.2
Spectra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.2.3
Dust Absorption Features . . . . . . . . . . . . . . . . . . . 24
2.2.4
Mid-IR Images . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.2.5
Intensity of the Radiation Field Striking the
Filaments . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.2.6
Adopted Error Bars . . . . . . . . . . . . . . . . . . . . . . 29
2.3
Cloudy Simulations of Knot 51 . . . . . . . . . . . . . . . . . . 30
2.3.1
General Methodology . . . . . . . . . . . . . . . . . . . . . 31
2.3.2
H2 Emission in Non-Equilibrium Environments . . . . . . . 32
2.3.3
A Basic X-ray Photoionization Model in Pressure Balance . 35
2.3.4
Models with Dense Cores . . . . . . . . . . . . . . . . . . . 53
2.3.5
Models with Additional H2 Mechanisms . . . . . . . . . . . 55
2.3.6
Models with Increased Dust Abundance . . . . . . . . . . . 65
2.4
Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
2.4.1
Understanding the H2 Formation Environment . . . . . . . 68
2.4.2
Timescale Considerations . . . . . . . . . . . . . . . . . . . 69
2.4.3
E↵ects of Geometrical Uncertainty . . . . . . . . . . . . . . 71
2.4.4
Are Models with Fully Molecular Cores Ruled Out? . . . . . 72
2.4.5
The Role of Dust . . . . . . . . . . . . . . . . . . . . . . . . 75
2.4.6
Mass Estimate . . . . . . . . . . . . . . . . . . . . . . . . . 78
2.4.7
Crucial Observational Questions . . . . . . . . . . . . . . . 80
2.5
Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
vii

3 Interpreting the Ionization Sequence in Active Galactic Nuclei Emission Line Spectra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
3.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
3.2
A Comparison Sample at Representative
Points Along the Locus . . . . . . . . . . . . . . . . . . . . . . 90
3.3
Modeling the Narrow Line Region . . . . . . . . . . . . . . . . 99
3.3.1
Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
3.3.2
Spectral Energy Distribution Optimization . . . . . . . . . . 102
3.3.3
X-ray E↵ects . . . . . . . . . . . . . . . . . . . . . . . . . . 109
3.3.4
Metallicity E↵ects . . . . . . . . . . . . . . . . . . . . . . . 112
3.3.5
Final Dust-Free Model . . . . . . . . . . . . . . . . . . . . . 113
3.3.6
Dusty Models . . . . . . . . . . . . . . . . . . . . . . . . . . 129
3.4
Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
3.4.1
Physical Meaning of the AGN Locus . . . . . . . . . . . . . 130
3.4.2
[O III] and [N II] Temperatures . . . . . . . . . . . . . . . . 137
3.4.3
The Possibility of a Dusty Narrow Line Region . . . . . . . 140
3.5
Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143
4 Theoretical Modeling of the Ionization Sequence in Star-Forming
Galaxy Emission-Line Spectra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145
4.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146
4.2
Observations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150
4.3
Preliminary Modeling Results . . . . . . . . . . . . . . . . . . . 161
4.3.1
Spectral Energy Distribution . . . . . . . . . . . . . . . . . 163
4.3.2
Excitation Mechanism Diagrams . . . . . . . . . . . . . . . 164
4.4
Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173
REFERENCES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 177

viii

LIST OF TABLES

Table 2.1

Archival HST Images . . . . . . . . . . . . . . . . . . . . . . . .

19

Table 2.2

K51 Observed and Dereddened Emission Line Fluxes . . . . . .

22

Table 2.3

K51 Chemical Abundances . . . . . . . . . . . . . . . . . . . . .

36

Table 2.4

Model Parameters and Predicted/Observed Surface Brightness
Ratios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

50

Table 3.1

Properties of aij subsets . . . . . . . . . . . . . . . . . . . . . .

92

Table 3.2a

The measured emission line strengths for the AGN Locus. Measurements are relative to H . . . . . . . . . . . . . . . . . . . . .

95

The dereddened emission line strengths for the AGN locus. Measurements are relative to H . . . . . . . . . . . . . . . . . . . . .

97

Table 3.2b

Table 3.3a

AGN Emission line ratio predictions - Dust Free Model . . . . . 117

Table 3.3b

AGN Emission line ratio predictions - Dusty Model . . . . . . . 120

Table 4.1

Properties of sij subsets . . . . . . . . . . . . . . . . . . . . . . 151

Table 4.2a

The measured emission line strengths for the SF locus. Measurements are relative to H . . . . . . . . . . . . . . . . . . . . . . . 154

Table 4.2b

The dereddened emission line strengths for the SF locus. Measurements are relative to H . . . . . . . . . . . . . . . . . . . . . 157

ix

Figure 4.4

Diagnostic diagrams from VO87 for a continuous SFH as a function of age. Each quadrant represents a di↵erent starburst age.
In each subpanel the squares connected by the dashed black line
are the dereddened measurements from the co-added spectra that
constitute the s41-s31-s21-s11-s01 sequence. We do not show the
other SF sequences due to their similarity to the central sequence.
The largest shape in the sequence indicates observations at the
top of SF locus (s41). The colored lines with circles represent
LOC integrations. The density-weighting indices in the LOC integrations, , have more negative values along the direction of the
arrow and are indicated by di↵erent colors. The radial weighting
indices, , become more negative as a function of increasing distance from the largest colored circle for each particular density
weighting. The hollow circle indicates our best fitting set of free
parameters that match the a41 subset (large black square) in the
0.1 Myr starburst. The cross in the lower right panel of each
section indicates the range of acceptable error. . . . . . . . . . . 164

Figure 4.5

Diagnostic diagrams from VO87 for a instantaneous SFH as a
function of age, displayed in the same manner as Figure 4.4. . . 168

xviii

1 A Brief Introduction to Nebular
Astrophysics
This chapter is meant to give a brief overview of the main physical concepts in nebular
astrophysics needed to understand subsequent chapters: photoionization equilibrium,
thermal equilibrium and emission line diagnostics. These concepts were developed
many decades before this work and therefore this summary draws upon the more
complete introductions presented in several textbooks including Spitzer (1978), Tielens (2005), Osterbrock & Ferland (2006). Certain physical processes, such as hydrodynamics, are omitted or briefly presented for simplicity and so the reader is urged to
seek these references and others for a more complete understanding. The astrophysical code known as Cloudy has long been the industry standard for self-consistently
incorporating all of these physical processes, and many more, to model gaseous clouds
over a vast range of ionization and physical conditions and predict the emitted spectrum. A brief summary of the assumptions, reliability and user input for Cloudy are
introduced.

1.1 Photoionization Equilibrium
Photoionization equilibrium or the balance between photoionization and recombination of ions with electrons determines the ionization structure of nebulae. The

1

generalized equation for photoionization equilibrium is
n(X i )

Z 1
4⇡J⌫
a⌫ (X i )d⌫ = n(X i+1 )ne ↵(X i , T )
h⌫
⌫i

[cm 3 s 1 ]

(1.1)

where X i is an element X with ionization i, ⌫i is the ionization threshold frequency,
J⌫ is the mean intensity, a⌫ is the ionization cross section and ↵ is the recombination
coe cient. The left side of the equation give the volumetric photoionization rate, or
the number density of X i atoms multiplied by the ionizing flux and photoionization
cross section. This is balanced by the right side of the equation, which represents the
volumetric recombination rate, or number density of X i+1 multiplied by the electron
density and the recombination coe cient. Note, the incident radiation field sets the
photoionization rate while the temperature sets the recombination rate. To a rough
approximation ↵ ⇠ T 1 , so as the temperature decreases the rate of recombination
increases.
Photons with

< 512 ˚ or h⌫ > 13.6 eV ionize hydrogen and originate from a
A

variety of objects including the surface of hot stars (e.g. O stars), pulsars generating
synchrotron emission, and radiation from active galaxies. Hydrogen is by far the most
abundant element in the universe and its simplicity allows a relatively basic physical
treatment that can provide insight about the typical scales in nebular astrophysics.
If we consider a pure hydrogen nebula, Equation 1.1 gives,
n(H0 )

Z 1
4⇡J⌫
a⌫ (H0 )d⌫ = np ne ↵(H0 , T )
h⌫
⌫i

(1.2)

where the proton (H0 ) number density np = ne . From the Einstein coe cients we
can see that the lifetime of an electron in an excited level of hydrogen is 1/Aij = tij ⇠
10 8

10 4 s. If we evaluate the left side of Equation 1.2 for a typical nebula, the

photoionization timescale, tphoto is ⇠ 108 s for the 1s2 S ground state, therefore, to
first approximation, all photoionizations of hydrogen are from the 1s2 S ground state.
2

The recombination timescale for the 1s2 S level is given by

trec =

1
ne ↵(H0 , T )

(1.3)

where ne ⇠ 10 cm 3 and ↵ ⇠ 2.0 ⇥ 10 13 cm 3 s 1 for a typical completely ionized
nebula with T ⇠ 104 K. This gives trec = 5 ⇥ 1011 s or approximately 104 years,
so a freshly ionized hydrogen atom will remain ionized for quite some time before
recombining with a free electron, only to quickly be ionized again. It follows under
these conditions that the neutral fraction of hydrogen will be very small and therefore
these regions are called H II regions, where H II stands for H+ . Electrons in H II
regions thermalize to a Boltzmann distributions very quickly since the cross section
for e + e

is orders of magnitude larger than that for photoionization. Typical

elastic e + e

timescales are

te,e ⇡

104

✓

◆
Ee 3/2 1
1eV
ne

[s]

(1.4)

where Ee is the energy of the electron. Taking the values above, ne ⇠ 10 cm 3 and
T ⇠ 104 K, we get te,e ⇠ 103 s which is much shorter than both the photoionization
and recombination time.
The mean free path of an ionizing photon is roughly the size of a typical H II region.
We can parameterize the neutral fraction by defining ⇠ such that n(H0 ) = ⇠n(H).
When the bulk of the ionized hydrogen begins to transition to neutral hydrogen, we
can take ⇠ = 0.5. The mean free path can then be evaluated as

l⇡

1
n(H0 )a⌫

(1.5)

where n(H0 ) = 5 cm 3 , since n(H) = 10 cm 3 . For the 1s2 S ground state ionization
cross section for hydrogen is ⇠ 6 ⇥ 10 18 cm 2 . Evaluating the mean free path
3

gives l ⇠ 3 ⇥ 1016 cm which is a typical size of an H II region. Thus after this
depth, ionizing photons will on average be absorbed and the hydrogen transitions
from ionized to neutral.

1.2 Thermal Equilibrium
The overwhelmingly dominant heating mechanism in nebulae and emission line galaxies occurs through photoionization of H and He. Photons with energy h⌫ ionize
atoms with an ionization potential h⌫0 , which produces a photoelectron with energy
h(v v0 ). As shown in the previous section, this photoelectron then elastically collides
with thermal electrons. Every photoionization is followed by a recombination, which
cools the gas. Therefore the energy di↵erence between an ionizing photoelectron and
the photons released during recombination with a thermalized electron, including the
subsequent cascade of the captured electron to the ground state, represents the net
energy injected into the gas by photoionization.
Assuming photoionization equilibrium, representing this numerically we have
3
Qnet = kTi ne n(X i )↵(X i , T )
2
where

ne n(X i )kT (X i , T )

[erg cm 3 s 1 ]

(1.6)

is recombination cooling coe cient, which depends on the recombination

cross section. The heating and recombination cooling rates are proportional to the
number densities of each ion; therefore heavy elements can be omitted from these
processes due to their low abundances relative to H and He, and Qnet ⇡ Qnet (H) +
Qnet (He).
However in spite of their low abundance, cooling from ion

electron collisions of

heavy elements is important to the overall thermal structure due to their low ionization potential. Other collisions, such as ion

ion and ion

proton, are ine cient due

to coulomb repulsion. Collisions between electrons and ions are e cient for metals

4

and make the largest contribution to overall cooling in nebulae.
It is simplest to consider a two level atom. Even with this simplification, the
details of collisionally excited line radiation are complex and therefore we will jump
straight to the final form of the cooling rate to consider the implications of this form
of cooling. The cooling rate from collisions is then,
!
C = ne n1 h⌫21 2 e h⌫21 /kT
!1

✓

ne q21
1+
A21

◆ 1

(1.7)

where h⌫21 is energy required raise the electron from level 1 to level 2, q21 is the
rate coe cient [cm3 s 1 ], !1 and !2 are the statistical weights for levels 1 and 2
respectively, and A21 is Einstein coe cient for spontaneous emission. As ne ! 0,
all collisions are followed by radiative de-excitations and collisional de-excitation is
unimportant. However as ne ! 1, collisional de-excitation is non-negligible and
suppresses the cooling rate. In this limit, the cooling rate is simply the thermal
equilibrium value.
The physical meaning from all of this is that electrons collide with metal ions,
which excites the bound electrons to higher energy states. If the electron density is
low enough then an excited electron de-excites to a lower energy state and the emitted
photon escapes the nebula, thereby cooling the gas. This is called the low-density
limit (LDL). But, if the electron density is high enough, then a second free electron
could collide with an excited ion and cause collisional de-excitation. In this case, the
second free electron carries the energy back into the gas and there is no net cooling
in the nebula. This is called the high-density limit (HDL).
For a two level atom, the transition from LDL to HDL is defined by,

ncrit =

A21
q21

(1.8)

where ncrit is the critical density. So, high densities are defined as any density above
5

the critical density, where collisionally excited line emission is suppressed. After this
limit, any further increases in ne will have a reduced impact on emission compared
to increases in ne below the critical density.
Finally, the highly simplified final form of thermal equilibrium in most nebulae
can be thought of as Qnet (H, He) ⇡ C(metals) where Qnet (H, He) is the net heating
rate from H and He after taking into account recombination cooling, and C(metals)
is cooling rate from metals due to electron

ion collisions.

1.3 Physical Conditions from Observations
In nebular astrophysics, understanding phenomena is often characteristic of the inverse problem in astronomy: we know the emitted spectrum and we want to know
the conditions that are able to produce it. Emission line ratios can probe several
properties of nebulae, but covering all of these properties in detail is beyond the
scope of this introduction. However, constraining electron temperatures and electron
densities using emission line ratios are two techniques used in the following chapters,
so we discuss them here.
As mentioned in Section 1.2, the ground states of metal ions are populated due to
collisional excitation. For a given ion, upper levels with di↵ering excitation energies
but with similar critical densities give rise to a temperature indicator. Specifically,
the upper levels of O III (i.e. O++ ) and N II are observed in the optical and are the
most common indicators of nebular temperature. For example, the [O III] 4363 line
populates an upper level while [O III] 4959 and [O III] 5007 (the ratio of which is
fixed by atomic physics) populate an intermediate level. Therefore, the ratio [O III]
4363 / [O III]

4959, 5007 constrains the electron temperature. Plugging in the

6

necessary atomic data for these lines we find,
4

j 4959 + j 5007
7.90e3.29⇥10 /T
=
j 4363
1 + 4.5 ⇥ 10 4 ne T 1/2

(1.9)

where j represents the emissivity [erg s 1 cm 3 sr 1 ] of each respective emission line. As expected, this ratio is highly sensitive to the temperature and weakly
dependent on the electron density.
Conversely, an ion with similar upper level excitation energies but di↵ering rate
coe cients or spontaneous emission coe cients gives rise to an electron density indicator. One example of this occurs with the emission lines [S II]
II]

6731. In LDL (⇠ 102 cm 3 ), the ratio [S II]

6716 and [S

6716/ 6731 is simply a con-

stant determined by the ratio of their statistical weights, !6716 /!6731 = 1.5. The
statistical weight signifies how many electrons a level can hold. The [S II] ratio is
greater than unity, since [S II] 6716 has a larger statistical weight and, in the case
of LDL, collisional de-excitations are unimportant. In HDL (⇠ 105 cm 3 ), the [S
II] ratio is also constant, but now collisional de-excitation does play role and thus
the radiative transition probabilities do as well. Now, the constant is determined by
(!6716 A6716 )/(!6731 A6731 ) = 0.44.
Therefore, between LDL and HDL, the [S II] ratio smoothly decreases from a
constant value to smaller constant value. In between these limits, the electron density
can be inferred from the observed [S II] ratio. As the electron density increases from
LDL, the [S II] 6731 emission line becomes stronger relative to [S II] 6716 since
the lifetime of an excited electron is smaller for the transition resulting in the [S II]
6731 emission line. Eventually however, the density becomes high enough that both
[S II] lines have Boltzmann populations and thus their ratio approaches a constant.

7

1.4 Plasma Simulations with Cloudy
Cloudy (Ferland et al. 2013)1 is an open source, plasma simulation code that selfconsistently calculates the ionization, thermal and chemical structure of an astrophysical cloud and, most importantly, predicts the emitted spectrum. Cloudy can reliably
simulate an astrophysical plasma with temperatures in the range of 100
densities in the range 10 9

108 K and

1015 cm 3 . For temperatures above 108 K, electrons

and ions decouple and this functionality is currently not included in version c10.00.
Cloudy assumes that the gas considered is in steady state. This is usually an
accurate assumption for temperatures above 104 K. The recombination timescale is
generally much longer than the ionization or electron thermalization timescale and
yet is still much shorter than the age of the nebula in question. However, in molecular
regions, reaction timescales such as those for associative detachment (Section 2.4.2)
can be longer than the age of the nebula indicating that a time dependent calculation
is needed.
Cloudy divides the cloud into many zones, in each of which the physical conditions are approximately constant. The equations governing conservation of energy,
conservation of charge and statistical equilibrium are solved self-consistently with a
minimal set of free parameters. Thus, a large set of observables is generated with
the parameters usually constrained by observations. Only the chemical composition, shape and intensity of the incident radiation field, geometry, equation of state
and boundary conditions must be specified to initiate a simulation. The boundary
condition criterion often entails a temperature, column density or cloud thickness
restriction.
Simulations with Cloudy can be performed in several seconds or can take several
days, depending on the complexity of the atoms and molecules included. The advent
of high performance computing allows the computation of massive grids of plasma
1 http://www.nublado.org/

8

simulations at the price of a single simulation. Cloudy comes equipped with an
“embarrassingly parallel” mode that allows a single model to run every computer
core supplied. When simulations finish running, new simulations still in the queue
run on newly available cores resulting in a very simple yet e cient algorithm. The
simulations presented in the following chapters relied heavily on this functionality.

9

2 The Nature of the H2 Emitting
Gas in the Crab Nebula
Understanding how molecules and dust might have formed within a rapidly expanding
young supernova remnant is important because of the obvious application to vigorous supernova activity at very high redshift. In previous papers we have mapped the
Crab Nebula in a roto-vibrational H2 emission line, and then measured the molecular
excitation temperature for a few of the brighter H2 -emitting knots that we have found
to be scattered throughout the Crab Nebula’s filaments. We found that H2 emission
is often quite strong, correlates with optical low-ionization emission lines, and has
a surprisingly high excitation temperature. Here we study Knot 51, a representative, bright example. It is a spatially isolated structure for which we have available
long slit optical and NIR spectra covering emission lines from ionized, neutral, and
molecular gas, as well as HST visible and SOAR Telescope NIR narrow-band images.
We present a series of Cloudy simulations to probe the excitation mechanisms, formation processes and dust content in environments that can produce the observed
H2 emission. There is still considerable ambiguity about the geometry of Knot 51,
so we do not try for an exact match between model and observations. Rather, we
aim to explain how the bright H2 emission lines can be formed from within a cloud
the size of Knot 51 that also produces the observed optical emission from ionized
and neutral gas. Our models that are powered only by the Crabs synchrotron radiation are ruled out because they are not able to reproduce the observed strong H2
emission coming from thermally populated levels. The simulations that come closest

10

to fitting the observations (although they still have conspicuous discrepancies) have
the core of Knot 51 almost entirely atomic with the H2 emission coming from just
a trace molecular component, and in which there is extra heating. In this unusual
environment, H2 forms primarily through H

by associative detachment rather than

by grain catalysis. In this picture, the 55 H2 -emitting cores that we have previously
catalogued in the Crab have a total mass of about 0.1 M , which is about 5% of the
total mass of the system of filaments. We also explore the e↵ect of varying the dust
abundance. We discuss possible future observations that could further elucidate the
nature of these H2 knots.1

2.1 Introduction
The Crab Nebula (hereafter “the Crab) presents a unique opportunity to study not
only the properties and evolution of a very young supernova remnant which is not yet
interacting with the ISM, but also in a more general way the details of the physical
processes that occur in filamentary condensations exposed to a harsh environment
of high-energy photons and particles. The relative proximity of the Crab, at a distance d ⇠2.0 kpc (Trimble 1968), and the well-known date of the explosion, 1054
A.D. (Duyvendak 1942; Mayall & Oort 1942), allow detailed observations and constraints. The Crab filaments have been the subject of many previous investigations
(see the review articles by Davidson & Fesen 1985 and by Hester 2008). The Crabs
synchrotron radiation field has been taken into account as the energy source for the
emitting regions, but not the e↵ects of the high-energy particles that must also be
present. Most of this previous work has concentrated on the ionized gas, generally
with the aim of using optical emission lines to measure the chemical abundances in
the aftermath of the supernova explosion. H2 molecules are also present in the Crab,
as was shown by the pioneering infrared observations and spectroscopy by Graham,
1 This chapter is almost entirely taken word for word from Richardson et al.

11

(2013a)

Wright & Longmore (1990; hereafter G90). Those measurements were made with
quite low (20” FWHM) angular resolution in just three locations.
We have been carrying out a much more detailed study of the molecular content
of the Crab, mapping most of the nebula at sub-arcsec angular resolution in the NIR
H2 2.12 µm emission line, and tying those results together with both existing and new
spectra and images over a wide range of wavelengths. In our first paper in this series
(Loh, Baldwin & Ferland 2010; hereafter Paper I), we investigated the brightest H2
feature, Knot 1, for which the molecular content was estimated to be Mmol 5 ⇥ 10 5
M . We arrived at this estimate under the assumption that the H2 emission is
collisionally excited in a predominately molecular core. Then, in Loh et al. (2011)
(hereafter Paper II) we presented our full near-infrared imaging survey made with the
Spartan Infrared Camera on the 4 m SOAR telescope. We detected 55 compact (1-2”
diameter) knots, spread across much of the face of the Crab, that radiate strongly in
the H2 2.12 µm line. We showed that the positions of these knots correlate strongly
with compact, low-ionization regions that are conspicuous in Hubble Space Telescope
(HST) [S II]

6720 images; [S II] and also [O I] are H2 tracers. The low-ionization

regions are found on or next to the bright filaments seen in higher ionization lines
such as [O III] 5007. We used optical spectra to measure radial velocities of most of
the [S II] features. Then in Paper III in this series (Loh et al. 2012), we used K-band
spectra to measure the H2 excitation temperature of 7 of the brightest knots. From
the ratio of the S(1) 2-1 and S(0) 1-0 lines we found Texc 2000

3000 K, which we

argue is also the kinetic gas temperature in the regions where H2 lines form under the
conditions present in the Crab (Ferland 2011). These temperatures are just below
the dissociation temperature of H2 .
There is also dust in the Crab, as can be seen both from the weak thermal dust
emission peak at

80 µm in the Crabs continuum spectrum (Marsden et al. 1984),

and from the presence of many small dust globules seen in silhouette against the

12

Crabs synchrotron continuum (Fesen & Blair 1990, Hester et al. 1990). The mass
of dust in the filaments is uncertain, ranging from 10 3 M
2012; Hester 2008) to 0.2 M

(Temim et al. 2006,

(Gomez et al. 2012). However, the gas to dust ratio

could vary considerably between the filament cores and the ionized gas, and di↵erent
regions could have both di↵erent abundances and di↵erent dust-to-gas ratios. Sankrit
et al. (1998) estimate a dust-to-gas mass ratio in the cores that is a factor of 10 above
the ISM value.
The ionized, neutral and molecular gas and the dust are distributed throughout
a complex system of filaments that extend down into the synchrotron plasma and
hence must be bombarded from all sides by intense synchrotron radiation and also
by the high-energy particles that emit the synchrotron radiation. A more complete
understanding of the nature of these filaments and condensations would shed light on a
number of open questions including: (1) What is the excitation process that produces
the strong molecular hydrogen emission? Photoionization, shocks, dissipative MHD
waves, and high-energy particles are all possible energy sources. (2) What are the
properties of the grains that formed in the cores of the filaments? The Crab serves
as a test bed for understanding grain formation in high-redshift galaxies where shortlived massive stars and supernovae must produce the observed dust. (3) What are the
molecular processes that occur in this hostile environment? In the ISM H2 commonly
forms by grain catalysis, a process dependent on temperature and composition of the
grains, but is that true in an environment like the Crab? Many H2 knots do not
show dust extinction. Could H2 form without dust? (4) What is the true mass of
the filaments? Molecular gas is less emissive than ionized gas and thus the molecular
regions of the filaments could contain a substantial amount of mass. (5) The Crab
filaments appear to be nearby, well-resolved prototypes of the physically much larger
filaments seen in the IGM of cool-core galaxy clusters (Fabian et al. 2008). Both show
a low-ionization spectrum, have similar geometries, and are surrounded by ionizing

13

particles. The Crab can serve as a test bed for understanding the physics of the
cluster filaments.
In this paper, we explore the physical conditions in the gas that produces the
H2 emission. We aim to build on the pioneering first look at this problem by G90.
Our approach is to fit photoionization models, supplemented by other sources of
heating and ionization, to our extensive observations of a single knot, using the plasma
simulation code Cloudy (Ferland et al. 1998) to explore the H2 excitation mechanisms,
formation processes, level populations and the role of dust in the molecular core.
While previous studies have addressed the structure of ionized skins on the filaments (Pequignot & Dennefeld 1983, Sankrit et al. 1998), none have included the
molecular cores. In addition, our observations provide us with more detailed constraints than have been used in previous modeling e↵orts. Many papers have explored a wide range of possible abundances and argued that point-to-point abundance
variations within the Crab filament system are quite plausible (e.g. Davidson 1978;
Davidson 1979; Pequignot & Dennefeld 1983; Henry and MacAlpine 1982; MacAlpine
and Satterfield 2008; Satterfield et al. 2012). Our study uses a wider variety of optical emission lines and focuses on a more spatially compact feature (a knot) than the
filaments considered in other papers.
The brightest H2 knot in our catalog, Knot 1, is not the best candidate for this
study. It su↵ers from lying in a confused region of the Crab. It is immediately
adjacent to, and probably physically associated with, the much more extensive bright
highly-ionized filament FK-10 (Fesen & Kirshner 1982), but FK-10 does not have
strong H2 emission while Knot 1 is not at all conspicuous in high ionization lines such
as [O III] 5007 (Paper I). It is unclear where the boundaries lie between these two
rather di↵erent types of emission regions. In addition, the radial velocities determined
from our long-slit optical spectra show that there is considerable other material at
di↵erent points along the same line of sight, and that the strong [S II] emission from

14

N"
E"
1’"="0.6"pc"

Pulsar"
Knot"51"

Figure 2.1 Positions of Knot 51 and the pulsar, shown on grey-scale rendering of the
well-known HST composite color image made with F502N, F631N and F673N filters
(courtesy of NASA, ESA and J. Hester, Arizona State University; StScI News Release
2005-37).
Knot 1 comes from material with a very complex velocity field.
We instead study a much simpler region, Knot 51 (K51) from our catalog in
Paper II. Figure 2.1 shows that K51 lies near the tip of a long tendril of filamentary
gas. Both K51 and the tendril have small radial velocities, indicating that they lie
approximately halfway through the nebula, where the direction of expansion is in the
plane of the sky. The tendril is actually a pair of elongated features extending back
down into the synchrotron plasma from one of the large-scale filaments that wraps
around the outer surface of the synchrotron plasma. The two join into a single tendril
before reaching K51. As we discuss below, the velocity structure in this region shows
that there is minimal contamination from foreground or background filaments. The

15

situation is still not ideal, because as is described below there is a second H2 knot
immediately adjacent to K51 that is part of the same filament, but this is the cleanest
case that we can find.
The following sections first present the available observations of K51, and then
describe a series of Cloudy models that explore a large parameter space with the goal
of reproducing the observed H2 2.12 µm emission along with the observed lines of
neutral and ionized atomic species. We do not aim to produce an accurate model
of the knot, in part because the existing observations do not uniquely specify the
geometry of the molecular region (higher spatial resolution data are needed) or reveal
the total molecular inventory (sub-mm observations are needed). Rather, this paper
is a plausibility study whose goal is to identify what general type of model can account
for the existing suite of observations. The observations of ionized and neutral species
in K51 serve as a guide to the overall structure of a typical condensation, but in spite
of being the simplest case, K51 is still a rather complicated structure for which we
do not attempt to fit all details of the observations.
In our previous papers we had assumed that the H2 emission comes from a predominantly molecular core of some type. The simulations presented here show that
in the unusual environment of the Crab it is instead possible that the observed H2
emission can be produced by a trace H2 component in what is basically a very extended H0 zone. Although dust is clearly present in the particular case of K51, grains
are not actually required to form H2 in this environment, because of the high electron
density where atomic hydrogen is present, which results in an enhanced e ciency of
the H

route. Finally we outline the types of observations that would be needed to

further constrain physical models of the atomic or molecular regions.

16

2.2 Observational Data
We assembled a wide variety of observations of K51 in order to constrain the nature of
its ionized, neutral and molecular gas, and its dust content. These include space- and
ground based images at visible, near-infrared (NIR) and mid-infrared wavelengths,
together with ground-based visible and NIR spectra.

2.2.1 Narrow-Band Emission Line Images
Figure 2.2 shows a series of visible and NIR images, all on the same scale, in boxes
covering the same 10” (0.097 pc) square patch on the Crab. These are magnified
images of the region shown in the insert on Figure 2.1. They include K51 as well as
Knot 50, which lies 3.4” away. The pair of small white crosses on each image mark
the position of the peak H2 emission from the two knots.
The H2 image is in the 2.12 µm line and is from our data described in Papers I
and II, taken with the Spartan Infrared Camera on the SOAR Telescope. It has 0.7”
FWHM angular resolution. The [S II]

6716+6731, [O I] 6300, [O III] 5007, H

and continuum images are from HST archival data in which the 0.1” pixels undersample the HST point spread function. Since the Crabs expansion causes significant
proper motion of the knots, for each filter we have only summed together HST data
taken at a single epoch (Table 2.1). Following the procedure described in Paper II,
the world coordinate system of each co-added HST image was then rescaled around
the center of the Crab’s expansion to adjust it to the same epoch (2009 Dec) as our
ground-based H2 images. This allows for an accurate comparison of the Knots morphology as seen at di↵erent wavelengths, smeared only by actual time evolution of its
internal structure.
The HST continuum image will be discussed below. The emission-line images
in Fig 2.2 show that in H2 emission, K51 appears slightly elliptical with a roughly
Gaussian surface brightness distribution. The [O III] emission has a conspicuous hole
17

Contin. 5470

H2
+
+

+

K51
+

K50

H

[O III] 5007
+

+
+

+

[O I] 6300

[S II] 6720

+

+
+

+

Figure 2.2 Images of Knots 50 and 51 through various filters. All boxes are 10”⇥10”
and cover the same area on the sky with N up and E to the left. The pair of small
white crosses on each image mark the position of the peak H2 emission from the two
knots. The H2 image is from our ground-based imaging, the others are HST archival
images. The [S II] 6720 and [O I] 6300 emission traces the H2 2.12 µm emission.
at the position of the H2 centroid and a pair of bright patches on opposite sides of that
hole, and then is strong along the filament that runs up to the north-west to join onto
the main system of high-ionization filaments seen in Figure 2.1. The H↵, [S II] and [O
I] images all not only show the same structure to either side of the molecular core that
is seen in the [O III] image, but also have emission from the position of the H2 central
core. This follows the general pattern reported by Sankrit et al. (1998) and others,
in which low-ionization cores are surrounded by more extended high-ionization halos.

18

Table 2.1: Archival HST Images
Line

Filter No. Images Total Expo. (s)

H

F487N

2

4000

1994 Feb 2

[O III] 5007 F502N

2

5200

2000 Jan 26

F631N

2

3900

2000 Jan 26

[S II] 6720 F673N

2

2600

2000 Jan 26

Continuum

2

11600

1994 Mar - 1995 Dec

[O I] 6300

F547N

Dates

In addition, all of the visible-wavelength lines are also emitted by a circular ring to
the SE (lower left in Figure 2.2) of K51 that includes H2 Knot 50 on the ring’s far
edge, and additional fainter H2 emission from along the ring.

2.2.2 Spectra
Our grid of long-slit optical spectra taken with the Kitt Peak National Observatory
(KPNO) 2.1 m telescope covers Knots 50 and 51 (Paper II). The spectrum covers the
wavelength range

3700-7400 ˚at 6.3 ˚resolution, with a 3.8” slit width and 2.7”
A
A

FWHM resolution along the slit. Figure 2.3 (left panel) shows the slit position of
the spectrum that crosses Knots 50 and 51, and next to it (in the center two panels)
on the same spatial scale are two segments of our KPNO 2D spectrum, showing the
various velocity systems in that region of the Crab. Our NIR spectrum across the
two knots (from Paper III) is shown in the right-hand box.
The emission lines from Knots 50 and 51 have the same radial velocities, vhelio 130
km s 1 . Regions along the slit close to the two knots also show emission, especially in
[O III], from the front and rear sides of the Crab’s expanding shell with vhelio roughly
±1000 km s 1 . The blueshifted component of this additional gas falls outside the
HST passbands for the emission-line images. The +1000 km s 1 component will
19

[O III]

[O III]

[S II]

1000
-1

km s

H2 2.12µm

K51
N

K50

+"
+"

5”

E

K51

K51

103 km s-1

Figure 2.3 HST [O III] image (left panel) with vertical lines showing the area covered
by our KPNO spectrum with its 3.8” wide slit. The two small white crosses mark
the positions of Knots 50 (lower) and 51 (upper). Portions of the KPNO spectrum
are shown in the middle two panels, aligned vertically with the HST image. The
11.2” region extracted along the slit is indicated by the horizontal lines, and the
solid-lined goal posts mark the wavelengths of the [O III]
4959, 5007 and [S II]
6716,6731 doublets from Knot 51 (blended along the slit with Knot 50). The
velocity scale is the same in the two central panels, with velocity increasing to the
right. Our NIR spectrum, taken through a 0.4” wide slit and with the scale on the
sky magnified relative to the other images, is shown in the right panel. The higher
spatial resolution of the NIR spectrum clearly separates Knots 50 and 51, and shows
that they have the same radial velocity.
contribute to the line emission seen in the HST images, but the line-of-sight to K51 is
relatively uncontaminated. The spectrum shows that Knots 50 and 51 are at nearly
the systemic velocity of the Crab. In addition, where adjacent KPNO slit positions
not shown here cross parts of the long double-tendril of emission reaching down to
K51 from the NW, the spectra show that the gas all along that tendril has the same
radial velocity as K51. From this we deduce that Knots 50 and 51 and the tendril
are all part of the same physical structure that is moving mostly in the plane of the
sky, and therefore must lie very close to halfway through the overall Crab Nebula.
Knots 50 and 51 are not spatially resolved in our optical spectrum, so we measured
them together by extracting a 1D spectrum that covers the 11.2 arcsec region along
the slit marked on Figure 2.3. We also extracted spectra over smaller lengths along
the slit trying to isolate Knots 50 and 51, and confirmed that those spectra do not
have any major di↵erences between them in their emission-line intensity ratios. Here
20

we have chosen to use the combined spectrum because we know that we have included
the total flux from the sum of the two knots.
The measured line strengths relative to H and the H flux are listed in Table 2.2.
The table also lists line strengths corrected for reddening assuming AV = 1.46 mag
(Davidson & Fesen 1985) and a standard Galactic extinction curve. This dereddening
is intended to remove the e↵ects of extinction from the ISM along our line of sight
to the Crab. The resulting H↵/H

ratio is 2.96, very close to the Case B value,

implying that there is little additional internal reddening of the Balmer lines. The last
column of Table 2 lists the estimated percentage uncertainties for the lines measured
from the KPNO spectrum, including the e↵ects of blending, sky subtraction and
systematic error in the flux calibration in addition to the pixel-to-pixel signal:noise.
Table 2.2 also lists line strengths relative to H for the total H2 2.12 µm and Br
fluxes for Knots 50 and 51 combined together (since the optical line strengths are for
the sum of the two knots). These H2 and Br line strengths were measured from
our NIR images. The H2 line strength is taken from Table 2 of Paper II and has an
estimated uncertainty of 10% including the random errors and the uncertainties in
the sky subtraction and flux calibration. The Br line strength was measured for the
aperture used for the H

spectrum, since the emission occurs over a larger region.

The dereddened Br /H intensity ratio derived from combining the optical spectra
and NIR images is I(Br )/I(H ) = 0.027, in good agreement with the value 0.028
for Case B. However, it should be noted that the Br line is very weak (see Figure 3
of Paper II) and so it is not as well measured as H2 or the optical lines.
Finally, we converted the reddeningcorrected fluxes to surface brightness. The line
emission comes from a much smaller area on the sky than is included in our extracted
spectrum (Figure 2.3, left panel). The HST images in Figure 2.2 show that the H
signal:noise ratio is too low to be able to tell exactly where that line comes from, and
that the [O III] line comes from a somewhat more extended area than [S II] or [O I].

21

Table 2.2: K51 Observed and Dereddened Emission Line Fluxes

Line

Observed

Dereddened

F (Line)/F (H ) I(Line)/I(H )

Estimated
Uncertainty(%)

[O II] 3727

7.38

12.48

25

[Ne III] 3869

1.59

2.55

10

[Ne III] 3967

0.38

0.58

50

[Fe V] 4071

1.34

1.96

50

H I 4340

0.42

0.54

25

[O III] 4363

<0.25

<0.32

50

He I 4471

0.16

0.19

50

He II 4686

0.28

0.31

25

H I 4861

1.00

1.00

-

[O III] 4959

2.45

2.35

5

[O III] 5007

7.31

6.87

5

[Fe II] 5159

0.08

0.07

50

[N I] 5198

0.17

0.15

25

He I 5876

0.65

0.46

25

[O I] 6300

5.06

3.29

10

[O I] 6363

1.78

1.14

10

[N II] 6548

3.06

1.88

20

H I 6563

4.82

2.96

10

[N II] 6584

9.35

5.71

10

[S II] 6716

8.62

5.10

10

[S II] 6731

11.25

6.64

10

He I 7065

0.25

0.14

50

[Ar III] 7136

1.27

0.69

25

[Fe II] 7155

0.26

0.14

50

22

Table'2.2'(cont0d)
Line

Observed

Dereddened

F (Line)/F (H ) I(Line)/I(H )

Estimated
Uncertainty(%)

[O II] 7320

0.35

0.18

50

[O II] 7330

0.28

0.15

50

[Ni II] 7377

1.24

0.63

25

H2 2.12 µm

0.62

0.14

10

Br 2.17µm

0.12

0.03

50

MIPS 24µm (21-26µm)

57.

11.

-

PACS blue (57-83µm)

260.

50.

-

PACS green (79 121µm)

120.

23.

-

Observed flux through full aperture: F (H ) = 1.04 ⇥ 10 14 erg cm 2 s 1 .

Dereddened flux through full aperture: I(H ) = 5.4 ⇥ 10 14 erg cm 2 s 1 .

Dereddened 15% of 3.8”⇥ 11.2”: S(H ) = 8.4 ⇥ 10 15 erg cm 2 s 1 arcsec 2 .

For purposes of comparison to the simplified models discussed below, we assume that
all lines come from the area within our spectroscopic aperture that emits [S II]

6716,

6731, which is 6.4 arcsec2 . This is the projected area on the sky, so for comparison
to our Cloudy simulations our measured surface brightness corresponds to a planeparallel sheet of material seen face on. The relative dereddened line strengths listed
in Table 2.2, scaled to the H surface brightness listed in the footnotes to the table,
are the primary constraints on the models that we will describe below. Although the
line strengths are derived from the sum of the two knots, the conversion to surface
brightness allows us to treat them as if they came just from Knot 51.
We also have a NIR long-slit spectrum across Knots 50 and 51, taken with the
Ohio State Infrared Imager/Spectrometer (OSIRIS) on the 4 m SOAR Telescope
(Paper III). For K51 it shows enough H2 lines for us to measure an H2 excitation
23

temperature TH2 = 2837±200 K. For comparison, previous studies have shown that
the ionized gas has a considerably higher kinetic temperatures, Tion ⇠ 15,000 K. The
1D-extracted spectrum is shown in Paper III, but here in the last panel of Figure
2.3 we show the 2D image of the region around the 2.12 µm line. The line has a
heliocentric radial velocity vhelio = 108±5 km s 1 , in reasonable agreement with the
value vhelio ⇠ 130 km s 1 quoted above which was measured from the [S II] lines.
This further establishes the spatial connection between the ionized and molecular gas.
The 2D NIR spectrum shown in Figure 2.3 has 0.7 arcsec spatial resolution, which is
high enough to show that Knots 50 and 51 clearly do have the same velocities and in
fact may be connected by a bridge of emission. We did not attempt to combine the
H2 line fluxes measured from this spectrum with the measurements from our optical
spectrum, because the NIR spectrum was taken through a slit only 0.4” wide with a
somewhat uncertain placement on the knot.

2.2.3 Dust Absorption Features
The HST intermediate-band continuum image shown in Figure 2.2 is of interest because it shows obvious dust absorption (which appears white in our reverse grey-scale
image) coincident with K51. We wished to measure the depths of these features in
order to find the dust extinction. However, the raw data through the F547M filter
include an appreciable contribution from the night sky, which must first be subtracted
o↵. We measured this by using the few existing F547M images that include the outer
parts of the Crab together with areas that clearly are outside of the Crab’s synchrotron emission region (as determined by comparison to Space Telescope Science
Institute Digitized Sky Survey visible and Two Micron All Sky Survey NIR images).
This let us produce several individual sky-corrected images, but they did not provide
high signal:noise coverage of the K51 region, so were not by themselves suitable for
this study. However, we were able to use them to calibrate the sky-subtracted surface

24

brightness at a number of points where the synchrotron continuum is not contaminated by emission from filaments. A further problem is that the sky level varied
over time. For each image used in our final co-added continuum image, we measured
the surface brightness at a series of fiducial points on the nebula that do not include
filaments, and worked out flux o↵sets for each image which then were subtracted
to remove relative di↵erences in the sky level. After making those corrections, we
combined the images and then finally subtracted a suitable o↵set from the co-added
image so that it has the same surface brightness at each calibration spot as is measured from our individual sky-corrected images. This produced the sky-subtracted
continuum image used here. At the position of K51, this sky correction amounts to
44% of the average level of the continuum at points not in the dust features.
The deepest part of the dust absorption feature at K51 has 69% of the flux level
of the adjacent (sky-subtracted) continuum regions. We need to account for the
fact that K51 is embedded in the synchrotron emitting gas. As noted above, the
radial velocity of K51 shows that it must lie approximately half-way through the
synchrotron-emitting region, so for our best estimate we subtracted o↵ half of the
continuum level to account for emission in our direction originating between us and
the knot. After this correction the K51 absorption feature has a central depth of 62%.
There is uncertainty in this measurement due to the patchy distribution of the
synchrotron emission. To estimate the e↵ect of this on our extinction measurement,
we integrated the continuum brightness along many 2”-wide lines that lie in the plane
of the sky and perpendicular to the long axis of the nebula, and formed ratios of the
light coming from either side of the midpoint of the nebula. The RMS scatter in these
ratios was 5%, which we used as an estimate of the uncertainty in the fraction of the
continuum emission that comes from in front of K51, along our line of sight through
the knot.
The extended source complicates the conversion from the measured 1.05 mag dip,

25

to the usual AV , which includes scattering as well as absorption. In the approximation
that K51 is surrounded by a uniform light source and the optical depth is small,
scattering into the line of sight compensates for scattering out of the line of sight,
and AV is larger by a factor of (1-albedo) 1 . This correction gives us the extinction
that is usually compared to the gas column density NH and which is what would be
observed for a sheet of absorbing material in front of a background point source, so
we will designate the corrected value as AV pnt . At 5470 ˚, the albedo is 0.77, 0.67,
A
and 0.81 for the LMC, RV = 3.1, and RV = 5.5 extinction curves given by Draine
(2003), corresponding to AV pnt = 4.6, 3.2, and 5.5 mag, respectively. The actual
grain composition and size distribution in the Crab is unknown. In the Cloudy models
described below, we have arbitrarily chosen the Orion Nebula grain distribution used
by Baldwin et al. (1991), for which the albedo at 5470 ˚is 0.6. This corresponds to
A
AV pnt = 2.6 mag, which for consistency we use as our observed value for comparison
to the models even though it is at the low end of the range of possible AV pnt . This
is about the depth where H2 initially forms in PDRs (Tielens & Hollenbach 1985).
However, the combination of the unknown grain properties and the uncertain and
possibly very large correction needed to the observed absorption depth means that
the extinction measurement is only a very marginal constraint.
There is also weaker dust absorption associated with the ring that includes Knot
50. Here the central depths indicate AV pnt ⇠ 0.6 1.3 mag for the same range of
assumptions about the albedo of the dust.

2.2.4 Mid-IR Images
There are also archival mid-IR images of the Crab available at 24 µm [Spitzer Multiband Imaging Spectrometer (MIPS)], 70 µm [Herschel Photodetector Array Camera
& Spectrometer (PACS) blue] and 100 µm (Herschel PACS green). These clearly detect the Knot 50+51 complex, but at such low angular resolution that it is unresolved

26

with no measurable structural details. The Knot 50+51 complex is blended with the
smeared-out images of surrounding structures, but we have crudely measured the total excess flux from the Knot 50+51 complex above what is our best estimate of the
underlying synchrotron continuum. The results for each of these passbands are listed
in Table 2. In principal, the wavelength covered by these images should include the
thermal emission from the dust in K51; a thermal spectrum peaking at about 80 µm
is seen in the integrated light of the Crab (e.g. Marsden et al. 1984) and the models
presented below predict dust emission with the correct temperature to produce that
emission. Knowing the strength of this feature for K51 together with the AV pnt
would place interesting constraints on the nature of the dust formed in the Crab (e.g.
Sankrit et al. 1998). However, the Cloudy models described below and the spectra
shown by Temim et al. (2012) indicate that the emission from the knots in these
bands should be completely dominated by the [O IV] 25.9 µm and [Fe II] 26.0 µm
emission lines for the Spitzer MIPS band, and by [O I] 63.2 µm for both Herschel
bands.

2.2.5 Intensity of the Radiation Field Striking the
Filaments
The incident ionizing flux over the 1 5 Ryd energy range is the crucial input parameter for our Cloudy models. Wu (1981) made careful measurements of the Crab’s
continuum flux down to 1550 ˚, from which it is only a modest extrapolation to
A
912 ˚(1 Ryd). Using Wu’s Figure 2, which shows the measured fluxes corrected
A
for extinction E(B V) = 0.5 (AV = 1.55 for a standard Galactic reddening curve)
and scaled to include the full Crab Nebula, we found log F⌫ (912 ˚) =
A

22.71 and

⌫F⌫ = 6.5 ⇥ 10 8 erg s 1 . Assuming a distance of 2 kpc, this gives ⌫L⌫ = 3.1 ⇥ 1037
ergs s 1 at 912 ˚and an integrated ionizing photon luminosity Q(H) = 1.7 ⇥ 1048
A
photons/sec.
27

We used that Q(H) value to calculate the incident flux on a slab of gas that is
facing the center of the Crab. For the radial position of K51 (approximately 80%
of the way from the center to the edge), and assuming a spherical nebula of radius
1.7 pc with uniform volume emissivity in the synchrotron continuum, this predicts
an incident ionizing photon flux on the slab log( (H)) = 9.7 (in units of photons
s 1 ) on the front face, and log( (H)) = 9.0 on the back face. The emission from
the photoionized region would be dominated by the front face, so we used the higher
value of (H) as a point of comparison with our models.
The actual flux incident on K51 is fairly uncertain. An ambiguity in our incident
flux calculation is that the synchrotron emission obviously is not uniformly distributed
throughout the Crab (as is clear from the HST visible-passband continuum image,
for example), plus there could be local beaming e↵ects. The face of the slab will
preferentially see the emitting regions directly out in front of it. Also, the incident
flux on the slab has a d 3 dependence on the distance d to the Crab, which according
to Trimble (1973) is probably about 2 kpc but really only is known to lie in the
approximate range 1.5 2.3 kpc.
In addition, there are factor-of-three di↵erences between di↵erent published luminosity values for the Crab at 1 Ryd. The broad-band SED constructed by Atoyan &
Aharonian (1996; and shown in the review article by Hester 2008) smoothly follows
the predicted synchrotron SED across the 1 4 Ryd region and corresponds to the
lower 1 Ryd luminosity. Atoyan & Aharonian state that this part of their composite
spectrum is from Wu (1981). However, the results actually shown in Wu’s paper
indicate the three times higher value that we have used here. This produces a local
bump above the synchrotron spectrum, as is shown in Figure 1 of the review article
by Davidson & Fesen (1985).
Wu’s result was substantiated by later near-UV observations by Hennessy et al.
(1992). We have checked, and the di↵erence in shapes between these two possible

28

SEDs does not a↵ect the model calculations we present below. It is really only the
luminosity at 1 Ryd that matters. In summary, we found that a reasonable guess at
the ionizing continuum flux is in the neighborhood of log( (H)) = 9.7 (in units of
photons s 1 below 912 ˚), but with considerable uncertainty.
A

2.2.6 Adopted Error Bars
This section has described our best estimates of a number of observed parameters
that we will next use to constrain simulations of K51, but it is important to keep in
mind the rather large uncertainties. Although K51 is the simplest case we could find,
it still has a complicated, poorly-known geometry, which limits our interpretation of
the measured emission lines strengths. For example, the H line strength is measured
to 5-10% accuracy, and so is the H2 2.12 µm line strength, so their intensity ratio in
Table 2.2 is listed as being known to about 10% accuracy. However, these lines come
(at least in part) from physically di↵erent regions of the K51 complex so it is very
hard to know how to compare the measured intensity ratio to the predictions from
the models discussed below in which K51 is described as a very simple plane-parallel
slab. The same is true for comparisons between emission lines from high-ionization
species to those from low ionization or neutral species; the actual spatial arrangement
is far more complicated than that in our simplified models.
The conversion to surface brightness is also uncertain. The deduced surface brightness uses the projected area of the Knot 50+51 complex as if it were a thin sheet seen
face-on. If instead the same observed flux in the lines from the ionized gas comes
from a thin ionized shell on the surface of a spherical core containing the H2 , there
would be a factor of 2 times more surface area. The measured emission-line surface
brightness also could easily include limb-brightening e↵ects, due to us sighting down
a part of the knot’s surface that is especially bright because it faces the strongest continuum source. This would a↵ect the lines that are powered by continuum radiation,

29

but it is unclear whether or not the H2 would have this same limb brightening.
In addition, the reddening correction is based in part on measurements of the [S
II] intensity ratio made by Miller (1973) at a di↵erent point on the Crab. We assume
that this corresponds to overlying reddening along the line of sight to the Crab, and
have not attempted to correct for internal reddening despite the clear presence of
dust within K51. The uncertainties in the geometry also a↵ect our estimates of the
incident flux.
All of this emphasizes that we are trying to hit a poorly defined target, which must
be kept in mind as we describe our Cloudy simulations in the following sections. For
our comparisons of the simulations to the observations involving surface brightness,
we will regard a factor of two agreement to be adequate, in order to account for the
geometrical ambiguities. In contrast, the actual observational error bars are appropriate to use for some key emission-line intensity ratios, such as the ratios between
di↵erent H2 lines that were simultaneously observed through the same spectrograph
slit.

2.3 Cloudy Simulations of Knot 51
We define a successful model of K51 as one that reproduces the observed emissionline surface brightness including the ionized, neutral and molecular species, from
within the observed volume of K51. This includes the thermally populated H2 levels
discussed in Paper III. Additionally, the line emissivities of [O I] and [S II] should
spatially correlate with H2 as is seen in the narrow-band images shown in Paper II.
The key issue is to identify the physics that leads to the strong H2 emission observed
from the Crab knots.

30

2.3.1 General Methodology
For our modeling of K51, we used the open-sourced plasma simulation code Cloudy
(Ferland et al. 1998), version 10.00. In one self-consistent calculation, Cloudy tracks
the microphysics and ionization structure through the H+ , H0 and H2 zones. This
is ideal for modeling these knots within the Crab, for which we have emission-line
measurements from all three regions. We assumed that K51 is in steady state, meaning that the timescale for any dynamical changes is much longer than that of any of
the atomic processes involved. The implications of this assumption are discussed in
Section 4.
The filaments are illuminated by the hard synchrotron radiation field, which we
refer to as the incident continuum. For simplicity, we used the synchrotron spectrum
averaged over the whole nebula, by digitizing the spectrum shown in Figure 1 of
Davidson & Fesen (1985). However we note that the observed SED depends on
position, being harder in the central regions and softer in the outer regions (Mori
et al. 2004). All models also include the e↵ects of a cosmic ray energy density
representing at least the Galactic background level. The “ionizing particles” model
discussed below has a far higher cosmic ray energy density representing the particles
that produce the Crab’s synchrotron continuum emission (see §3.5).
We assumed a plane parallel geometry for all models, because it is the simplest
case and the actual geometry is unknown. Since the actual knot is illuminated from
all sides, all models are constrained to have a thickness of 1016.5 cm, which is half the
observed width on the sky of K51. The synchrotron continuum penetrates successive
H+ , H0 , and H2 layers. The outer, ionized regions attenuate the incident continuum,
eventually allowing H0 and H2 to form. A similar analysis has been conducted by
Sankrit et al. (1998), but using a cylindrical geometry, a smaller suite of observational
constraints, and, most importantly, not considering molecular gas.
We included the H2 molecule described by Shaw et al. (2005) in our calculations.
31

This model includes all rotational/vibrational levels within the ground electronic
state, and the six lowest electronic excited states that are coupled to the ground state.
Photoexcitations, collisional excitation and deexcitation, bound-bound transitions,
ortho-para conversion, collisional dissociation, high energy e↵ects and the chemistry
network are all included, self consistently with other properties of the gas. These
calculations also include the model Fe II atom described by Verner et al. (1999).
The Crab contains dust, as is evident from the thermal emission feature peaking
at ⇠80 µm in the SED of the full Crab Nebula (Glaccum et al. 1982; Marsden et
al. 1984) and from the dust absorption features seen scattered across the face of the
Crab in visible and IR images (Fesen & Blair 1990; Hester et al. 1990). Sankrit et al.
(1998) estimated that the Crab filaments might contain up to an order of magnitude
greater dust/gas ratio than the ISM, and recent Herschel results (Gomez et al. 2012)
appear to confirm this when averaged over a large area on the sky. However, we do
not have a strong constraint on the actual dust/gas ratio in the particular knot K51,
so we adopted as our default value the ISM dust/gas ratio (although we did check
the e↵ect of varying this ratio in our final models; see §3.6).
Dust provides extinction that shields the molecular regions from the incident continuum. Furthermore, dust acts as a catalyst for H2 formation. We adopted the rate
given in Cazaux & Tielens (2002) for formation of H2 through grain catalysis. We
assumed the same mix of silicate and graphite grains and grain size distribution used
by Baldwin et al. (1991) to represent the Orion Nebula dust.

2.3.2 H2 Emission in Non-Equilibrium Environments
One major goal of our models is to reproduce the observed H2 emission strength and
thermal level populations presented in Paper III. The H2 lines are very strong relative
to nearby H I lines, and the excitation diagrams indicate a temperature of around
3000 K.

32

H2 emission can be produced by a variety of processes (see Chapter 8 and Appendix 6 in Osterbrock & Ferland 2006; hereafter AGN3). Starlight fluorescent excitation is important in H II regions like the Orion Nebula. The strength of the H2
emission relative to H I recombination lines is determined by the shape of the SED,
as discussed in Paper III, with harder SEDs producing weaker H2 /H I ratios. The
SED in the Crab is exceptionally hard, but the measured (H2 2.12 µm)/Br intensity
ratio is also large. Paper III shows that with the Crab’s SED, UV fluorescence photoexcitation and H I recombination will produce a ⇠3 dex smaller H2 /H I intensity
ratio than is seen in the Crab. Fluorescent excitation is negligible in Crab filaments.
X-rays and ionizing particles penetrating into knots’ cores are possibilities. These
actually behave in a very similar manner (AGN3 Chapter 8). An ionizing particle
entering cool atomic or molecular regions will dissociate and ionize the gas. A newly
created free electron has a large amount of energy and creates a shower of secondary
suprathermal electrons with a typical energy of about 40 eV. These can excite emission
lines, heat the gas, and create further ionization / dissociation. X-rays would produce
very similar results. In the central regions of the knots the lower-energy ionizing
photons have been removed by the photoelectric opacity of the outer layers, and
the penetrating photons have high energy, roughly 10 keV (Figure 2.4). Inner shells
of the heavy elements absorb such high-energy photons. The photoelectron has a
large kinetic energy, and the inner-shell hole produces additional high-energy Auger
electrons. The result is a shower of secondary electrons producing heating, ionization,
and excitation.
The net e↵ect of the secondary electrons produced by X-rays or ionizing particles
depends on the degree of ionization of the gas. Collisions between the secondary
suprathermal electrons and thermal electrons will change kinetic energy into heat
since electron electron collisions are elastic. In this case the H2 emission will be
the result of collisional excitations by thermal electrons and the level populations

33

Figure 2.4 The spectral energy distribution in the center of the knot, including the
transmitted synchrotron continuum, for the dense core model (solid line). The dashed
line is the incident continuum. X-rays with energies E = 1 10 keV, corresponding to
= 10 3 10 4 µm, penetrate deep into the core of the knot and are preferentially
absorbed by the inner shells of heavy elements. This produces high-energy Auger
electrons and a shower of secondaries which then heat, ionize and excite the gas. This
is the mechanism responsible for producing H2 emission in the dense core model.
will appear thermal. When there are few free electrons the energy of the secondary
electrons is used to heat, ionize, and excite the gas. Direct excitation of the H2
electronic Lyman-Werner bands will produce non-thermal emission that is nearly
identical to that produced by fluorescent excitation. The result in this case will be a
highly non-thermal excitation diagram.
Other sources of heat such as shocks or dissipative MHD waves are possible. These
generate heat that raises the gas kinetic temperature. The resulting H2 emission
will be thermal if the density is high enough the H2 level populations will have an
34

excitation temperature equal to the gas kinetic temperature. At lower densities the
level populations will appear thermal but at a temperature that is smaller than the
kinetic temperature.
In the following subsections we describe a series of Cloudy models that systematically test for these di↵erent possible H2 excitation mechanisms.

2.3.3 A Basic X-ray Photoionization Model in Pressure Balance
We first computed a series of Cloudy models with constant total pressure, photoionized by the Crab’s hard synchrotron continuum. The inward pressure of the absorbed
synchrotron radiation and the gas pressure at the illuminated face are balanced by
the gas pressure within the cloud. The gas is basically in hydrostatic equilibrium
(Baldwin et al. 1991). The goal here was to establish a base-line model that would
at least properly represent the H+ zone for a simple case.
Setting the Ionization Parameter
First, the ionizing photon flux,

(H), and the initial hydrogen density, nH , were

adjusted so that the model reproduces the observed H surface brightness and the
density-sensitive [S II] 6716/ 6731 ratio (AGN3). The Crab’s filaments are made
of material expelled from the interior of the progenitor star, and therefore can have
very unusual abundances which are known to vary from point to point. As our
starting point we adopted the chemical abundances from Model 3 of Pequignot &
Dennefeld (1983), which represents gas in the helium rich band, along with a solar
Fe abundance and a solar Ar abundance. For reference, these initial abundances are
listed here in Column 2 of Table 3.
For a solar composition H or He absorb most hydrogen-ionizing photons, and there
is a simple linear relationship between (H) and the H surface brightness (AGN3).
For highly enhanced abundances, such as those found in the Crab, this is not the case

35

Table 2.3: K51 Chemical Abundances

Element

Initial Abundances Final Abundances Final Abundances Final Abundances2
by Number1

by Number1

Mass Fraction2

Relative to Solar

He

-0.16

-0.53

-0.27

0.28

C

-2.80

-3.40

-2.66

0.02

N

-3.84

-4.25

-3.45

-0.37

O

-2.62

-3.28

-2.42

-0.16

Ne

-3.11

-3.74

-2.78

0.07

Mg

-4.10

-4.70

-3.66

-0.43

Si

-

-5.05

-3.95

-0.79

S

-4.38

-4.71

-3.55

-0.16

Cl

-

-7.33

-6.11

-0.79

Ar

-5.60

-5.32

-4.06

0.09

Fe

-4.49

-4.61

-3.20

-0.25

1

Relative to H.

2

Fraction of total mass.

36

since heavy elements also remove ionizing photons. The (H) required to reproduce
the H surface brightness does depend on the heavy-element composition. As noted
previously in §2.5, the observations suggest that log (H) is roughly 9.7, but with
considerable uncertainty. This can be varied to some extent to match observations.
This set of abundances gave an optimal value of log( (H)) = 10.06 and log(nH ) = 2.97
cm 3 .
Abundances in Ionized Gas
Determining abundances is not the primary motivation for this work, rather we are
trying to understand why there is strong H2 emission. There is an extensive literature
about the abundances in the ionized gas, and we do not have useful measurements
of the [O III] temperature indicator that would be needed to accurately determine
abundances for K51. However, there are known to be point-to-point abundance variations within the Crab (see the introduction), so we did need to carry out a rough
determination of the overall metal abundances.
A preliminary adjustment was made to the abundances of O, N, Ne, S, and Ar
using the optimizer included within Cloudy to match to match the strengths of [O
II] 3727, [Ne III] 3869, [O III] 4959, [O III] 5007, [N II] 6584, [S II] 6716, [S
II] 6731 and [Ar III] 7135 to within a factor of two. The Fe abundance was not
allowed to vary because doing so would have greatly slowed the calculations, so we
varied the Fe abundance to match [Fe II] 7155 after using the optimizer. Silicon and
chlorine, for which no lines are observed, were included because of their importance
to the chemistry network but with solar abundances. The helium abundance was
scaled linearly to match the line strengths of He II

4686 and He I

5876. The

derived He/H is a bit lower than found in some other regions of the Crab such as
the helium rich belt, but within the range of values reported elsewhere. The model
reproduces the observed He II / He I ratio. This ratio is mainly sensitive to the shape

37

of the SED (AGN3; Chapter 2), suggesting that the Davidson SED is a reasonable
fit to reality. The values of (H) and nH were iterated to maintain the observed H
surface brightness and [S II] intensity ratio.
These preliminary models still produced [O III], [Ne III] and [Ar III] forbidden
lines that were systematically too weak relative to H , so further abundance adjustments were required. Increasing the heavy element abundances will not produce
stronger forbidden lines in photoionization equilibrium. These lines are among the
strongest gas coolants so raising their abundances will only make the gas cooler, the
so-called “thermostat e↵ect. However, decreasing all metal abundances raises the gas
temperature and shifts the cooling to the optical. This produces stronger optical
forbidden lines. Figure 2.5 shows how surface brightness predictions, along with the
H2 temperature and S II ratio, vary as a function of scaling all metal abundances by
a given factor.
As the metal abundances decrease, the forbidden line strengths increase until the
cooling shifts to collisional excitation of Lyman hydrogen lines. The strengths of these
forbidden lines are quite sensitive to the gas temperature. The best abundance set
produces Te =12500 K in the O++ -emitting zone of the ionized gas, very similar to
previous modeling results (Davidson 1979, Henry & MacAlpine 1982, Sankrit et al.
1998). Unfortunately, our observations do not detect the [O III] 4363 line, which
could directly measure the temperature. The continuum in this region is both noisy
and contaminated by emission from other velocity systems, so only a poor upper limit
can be obtained. The models predict that the line should be a factor of two or more
below that upper limit.
The best abundance set obtained from this process and used throughout the remainder of this paper is listed in Column 3 of Table 2.3. It corresponds to scaling
all metal abundances down by a factor of 4 relative to those obtained in the preliminary optimization. These abundances lie between the most metal poor (Henry &

38

Figure 2.5 The e↵ect of scaling metallicity in our constant pressure model, by number relative to the abundances obtained from our preliminary round of abundance
optimization. As the metallicity decreases, the electron temperature rises and cooling shifts to the optical. The net e↵ect is an increase in the optical forbidden line
strengths. We adopt a scaling of 0.25, which gives the best fit to the higher ionization
lines. For interpretation of the references to color in this and all other figures, the
reader is referred to the electronic version of this dissertation.
MacAlpine 1982) and metal rich (Pequignot & Dennefeld 1983) abundance sets found
in previous studies.
Sensitivity to the Ionizing Flux
A change in (H) a↵ects the ionization structure of the knot and also the strengths
of higher ionization lines, in particular, [O III] 5007 and He II 4686. Larger values
of (H) cause regions that form these lines to extend deeper into the knot. In Figure
2.6a we show the sensitivity of various quantities with respect to the ionizing flux for
39

the constant pressure model. The ionizing flux cannot increase since it leads to the
H+ region becoming substantially larger thereby overestimating [O III] 5007 and He
II 4686. Decreasing (H) doesn’t provide enough ionizing photons and all lines get
weaker, except for the H2 2.12 µm line which becomes bright.
Despite producing strong H2 emission, this low (H) model is not able to predict
the observed H2 level populations. This is illustrated in Figure 2.7b, which compares
the observed H2 level population diagram to the predicted diagram for a constantpressure model with log( (H)) =9.66. The observed populations in Knot 51, and also
in the other 6 knots observed for Paper III appear thermal, but the higher-energy

Figure 2.6 The predicted/observed ratios of surface brightness or intensity ratios of
key emission lines and of TH2 , as a function of the ionizing flux (H), for the constantpressure case (a), dense core case (b), temperature floor case (c), and ionizing particles
case (d). In the ionizing particles case, H and He I depend only weakly on (H)
because ionizing particle heating accounts for over 50% of the surface brightness.
In both cases the H2 weighted core temperature is set almost entirely by the extra
heating in the core rather than by the ionizing flux.
40

Figure'2.6'(cont0d)

41

Figure'2.6'(cont0d)

42

Figure'2.6'(cont0d)

levels in this series of models have a decidedly non-thermal population. This produces
a predicted H2 spectrum that is very di↵erent from what is observed. The Cloudy
simulation includes all known resonance and fluorescence processes, and reproduces
the relative intensities given by Le Bourlot et al. (1995), Sternberg & Neufeld (1999),
and Black & van Dishoeck (1987) (see figure 21 of Shaw et al. (2005)) when the
same input parameters are used, so Cloudy is able to correctly handle the physical
conditions found here.
The strong 2.12 µm emission in the low

(H) constant pressure model is obtained

by combining a very high column density of H2 (nearly 1021 cm 2 ) with a core that
is very dense. The lower temperature is the result of the high density since the cooling
per unit volume depends on the square of the density. The observed H2 emission in
this case is the result of secondary electron excitation of electronic H2 levels that
43

Figure 2.7 Observed and predicted H2 level populations, normalized to the 1-0 S(1)
2.12 µm line. Panel (a): observed values for Knot 51 (from Paper III), showing
four detected H2 lines with error bars, and upper limits for two undetected lines.
Panels (b), (c), (d), (e) and (f): predicted values for five models, showing that the
constant-pressure model (here with log( (H)) = 9.66) and dense core model (here
with log(ncore ) = 6.0) do not fit the observations, while the temperature floor, ionizing
particles, and 5x Dust models do fit.
produces a fluorescence spectrum very similar to starlight pumping. This results
in a highly nonthermal excitation diagram, in conflict with the observations, which
together with the weakness of the lines from the H+ zone rules out this low (H)
model.
Final Constant Pressure Model
The constant pressure model that best fits the H+ zone has log( (H)) = 10.06 and
metal abundances decreased by a factor of 4 from our initial values. We explore it
44

further here. The physical conditions, line emissivities, and hydrogen structure for
this model as a function of depth into the cloud are shown in Figure 8a. The top
panel shows the number density of the indicated hydrogen state as a fraction of the
total hydrogen density including all states. The middle panel provides the electron
temperature as well as the electron and hydrogen number densities.
The depth dependence of the emissivity, or net emission produced per unit volume
at each point and escaping the cloud, is given in the bottom panel of Figure 8a. The
surface brightness in each of these lines is not shown in the figure, but is found by
integrating the line’s emissivity over all depths through the cloud. Examination of
Fig. 8a shows that in the constant-pressure model the gas is mostly ionized with
very little H0 and no fully molecular zone. The remaining parts of Figure 8 show
additional models that will be discussed below.
To quantify the comparison to observations, Table 2.4 shows the ratio of the predicted/observed surface brightness of various emission lines, and also some key Cloudy
input parameters, for this and the series of additional models described below. The
models are compared to the observations by considering predicted/observed surface
brightness ratios for each separate line in order to give a clearer picture of which lines
are properly reproduced by which models in a situation where the strength of H
(typically used as the denominator in tables of emission-line strengths) also varies
from model to model. A value of twice the predicted AV pnt is used in Table 2.4
because our models only are computed to the center of the knot, while the observed
extinction is for all the way through the knot.
The constant-pressure model matches the observed surface brightness of the higherionization forbidden lines and of the recombination lines to within ±30%. However,
the [O III], [Ne III] and [Ar III] lines are still systematically 10-20% too weak relative
to H , suggesting a residual energy balance problem.

45

Figure 2.8 Each block of three panels shows the hydrogen ionization structure (top
panel), Te [K] and densities [cm 3 ] (middle panel), and line emissivities [erg s 1
cm 3 ] (bottom panel), as a function of depth. The di↵erent blocks show results for
the constant pressure model (a), the dense core model (b), the temperature floor case
(c), and the ionizing particles case (d). Comparison of the top panels in each block
shows that of the three models which produce H2 emission, the dense core model
has a substantial molecular fraction in its center, while in the temperature floor and
ionizing particles cases the H2 is only a small trace constituent in an extended H0
zone which emits H2 as well as low-ionization and neutral forbidden lines.

46

Figure'2.8'(cont0d)

47

Figure)2.8)(cont1d)

48

Figure)2.8)(cont1d)

49

Table 2.4: Model Parameters and Predicted/Observed Surface Brightness Ratios

Parameter

Constant

Dense Core

Pressure log(ncore ) = 6.0

Temperature

Ionizing

Floor

Particles

5⇥ Dust

log (H) (cm 2 s 1 )

10.06

10.06

10.06

10.06

10.06

log Thickness (cm)

16.5

16.5

16.5

16.5

16.5

Additional Heating

None

None

log ncore (cm 3 )

3.6

6.0

5.1

5.25

5.25

Tcore (T in last zone)

6890

108

2837

2890

2880

Predicted AV pnt

0.1

1.7

0.3

0.4

2.0

log NH2 (cm 2 )

10.6

20.5

16.9

17.0

17.2

< n(H+ )/n(Htot ) > 5.8⇥10 1

1.9⇥10 2

1.0⇥10 1

1.0⇥10 1

1.0⇥10 1

< n(H0 )/n(Htot ) >

4.2⇥10 1

5.6⇥10 1

9.0⇥10 1

9.0⇥10 1

9.0⇥10 1

< n(H2 )/n(Htot ) >

1.7⇥10 9

4.2⇥10 1

6.3⇥10 4

5.5⇥10 4

8.6⇥10 4

ne /nH in core

1.0⇥10 1

7.8⇥10 5

5.7⇥10 3

1.6x⇥10 2

1.38⇥10 2

Line

Temp. floor ⇠H /⇠0 = 105.3 ⇠H /⇠0 = 105.3

Predicted/Observed Surface Brightness

H2 2.12 µm

0.0

0.9

0.9

1.1

1.5

[O II] 3727

0.7

0.6

0.7

1.0

0.8

[Ne III] 3869

0.8

0.6

0.7

1.1

0.9

50

Table)2.4)(cont1d)
Parameter

Constant

Dense Core

Pressure log(ncore ) = 6.0

Temperature

Ionizing

Floor

Particles

5⇥ Dust

H I 4340

0.8

0.8

0.9

2.0

1.3

He I 4471

0.7

0.7

0.8

1.2

0.9

He II 4686

1.3

1.0

1.2

1.2

1.0

H I 4861

1.0

1.0

1.1

2.3

1.6

[O III] 5007

0.9

0.7

0.8

1.2

0.9

[N I] 5198

0.4

0.2

0.3

1.6

1.3

He I 5876

0.9

0.9

1.0

1.6

1.2

[O I] 6300

0.2

0.1

0.2

0.7

0.5

H I 6563

1.0

1.0

1.1

2.5

1.8

[N II] 6584

0.4

0.3

0.3

0.5

0.4

[S II] 6716

0.3

0.2

0.3

0.8

0.7

[S II] 6731

0.3

0.2

0.3

0.9

0.7

He I 7065

0.8

0.8

0.8

1.3

1.0

[Ar III] 7136

0.8

0.7

0.8

1.0

0.8

[Fe II] 7155

0.6

0.5

0.6

2.0

1.5

51

Table)2.4)(cont1d)
Parameter

Constant

Dense Core

Pressure log(ncore ) = 6.0
[O II] 7320
H I Br

2.17 µm

Temperature

Ionizing

Floor

Particles

5⇥ Dust

1.2

1.2

1.2

2.3

2.2

0.9

1.1

1.1

2.1

1.7

52

There has to be more heating per unit volume that is provided by the observed
SED, (H), and nH . Previous modeling of the filaments was not so observationally
constrained and thus did not encounter this problem. Throughout the process of
fitting these models we have accepted cases in which the low-ionization forbidden
lines ([O I], [O II], [N II], [S II]) are underpredicted if H2 is also underpredicted,
because these forbidden lines are observed to be H2 tracers (Paper II). Our aim is
to use the constant pressure model as a baseline to which we will then add an H2
production mechanism which we hope will also contribute lines such as [S II] and [O
I]. Indeed, this best-fitting constant pressure model makes a negligible amount of H2
2.12 µm emission (Table 2.4), ruling it out as a complete model.

2.3.4 Models with Dense Cores
The results from the constant-pressure model suggested that a denser core is needed
to have enough column to produce the observed H2 emission from within the size
of K51 while still matching the lines from the H+ zone. We provided this dense
core by introducing a density profile that followed the constant pressure law through
the ionized zone, but then quickly ramped up to a higher, constant central core
density ncore H atoms cm 3 . The ramp-up started at a depth of 1016.37 cm which
is the point at which most of the [S II] emission has occurred (so that the density
determined from the [S II] doublet ratio would not exceed the measured value), and
reached the constant value at ncore at a depth of 1016.50 cm.
We experimented with a range of ncore values. The results are given in Figure
2.9, where the solid line shows the predicted/observed H2 intensity ratio. We found
that log(ncore ) = 6.0 reproduces the observed H2 2.12 µm line. This is our best-fitting
model of a dense knot that is powered only by the incident synchrotron continuum
from the Crab. The physical conditions, line emissivities, and hydrogen structure
for this model are shown in Figure 2.8b, and predicted/observed line intensities are

53

Figure 2.9 The predicted/observed ratio of the surface brightness of the 1-0 S(1) line
of H2 (solid line) and the predicted/observed temperature (dashed line), as a function
of core density ncore , for the dense-core models. The factor-of-two error bars shown
for the H2 surface brightness reflect the uncertainty in the geometry of the knot. TH2
is the H2 -weighted kinetic temperature. The model with the best-fitting H2 surface
brightness is for log(ncore ) = 6.0
given in column 3 of Table 4. Note that the line emissivities of H and He I, shown
in the bottom panel of Figure 2.8b, actually increase in the core, where the electron
temperature Te ⇠ 500 K. The predicted line strengths only vary 20% from Case B
showing that these lines are mainly due to recombinations.
This best-fitting dense-core model reproduces the observed H I, He I, He II and
higher ionization forbidden lines as well as, and in some cases better, than our best
constant pressure model, but now also reproduces the observed H2 2.12 µm emission.
However, like the low (H) constant-pressure case discussed above, it conspicuously
54

fails to reproduce the relative strengths of other H2 lines for which measurements
or observational upper limits are given in Paper III. This is shown in Figure 2.10,
which displays the predicted/observed ratios for the other H2 lines relative to the
2.12 µm line, as a function of ncore . The 1-0 S(0) line is predicted to be excited by
non-thermal mechanisms and to be several times stronger relative to the 1-0 S(1) 2.12
µm line than is actually observed. This is further illustrated in Figure 2.7c, which
compares the observed H2 level population diagram to the predicted diagram for the
log(ncore ) = 6.0 dense core model. The left most set of points on Figure 2.10 show
that the dense-core case which comes the closest to fitting the observed H2 relative
intensities has a lower core density, log(ncore ) = 4.2. However, Figure 2.9 shows that
this model underpredicts the 2.12 µm line by more than two orders of magnitude so
it is not viable.

2.3.5 Models with Additional H2 Mechanisms
None of the above models get the H2 -emitting core hot enough to produce the observed
bright, thermal H2 spectrum. Additional heating is needed in the core. As was
suggested by G90, the relativistic particles that fill the Crab’s synchrotron-emitting
plasma are a likely source. We proceeded as in Ferland et al. (2009) and considered
two possibilities. The first is a generic heating source corresponding to mechanical
energy input, possibly due to shocks or the presence of non-dissipative MHD waves.
We will refer to it as the “temperature floor” case. In this model any collisional
ionization would be produced by thermal particles so warm temperatures are needed
to produce, for instance, H+ . The second possibility that we considered is heating by
relativistic particles, which we will refer to as the “ionizing particle” case. Ionizing
particles, as suggested by their name, can produce ions in gas that has a low kinetic
temperature. For more details see the study by Ferland et al. (2009) of the filaments
seen in the intergalactic medium (IGM) of cool-core galaxy clusters.

55

Figure 2.10 The predicted H2 line ratios divided by the observed H2 line ratios from
Paper III, for the series of dense core models. In a perfect fit, all predicted/observed
ratios would be 1. The two curves with upward arrows are cases where the observed value for the line in the numerator is only an upper limit, so that any predicted/observed ratio smaller than 100 represents an acceptable fit. For each of the
other curves, the error bars shown at a single point are the same size for all points.
The model at log(ncore ) = 4.2, fits the observations the best but it is clearly not ideal,
and Figure 2.9 shows that the 2.12 µm emission is too weak.
Temperature Floor Case
In the temperature-floor case, the electron temperature is artificially fixed at the observed H2 temperature, TH2 = 2800 K, corresponding to a heating rate of 2.5 ⇥ 10 16
ergs s 1 cm 3 . The only source of ionization in the deeper zones of the knot is
through thermal collisions, and the dissociation of H2 occurs only when the temperature is quite high. We included ionization by cosmic rays at the Galactic background
rate of 2 ⇥ 10 16 s 1 (Indriolo et al. 2007), but not the e↵ects of the additional
56

high-energy particles localized in the Crab Nebula. We iterated over a range of core
densities until we reproduced the observed H2 2.12 µm surface brightness. From this
procedure, we found the optimal core density to be ncore = 105.1 cm 3 .
Ionizing Particle Energy Density
In the ionizing particles case, we included the e↵ects of the high-energy particles
that must pervade the Crab’s synchrotron plasma. We assumed that these particles
penetrate into the dense knots, heating, ionizing and exciting the gas. In regions
with high electron fractions, energetic particles predominately heat the gas, while
in regions with low electron fractions, they predominately ionize and excite the gas
(AGN3; Ferland et al. 2009). We add ionizing particles as in Ferland et al. (2009).
Cloudy uses the results of Dalgarno, Yan & Liu (1999) to include the e↵ects of cosmic
rays. The key input parameter for Cloudy in this process is the ionization rate per
H2 . We parameterize this rate as a ratio relative to the ionization rate in the ISM,
⇠0 = 2 ⇥ 10 16 s 1 .
We ran Cloudy models for a large grid of ionization rates per particle ⇠H and core
densities ncore to find the models which reproduce the observed H2 2.12 µm surface
brightness and the observed H2 temperature. Figure 2.11 compares the constraints
placed by the observed H2 temperature and by the H and H2 surface brightnesses in
the ⇠H /⇠0

ncore parameter space. The slanting dashed line marks models which re-

produce the observed temperature TH2 =2800 K, with the gray shaded area indicating
the range due to the observational uncertainty in the temperature (as given in Paper III). The solid contours show the logarithm of the predicted/observed H2 surface
brightness, while the dotted contours show the logarithm of the predicted/observed
H surface brightness. A perfect model would have the dashed line simultaneously intersecting with the 0.0 contours for both H2 and H , but in fact all models overpredict
H . We adopted a scale factor of ⇠H /⇠0 = 105.3 and a core density of ncore = 105.25

57

Figure 2.11 The optimization of core density and ionization rate scale factor for the
ionizing particle case. Models reproducing the measured H2 kinetic temperature fall
along the dashed line while the gray band indicates the range due to the error bars
of the measurement. The log contours show the ratio of the predicted to observed
surface brightness of H2 (solid line) and H (dotted line). We adopt ⇠H /⇠0 = 105.3
and ncore = 105.25 cm 3 as the optimal parameters despite the overestimated H
emission.
cm 3 to be the best parameters for the ionizing particles case. This model overpredicts the H surface brightness by approximately a factor of 2.
We next check to see if the deduced value ⇠H /⇠0 = 105.3 is reasonable. Positively
charged nuclei may be an important component of the ionizing particles (e.g. Hoshino
et al. 1992; Gallant & Arons 1994), but there is no way to directly measure their presence. For the electrons and positrons, on the other hand, the energy distribution can
be determined from the observed synchrotron radiation, so in the following discussion
we will consider just the contribution of these particles. In the case of the ISM this
underestimates the cosmic ray energy density by about 1 dex (Webber 1998).
58

The number of ionizations that the electrons and positrons will cause within the
knot depends on the degree to which they are trapped by tangled magnetic fields.
One limit is to assume that the magnetic field threads though the knot simply so
that particles make only one pass through the knot. Then the flux of particles inside
the knot is the same as the flux on the outside where it can be deduced from the
synchrotron radiation. Using the number distribution of electrons as a function of
energy found by Atoyan & Aharonian (1996) and the ionization cross section of atomic
H (Indriolo, Fields & McCall 2009) scaled up by a factor 2.3 to account for ionization
by secondary electrons and the di↵erence between H and H2 (Glassgold & Langer
1974), and integrating over all energies above 50 keV, the minimum energy for which
a particle reaches the center of the knot, we found ⇠H /⇠0

103.2 . This is 100 times

smaller than is inferred from the Cloudy models.
An upper limit from the electron/positron contribution is to assume complete
trapping so that all of the energy carried by ionizing particles entering the face of
the knot is used up in heating or ionizing the large column density of neutral gas.
Assuming energy equipartition with a 300 µG magnetic field (Marsden et al 1984)2
the ionizing particle energy density in the synchrotron plasma is UIP = 2000 eV
cm 3 . The energy flux into the knot is cUIP , and assuming that there is about one
ionization per 40 eV deposited, the total ionization rate for a knot of half-thickness
lknot = 1016.5 cm will be ⇠H (cUIP )/(40lknot ncore ) ⇠ 1 ⇥ 10 9 s 1 . This corresponds
to ⇠H /⇠0  106.8 , about 30 times larger than the model requirements. The value
⇠H /⇠0 = 105.3 found for our ionizing particle model lies within the range set by these
two limiting cases.
For energy equipartition with the magnetic field, the energetic electron energy
density corresponds to a pressure of 8 ⇥ 106 K cm 3 . This is almost two orders of
2 Recent

estimates for the magnetic field at the radial distance of Knot 51 range
from 240 µG to produce the synchrotron emission (Atoyan & Aharonian 1996) to 540
µG in the MHD simulations by Hester et al. (1996).

59

magnitude less than the gas pressure in the core, so even for the case of maximum
trapping the ionizing particles are not dynamically important and hence do not a↵ect
the geometry of the gas.
Extended H0 Zones
A crucial result for both cases is that hydrogen in the core of K51 is almost entirely
atomic; only a trace amount of H2 is required to reproduce the observed 2.12 µm
line. The lower panels of Figures 2.8a-2.8d show the emissivity as a function of depth
for several key lines. In the temperature floor and ionizing particle cases, the H
and [S II] line emissivities increase in the H0 /H2 zone. The middle panels in Figure
2.8a-2.8d displays the temperature and density as a function of depth.
We also found that less than half of the H2 formation proceeds on grains. H2 is
mainly formed through associative detachment: H + H0 ! H2 + e

(see Ferland,

Fabian & Johnstone 1994). The main H2 destruction paths include UV fluorescence
and reverse molecular reactions involved in the formation of H2 . The dependence of
these formation and destruction mechanisms on depth is shown in the middle panels
and bottom panels, respectively, of Figure 2.12.
Heating and Cooling
Figure 2.13 identifies contributors to the total heating and cooling for both cases.
Photoionization of H or He is the dominant heating process in the knot until a depth
of 2.7⇥1016 cm where the temperature floor is imposed in the temperature floor case,
and to a depth of 6.3 ⇥ 1015 cm in the ionizing particles case. The heating sources in
the core for the temperature floor model are not meaningful, because the temperature
has been artificially set to a fixed value. In the ionizing particles case, the H0 /H2
zone is heated primarily by the ionizing particles, followed by a 10% contribution by
ionization of Fe0 . Cooling proceeds mainly though electron collisional excitation of

60

Figure 2.12 The H2 formation mechanisms (middle panels) and H2 destruction mechanisms (bottom panels), as a function of depth for the temperature floor case (a) and
the ionizing particles case (b). A large fraction of H2 forms through H + H0 associative detachment (solid line in middle panels) rather than grain catalysis (“grains”).
This is due to the relatively large ne /nH ratio in the core (see middle panels, Figure
2.8). The hydrogen ionization structure is shown in the top panels for reference.

61

[O III] lines in the H+ zone and of [Fe II] lines in the H0 /H2 zone for the ionizing
particles case.
E↵ect of the Ionizing Flux in Extra Heating Models
As explained in §3.3.3, the ionizing photon flux,

(H), most greatly a↵ects higher

ionization lines, in particular, [O III] 5007 and He II 4686. Larger values of (H)
allow regions that form these lines to extend deeper into the knot. Figure 2.6 shows
the dependence of various lines and physically relevant parameters on (H) for the
temperature floor and ionizing particle cases. We used log ((H)) =10.06 in the models
presented in Table 2.4. We note that H and He I are not strongly correlated with
(H) in the ionizing particles case since cosmic-ray heating accounts for over 50%
of the surface brightness. Similarly, the gas kinetic temperature weighted by the H2
abundance at each point, TH2 , is almost entirely set by ionizing particle heating with
the ionizing photon flux having only a small e↵ect.
Results for Models with Additional H2 Excitation
Table 2.4 shows that the core temperatures, H2 column densities and AV pnt for the
temperature floor and ionizing particle models are very similar. Both models produce
a satisfactory H2 emission spectrum (as shown in panels (d) and (e) of Figure 2.7).
The models are largely equivalent: the ionizing particles have provided a physical
reason for the temperature floor.
The di↵erences in the results from these two models come from their ability to reproduce the observed neutral and recombination line emission. In the ionizing particle
case, the [O I] emission is correct to within 30% and [N I] is overpredicted by 60%.
Ionizing particles ionize and heat the gas, producing stronger lines in the neutral zone.
In this regard, the ionizing particles case is more successful than the temperature floor
case which produces only about 20% of the [O I] and [N I] emission. The reason for

62

Figure 2.13 The heating fractions and cooling fractions as a function of depth, for
the temperature floor case (a) and the ionizing particles case (b). The knot is heated
by photoionizations followed by recombinations of the di↵erent ions indicated in the
legend in the H+ zone. But, heating by ionizing particles dominates in the H0 /H2
zone for the ionizing particles case.

63

this stems from the temperature at which the H2 is observed. Neutral oxygen requires
temperatures of about 7000-8000 K to thermally excite substantial emission, while
the measured H2 temperature is only about 2800 K. In the temperature floor case, the
[S II] doublet is underpredicted by a factor of 3, however in the ionizing particles it is
correct to within ⇠15%. The extra [S II] emission in the ionizing particles case comes
from the core, which we explain in detail below. Finally, the temperature floor case
correctly reproduces the [O II] 7320 line but the ionizing particles case overpredicts
it by a factor of two. However, there is considerable uncertainty in this particular
measurement.
Conversely, the H0 and He0 recombination lines are overpredicted by approximately a factor of two in the ionizing particles case. This is a failure of the ionizing
particles case, but the temperature floor case correctly predicts all recombination
lines to within 20%. The discrepancy here arises from the physics in the core. Figure
2.8 displays the emissivity of H for both cases (shown as the dashed black line in
the bottom panels). The ionizing particles case features an increase in H emission
in the core due to ionization of H0 , the majority of the hydrogen, followed by recombinations. The temperature-floor case does not display this feature since the gas is
too cool to collisionally ionize H. In this aspect, the temperature floor case is more
successful than the ionizing particles case. The existing HST H images have too low
signal/noise ratio to be able to tell whether or not there is emission from the core.
Future high-resolution H imaging with better signal/noise ratio will help to decide
whether the ionizing particle case is a viable model.
A successful feature of both the temperature-floor case and the ionizing-particles
case is that part of the [S II]

6720 doublet is emitted from the same region as

molecular hydrogen (bottom panel of Figure 2.8c and 2.8d). This is observed to be
the case in Knot 51 (Fig 2), and in Paper II we showed that many other, fainter
regions of the Crab’s system of filaments produce H2 emission that correlates with [S

64

II]. Our two viable models are consistent with this observation. Conversely, Figure
2.8a shows that the constant pressure model does not emit any H2 and Figure 2.8b
shows that the dense core model emits H2 only from a region largely separate from
that which emits [S II].
Both the temperature floor and ionizing particle cases are able to reproduce the
major observational signatures of the H2 knots. However, each of these models has
failures as well as successes, making it unclear which model is to be preferred. In
both cases, we are able to successfully reproduce the forbidden lines from the ionized
gas to within a factor of two, the H2 2.12 µm line to within 10%, and the correct
relative strengths of the observed H2 lines.
A variation on the ionizing particles model is to accept a worse fit to the H2 2.12
µm surface brightness, but one that is still within the uncertainty of the measurement,
in order to lower the recombination line strengths enough so that they fit the observations to within a factor of two. A model with log(⇠H /⇠0 ) = 5.2 and log(ncore ) = 5.1
produces only 70% of the observed H2 emission, but fits all the optical lines to within
a factor two except for [O II] 7320. However, the purpose of this paper is to explore
models that match the H2 emission. Therefore, this model which purposefully did
not attempt to fit the H2 emission will not be considered any further.

2.3.6 Models with Increased Dust Abundance
The models with extra H2 excitation mechanisms reproduce the observed H2 spectrum along with the surface brightness of many optical lines. However, these models
underpredict the observed AV pnt ⇠3 mag. Recent estimates of the Crab’s total
dust content range from 0.002 M

(Temim et al. 2012) to 0.25 M

Gomez et al.

2012), and the dust/gas ratio is known to vary from place to place within the filamentary structure (Sankrit et al. 1998). Sankrit et al. (1998) estimate the mass
ratio Mdust /Mgas to be as large as ten times the ISM value for one filament. We
65

do not have additional information that allows us to constrain the dust abundance
in the Crab. However, we tested the robustness of our results by varying the dust
abundance, ranging from dust-free to an order of magnitude above the ISM value. For
simplicity, in the ionizing particles case, we did not reoptimize the heating rate and
core density as this is a computationally expensive task. Instead, we simply explore
the parameter space here to present the sensitivity of various relevant values.
The top panel of Figure 2.14 shows the fraction of H2 formation due to grain
catalysis (“grains”) and associative detachment (“H +H0 ”) at the center of the knot
and the predicted/observed AV pnt ratio measured clear through the knot, plotted
as a function of dust abundance for the temperature floor and ionizing particle cases.
In both cases, the H2 formation rate by grain catalysis eventually becomes faster than
associative detachment. The observed AV pnt is found at approximately 10 times
the ISM dust abundance. The bottom panel of Figure 2.14 shows a slow decrease in
the optical line strengths as function of dust abundance due to increased extinction.
In contrast, H2 emission increases with dust
abundance due to an increase in the H2 abundance due to grain catalysis and shielding.
An ionizing particles model with a dust abundance a factor of 5 above the ISM
value gives a reasonable fit to observed surface brightnesses as well as matching the
observed AV pnt to within about a factor of three (Table 2.4). Internal extinction
removes many of the H photons generated in the central core and brings this line into
better agreement with the observations. This model matches all the higher ionization
forbidden lines to within 20%, all recombination lines to within at least a factor of 2,
and the thermal H2 spectrum (Figure 7f). However, it does fail to produce enough
[N II] 6584 and [O II] 7320.
Table 2.4 does not include the Spitzer or Herschel broad-band measurements because these are extremely crude measurements which combine continuum emission

66

Figure 2.14 The sensitivity of several key parameters as a function of the grain abundance relative to the ISM grain abundance. The top panels show the formation fraction at the center of the knot due to the grain catalysis (“grains”) and the H + H0
formation routes and twice the predicted visual extinction AV pnt relative to the
observed values. At a grain abundance of ⇠3.0 times the ISM value, grain catalysis
overtakes associative detachment in the temperature floor case. The bottom panels
show the predicted over observed values for various quantities as a function of grain
abundance.

67

with a variety of emission lines, and therefore do not deserve the same weight as the
individual lines listed in the table. However, the predictions from the extra heating
models are 25 times too low for the temperature floor model, within a factor of two
for the ionizing particles model, and 14 times too high for the 5⇥ dust model. We
regard all of these as being satisfactory matches.
We briefly explored a model with no dust at all, to determine whether su cient
H2 would be formed by routes not involving dust. Although it is clear that K51 must
contain some amount of dust because we see dust absorption in the HST continuum
image (§2.3), there are other H2 knots in the Crab that do not have detectable dust
absorption and which therefore might be described by such a model. As Figure 14
(bottom panels) shows, models without dust can still produce strong H2 emission.
For example, a dust-free model with density ncore = 105.2 cm 3 at the center of the
knot and a temperature floor is able to correctly reproduce the observed H2 2.12 µm
surface brightness and thermal H2 spectrum. It underpredicts [N I] 5198, [N II]
6584, [O I] 6300, and [S II] 6720 by about the same amount as in the temperature
floor model described above (in which dust is present at the ISM abundance). This
dust-free model is not intended to describe K51, but it does show that in principle
dust-free knots could form H2 and produce strong H2 emission.

2.4 Discussion
2.4.1 Understanding the H2 Formation Environment
For both the temperature floor and the ionizing particle models, we find the surprising
result that the strong H2 emission comes from a region where hydrogen is predominantly atomic. The overwhelming reason is that the core is insu ciently shielded
from the Balmer continuum to allow H2 to fully form. The energetic photons penetrate into the core and produce a region where hydrogen is mainly atomic but there

68

is a significant electron density.
The H2 formation mechanism in this environment is primarily associative detachment (H + H ! H2 + e). This is because the product of the electron density
ne and the hydrogen density nH is large enough to allow H2 formation to proceed
faster by associative detachment rather than by the usual grain catalysis method
(2H + grain ! H2 + grain). For typical conditions present in the ISM at the H0 /H2
interface,
rH
rgrain

=

ne ↵ H
ne
⇡ 250
nH ↵grain
nH

(2.1)

(Ferland et al. 1994 equation A10), where rgrain and rH are the rates, and ↵grain
and alphaH are the rate coe cients for grain catalysis and associative detachment,
respectively. For ne / nH > 4 ⇥ 10 3 associative detachment is faster than grain
catalysis at forming H2 . In the cases we have considered in this paper, ne /nH ⇠ 10 2
in the core of the knot, showing that associative detachment is dominant. This result
will hold regardless of any changes to the abundances or (H). If we are overpredicting
the H2 emission due to uncertainties in the geometry, then the core density ncore
would need to be lowered to compensate, producing a larger ne /nH and still favoring
associative detachment.

2.4.2 Timescale Considerations
We have assumed that the age of the SNR is long enough for physical processes
occurring in K51 to be in steady state. Here we discuss processes for which this is not
true: (1) the known variability in the incident continuum spectrum due to high-energy
gamma ray bursts; (2) the H2 formation timescale; (3) the dynamical timescale; (4)
the evaporation timescale.
The gamma-ray activity of the Crab is the quintessential example for high-energy

69

plasma physics, and as such, has been the subject of several recent investigations (e.g.
Abdo et al. 2011; Tavani et al 2011). In order to consider the most extreme case
in which a gamma ray flare would alter the physical conditions in which H2 forms,
we have chosen an exceptionally energetic flare to compare with the total luminosity
produced by the pulsar. Buehler et al. (2012) detected a gamma-ray burst that
corresponds to a peak luminosity of ⇠ 4 ⇥ 1036 erg s 1 . However, that gamma-ray
peak only corresponds to ⇠5% of the total luminosity currently being injected into the
nebula, after the luminosity responsible for acceleration of the filaments is removed.
Such slight variations are negligible when considering the e↵ects on the molecular
knots.
Of greater concern is the H2 formation timescale. In K51, this is set by the
associative detachment reaction. The rate for this reaction is limited by the rate for
radiative attachment, H0 +e ! H + . Therefore, from Equation 1 the rate is given
0.663 is the rate coe cient for radiative
by rH = ne ↵H where ↵H = 8.861 ⇥ 10 18 Te

attachment at the temperature in core of the Crab (Ferland et al. 1994). The electron
(or kinetic) temperature in the core is very close to the H2 temperature of 2800 K. The
smallest electron density in the core is ⇠ 7 ⇥ 102 cm 3 in the temperature floor case,
and this will set the most extreme timescale. In this situation, ⌧H = 1/rH ⇠ 2.5 ⇥ 104
yr, which is much longer than the age of the Crab Nebula. This indicates that future
calculations must include time dependent chemistry. We hope to explore the e↵ects
of a time-dependent calculation in a future paper.
Our extra-heating models are highly over-pressured in the core relative to the
envelope of ionized gas, raising the question of how long such knots could survive.
The sound crossing time is set by the sound speed in the core, cs = (kT /P0 mH )1/2 ⇠
7.0 ⇥ 105 cm s 1 , and is ts ⇠ 1400 yr, roughly the age of the Crab. However, the
evaporation (or mass-loss) timescale is more physically relevant. Assuming that the
heated gas in the ionized layer flows outward at the speed of sound in the ionized gas,

70

2
cs ion ⇠ 1.4 ⇥ 106 cm s 1 , the mass loss rate is dm/dt ⇠ 4⇡lknot cs ion ⇢ion where

⇢ion is the total mass density in the ionized skin. Using a total K51 mass of 2 ⇥ 10 3
M

(as found below in §4.6) gives tpe ⇠ m/(dm/dt) ⇠ 1800 yr. Therefore, even if

K51 formed instantaneously after the supernova explosion, it is reasonable that it has
marginally survived.
A related possibility is a dynamical environment where dense seeds are ablated to
produce the emission in a non-equilibrium flow. This is highly similar to the knots
in the Helix planetary nebula (Henney et al. 2007).

2.4.3 E↵ects of Geometrical Uncertainty
We assumed a plane-parallel geometry for our modeling of K51, and compared our
results to emergent line intensities deduced as if we were in fact measuring a plane
parallel slab face-on. In reality, although K51 is the simplest case available, it has
a complicated and poorly defined morphology. Our modeling does not intend to
accurately describe all of the features present in K51. Instead, we o↵er a proof of
concept model in which our main results hold regardless of the actual geometry. As
an example of the uncertainties due to the poorly-known geometry, we consider the
e↵ects if the actual geometry were spherical as opposed to the assumed plane-parallel
situation.
As described in §2.2, we determined an observed surface brightness for each emission line by taking the total flux measured from the knot and dividing by the knot’s
projected surface area on the sky. But a sphere of a given radius has two times more
surface area than the front and back surfaces of a thin disk of the same radius, so the
emergent energy per unit area would be two times smaller.
In our plane-parallel models, a simulation stopping halfway through the knot has
an H0 /H2 core that extends over a range of depth that is one-third the depth of
the ionized zone, producing a core:ionized volume ratio of 1:3. But if this same

71

dependence on depth into the cloud were true for a spherical cloud, the core:ionized
volume ratio would become about 0.02. Thus, if our Cloudy results were adjusted to
account for a spherical geometry, the predicted ratio of H2 :ionized surface brightness
should be decreased by approximately a factor of 16.
An additional uncertainty comes from the ⇠1 mag of internal extinction at visible
wavelengths. If the knot is opaque to optical lines then emission lines coming from
the rear half of the knot (which we assume to be illuminated from the rear side) would
need a greater reddening correction. However, the knot is transparent to IR lines,
thus, only half of the volume emitting optical lines from the ionized outer skin would
be observed, while all of the H2 -emitting volume would be seen at NIR wavelengths.
This would improve the fit for the ionizing particles case, since the overpredicted
H0 and He0 visible wavelength lines originate mostly from the core, and on average
would experience greater extinction than has been accounted for and hence would
become weaker. Although in that case the Br /H ratio would no longer fit the Case
B values, the measurement of the Br line strength has considerable uncertainty and
is not a tight constraint on our models.
These factors may or may not partially cancel each other out, but taking them
into account would change our optimization of the abundances, of (H), and of the
density law. However, we stress that such a reoptimization is unlikely to change our
basic results about H2 formation and the state of hydrogen in the core.

2.4.4 Are Models with Fully Molecular Cores Ruled Out?
In Papers I and III we presented estimates of the H2 emissivity per H baryon based
on the assumption that the emitting cores are fully molecular and are excited solely
by collisions. We suggested that our H2 survey selectively found molecular regions
in which the temperature and density are in the range where H2 has its highest
emissivity. To create a fully molecular core, the internal extinction must be large to

72

shield UV radiation from inducing dissociation (Tielens & Hollenbach 1985). As is
discussed in Paper III, given the observed 2800 K H2 temperature, the measured H2
emission would need to come from a much smaller volume than is contained within
the ⇠1” diameter structures seen in our H2 images, suggesting that the emission
might come from either a thin shell around a completely obscured inner core or from
dense sub-condensations with a very small volume filling factor.
Here, we have arrived at a very di↵erent picture. At least in the case of K51,
the Crab SED is insu ciently shielded to neglect photodissociation. The two models
which include additional excitation mechanisms produce an extended H0 zone in
which hydrogen is almost entirely atomic, not molecular, but has a su ciently large
column density so that the observed H2 emission still can be produced within the
observed volume of the knot. The measured H2 emission comes from an atomic
region that fills the volume of the knot, rather than from the fully molecular zone
posited in Paper III.
However, the statistics suggest that other H2 knots may well have fully molecular
regions. In Paper II we found that the H2 knots are predominately redshifted (34
have v > 100 km s 1 while only 7 have v <

100 km s 1 ). This means that the knots

we have detected in the H2 line are on the back side of the expanding nebula, which
has a systemic velocity near 0 km s 1 . Our translucent model does not account for
this observation. A possible explanation is that in most H2 knots the excitation is
stronger on the side facing the pulsar and that a fully molecular zone with significant
2.12 µm extinction has formed on the side away from the pulsar, which for the missing
blueshifted knots would be the side nearer to us.
In addition, some of the other knots found in our H2 survey could turn out to
be fully molecular gas that is optimally emitting in H2 as suggested in Paper III.
The extended very faint H2 emission that seems to follow many of the large-scale
filament structures (Paper II) could represent fully molecular regions with less than

73

the maximum H2 emissivity. There could also be additional discrete condensations
of cooler molecular gas that do not produce detectable H2 emission.
An important check on the possible existence of fully molecular gas will come from
future mm-wavelength observations. Figure 15 compares the very di↵erent appearance of the mid-IR spectra expected from an extended H0 zone of the type we propose
for K51, and from a fully molecular core. For both cases we show the predicted spectrum from an H2 knot, including the line emission, bound-free emission, free-free
emission and dust continuum emission produced within the knot, together with the
transmitted part of the incident synchrotron continuum. The mm-wavelength range
in the fully molecular model (lower panel) has a spectrum rich with lines from CO,
CS, HCN, and a host of other molecules, which are mostly absent in the temperature
floor model (upper panel). As an example of the di↵erence between the models, for
the temperature floor and ionizing particle models the predicted strength of the CO
j=1

0 line at 2.6 mm is ⇠ 107 times fainter than H2 2.12 µm. This corresponds to

a flux of about 0.6 mJy km s 1 . For a line width of 20 km s 1 the peak brightness
would be 0.05-0.24 mJy, which is very faint but probably just detectable with Atacama Large Millimeter Array (ALMA) if the line is narrow and all of the flux falls
in one beam. In the fully molecular core case, the CO emission would instead come
from a region where H2 is not emissive because it is much cooler than the 2800 K
that we found for K51, so for the same amount of H2 emission the CO lines would be
extremely bright.
In Paper III we suggested a model in which the H2 emission comes from a thin,
heated outer layer on a cool core of molecular gas. In that sort of picture, the column
density through K51 of warmer (emitting) H2 is 2 ⇥ 1018 cm 2 , and (based on the
ratio of the knot diameter to the thickness of the H2 -emitting layer) the total column
density including the cold (non-emitting) H2 is roughly 2 ⇥ 1021 cm 2 . Using a
conversion factor X = 3 ⇥ 1020 H2 cm 2 (K km s 1 ) 1 (Young & Scoville 1991),
74

the surface brightness in the CO line should be about 7 K km s 1 . This is 104
times brighter than the previous case, so CO would be relatively easy to detect with
current interferometers. However, this latter result assumes that the CO-emitting gas
fills the entire volume corresponding to the observed 2” diameter of K51; a smaller
filling factor (as is likely to be the case) would make the line much harder to detect.
Also, the X-factor value is based on measurements of optically thick gas with typical
ISM metallicities, so it is unclear how accurate it is for the case of the Crab. In spite
of these two caveats, a deep CO measurement is likely to provide a critical test of our
model.

2.4.5 The Role of Dust
Thermal emission from dust is clearly seen in the integrated continuum spectrum from
the full Crab Nebula (Marsden et al. 1984). Temim et al. (2006, 2012) estimated
the total dust mass associated with the filament system to be 10 3

10 2 M while

Gomez et al. (2012) estimated the range to be 0.26-0.68 M . Well-resolved dust
absorption features are known to be present at a number of points across the Crab
(Fesen & Blair 1990, Hester et al. 1990). Previous photoionization models of the
filaments have yielded a dust to gas ratio estimate an order of magnitude above the
ISM (Sankrit et al. 1998).
Our extra heating models correctly predict the shape of the observed thermal
emission bump from dust in the Crab. The filled circles in the upper panel of Figure
2.15 (plotted on top of the temperature floor model) show observed points derived
from Herschel, Spitzer, Wise and Planck images, taken from Table A1 of Gomez et al.
(2012). The Gomez et al. (2012) results are integrated over an aperture that covers
most of the Crab. They are corrected for the average line emission, but include a st-rong contribution from the unfiltered synchrotron emission. We adjusted their observed flux values by subtracting o↵ the synchrotron component, scaling the remaining

75

Figure 2.15 Predicted IR-mm wavelength spectra for the Temperature Floor model
and for a typical fully molecular core, computed with Cloudy. Examples of the
stronger atomic lines are marked in the upper panel, while the CO ⌫ = 0 lines and
the region dominated by H2 emission lines are marked in the lower panel. The region
longward of 100 µm in the molecular core model includes strong lines of CO, CS, HCN,
SiO, NH3 , H2 O, SO and NO. The red dots show the continuum emission summed
over a wide area on the Crab as measured by Gomez et al. (2012), after adjusting
their observations to match the underlying continuum shape not due to dust emission
and also the height of the dust emission feature predicted for a single H2 knot. This
shows that our model correctly predicts the peak wavelength and width of the dust
emission feature.

76

dust emission to match the height of our predicted dust emission, and finally adding
back on the underlying continuum component predicted by our models. This means
that the observed points at

< 10 µm and

> 500 µm have been forced to match our

simulated spectrum, but that the points between those two wavelength limits show
that the shape and peak wavelength of our predicted dust emission closely match the
observations. This tests the dust temperature that is predicted by our models, but
not the dust abundance. Our models predict that as a result of the assumed grain-size
distribution and compositions, the dust will have a range of temperatures lying between 38 and 54 K. This is in reasonable agreement with the cold and hot components
found by Gomez et al (2012), which had Tdust ⇠28 and 63 K, respectively.
While we do not have observations that allow us to directly constrain the dust
abundance in K51, we explored the e↵ects of varying this quantity in §3.6. The best
model from this analysis arises from adopting a dust abundance five times greater
than the ISM. This creates an environment in which H2 forms at equal rates via grain
catalysis and associative detachment, emphasizing the importance of less common H2
formation routes even in dustier environments.
As was noted in §2.3, we assumed dust with optical properties that produce a
rather low albedo. If we were to repeat the entire analysis using dust properties at
the other end of the albedo range, we would correct the observed absorption dip to
a significantly larger AV pnt value but then the models would also predict a larger
AV pnt . The two AV pnt values should match about as well as they do for the
dust properties assumed here, but there would be changes in the dusts e↵ect on the
internal structure of the knot and the emission line spectrum which we would expect
to be fairly modest.
Are there dust-free knots in the Crab? Although K51 clearly contains dust,
archival HST continuum images show that there are other H2 knots that do not
appear to be associated with specific dust absorption features. Therefore, it is possi-

77

ble that the observed dust emission originates only in some parts of the filamentary
system, while other parts have no dust. The dust-free model briefly described in §3.6
represents a first pass at describing dense condensations in such regions. It is able to
produce the observed 2.12 µm line in a unique environment in which grain catalysis,
typically a primary mechanism for H2 formation in the ISM, is absent and associative detachment is the only major H2 formation mechanism. Our dust-free model
is roughly as successful as the temperature floor model in describing the observed
emission line intensities from K51.

2.4.6 Mass Estimate
Previous sections have developed a type of model that largely reproduces the observed emission from the ionized and molecular regions of K51. Here we use this
model to estimate the total mass in the knot. Previous studies have used hydrogen
recombination lines to estimate the mass of ionized gas, and H2 lines to measure the
mass in molecular regions. The latter is especially uncertain because we do not have
any direct measurement of the density of the molecular core, and the derived mass
depends on the density when the density is low.
We use our model to convert line intensities into a mass. In this approach, the
model self consistently accounts for regions where hydrogen is ionized, molecular, and
atomic. The latter can only be estimated with a model because there are no existing
21 cm detections of the knots3 . The model reports the total column density required
to produce the surface brightness in H I or H2 lines, which we then rescale to subtract
o↵ the H+ zone. A surface brightness can then be converted directly into a column
density and, if the size of the knot is known, into the total mass, for the atomic and
molecular regions.
3 The

predicted 21 cm emission from Knot 51 would only be ⇠ 1017 cm 2 K 1 ,
which is much too faint to be detectable with an instrument like the Expanded Very
Large Array (EVLA).
78

We can use the temperature floor and/or ionizing particle model to estimate the
total mass that is needed to produce the observed H2 emission. To within 5%, the
two models give the same result. Only a modest amount of H2 is required, because
the density and kinetic temperature in the extended H0 / H2 core are in the range
over which H2 will emit with nearly maximum e ciency (see Figure 5 of Paper III).
However, both models predict that there will be a factor of 2000 times more H atoms
present in the form of H0 than in the form of H2 , which is a significant amount of
mass that would not otherwise be detectable.
For K51 alone, using the ionizing particle model, the mass of H2 is mH2 = 8⇥10 7
M

while the total mass in the core including H0 and the He and heavy elements

in the H0 / H2 zone is mcore = 2 ⇥ 10 3 M . For the sum of all 55 H2 knots
reported in Paper II, scaling by the observed H2 flux gives mH2,total = 4 ⇥ 10 5
M

and mcore,total = 0.1 M . This assumes that the surface brightness to mass

conversion factor derived for K51 applies to all knots. The latter value is an order of
magnitude smaller than the estimate made by Fesen, Shull & Hurford (1997) based on
the assumption that fully molecular cores were needed to explain the dust absorption
features and H2 emission known at that time. Combining our results with Fesen et
al.(1997)’s estimate of 3.4 M

in the form of ionized and previously-known neutral

gas, the H0 / H2 cores contribute about 5% of the total mass in the system of filaments.
In Paper II we noted very weak H2 emission coming from many parts of the filaments outside of the catalogued H2 knots, which could trace additional similar cores
that do not e ciently emit H2 because they do not have the optimum temperature
and density for H2 emission. These would not be counted in the above mass estimate.
Figure 8 suggests that such cores would emit low surface-brightness [O I] and [S II]
lines, but we have not run specific models of them.
There could also be additional cores that in fact are fully molecular. We ruled out
such a possibility for the particular case of K51 because the H2 spectrum predicted

79

by our dense core model does not match the observations. But there could be other
knots for which we do not have a NIR spectrum that in fact would be described by
our dense core model, or there could be cores that are too cold to strongly emit H2
but which could be detected through CO observations.

2.4.7 Crucial Observational Questions
Further progress on understanding the H2 emitting knots rests in part on filling in
some missing key pieces of information. Here we describe four open questions that
can be addressed in a practical way by new observations.
Where is the dust and what is the dust-to-gas ratio? As is discussed above, many
H2 emitting knots do not show accompanying dust absorption. Are they dust free?
Dust absorption is seen in many locations where we have not been able to detect
strong H2 emission. What are the detailed correlations between the absorbing dust
features and the H2 emitting cores? We are currently using archival HST images
to investigate this question, but the available continuum images do not fully sample
HST’s point spread function and have useful signal:noise ratio over only about 2/3
of the Crab. A survey to uniform depth over the full Crab at the higher spatial
resolution o↵ered by WFC3 would be very helpful.
What is the geometry? The existing H2 images have only about 0.7” FWHM
spatial resolution, which is insu cient to discern the exact relationship between the
ionized and molecular regions. Observing the 2.12 µm line in K51 and a few other
knots with an adaptive optics system on an 8 m-class telescope could produce images
with spatial resolution similar to or better than the HST visible-wavelength narrowband images, which would be a major improvement.
Does significant H

emission come from the same zone that emits H2 ? Both

of our models with extra heating say that it does. The available HST H

image

(Fig. 2) is not nearly good enough to answer this question. A new HST H image

80

with high signal:noise ratio would directly test these models, which predict that the
H

emission should largely trace the H2 morphology, rather than coming from a

surrounding sheath. It is important that K51 has a radial velocity that is measured
to be close to zero, because this means that its emission lines fall within the velocity
range covered by the HST narrow-band filters.
Finally, are there any fully molecular cores? CO observations are needed to search
for these. Conversely, detection of strong CO emission from K51 would rule out
our extra heating models (§4.4). We are currently carrying out a program using
the Institut de Radioastronomie Millimetrique (IRAM) 30 m and Plateau de Bure
telescopes for a preliminary search for CO emission from the Crab, but ALMA will
be able to carry out a much deeper survey.

2.5 Conclusions
We have explored four di↵erent types of models of the regions in small knots within the
Crab Nebula. Our goal is to simultaneously reproduce the classical optical emission
lines and the exceptionally strong emission from H2 that has a thermal distribution
of populations. The first two of these models are powered only by the synchrotron
radiation field and can be ruled out. This includes the pressure-balance model, which
does not produce su ciently strong H2 emission from within the size scale of Knot 51,
and the dense-core model, which can produce the necessary H2 2.12 µm emission by
secondary electron excitation but predicts a highly non-thermal H2 spectrum. From
those results we deduced that extra heating is needed in the core of the knot in order
to excite the observed H2 emission.
We considered two cases for this additional heating: generic heating by any mechanism (such as shocks or dissipative MHD waves) that is able to supply a heating
rate of 2.5 ⇥ 10 16 ergs s 1 cm 3 (the temperature floor case), and as a specific
heating mechanism, ionizing particles. As was discussed above, both of these models
81

successfully reproduce the observed H2 emission including both the surface brightness
and the H2 line strengths relative to each other. Both also successfully reproduce the
forbidden lines from the more highly ionized ionized gas. The temperature floor case
underpredicts the strengths of forbidden lines from neutral and low-ionization species
while getting the recombination lines approximately right, while the ionizing particle case comes closer to correctly predicting the neutral and low-ionization forbidden
lines but overpredicts the H I and He I recombination lines. The 5⇥ dust model is
able to also fit the observed AV pnt to within a factor of two while still tting the
emission line surface brightnesses to the same accuracy that is achieved by the other
two models, but at the price of adding the dust abundance as an additional free parameter. We conclude that all three of these models t the observations to within the
error bars that are set by the uncertainties in the geometry
We have found that the dust content in our simulations can be varied appreciably,
from dust free to a factor of 5 above ISM abundances, while still maintaining overall
reasonable fits to the ionized, neutral and molecular line strengths. Thus, we can
conclude that dust might play an important role in understanding certain knots but
it is not necessary for all knots.
The goal of this paper has been a general exploration of the types of models
that might be able to explain the strong, thermal H2 emission seen coming from the
molecular knots. This paper presents a progress report, not a fully accurate model.
In that spirit, we conclude that the models with extra heating do in fact represent
progress and are likely to be a guide to what is really happening in the H2 emitting
zone. We have ruled out one of the two cases (dense cores) proposed by G90, while
the second, ionizing particles, remains a possibility.
We have found an unusual astrophysical environment, independent of the exact
H2 excitation mechanism, in which the relatively large ne /nH ratio allows H2 to form
most quickly through associative detachment rather than through grain catalysis.

82

There are only a few published descriptions of such environments (e.g. Ferland et al.
1994), and this paper presents a di↵erent look at the possible conditions in supernova
remnants.
We have assumed that steady state holds in all of our models. In §4.6 we showed
that in order to accurately understand the expected emission, future simulations
including time dependent chemistry must be performed. Time dependent plasma
/ chemistry simulations in an environment like the Crab filaments have never been
performed and this paper provides the basis for such an investigation. The existence
of H2 requires a dense core to provide shielding. The observed H2 temperature shows
that this would be highly overpressured relative to the envelope of ionized gas.
The mass contained in these extended H0 /H2 zones is significant, amounting to
about 0.1 M

for the sum of the 55 H2 -emitting knots reported in Paper II. This

is about 5% of the total mass estimated to be in the form of ionized gas or in more
conventional H0 zones. The observed H2 is producing the 2.12 µm line at essentially
its peak emissivity per unit mass, so the Crab could easily include many additional
extended H0 / H2 zones or fully molecular zones that are too cool to produce detectable
H2 lines; CO observations are needed to look for such gas.
Improved high-angular resolution, high signal-to-noise images are needed, especially in the visible continuum, in H and in the H2 2.12 µm line. These would better
trace the morphological connections between the dust, the ionized gas, and the H2 .
This would greatly improve our ability to discriminate between possible heating and
H2 excitation mechanisms. Finally, a survey of CO would provide the means to distinguish between models with dense molecular cores and our favored models, which
include extra heating.

83

3 Interpreting the Ionization
Sequence in Active Galactic
Nuclei Emission Line Spectra
We investigate the physical cause of the great range in the ionization level seen in
the spectra of narrow lined active galactic nuclei (AGN). We have used a recently
developed technique called mean field independent component analysis (MFICA) to
identify examples of individual SDSS galaxies whose spectra are not dominated by
emission due to star formation. From these objects, we assembled high S/N ratio composite spectra into large subsets of similar objects, defining an ionization sequence
that includes exceptionally low-ionization cases. This is our AGN sequence. We then
used a local optimally emitting (LOC) cloud model to fit emission-line ratios in this
sequence, including the weak lines that can be measured only in the co-added spectra. These weak line ratios provide consistency checks on the density, temperature,
and ionizing continuum of Seyfert galaxies determined from strong-line ratios. After
integrating over a wide range of clouds at di↵erent radii and densities, we find in
our models that the radial extent of the NLR is the major parameter in determining
the position of higher to moderate ionization AGN along our sequence, providing a
physical interpretation for their systematic variation. Higher ionization AGN contain
optimally emitting clouds that are more concentrated towards the central continuum
source than in lower ionization AGN. We also find that the ionizing luminosity is probably anti-correlated with the NLR ionization level, and hence anticorrelated with the
radial concentration of the NLR. This suggests that the ionization sequence might be
84

an age sequence, where the low-ionization objects are older and have systematically
cleared out their central regions by radiation pressure. The success with which we
have fitted classical excitation mechanism diagnostics demonstrates the e↵ectiveness
of the MFICA technique in isolating AGN contributions to emission lines, which future studies can utilize to further investigate the structure of the NLR and the origin
of emission.1

3.1 Introduction
The narrow line region (NLR) of active galactic nuclei (AGN) is a region of ionized and
neutral gas generally classified by strong [N II] 6584 and [O III] 5007 emission. In
contrast to the very small broad line region (BLR), the NLR is on the order of 102 pc
in size and contains relatively low-density gas, with electron densities ne ⇠ 104 cm 3 .
However, the physical picture of the NLR is unclear as there are few resolved NLRs.
An important route for advancing our understanding of this region is to compare
observed emission-line strengths to the ones predicted by models simulating an entire
NLR.
Diagrams constructed from emission-line intensity ratios are widely used for this
sort of comparison. Early empirical work by Baldwin, Phillips and Terlevich (1981)
pointed out the importance of the [O III] 5007/ H vs. [N II] 6584/ H↵ diagram
(hereafter the BPT diagram) for classifying emission-line galaxies as either active
galactic nuclei (AGN) or star forming (SF). The optical diagnostic diagrams presented
in Veilleux and Osterbrock (1987; hereafter VO87) added further diagnostic intensity
ratios that are minimally a↵ected by reddening. More recent work has extended this
method into the ultraviolet (UV) and infrared (IR) (Spinoglio & Malkan 1992; Allen,
Dopita & Tsvetanov 1998; Sturm et al. 2002; Groves, Dopita & Sutherland 2004a,
1 This

chapter is almost entirely taken word for word from Richardson et al.

(2013b)

85

2004b (hereafter G04a, G04b)). The class of low-ionization galaxies called LINERs
(Heckmann 1980), which has also been studied in many papers, is characterized by
very strong emission lines from neutral species and includes many objects which are
well separated from both AGN and SF galaxies in the BPT and similar diagrams.
However, there are other low-ionization galaxies called “transition” objects which
merge down into the SF galaxies on such diagrams (Ho, Filippenko & Sargent 1993;
Kau↵mann et al 2003; Kewley et al 2006). Here, we study the relationship between
this latter class of low-ionization AGN and typical higher-ionization Type II AGN.
There have been many studies using the observed intensity ratios to identify the
ionization source in narrow line Seyferts and SF galaxies. Photoionization models
have been remarkably successful in reproducing the observed emission-line ratios in
high-ionization AGN (Ferguson et al. 1997 (hereafter F97); Komossa & Schulz 1997;
G04b). The excitation mechanism for low-ionization narrow line regions in the region
overlapping with SF galaxies remains uncertain and it is not even clear whether all
such objects contain AGN. Photoionization models of AGN have focused on matching
the emission from either high-excitation and low-excitation galaxies but none have
attempted to simultaneously match both. Models that incorporate shock excitation
(Dopita & Sutherland 1995, 1996) can account for LINERS but fail to fit high ionization AGN.
The structure of the NLR is another open question. Three di↵erent models that
can reproduce the observed emission from the NLR of Seyferts are: (1) a combination
of matter and ionization bounded clouds (Binette et al. 1996); (2) radiation pressure
dominated dusty clouds (G04a); and (3) locally optimally-emitting clouds (F97; see
also Komossa & Schulz 1997).
The model that combines matter- and ionization-bounded clouds is specified by
the ratio of the solid angles subtended by matter-bounded and ionization-bounded
clouds.

The matter-bounded component is responsible for high ionization lines,

86

namely He II 4686, and the ionization-bounded component is responsible for low
to intermediate ionization lines. The dusty, radiation-pressure-dominated model is
based on the principle that radiation pressure incident on grains is su cient to moderate the density, excitation and surface brightness of high ionization regions. Finally,
the locally optimally emitting cloud (LOC) model assumes that the emission lines
come from a vast sea of clouds distributed over a wide range of densities and radial
distances from a central source of ionizing radiation and that the observed emissionline intensity ratios are the result of a powerful selection e↵ect: emission lines are
produced in clouds that optimally emit them.
Each of these models can successfully represent the NLR in at least some galaxies,
but it is not clear whether any of them can explain all NLRs. The low-ionization AGN,
which we include in this present study, fall in a region of the BPT diagram that is
also occupied by low ionization SF galaxies (cf. Tanaka 2012a, 2012b). This overlap
with SF cases greatly confuses the issue of whether or not one general type of model
can explain the full sequence from low-ionization non-SF objects to the much more
luminous high-ionization objects such as classical Seyfert galaxies. The low-ionization
AGN may or may not include cases where the ionization is due to a composite nonthermal and SF spectrum. What is needed is an unbiased set of AGN that span the
full ionization range.
In Allen et al. (2013; hereafter Paper I) we developed a new technique that is able
to separate AGN from SF galaxies with much greater purity than is usually achieved
with older methods such as those based on BPT diagrams or principal component
analysis. This mean field independent component analysis (MFICA) technique allowed us to separate galaxies along a star-forming (SF) locus from those along an
active galaxy (AGN) locus extending down to even quite low ionization levels. We
used a sample of about 104 low redshift (0.1 < z < 0.12) SDSS emission-line galaxies
to isolate sequences, or “loci,” of pure SF and of AGN cases.

87

First, the MFICA technique was used on small subsets of the full sample to generate five “continuum” and five “emission-line” components. We found that three of
the emission-line components could be combined with various weightings to reproduce the spectra of any SF galaxy, while all five emission-line components were also
needed in order to reproduce the AGN spectra. The relative strengths of the di↵erent
emission-line components therefore specified the relative SF and AGN contributions
to each galaxy’s spectrum.
Next, the components were fitted to each member of the full sample of ⇠ 104
emission-line galaxies in a two-step process, first fitting and subtracting the continuum
and then fitting the emission-line components. This produced a set of five emissionline component strengths for each galaxy in the full sample, which defined a fivedimensional parameter space. We then used an existing technique (Newberg & Yanni
1997) to identify two loci of galaxy points running through this parameter space, one
corresponding to pure SF galaxies and the other to AGN.
Figure 3.1a shows the BPT diagram for the full sample of galaxies used in Paper
I. The dashed curve is the boundary between SF and AGN defined by Kau↵mann et
al. (2003), while the solid line is the theoretical upper limit on SF line ratios found
by Kewley et al. (2001). The diagonal line is the boundary between Seyferts and
LINERs suggested by Kau↵mann et al. (2003), although this definition was later
refined (Kewley et al. 2006). Objects between the dashed and dotted lines may be
a mix of SF, AGN and composite cases, introducing considerable confusion about
the ionization mechanism over an extensive part of the full diagram. In contrast,
Figure 3.1b displays a diagram that can be interpreted as a more abstract version
of the BPT diagram. The three axes represent the contribution of the first three
MFICA emission-line components (not including the AGN-only components) to each
individual galaxy’s spectrum, and shows how the SF and AGN loci from Paper I
project into this particular parameter space. There is considerable separation between

88

them except at the very lowest-ionization end of the AGN locus. This degree of
separation shows that AGN galaxies fall along a single locus implying that their
variation can be represented by a single free parameter.
The triangles on Fig 3.1a (discussed in more detail below) represent the spectra
that we will analyze below, and show how the AGN sequence from Paper I maps
back onto the BPT diagram. On this and other line-ratio diagrams, the AGN locus
does not pass through the regions on line-ratio diagrams usually considered to be the
realm of LINERs. Rather it starts out below the LINER region in the area occupied
by transition objects, and then angles up into the AGN region.
The MFICA technique has picked out a sequence of AGN and has followed this
sequence down into the transition object region, but in a way that rejects SF galaxies
with similar emission-line spectra. Our hope is that we have separated out the “wheat
from the cha↵” in the transition zone. We can use this fact to gain a better understanding of how NLR properties change in response to changes in more basic AGN
properties, over a wider range of parameter space than has previously been studied.
Since the AGN sequence is mainly a monotonic progression along the AGN locus, a
tantalizing possibility is that we may be able to find a single physical parameter that
can be varied to determine the position of an individual object along sequence. As a
check on this, we will also study objects that represent the observed scatter orthogonal to the AGN sequence. In §2 we describe composite observed spectra formed at
five roughly equally spaced points along the MFICA AGN locus, and at two additional points to either side of the AGN sequence at each of those positions on the
sequence. These composites range from very low ionization objects up through very
highly ionized Type II AGN, and have su ciently high signal:noise (S/N) ratios that
many additional observed line ratios involving weak lines can be measured as consistency checks in addition to using the usual strong line ratios which have been used
to constrain previous models. None of these composite spectra show any evidence

89

of a non-thermal continuum, except for their emission lines. Then in §3, following
F97 we adopt an LOC model and use it to reproduce the observational line-ratio diagrams. We investigate the sensitivity of the line ratios to various physical parameters.
The physical interpretation of our results is discussed in §4, while §5 summarizes our
conclusions. We will address the SF locus in a future paper.

3.2 A Comparison Sample at Representative
Points Along the Locus
The clean separation of AGN from SF galaxies enabled by the technique developed in
Paper I enables a deeper exploration of the NLR than has been possible in previous
work. Our approach here is to use the MFICA technique to identify subsets of galaxies
which are essentially AGN cases at five di↵erent, equally spaced points along central
spine of the AGN locus, plus additional point to either side of the AGN locus. The
continuum was fitted to and subtracted from each individual galaxy spectrum using
the MFICA components, and then the continuum-subtracted AGN spectra within
each subset were co-added in their rest frames. The subsets for these AGN are named
aij, where the first index indicates the position along the AGN locus ranging from
i = 0 at the low-ionization end to i = 4 at the high-ionization end. The second index
is the position in a direction orthogonal to the AGN locus, with j = 0 corresponding
to the “wing” closest to the SF sequence and j = 2 corresponding to “wing” further
from the SF sequence. Figure 3.1 shows where these points (shown as triangles) fall
on the conventional BPT diagram and also along the particular 3D projection of the
AGN locus used in Figure 3.1b. Note that on the BPT diagram the lowest ionization
subsets (a0j) fall very close to the line representing the Kau↵mann et al. (2003) lower
boundary for finding AGN and are well below the line representing the Kewley et al.

90

Figure 3.1 (Upper panel) The BPT diagram for our sample of galaxies from Paper
I with classification curves Kau↵mann et al. (2002; solid line), Kewley et al. (2001;
dashed line) and Kau↵mann et al. (2003; dotted line). The red triangles show our
sequences of AGN picked out by MFICA with a4j representing the highest ionization observations. (Lower panel) The 3-D MFICA classification diagram with axes
representing the strengths of the first three MFICA components. The ovals along
each locus represent the range in galaxies still identified as AGN or pure SF, and
the sequences described in the text as “wings” (ai1 and ai2) would follow the outer
boundaries of these ovals for the AGN locus.
91

Table 3.1: Properties of aij subsets
Subset
E(B

a00 a01 a02 a10 a11 a12 a20 a21 a22 a30 a31 a32 a40 a41 a42
V)

0.80 0.83 0.68 0.43 0.49 0.56 0.40 0.43 0.45 0.30 0.33 0.39 0.20 0.20 0.21

Observed L(H )1

1.8

2.0

2.2

2.6

3.1

2.9

3.2

3.1

2.4

2.8

2.9

2.9

2.8

Observed L(5007)1

0.95 0.85 0.92 3.1

3.6

3.7

9.2

9.6

9.1

16.

12.

13.

24.

25.

24.

Dereddened L(H )1

29.

38.

25.

9.0

12.

18.

12.

13.

15.

8.8

7.7

11.

5.9

5.9

5.7

Dereddened L(5007)1 14.

14.

8.8

13.

18.

24.

35.

40.

41.

43.

35.

47

47.

49.

48.

2.1

2.4

Component weights
Component 1

0.03 0.02 0.01 0.14 0.13 0.11 0.25 0.25 0.25 0.40 0.40 0.40 0.54 0.57 0.59

Component 2

0.02 0.01 0.01 0.01 0.01 0.02 0.02 0.01 0.01 0.00 0.01 0.01 0.01 0.00 0.00

Component 3

0.73 0.74 0.75 0.52 0.54 0.55 0.34 0.34 0.34 0.13 0.16 0.18 0.01 0.01 0.01

Component 4

0.09 0.17 0.24 0.16 0.34 0.32 0.17 0.26 0.35 0.22 0.30 0.38 0.22 0.26 0.29

Component 5

0.13 0.07 0.00 0.17 0.09 0.00 0.22 0.14 0.05 0.24 0.14 0.04 0.23 0.16 0.11

1

Luminosities are in units 1040 erg s 1 .

92

(2001) upper limit for SF galaxies. Table 3.1 lists some general properties of these 15
aij subsets, each of which contains a sample of the 50 galaxies whose MFICA weights
most closely matched the selected points along the AGN locus. The E(B
were determined from the H↵/H

V ) values

intensity ratio as described below. We then list

both the observed and dereddened luminosities of H and [O III] 5007. Finally, the
table lists the relative weightings of each of the MFICA emission-line components
when fitted to the coadded spectrum representing each individual subset along the
AGN locus. Selected regions of the co-added spectra of the 5 central (j = 1) subsets
are shown in Figure 3.2. The additional subsets (j = 0, 2) are not shown in Figure
3.2 due to their close similarity to the central subset.
We then measured emission-line strengths from the co-added spectrum for each
subset. For most lines, we just integrated the flux within the line profile above a
locally fitted continuum. The accuracy was about +/-5 per cent for the stronger
lines ranging to +/-20 per cent for the weakest lines based on the S/N ratio in the
adjacent continuum. We separated the [S II] doublet by simply dividing the blended
profile at the lowest point between the two lines. He I

5876 is on the wing of

Na D, which is the one absorption feature in the underlying galaxy spectrum that
obviously did not subtract o↵ properly along with the rest of the continuum. Our
He I 5876 measurements are based on fitting the He I line with the profile of the
H

emission line, with a typical uncertainty of about 10 per cent. Table 2a lists

the observed (reddened) value, while Table 2b lists the dereddened values assuming
a standard Galactic reddening curve with E(B

V ) chosen to produce dereddened

I(H↵)/I(H ) = 2.86, appropriate for Case B recombination at ne = 102 cm 3 and
Te = 104 K (Osterbrock & Ferland 2006; hereafter AGN3). The dereddened values
are the ones used in the following analysis. The S/N obtained from co-adding spectra
allows weak lines such as [S II] 4070 (the sum of the [S II]

4068,4076 doublet),

[Ar IV] 4711, [N I] 5200, [N II] 5755, He I 6678 and [O II] 7325 to be measured

93

Figure 3.2 Observed, coadded spectra for the five subsets that fall directly along
central the AGN locus. Panel (a) shows the full spectra. The other panels show
enlargements that slightly overlap in wavelength. All flux values are shown in units
of the peak H intensity in the particular spectrum.
94

Table 3.2a: The measured emission line strengths for the AGN Locus. Measurements
are relative to H .
Ion
[O II]

air

a00

a01

a02

a10

a11

a12

a20 a21 a22 a30 a31 a32 a40 a41 a42

3727

1.47

1.43

1.18

1.48

1.30

1.31

1.78 1.58 1.48 2.05 1.82 1.71 2.48 2.37 2.35

HI

3798 <0.03 <0.04 <0.04

0.01

0.03

0.04

0.01 0.03 0.05 0.03 0.01 0.01 0.05 0.04 0.04

HI

3835 <0.01 <0.02 <0.05

0.02

0.02

0.03

0.01 0.03 0.04 0.03 0.01 0.02 0.08 0.07 0.08

[Ne III] 3869

0.10

0.02

0.08

0.16

0.17

0.18

0.41 0.42 0.38 0.57 0.54 0.56 0.78 0.79 0.81

HI

3889

0.04

0.03

0.08

0.11

0.09

0.11

0.14 0.14 0.14 0.13 0.13 0.14 0.16 0.16 0.16

[S II]

4070

0.12

0.06

0.04

0.07

0.03

0.01

0.10 0.13 0.11 0.11 0.10 0.08 0.12 0.11 0.11

HI

4102

0.17

0.14

0.18

0.18

0.17

0.20

0.21 0.20 0.20 0.19 0.20 0.19 0.22 0.21 0.21

[Fe V]? 4229

0.04

0.10

0.09

0.05

0.02

0.04

0.04 0.04 0.04 0.04 0.05 0.06 0.03 0.01 0.01

0.33

0.32

0.35

0.42

0.40

0.44

0.42 0.44 0.44 0.38 0.39 0.41 0.42 0.42 0.43

[O III] 4363 <0.07 <0.04 <0.04

0.05

0.04

0.08

0.07 0.07 0.07 0.10 0.11 0.13 0.11 0.12 0.12

He II

4686 <0.05 <0.01 <0.03

0.05

0.04

0.03

0.09 0.09 0.09 0.14 0.15 0.14 0.23 0.23 0.24

[Ar IV] 4711 <0.02 <0.02 <0.04

0.03

HI

HI

4340

<0.02 <0.03 0.02 0.02 0.01 0.04 0.04 0.05 0.06 0.06 0.06

4861

1.00

1.00

1.00

1.00

1.00

1.00

1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00

[O III] 4959

0.17

0.18

0.14

0.60

0.60

0.50

1.31 1.23 1.11 1.85 1.81 1.71 2.76 2.86 2.92

[O III] 5007

0.53

0.40

0.39

1.52

1.60

1.39

2.94 3.25 2.88 5.11 4.76 4.62 8.34 8.64 8.76

95

Figure)3.2a)(cont1d)
Ion

air

a00

a01

a02

a10

a11

a12

a20 a21 a22 a30 a31 a32 a40 a41 a42

[N I]

5200

0.05

0.07

0.08

0.04

0.08

0.08

0.06 0.08 0.07 0.06 0.07 0.07 0.07 0.07 0.07

[N II]

5755 <0.07 <0.07 <0.04 <0.06 <0.07 <0.06 0.03 0.05 0.04 0.03 0.03 0.04 0.02 0.02 0.02

He I

5876

0.14

0.08

0.07

0.11

0.09

0.09

0.12 0.11 0.11 0.13 0.12 0.11 0.12 0.12 0.12

Fe VII 6087 <0.00 <0.01 <0.00 <0.01 <0.01 <0.01 0.02 0.03 0.03 0.04 0.05 0.07 0.08 0.08 0.08
[O I]

6300

0.24

0.25

0.16

0.23

0.23

0.22

0.33 0.35 0.31 0.35 0.31 0.33 0.40 0.40 0.40

[O I]

6363

0.10

0.06

0.04

0.07

0.06

0.06

0.10 0.10 0.10 0.11 0.13 0.14 0.11 0.11 0.11

[N II]

6548

1.05

1.35

1.30

0.91

1.15

1.40

0.97 1.20 1.35 1.03 1.19 1.34 0.94 0.95 0.98

HI

6563

6.58

6.79

5.82

4.49

4.73

5.11

4.35 4.45 4.58 3.92 4.03 4.29 3.51 3.52 3.54

[N II]

6584

3.36

4.13

3.78

2.85

3.34

3.92

2.92 3.44 3.55 2.88 3.22 3.45 2.58 2.54 2.64

He I

6678

0.04

0.03

0.02

0.02

0.02

0.04

0.04 0.03 0.03 0.04 0.04 0.03 0.04 0.04 0.04

[S II]

6716

1.25

1.21

0.89

0.90

0.86

0.86

0.93 0.93 0.88 0.97 0.94 0.92 1.01 0.96 0.96

[S II]

6731

0.98

1.01

0.76

0.76

0.76

0.75

0.80 0.83 0.78 0.85 0.84 0.82 0.85 0.81 0.81

[Ar III] 7135

0.08

0.10

0.06

0.09

0.09

0.09

0.14 0.14 0.14 0.19 0.18 0.18 0.25 0.25 0.25

0.13

0.12

0.07

0.08

0.13

0.17

0.16 0.19 0.13 0.16 0.15 0.16 0.15 0.16 0.16

[O II]

7325

Table 3.2b: The dereddened emission line strengths for the AGN locus. Measurements
are relative to H .
Ion
[O II]

air

a00

a01

a02

a10

a11

a12

a20 a21 a22 a30 a31 a32 a40 a41 a42

3727

3.60

3.62

2.54

2.40

2.24

2.45

2.79 2.54 2.46 2.88 2.63 2.64 3.10 2.97 2.96

HI

3798 <0.07 <0.09 <0.07

0.01

0.05

0.07

0.02 0.05 0.08 0.05 0.02 0.02 0.06 0.06 0.06

HI

3835 <0.03 <0.04 <0.10

0.03

0.04

0.06

0.01 0.05 0.06 0.05 0.02 0.04 0.09 0.09 0.10

[Ne III] 3869

0.22

0.04

0.17

0.25

0.28

0.31

0.62 0.64 0.60 0.77 0.76 0.83 0.96 0.97 1.00

HI

3889

0.08

0.07

0.15

0.17

0.15

0.19

0.22 0.22 0.22 0.17 0.18 0.21 0.20 0.19 0.20

[S II]

4070

0.23

0.12

0.07

0.10

0.05

0.02

0.14 0.18 0.16 0.14 0.13 0.11 0.14 0.13 0.13

HI

4102

0.32

0.27

0.31

0.25

0.25

0.31

0.29 0.28 0.29 0.24 0.25 0.26 0.25 0.25 0.25

[Fe V]? 4229

0.06

0.18

0.15

0.06

0.03

0.06

0.05 0.05 0.06 0.04 0.07 0.07 0.03 0.02 0.01

0.50

0.51

0.51

0.53

0.52

0.59

0.53 0.56 0.56 0.45 0.46 0.51 0.47 0.47 0.48

[O III] 4363 <0.10 <0.06 <0.06

0.06

0.05

0.10

0.09 0.09 0.09 0.12 0.13 0.16 0.12 0.13 0.13

He II

4686 <0.06 <0.01 <0.03

0.05

0.04

0.03

0.10 0.10 0.10 0.15 0.16 0.15 0.24 0.24 0.24

[Ar IV] 4711 <0.02 <0.02 <0.05

0.04

HI

HI

4340

<0.03 <0.04 0.03 0.02 0.01 0.04 0.04 0.05 0.06 0.06 0.06

4861

1.00

1.00

1.00

1.00

1.00

1.00

1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00

[O III] 4959

0.16

0.16

0.13

0.58

0.57

0.47

1.26 1.18 1.07 1.80 1.76 1.65 2.71 2.81 2.87

[O III] 5007

0.48

0.36

0.35

1.43

1.50

1.29

2.78 3.07 2.71 4.91 4.56 4.39 8.12 8.41 8.53

Figure)3.2b)(cont1d)
Ion

air

a00

a01

a02

a10

a11

a12

a20 a21 a22 a30 a31 a32 a40 a41 a42

[N I]

5200

0.04

0.05

0.06

0.04

0.07

0.07

0.06 0.07 0.06 0.05 0.07 0.06 0.06 0.06 0.07

[N II]

5755 <0.04 <0.04 <0.02 <0.04 <0.05

0.04

0.03 0.04 0.03 0.03 0.03 0.03 0.02 0.02 0.02

He I

5876

0.06

0.09 0.08 0.08 0.10 0.10 0.09 0.11 0.10 0.10

0.08

Fe VII 6087 <0.01

0.04

0.04

0.08

0.01

<0.01

0.01

0.06

<0.01 <0.01 0.02 0.02 0.02 0.03 0.05 0.07 0.07 0.07 0.08

[O I]

6300

0.11

0.12

0.09

0.15

0.15

0.13

0.23 0.23 0.21 0.26 0.23 0.23 0.33 0.34 0.33

[O I]

6363

0.05

0.03

0.02

0.05

0.04

0.03

0.07 0.07 0.07 0.08 0.09 0.10 0.09 0.09 0.09

[N II]

6548

0.46

0.57

0.64

0.58

0.70

0.79

0.64 0.77 0.85 0.75 0.85 0.90 0.77 0.77 0.79

HI

6563

2.86

2.86

2.86

2.86

2.86

2.86

2.86 2.86 2.86 2.86 2.86 2.86 2.86 2.86 2.86

[N II]

6584

1.45

1.72

1.84

1.81

2.01

2.18

1.91 2.20 2.21 2.10 2.28 2.29 2.09 2.05 2.13

He I

6678

0.02

0.01

0.01

0.01

0.01

0.02

0.02 0.02 0.02 0.03 0.03 0.02 0.03 0.03 0.03

[S II]

6716

0.51

0.48

0.41

0.56

0.50

0.46

0.60 0.58 0.53 0.69 0.65 0.59 0.81 0.76 0.77

[S II]

6731

0.40

0.40

0.35

0.46

0.44

0.40

0.51 0.52 0.47 0.60 0.58 0.53 0.68 0.64 0.65

[Ar III] 7135

0.03

0.03

0.02

0.05

0.05

0.05

0.08 0.08 0.08 0.13 0.12 0.11 0.19 0.19 0.19

0.04

0.04

0.03

0.05

0.06

0.08

0.09 0.10 0.07 0.10 0.10 0.09 0.11 0.11 0.11

[O II]

7325

in most cases, providing important consistency checks on our models.

3.3 Modeling the Narrow Line Region
3.3.1 Method
We used the LOC model, which treats the NLR as the sum of a large number of
separate gas clouds spread out around a central ionizing source (Baldwin et al 1995).
Once the incident spectral energy distribution (SED) and the chemical abundances
are set, the major free parameters in the LOC model are the distributions of the
individual clouds in their radial distances r from the source of ionizing radiation, and
in their gas densities nH .
The individual clouds were modeled using version 10.0 of the plasma simulation
code Cloudy (Ferland et al. 1998). For simplicity, we assumed that there is no ISM
attenuation of the AGN continuum radiation incident on the clouds. The consequences of this are described in F97. We used constant density models; the results
of Pellegrini et al. (2007) show that constant-density and constant-pressure models
give very nearly the same result for optical lines, which form in the warm gas within
the H+ zone. The distinction is mainly important for infrared lines which form in
cool atomic gas and which are not considered here. We used the solar abundances
summarized in G04a, taken from a series of papers by Asplund and his collaborators (Asplund et al. 2000; Asplund 2000, 2003; Allende Prieto et al. 2001, 2002).
Abundances for additional elements not given by G04a are taken from F97. Following
F97, we used an ionizing luminosity Lion = 1043.5 erg s 1 , typical of Seyfert galaxies.
Our simulations proceeded until the electron temperature fell below 4000 K or above
105 K. Gas below 4000 K does not significantly contribute to optical emission lines,
although it does produce infrared lines, and gas above 105 K contributes principally
to X-Rays. These models did not include any cosmic rays. The Cloudy models all

99

terminate before reaching molecular regions where cosmic ray ionization from the
Galactic background becomes important. Komossa & Schulz (1997) studied the effect of adding the standard Galactic cosmic ray background to similar models, and
found that the emission line ratios changed by only 1 per cent.
A series of grids of individual clouds were computed over a range of densities and
radii that represent plausible values for gas in the NLR. Each grid consisted of 7171
cloud models. The radial distance from the ionizing source to the individual cloud,
r, was varied in 0.1 dex steps from 1016 to 1023 cm and the hydrogen density, nH ,
was varied in 0.1 dex steps from 100 to 1010 cm 3 . We assumed that the ionizing
flux was isotropic so that a di↵erence in radii reflects a di↵erence in flux proportional
to r 2 . Di↵erences in Lion would appear as a rescaling of the radial distances by
(Lion /1043.5 )1/2 . The total hydrogen density was kept constant for each individual
cloud but the molecular and electron densities were determined self-consistently and
so vary with depth. The results from these grids are shown below in terms of the
equivalent width of various emission lines, relative to the continuum level at 4860˚,
A
which indicates the e ciency of each cloud in reprocessing the continuum into a
particular emission line.
We computed both dust-free and dusty grids of clouds. By “dusty”, we mean LOC
grids that include dust in the individual clouds in parts of the NLR where it is not
expected to have been sublimated, while “dust free” refers to LOC grids with no dust
anywhere. Dust is di cult to accurately model in a physically consistent way, since
sublimation occurs at small distances from the central object, so following F97 we
adopted a step function in our dusty grids to artificially account for this process. We
used a mix of graphite and carbonaceous grains and a grain size distribution typical
of an H II region (Baldwin et al. 1991). At r  1016.9 cm the grain temperature
is high enough to sublimate both types of grains so our models were dust-free with
solar abundances. For 1016.9 < r  1017.6 cm, the silicates are too hot to exist,
100

so we included only graphite grains and used a solar composition, except for carbon
which was depleted by 15 per cent. Finally, for r

1017.6 cm, we include both

graphite and carbonaceous grains and adopt Orion Nebula abundances (Baldwin et
al. 1996) that are typical of an H II region. This is the procedure used by F97.
In reality, di↵erent sized grains sublimate at di↵erent temperatures. However, we
followed F97 and ignored this e↵ect, assuming that there is threshold radius for each
grain composition after which those particular grains are allowed to form.
The LOC model is based on the fact that di↵erent emission lines are optimally
emitted in di↵erent individual clouds according to the correct density and incident flux
for e cient production of each particular line, while we observe only the integrated
spectrum of the ensemble of clouds. This gives rise to a powerful selection e↵ect in
which the overall radial and density distributions of the clouds largely determines the
measured spectrum. As in F97 (see also Baldwin et al. 1995), we modeled the total
emitted spectrum integrated over radial distance and density distributions defined
as f (r) / r and g(n) / nH , respectively, where

and

are free parameters. An

average spectrum over many galaxies as we have produced here is more likely to be
correctly described by a broad power-law distribution than might be the case for an
individual galaxy. The total luminosity of the line is then,

Lline /

Z Z

r2 F (r, nH )r nH dnH dr

(3.1)

where F (r, nH ) is the flux of the line and r nH is the spatial distribution function
(Bottor↵ et al. 2002). Throughout the remainder of the paper, when the log of
densities and radii are given, the units are cm 3 and cm, respectively. We chose the
integration limits to include the range of parameter space that is physically relevant
for the NLR. Gas with densities log(nH ) > 8 is above the maximum critical density
for optical forbidden lines. Gas at radii log(r) < 17.48 would either be much too
hot to produce strong optical emission lines, or would have high density which would
101

suppress the forbidden lines. Gas with densities log(nH ) < 2.0 does not optimally
emit many lines, the notable exception being [S II] 6720. This is the justification
for the integration limits used by F97.
We experimented with changing the integration limits to include a larger range,
but this only produces small changes to the integrated spectrum. Changes in the
abundances, the SED,

or

all produced much more noticeable e↵ects. For this

reason, we will only present integrations within the limits 2< log(nH < 8 and 17.48
<log(r) <23.
We compared the observed integrated predictions to the predictions by using lineratio diagrams as in Baldwin, Phillips & Terlevich (1981), Villeaux & Osterbrock
(1987), and elsewhere. We discuss these diagrams below. We will exclusively describe
dust-free LOC grids until §3.6 where we address the fact that our dusty models do
not fit some of the AGN subsets as well as our dust free models.

3.3.2 Spectral Energy Distribution Optimization
The SED in AGN is uncertain. The various forms of the incident continuum assumed
in previous work include simple power laws, f⌫ / ⌫ ↵ ,(Binette et al. 1996, G04a),
a combination of power law and blackbody continua (Komossa & Schulz 1997) and
multi-component distributions (F97). In addition, the SED viewed by the BLR (Korista et al. 1997) may be di↵erent than the SED viewed by the NLR, and parts of
the SED, such as emission originating from a dusty torus are only visible for AGN at
certain viewing angles. We adopted a multi-component approach as in F97 and used
He I/H I line intensity ratios to infer the SED.
We considered the three di↵erent spectral energy distributions (SEDs) shown in
Fig. 3. They all have the general form,
f⌫ / ⌫ ↵uv e h⌫/kTcut kTIR /h⌫ + C⌫ ↵x
102

(3.2)

Figure 3.3 The SEDs for Ferguson et al. (1997), Mathews and Ferland (1987), our
best model (dust-free) and our dusty model. Our optimized models result from an
SED that has a shifted UV “Big Bump” in order to match He II 4686, a well known
indicator of the SED.
where ↵uv is the low-energy slope of the “big bump,” Tcut is the UV temperature
cuto↵, TIR is the IR temperature cuto↵, C is a normalizing constant, and ↵x is the
slope of the X-ray component. The factor C sets the scaling between the big bump
and the X-ray power law, and in Cloudy is determined by specifying ↵ox , which is
given by,
f⌫ (2keV)
= 403.3↵ox
˚)
f⌫ (2500A

(3.3)

(Korista et al. 1997). The last term of Equation 3.2 was assumed to fall o↵ as ⌫ 2
above 100 keV.

103

A commonly used AGN SED is given by Mathews & Ferland (1987; hereafter
MF87), and is shown as a dotted line in Figure 3.3. A combination of direct observations, with inferences based on the emission-line spectrum, determined the components of this SED. The X-ray slope of the SED, ↵x ⇠

0.7, is typical for radio-loud

AGN (Zamorani et al. 1981). As we describe below, this has a significant e↵ect on
the equivalent width contours but minimally a↵ects the LOC integrations. F97 used
a modified SED, shown on Figure 3.3 as the dashed line. This SED di↵ers from that
of MF87 in that it lacks infrared emission from a dusty torus, includes a soft X-ray
component, and has a softer X-ray slope ↵x ⇠

1.0 that is typical of radio quiet

quasars and AGN (Elvis et al. 1994).
The He II 4686/H intensity ratio is an important test of these possible SEDs,
because (1) He II and H

are recombination lines and therefore less sensitive to

Te than collisionally excited lines (AGN3); (2) the ratio is relatively insensitive to
large changes in the ionization parameter (Binette, Courvoisier & Robinson 1988).
Photoionization models of AGN often underpredict this ratio (Ferland & Osterbrock
1986).
Using the F97 and MF87 SEDs, we computed two separate grids of Cloudy models
as outlined in §3.1, and then ran LOC integrations. The predicted [O III] 5007/H vs
He II 4686/H line ratio diagrams for our dust-free LOC models using each of these
SEDs are shown in the upper two panels of Figure 3.4 as colored lines connecting a
series of circles. Each line represents the LOC results for a di↵erent density weighting
in the range
in

1.8 



0.6 with the di↵erent lines corresponding to increments

of -0.4. Along each individual line the radial weighting is varied over the range

2.0 

 2.0, and the circles along the line represent di↵erent LOC models run at

increments in

of 0.25. The circles with white centers represent the best-fitting LOC

models described below in § 3.5. The black diamonds, triangles and squares represe-

104

Figure 3.4 The He II 4686 /H line ratio for various SEDs: Ferguson et al. (top
left panel), Mathews and Ferland (top right panel), and our optimized SED (lower
panel). In each panel the triangles connected by the solid black line are the dereddened measurements from the co-added spectra that constitute the a40-a30-a20-a10a00 sequence. Similarly, the squares connected by the dashed black line show the
a41-a31-a21-a11-a01 sequence, and the diamonds connected by the dotted black line
show the a42-a32-a22-a12-a02 sequence. The largest shape in each sequence indicates
observations at the top of AGN locus (a40, a41, a42). The colored lines with circles
represent LOC integrations. The density-weighting indices in the LOC integrations,
, have more negative values along the direction of the arrow and are indicated by
di↵erent colors. The radial weighting indices, , become more negative as a function of increasing distance from the largest colored circle for each particular density
weighting. The hollow circles indicate our best fitting set of free parameters, with
the largest hollow circle indicating our best fit to the a41 subset (large black square).
The cross in the corner of the plot located in the lower right corner of the bottom row
indicates the range of acceptable error. The He II 4686 /H line ratio successfully
matches the higher to moderate ionization AGN subsets within a factor of two for
our optimized SED and the Mathews and Ferland SED.

105

Figure)3.4)(cont1d)

106

Figure)3.4)(cont1d)

-nt the ratios measured from the AGN observations, with each sequence in the observations indicated by a di↵erent line style and symbol. The largest of these symbols
represents the top of the AGN locus for that sequence. The largest symbol should be
compared to the largest white-filled circle, which represents our best-matching LOC
result for the high ionization end of the observed sequence.
While in Figure 3.4 and in other similar figures discussed below we have shown
all three observed AGN sequences (the central ai1 locus and the two “wings” corresponding to ai0 and ai2), all of the sequences are very similar. In fact, even in the
full set of line ratio diagrams described below the di↵erence between each sequence
is so slight that we cannot determine a physically meaningful interpretation for the
wings; they appear to just represent object-to-object scatter about the central locus.
For this reason, we compare our models only to the central locus.
107

It can be seen that the higher-ionization observed points in the upper two panels
of Fig. 4 are in good agreement with the calculated curves. However for the F97
SED, the models producing integrated spectra with moderate [O III] 5007/H ratios
(moderate ionization) overpredict the He II/H

ratios by almost a factor of three.

In contrast, the MF87 SED produces a He II/H ratio within a factor of 2. This is
closely related to the “Stoy method” of determining the SED (AGN3) and we adjust
its hardness to reproduce the ratio. Here and throughout the paper, we adopt a
factor of 2 as our criterion for an acceptable fit between the predicted and observed
line ratios, taking into account the uncertainty in the measurements of the weaker
lines, the uncertainties (and object to object scatter) in the SED and abundances,
and the fact that the LOC models just use simplified power laws to describe the
distributions of r and nH .
We optimized a third SED, shown on Figure 3.3 as the solid line, to try to produce
a better fit to the measurements from our galaxy sample. In addition to a SED similar
to radio-loud quasars (MF87), we wanted one that was typical of AGN (F97) but that
still matched the moderate ionization He II/H ratio. To achieve this, we modified the
F97 SED by shifting the “big bump” until the He II/H ratio matched the moderate
to high ionization subsets. This modified “bump” turned out to match the one in the
MF97 SED. We also adjusted ↵ox to give closer agreement with the mean value of
-1.4 found in AGN (Zamorani et al. 1981; Ste↵en et al. 2006).
Our optimized SED is characterized by Tcut ⇠ 4.2 ⇥ 105 K, ↵UV =
1.35, ↵x =

0.3, ↵ox =

1.0, and kTIR = 0.14 eV. The bottom panel of Figure 3.4 verifies the fit

to the moderate and high-ionization He II/H ratio. Further adjustments were made
to match the lower ionization He II/H ratios but we were unable to simultaneously
account for both the lower and upper end of the locus. By matching the higher and
moderate ionization He II 4686/H ratios to within a factor of two, we were able
to vast majority of our SED-indicating line-ratios ([O II] 3727/ [O III] 5007, He II

108

4686/H , [O I] 6300 /[O III] 5007, He II 4686/ He I 5876 and [N I] 5200/ [N
II] 6584) to within a factor of two without further tweaking, including many ratios
at lower ionization. The MF87 and our optimized SED give essentially identical fits
to the observed points on Figure 3.4, which are slightly better than the fit for the
F97 SED. We chose to use the optimized SED for the remainder of our analysis.

3.3.3 X-ray E↵ects
Figure 3.5 shows that the peak emissivity of some strong lines occurs at two locations on the LOC plane for the Mathews & Ferland (1987) AGN continuum, but
the Ferguson et al. (1997) continuum only creates a single maximum. We found
that this is caused by the di↵erent assumptions in the form of the X-ray continuum. The emission-line equivalent widths indicate the e ciency with which individual clouds convert continuum radiation into emission lines, and, as Figure 3.5 shows,
are generally similar for each of the SEDs shown in Figure 3.3. In high density, high
incident-flux (low-r) environments, the equivalent width contours of these lines become smeared out over a larger range of parameter space. In the case of the Mathews
and Ferland SED, these regions become isolated and form a second optimally emitting
region.
By introducing gradual changes to the MF87 SED until it resembled the F97 SED,
we were able to determine that the shape of the continuum at energies greater than
1 keV is responsible for creating the second optimally emitting region. The infraredwavelength “shelf” in the MF87 SED represents emission from the torus and does not
a↵ect the equivalent widths. The densities at which the second optimally emitting
regions occur are most likely larger than those physically relevant to the NLR, however
X-ray e↵ects are important in other contexts of AGN (e.g. in the BLR). It is important
to note that the thicknesses through some of the clouds in these optimally emitting
regions became comparable to their distance from the source. While F97 did not

109

Figure 3.5 Contours for log equivalent widths as a function of radius, r, and hydrogen
density, nH . All equivalent widths are relative to the continuum at 4860 ˚. The
A
contour increments represent steps of 0.2 dex with a lower cuto↵ at 1 ˚. The squares
A
are labeled with log equivalent width and indicate the local maxima for each grid. We
show equivalent width diagrams for three di↵erent SEDs: (a) Ferguson et al. (1997),
(b) Mathews and Ferland (1987), (c) our optimized SED. Panel (b) clearly shows a
second region in which the gas is optimally reprocessed into emission lines. This can
be seen to a lesser degree in panels (a) and (c) where a “shelf” starts to form in high
density, high flux regions.

110

Figure)3.5)(cont1d)

111

Figure)3.5)(cont1d)

allow clouds to have thickness greater than 10 per cent of their distance to the source,
they note that this criterion did not have a real impact on the resulting contours.

3.3.4 Metallicity E↵ects
We found that [N II]

6584 was underpredicted with our solar abundance set, so

we adjusted the nitrogen abundance. Nitrogen abundances have long been known to
deviate from solar abundances in AGN (e.g. Storchi-Bergmann & Pastoriza 1990; Dopita & Sutherland 1995; Hamann & Ferland 1999). Baldwin et al. (2003) found that
nitrogen abundances a factor of 15 above solar are needed to accurately represent the
QSO Q0353-383, and Bradley, Kaiser & Baan (2004) found a factor of 34 enhancement in the weak NLR in M51. Nitrogen enhancement in AGN is thought to stem

112

from secondary CNO nucleosynthesis, where carbon and oxygen are pre-existing, and
then distributed by stellar winds from massive stars. This creates a unique scaling
relationship with metallicity, [N/H] / Z 2 while other elements scale linearly with
metallicity (Hamann & Ferland 1999).
We used that scaling relation to tune the abundances to fit [N II]

6584/ [O

II] 3727, a metallicity indicator (Groves, Heckman & Kau↵mann 2006) and [N II]
6584/H↵. After iterating over a range of metallicities2 , we found that a metallicity
of 1.4 Z provides the best agreement with these two line ratios. Enhancing only nitrogen by a factor of 1.5 above solar would give an equally good fit to the observations,
but we used the 1.4 Z

value appropriate for secondary enrichment.

Our final abundance set for the dust-free models, by number relative to hydrogen,
is given in the log as follows: He: -0.987; Li: -8.54; Be: -10.4; B: -8.97; C: -3.46; N:
-3.90; O: -3.16; F: -7.37; Ne: -3.77; Na: -5.53; Mg: -4.27; Al: -5.36; Si: -4.34; P: -6.28;
S: -4.65; Cl: -6.57; Ar: -5.45; K: -6.72; Ca: -5.50; Sc: -8.65; Ti: -6.81; V: -7.83; Cr:
-6.16; Mn: -6.32; Fe: -4.39; Co: -6.93; Ni: -5.60; Cu: -7.58; Zn: -7.19.

3.3.5 Final Dust-Free Model
Our goal in this paper is find a single tunable LOC parameter responsible for the
variation in AGN ranging from Seyferts all the way down to the lowest-ionization
objects as evident in Figure 3.1. The best case for achieving that goal is to use our
final dust-free models for the individual clouds.
In order to get a fuller picture of the physical conditions in the individual clouds as
a function of their position on the log(nH )log(r) plane, Figure 3.6 shows the ionization
parameter, U (dotted lines), and the gas temperature in the [O III]-emitting zone Te
(dashed and solid lines) at the illuminated face of the cloud. The cuto↵ for detectable
2 Here

we adopted the common usage of the term metallicity as the metals abundance relative to hydrogen, by number, normalized to solar values. However, we note
that metallicity is more properly defined as a mass fraction.
113

Figure 3.6 The ionization parameter and temperature at the illuminated face of each
cloud. The ionization parameter is shown by a representative sample of evenly spaced
blue, dotted contours across the LOC plane. The temperature is given in 1.0 dex
(black lines) and 0.2 dex (red lines) increments.
lines is taken to be EW = 1 ˚, which roughly corresponds to log(U ) = 3.0. Each
A
U contour has a slope of (nH dr)/(rdnH ) =

0.5. Figure 3.7 displays the equivalent

width contours for 20 important emission lines, as a function of log(nH ) and log(r).
Next, Figs. 3.8-3.12 show 24 di↵erent diagnostic diagrams predicted for the LOC
integrations using the dust-free model, and compare them to the observations in the
same way as was done in Figure 3.4. The lines and symbols in Figs. 3.8-3.12 are
keyed the same as in Figure 3.4, with black lines and symbols showing the observed
line ratios, and the colored lines and symbols showing predicted line ratios, with a

114

separate predicted line for each sequence of LOC models in which
but

is held constant

is varied. Each panel is on the same scale to emphasize at a glance which line

ratios are better indicators. The error bars in the bottom right corner of the last
panel in each figure represent a factor of two error.
We discuss the individual diagrams from Figure 3.8-3.12 below, but a first glance
over the full set of diagrams shows that our LOC models can fit almost all line ratios
for the upper half of the observed AGN locus (the black points starting from the
larger black symbols and working back down the observed sequence). The larger
black square represents the observed extreme AGN case, a41. The LOC model which
provides the overall best fit to those observations has

=

1.4,

=

0.75 and is in-

-dicated by the larger open (filled in white) green circle. The observed a41-a01 sequence is matched using a density weighting of

=

1.4 (the green line on Figure

3.8-3.12) while changing the radial weighting as parameterized by

in Equation 3.1.

The a41-a31-a21-a11-a01 observations are fitted to reasonable accuracy in most panels on Fig. 3.8-3.12 by the LOC models with values
= 0.0,

=

0.75,

=

0.5,

=

0.25,

= 1.0 respectively, each of which is shown as an open circle. A perfect fit

would have the larger black square fall exactly on top of the larger open circle, and
then to have each successive black square fall on top of each successive open circle.
At high - moderate ionization, the fits are usually within a factor of two, which
demonstrates the e↵ectiveness of the MFICA method at picking out an AGN ionization sequence that is largely described by the variation of a single physically meaningful parameter the radial distribution of clouds. Table 3.3a shows the line ratio
predictions for our model at each of these locations as well as comparison to the observed value. Next, we discuss the details of these comparisons between the models
and the observations.

115

Figure 3.7 Contours of log equivalent width as a function of radius, r, and hydrogen
density, nH , in the same manner as Figure 3.5 but for our best fitting model.
Excitation Mechanism
Panels (a)-(d) of Figure 3.8 show classic diagnostic diagrams from BPT and VO87,
which are particularly sensitive to whether the excitation source is a power law (AGN)
or starburst. The degree to which our model fits the AGN observations emphasizes
its consistency with line ratios widely used to identify the excitation source and

116

Table 3.3a: AGN Emission line ratio predictions - Dust Free Model
Model Predictions (log of line ratio)
Line Ratio

a01

[O II] 3727 / [O III] 5007

a11

a21

a31

Log(Predicted / Observed)

a41

a01

a11

a21

0.45 -0.05 -0.22 -0.41

-0.63

-0.55

-0.23

-0.14 -0.17 -0.18

[O II] 3727 / [O I] 6300

0.70 0.73 0.74

0.75

0.74

-0.80

-0.46

-0.29 -0.31 -0.20

[O II] 3727 / [O II] 7325

2.10 1.90 1.81

1.69

1.53

0.10

0.36

0.41 0.27 0.13

[Ne III] 3869 / He II 4686

0.24 0.37 0.43

0.50

0.57

<-0.20

-0.47

-0.38 -0.19 -0.04

[S II] 4070 / [S II] 6720

-1.38 -1.27 -1.21 -1.13

-1.02

-0.53

[O III] 4363 / [O III] 5007

-2.52 -2.38 -2.34 -2.29

-2.22

>-1.74

-0.87

-0.78 -0.74 -0.41

He II 4686 / H

-0.90 -0.84 -0.81 -0.77

-0.74

>0.92

0.55

0.20 0.03 -0.11

He II 4686 / He I 5876

-0.18 -0.10 -0.05 -0.002

0.06

>0.28

0.10

-0.14 -0.21 -0.30

[Ar IV] 4711 / [Ar III] 7135 -2.15 -1.41 -1.20 -0.98

-0.76

>-2.00

-1.11

-0.59 -0.49 -0.25

[O III] 5007 / H

-0.24 0.28 0.43

0.60

0.75

0.20

0.10

-0.05 -0.06 -0.17

[O III] 5007 / [Ar III] 7135

0.69 1.16 1.30

1.46

1.64

-0.33

-0.32

-0.27 -0.12 0.00

[N I] 5200 / H↵

-1.68 -1.75 -1.80 -1.87

-2.00

0.04

-0.12

-0.19 -0.24 -0.33

[N I] 5200 / [N II] 6584

-1.43 -1.48 -1.51 -1.54

-1.59

0.07

-0.01

-0.02 -0.01 -0.07

[N I] 5200 / [O II] 7325

0.69 0.40 0.27

0.10

-0.11

0.51

0.39

0.43 0.29 0.16

[N II] 5755 / [N II] 6584

-2.38 -2.28 -2.23 -2.16

-2.05

117

a31

a41

-0.002 -0.42 -0.16 0.02

>-0.75 >-0.68 -0.50 -0.20 -0.02

Table)3.3a)(cont1d)
Model Predictions (log of line ratio)

Log(Predicted / Observed)

He I 5876 / H

-0.72 -0.74 -0.75 -0.77

-0.79

0.64

0.45

0.33 0.25 0.19

[O I] 6300 / [O III] 5007

-0.25 -0.79 -0.96 -1.16

-1.37

0.25

0.23

0.16 0.14 0.02

[O I] 6300 / H↵

-0.97 -0.99 -1.00 -1.03

-1.09

0.42

0.31

0.08 0.06 -0.16

[O I] 6300 / [N II] 6584

-0.72 -0.72 -0.71 -0.70

-0.68

0.45

0.42

0.26 0.29 0.10

H↵ / H

0.48 0.48 0.48

0.48

0.47

0.03

0.02

0.02 0.02 0.02

[N II] 6584 / [O II] 3727

0.02 -0.01 -0.03 -0.04

-0.06

0.35

0.03

0.04 0.02 0.01

[N II] 6584 / H↵

-0.25 -0.27 -0.29 -0.33

-0.41

-0.03

-0.12

-0.18 -0.23 -0.26

[S II] 6716 / [S II] 6731

0.08 0.06 0.05

0.04

0.03

-0.01

0.00

0.00 -0.01 -0.05

[S II] 6720 / H↵

-0.09 -0.15 -0.19 -0.25

-0.36

0.42

0.33

0.23 0.11 -0.05

[O II] 7325 / [O I] 6300

-1.40 -1.17 -1.07 -0.94

-0.79

-0.89

-0.82

-0.71 -0.58 -0.33

[O II] 7325 / H↵

-2.37 -2.16 -2.07 -1.98

-1.88

-0.47

-0.51

-0.62 -0.52 -0.49

-

0.50

0.46

0.62 0.81 0.85

Fraction of Ratios Fitted

-

-

-

-

Free Parameters

Free Parameters

Density Weighting -

-1.4

-1.4

-1.4

-1.4

-1.4

-1.4

-1.4

-1.4

Radial Weighting -

1.0

0.0

-0.25

-0.5

-0.75

1.0

0.0

-0.25 -0.5 -0.75

118

-1.4

-1.4

Table 3.3b: AGN Emission line ratio predictions - Dusty Model
Model Predictions (log of line ratio)
Line Ratio

a01

[O II] 3727 / [O III] 5007

a11

a21

a31

Log(Predicted / Observed)

a41

a01

a11

a21

a31

a41

0.79 -0.08 -0.26 -0.45

-0.65

-0.21

-0.26

-0.18

-0.22

-0.20

[O II] 3727 / [O I] 6300

0.77 0.78 0.76 0.72

0.64

-0.49

-0.41

-0.27

-0.34

-0.30

[O II] 3727 / [O II] 7325

1.94 1.76 1.68 1.57

1.41

-0.06

0.23

0.28

0.15

-0.002

[Ne III] 3869 / He II 4686

0.63 0.57 0.54 0.51

0.48

<0.18

-0.27

-0.27

-0.18

-0.14

[S II] 4070 / [S II] 6720

-1.24 -1.18 -1.15 -1.12

-1.06

-0.38

0.09

-0.36

-0.15

0.02

[O III] 4363 / [O III] 5007

-2.32 -1.73 -1.60 -1.46

-1.31

>-1.54

-0.22

-0.05

0.08

0.50

He II 4686 / H

-1.37 -1.07 -0.98 -0.87

-0.77

>0.45

0.32

0.03

-0.07

-0.14

He II 4686 / He I 5876

-0.53 -0.21 -0.12 0.01

0.11

>-0.07

0.02

-0.20

-0.22

-0.25

[Ar IV] 4711 / [Ar III] 7135 -2.35 -1.25 -1.09 -0.95

-0.84

>-2.11

-0.95

-0.49

-0.47

-0.33

[O III] 5007 / H

-0.18 0.57 0.69 0.81

0.90

0.26

0.39

0.21

0.15

-0.02

[O III] 5007 / [Ar III] 7135

0.47 1.19 1.32 1.43

1.55

-0.54

-0.28

-0.26

-0.14

-0.10

[N I] 5200 / H↵

-1.13 -1.30 -1.36 -1.43

-1.51

0.59

0.33

0.24

0.21

0.15

[N I] 5200 / [N II] 6584

-1.21 -1.24 -1.24 -1.23

-1.22

0.29

0.24

0.25

0.31

0.30

[N I] 5200 / [O II] 7325

0.68 0.46 0.38 0.28

0.15

0.50

0.45

0.54

0.47

0.43

[N II] 5755 / [N II] 6584

-2.11 -1.97 -1.89 -1.75

-1.55

>-0.49 >-0.37 -0.16

0.21

0.49

119

Table)3.3b)(cont1d)
Model Predictions (log of line ratio)

Log(Predicted / Observed)

He I 5876 / H

-0.84 -0.85 -0.86 -0.87

-0.88

0.52

0.34

0.23

0.15

0.10

[O I] 6300 / [O III] 5007

0.02 -0.86 -1.03 -1.17

-1.30

0.52

0.15

0.09

0.12

0.10

[O I] 6300 / H↵

-0.64 -0.78 -0.83 -0.87

-0.90

0.75

0.51

0.26

0.22

0.03

[O I] 6300 / [N II] 6584

-0.73 -0.72 -0.70 -0.67

-0.62

0.45

0.42

0.27

0.32

0.17

H↵ / H

0.48 0.49 0.49 0.50

0.51

0.02

0.03

0.04

0.04

0.05

[N II] 6584 / [O II] 3727

0.05 -0.06 -0.06 -0.05

-0.03

0.27

-0.02

0.003

0.01

0.13

[N II] 6584 / H↵

0.08 -0.06 -0.12 -0.20

-0.29

0.30

0.09

-0.01

-0.10

-0.14

[S II] 6716 / [S II] 6731

0.09 0.08 0.07 0.06

0.06

0.01

0.02

0.02

0.02

-0.02

[S II] 6720 / H↵

0.01 -0.20 -0.28 -0.38

-0.49

0.52

0.28

0.14 <0.001 -0.18

[O II] 7325 / [O I] 6300

-1.17 -0.98 -0.92 -0.85

-0.76

-0.66

-0.63

-0.56

-0.48

-0.30

[O II] 7325 / H↵

-1.81 -1.77 -1.74 -1.71

-1.67

0.09

-0.12

-0.30

-0.26

-0.27

-

0.54

0.62

0.85

0.77

0.85

Fraction of Ratios Fitted

-

-

-

-

Free Parameters
-1.8

Free Parameters

Density Weighting -

-1.8

-1.8

-1.8

-1.8

-1.8

-1.8

-1.8

-1.8

-1.8

Radial Weighting -

0.5

-1.0 -1.25 -1.5

-1.75

0.5

-1.0

-1.25

-1.5

-1.75

120

Figure 3.8 Line ratio diagrams constraining the excitation mechanism, displayed in
the same manner as Figure 3.4. Our model successfully fits the high to moderation
ionization AGN subsets. In the cases of [N II] 6584/ H↵ and [O III] 5007/ H we
have successfully matched the entire AGN sequence to within a factor of two.
classification of emission line galaxies. Comparing simulations with power law spectra
and starlight to observations resulted in the development of these diagrams.
We used these diagrams to make a baseline assessment of our model, since we
must be able to reproduce typical emission line ratios before exploring other diagnostic diagrams. Using the criterion of a factor of two representing an adequate
fit, we successfully matched [N II] 6584/H↵ over the entire AGN sequence, [O III]
5007/H , [O II] 3727/ [O III] 5007 and [O I]

6300/ H↵ for the a41-a31-a21-a11

sequence (higher to low ionization subsets), and [S II] 6720/H↵ for the a41-a31-a21
121

Figure 3.9 Line ratio diagrams that constrain the SED, displayed in the same manner as Figure 3.4. The extreme AGN subsets
(large symbols) are successfully fitted in essentially every diagram. Our model reproduces the high to moderate ionization
observations for all but [Ar IV] 4711 / [Ar III] 7135 and the entire sequence for all but He II 4686/ H and [O II] 7325 /
[O I] 6300.

122

sequence (higher to moderate ionization subsets). These results show that we have
su cient agreement with established excitation diagnostics to proceed on to other
line ratio diagrams.
Spectral Energy Distribution
In Panels (a)-(f) of Figure 3.9 we show less common line ratio diagrams, which draw
upon weaker lines to constrain the SED. This adds important cross-checks to the
results from the usual strong-line diagrams in previous investigations. In order to
limit the e↵ects of di↵ering abundances, these diagrams mostly use ratios of lines
from di↵erent ionization states of the same element. The S/N obtained from coadding spectra instead of using emission lines from individual AGN makes these
diagrams more trustworthy indicators. Panel (e) combines two SED indicators, [N
I] 5200/ [N II] 6584 and [O II] 7325 / [O I] 6300, both of which are minimally
a↵ected by reddening. The results show that although the observations do not show
the predicted modest changes of these line ratios as the radial weighting is changed in
the LOC model, the entire [N I] 5200/ [N II] 6584 sequence and higher ionization
observation of [O II]

7325 / [O I]

6300 still does agree with a change in radial

weighting to within approximately a factor of two.
Panel (f) on Figure 3.9 combines two SED indicating ratios, [S III] 9069/ [S II]
6720 vs. [O II] 3727/ [O I] 6300. The observations used here do not extend out
to the [S III] lines, but these lines are commonly measured in low-z galaxies and we
include them here so that other investigations can take of advantage of these results.
The ratio [O II] 3727/ [O I] 6300 is included in our observations, and ranges from
1.0 (in the log) at the high ionization end to 1.5 for the lowest ionization cases. The
adopted LOC models all predict that this ratio should be about 0.7 in the log, which
agrees with the observations to within a factor of two except for the lower ionization
cases.

123

Figure 3.10 Line ratio diagrams that constrain the physical conditions, displayed in the same manner as Figure 3.4. The excellent
agreement between our model and the [S II] 6716 / [S II] 6731 ratio over the entire range of ionization shows that low density
regions are correctly predicted by our model. High density regions, probed by the [S II] 4070 / [S II] 6720 ratio, also agree
with the high to moderate ionization observations. The [O II] 3727 / [O III] 5007 ratio, sensitive to ionization parameter, fits
the observations at all but the lowest ionization points. The temperature, however, given by the [O III] 4363 / [O III] 5008
ratio and [N II] 5755 / [N II] 6584 ratio contradicts observations in all but the highest ionization cases.
124

Physical Conditions
Figure 3.10 probes the average temperature and density. Except as noted, these
figures use standard temperature- and density-sensitive intensity ratios (see AGN3,
chapter 5). In the case of the temperature sensitive [O III] 4363/ [O III] 5007 ratio
(Panels (a)-(b)), the models come close to matching the higher ionization subsamples
but the lower ionization subsamples give higher temperatures than those given by
our models. Panel (d) confirms this finding using another temperature indicator, [N
II] 5755/ 6584. The observations agree with our models at the higher ionization
end of the sequence but at lower ionization the ratio disagrees with our models by
indicating higher temperatures than predicted. This is an important discrepancy that
we will explore further in §4.2.
Panel (a) shows a typical density sensitive ratio, [S II] 6716/ 6731, which is
in excellent agreement with our models. This ratio is very insensitive to our free
parameters, which is not surprising in these LOC models because the density indicated
by the [S II] line ratio is close to the density for optimal emission of these lines.
Panel (c) combines a low density indicator, [O II] 3727/ 7325, with a high density
indicator, [S II] 4070/ 6720 (Holt et al. 2011). Our model agrees with the higher
and lower ionization subsamples, but fails to match the moderate ionization AGN
subsamples. Panel (e) shows a diagram originally proposed by BPT, and therefore
also appears in Figure 3.8, but which is a good indicator of the ionization parameter
(Komossa & Schulz 1997; G04b). The observed [O II] 3727/ [O III] 5007 agrees
nicely with all but the lower end of the sequence and even in this region only disagrees
by a factor of 3-4.
Grains
The diagnostics in Figure 3.11 are taken from G04b, who found these line ratios to
be sensitive to the radiation pressure from dust grains and to emphasize the success
125

Figure 3.11 Line ratio diagrams sensitive to radiation pressure from grains taken
from G04b, displayed in the same manner as Figure 3.4. These dust-free models
successfully fit all of the high to moderate ionization observations except for [O III]
4363 / [O III] 5007, which also failed in G04b.
of their dusty model. However, we find that our dust-free models match the observations just as well as the G04b dusty models. For example, Panel (a) shows [N I]
5200 and [O II] 7325, which are typically weak in Seyfert galaxies, and there is decent agreement between the AGN observations and our dust free model. G04b found
dramatic di↵erences between the observations and dust-free models on this diagram.
G04a argued that the diagram in Panel (b) ([O II] 3727/ [O III] 5007 vs. He II
4686/H ) points out the success of their dusty models over their dust-free models
at high ionization parameters. In contrast, our dust-free models fit quite nicely over

126

Figure 3.12 Line ratio diagrams that constrain abundances, displayed in the same manner as Figure 3.4. The best abundance
sensitive ratio, [N II] 6584 / [O II] 3727, was used to optimize our abundance set and fits observations over the entire range
of ionization except for the lowest subsample. The models and observations agree to within the factor two error bars except for
the [O II] 7325/ H↵ ratio, which agrees with the highest ionization observation to within a factor of three.

127

a large range of ionization. The integrations in Panel (c) match our high to moderate
ionization region of the sequence. Panel (d) was used by G04b to emphasize that
their dusty model is not devoid of problems because of the inability to reproduce
high temperatures at moderate ionization; our model shows that this problem is also
present in our dust-free models. It should be noted that the classical Seyferts to
which G04b compared their models almost all lie above the dotted and solid lines in
Figure 3.1a, so that in terms of the BPT and VO87 diagrams their results should be
compared only to our 3 highest-ionization levels. These are the cases for which our
dust-free models give the best fits to the observations.
Abundances
Figure 3.12 shows ratios of lines from ions with similar ionization potentials, but from
di↵erent elements, leading to abundance sensitivity. Panel (a) stems from a discussion
in Storchi-Bergmann & Pastoriza (1990) where only nitrogen and sulfur abundances
needed enhancements to account for emission in di↵erent AGN covering a wide range
in ionization. However, our 1.4 Z

solar abundance set nicely fits this diagram for

the high - moderate ionization subsamples. Panel (b) incorporates [O II] 7325/H↵,
a ratio involving a weak line only measurable because of our large sample, which fits
our highest ionization subsamples within a factor of 3. He I 5876 depends linearly
on the helium abundance; therefore Panel (c) shows that our default solar helium
abundance is correct for all but the lower ionization parts of the sequence.
Panel (d) displays the main determinant in our metallicity optimization, [N II]
6584/ [O II] 3727, and our model fits the entire sequence to within 10 per cent
except for the lowest ionization observation. The diagram in Panel (e), [Ne III] 3869/
He II 4686 vs. [O III] 5007/ [Ar III] 7135 uses lines with fairly high ionization
potentials, and shows agreement between the models and observations over almost
the entire length of the AGN locus. The diagram in Panel (f), [O I] 6300/H↵ vs.

128

[N I] 5200/H↵ uses lines with lower ionization potentials. In this case the predicted
ratios for our adopted sequence of models all cluster at one point on this diagram, and
the observed ratios fall at that same location except for the case of [O I] 6300/H↵
for the very lowest ionization subset.

3.3.6 Dusty Models
As described in §3.1, in addition to the dust-free models we also ran grids of dusty
models, to test whether or not a dusty NLR could do an equally good job of reproducing the observations along the entire length of the AGN locus, again by tuning only
a single free parameter. A dusty NLR model was also explored by F97 and preferred
by G04a, making it natural to see if such a model is a viable alternative to a dust-free
interpretation.
The abundances we adopted for modeling a dusty NLR are given in §3.1; they
allow for depletion of refractory elements into grains at large radial distances from
the central engine. As with our dust-free models, the He II

4686/H

ratio was

overpredicted using the F97 SED. We optimized the SED using the same process
as discussed in §3.2 and arrived at an incident continuum characterized by Tcut ⇠
1.5 ⇥ 105 K, ↵UV =

0.5, ↵ox =

1.4, ↵x =

1.0, and kTIR = 0.14 eV (Figure 3.3).

Depending on the radius, our baseline abundance set changes due to the formation
of grains, therefore there are actually slightly di↵erent metallicities within a single
dusty grid. As with the dust-free models, there was the need for an increase in
metallicity or selective nitrogen enhancement. We scaled the abundances using the
empirically determined scaling relationships given in §3.4 and found that a metallicity
of 1.4 Z

in the dust-free portion of the grid gave a satisfactory fit.

Table 3.3b provides the emission-line ratio predictions and their comparison to
our subsets for our best dusty model. Our dusty model fits a few line ratios slightly
better in the moderate ionization part of the sequence. In particular, the entire [O

129

II] 7325/H↵ sequence and the moderate ionization [O III] 4363 / [O III] 5007
subsets are fit within a factor of two.
However, our dust-free model fits all but the moderate ionization subsets of [O II]
7325/H↵ to within a factor of three. Furthermore, our dusty models still have the
same problem as the dust-free models in that they predict that the [O III] 4363 /
[O III] 5007 line ratio should get larger as the ionization level increases, while the
opposite is observed to happen. The dusty models predict a systematically higher
[O III] 4363 / [O III] 5007 emission line ratio than do the dust-free models, as a
result of the gas being hotter due to photoelectric heating by dust (see Baldwin et
al 1991), so the subsets are matched at a lower observed ionization level. Our dusty
model similarly fails to reproduce the observed [N II] temperature at high ionization.
Our dusty models do not reproduce the lower ionization portion of some of the
VO87 diagrams as well as our dust-free models and fail to account for the high
ionization observations in a few ratios. Figure 3.13 shows some line ratios, which
are not fit as well in the dusty case as in our dust-free model. Fig. 3.13a should be
compared to Figure 3.8b and Figure 3.13b to Figure 3.11a. The dusty case is not ru-led out by these results, but the dust-free models give a marginally better fit at
higher ionization.

3.4 Discussion
3.4.1 Physical Meaning of the AGN Locus
Analysis of the AGN region of the BPT diagram has historically been limited to wellidentified emission-line galaxies displaying strong AGN properties. F97 studied these
galaxies using LOC models, and Komassa & Schulz (1997) also studied luminous,
high ionization Seyfert 2 galaxies using similar composite models having clouds spread
over wide ranges in radial distance and density. In this paper, we have for the first

130

Figure 3.13 Line ratio diagrams for our dusty NLR models, displayed in the same
manner as Figure 3.4. The AGN observations of some key line ratios, which are
matched by our dust-free models, are not fit as well by the dusty models shown here.
The hollow marker indicates our best set of free parameters in the dusty case.
time broadened this LOC-type analysis to include low ionization AGN, where the
spectra become very similar to those of SF galaxies. Diagnostic diagrams that include
weak lines can then provide valuable consistency checks for highly ionized AGN when
developing physically meaningful interpretations of the AGN locus.
An ideal interpretation would be to find a single, tunable physical parameter which
recreates the variation from high to low ionization along the AGN locus (first index
i designating the aij subsets), as well as a second physical parameter responsible for
moving galaxies o↵ at right angles to the AGN locus (second index j). As we discuss
below, we have found that variations in the radial weighting in the LOC model can
account for the di↵erences in the high and intermediate ionization AGN.
We note that we used a rather coarse interval of radial and density weightings
to match the AGN observations so that their impact on the LOC integrations would
be readily apparent from looking at line ratio diagrams. A finer spacing would yield
slightly better fitting to all of the observations along the length of the AGN locus,
but should not change any of our conclusions.
131

Higher - Moderate Ionization AGN
Figures 3.8-3.12 show that the higher ionization AGN are matched to within a factor
of 2 by

=

1.4 models with either

=

0.75 or

=

0.5. The only emission-

line ratios completely contrary to this statement are [O III] 4364 / [O III] 5007
and [O II] 7325/H↵ but even these ratios are in agreement to within a factor of 3
or 4. The ratio [Ar IV] 4711/ [Ar III] 7135 agrees with highest ionization locus
point and [O II] 7325/ [O I] 6300 is approximately within a factor 2 of the highest
ionization subsample. Although they do not agree for sequentially lower ionization
observations, they again are within a factor of 3 or 4. The moderate ionization AGN,
which correspond to the middle point along the AGN locus, are best fit by
and

=

=

0.5

1.4. Indeed, every strong-line ratio in our sample agrees with this set of

free parameters except for He I 5876/H . In addition to fitting ratios with strong
lines, we have also incorporated many weak lines as consistency checks, and these
also fit the higher to moderate ionization spectra. These results indicate that the
progression along the AGN locus from high ionization to moderate ionization can be
represented by a positive change in the radial weighting. Indeed, Figures 3.8-3.12
show this trend to be true for the majority of line-ratio diagrams presented. The
change in radial weighting is in the sense that the galaxies with more luminous AGN
(the ones higher up the AGN sequence) have more centrally concentrated NLRs.
However, there are still obvious discrepancies in [O III]
[N II]

5755/ 6584 and [Ar IV]

4711/ [Ar III]

an explanation. In particular, the [S II]

4363/ [O III]

5007,

7135 for which we do not have

6716/ [S II]

6731 and [O III]

4363/

[O III] 5007 ratios indicate low density, high ionization gas defining the so called
“temperature problem” in photoionization models (Komossa & Schulz 1997). As in
F97, this problem has been solved for higher ionization AGN but there does not seem
to be a relationship between a decrease in ionization and our free parameters. This
is discussed in more detail below.
132

Low Ionization AGN
Several line ratios are successfully fitted over the entire length of the AGN sequence
including the low-ionization cases. Specifically, [O III] 5007/H , [N II] 6584/H↵,
[O I] 6300/ [O III] 5007, [S II] 6716/ [S II] 6731, He II 4686/ He I 5876, [N I]
5200/ [N II] 6584, fit every observational point over the full range in ionization, and
[O III] 5007/ [Ar III] 7135 and [N I] 5200/H↵ are very nearly within a factor of 2
for every subsample. However, there are still many line ratios that are not fit at low
ionization. In addition to the [O III] 4363/ [O III] 5007, [N II] 5755/ 6584 and [Ar
IV] 4711/ [Ar III] 7135 ratios mentioned in the previous section, [O I] 6300/H↵
also deviates from the trend that a more positive radial weighting corresponds to a
lower ionization point on the locus.
We chose to vary the SED, metallicity, integration limits, density weighting and
radial weighting, as these were the most likely parameters to provide a simple understanding of the AGN locus. Not counting the integration limits, we have used four
free parameters to fit 24 di↵erent line ratios. The bottom rows in Tables 3.3a and
3.3b show the fraction of line ratios that are fitted to within a factor of two by our
models. We fitted about 80 per cent of the observed ratios at the higher ionization
end of the AGN sequence, but this drops to about 50 per cent (1113 ratios fitted)
in the two lowest-ionization AGN subsets. The variation in the free parameters cannot completely account for the lowest ionization AGN, but varying other parameters
could possibly produce an explanation still with the context of a photoionized, LOCtype model. The simulations in this paper ended when temperatures fell outside the
range in which optical emission is strong. Other stopping conditions, such as cloud
thickness, could be implemented to test their e↵ects. We have also assumed a simple
constant density model for our simulations. Constant density and constant pressure
conditions have been implemented in other models (i.e. G04a, G04b) and using similar assumptions in our individual cloud models could produce better fits at the low
133

ionization end of the AGN locus.
Another consideration is whether photoionization or shocks are the prevailing
excitation mechanism for low ionization AGN. It is generally accepted that photoionization is responsible for exciting the NLR in luminous Seyfert galaxies. Early work
suggested that shock heating was applicable in LINERS or radio galaxies (Dopita &
Sutherland 1995, 1996) but this fails to reproduce the emission found in Seyferts and
our lowest-ionization subsets. Shock models are not likely to explain the entire AGN
locus.
Instead of non-stellar continua, a starburst could photoionize gas in the low ionization AGN. Alonso-Herrero et al. (2000) found that many LINERs can be explained
by an aging starburst. While our low ionization subsets do not occupy the same region as LINERs in the VO87 diagrams, several starburst models in Alonso-Herrero et
al. (2000) overlap with these subsets. Indeed, the low ionization AGN observations
descend into a region of the BPT diagram that is also occupied by low ionization
SF galaxies. We have assumed that adjusting a small number of parameters could
account for the variation in both AGN and pure SF sequences; however, it is possible
that there is considerable overlap in the physical properties in low ionization AGN
and SF galaxies. A follow-up paper will address LOC models with multiple excitation
mechanisms, to see if these fit the low-ionization end of the AGN sequence better than
our current models. We also plan to address the SF locus in a future paper, which
may also provide clues about the deeper nature of low ionization AGN.
Interpreting the Radial Distribution of Clouds
We studied the relationship between AGN spanning a wide range of ionization. We
found that the di↵ering properties of these AGN relate to the radial distribution of
their NLR clouds, or more specifically, to the concentration of the clouds towards the
central ionizing continuum source. Higher ionization AGN, such as Seyferts, have

134

more centrally concentrated gas, while lower ionization “transition” AGN have their
NLR spread out over more extended regions of the host galaxies. We are unable to
comment on the nature of LINERs since our AGN subsets narrowly miss this region
of the BPT diagram (Figure 3.1a), and clearly miss the region of LINERs identified
by Kewley et al. (2006) on other diagrams including [S II]/H↵ and [O I]/H↵. The
above discussion is in terms of our results for dust-free models, but the dusty models
lead to the same conclusions.
Table 3.1 shows that L([O III]) increases fairly strongly as the ionization level
increases. However, this is just an ionization e↵ect caused by the clouds being able
on average to more e ciently produce [O III] emission. Kau↵mann et al. (2003)
found in their sample of SDSS AGN that L([O III]) is fairly constant as a function
of distance from the star-forming locus, which they interpreted as support for the
idea that in their sample the increase in ionization level at greater distances from
the star-forming sequence represents a “mixing sequence” involving ionization both
by starlight and by non-thermal AGN emission. Here we have found that L([O III])
does systematically depend on the position along the MFICS AGN sequence, which
is roughly equivalent to the distance from the SF locus as measured by Kau↵mann
et al. (2003). We view this as an encouraging sign that we have cleaned up the
Kau↵mann et al. (2003) sample to produce one more dominated by just a single type
of ionization and therefore better suited to studying the e↵ects of additional variable
parameters such as the size of the NLR.
A more fundamental result from Table 3.1 is that although the observed L(H ) is
roughly constant along the entire length of the AGN sequence, the dereddened L(H )
values are higher for the lower-ionization objects and then decrease by a factor of
about 5 between the low-ionization end of the AGN sequence and the high-ionization
end. Since the NLR is in photoionization equilibrium, L(H ) / Lion , assuming that
the covering factor is the same for all AGN subsets. Therefore, we can infer from

135

L(H ) that the ionizing luminosity drops as the ionization level increases along the
AGN locus. We could not measure the ionizing luminosity more directly, since even
the high ionization subsets do not show any evidence of a non-thermal continuum,
and X-ray measurements have not been made for the AGN in our subsets.
We deduce that the underlying physical trend is that the objects with more luminous central continuum sources have more extended NLRs. This suggests that the
NLR might be a wind driven by radiation pressure. The decrease in ionization level
in the more luminous objects is just an incidental result of their NLRs being more
extended. The change in the radius weighting might represent a “clearing out” of the
central regions as the AGN starts influencing the structure of the NLR. The usual
scenario is that the AGN epoch begins after the ULIRG epoch and results in winds.
The change in radius weighting from central to more distant might represent blowing
out material from closer in, leaving distant material in place.
The Possible Role of Selection E↵ects
The above interpretation assumes that (1) the NLRs in all objects are photoionized
by a central continuum source, and (2) the covering factors are the same from subset
to subset. This interpretation does not explain the other strong trend in Table 1: the
E(B

V ) values indicate that the low-ionization objects (the more extended NLRs)

are more strongly reddened than are the high-ionization objects.
The luminosity correlation does not appear until the reddening correction has
been applied. Could the luminosity correlation could be an artifact due to an incorrect reddening correction? We used H↵/H as our primary reddening indicator
and assumed that the lines have their Case B intensities. Cloudy explicitly solves
a many-level model of H I emission, including radiative transfer and collision processes (Luridiana et al. 2009), so our calculations result in explicit predictions of the
H↵/H ratio. These are listed in Table 3. While the predicted Balmer decrements

136

are slightly steeper than the Case B value that we used (log(H↵/H ) = 0.46), there
is little variation in the predicted H↵/H ratio, so that the physics of the emission
line formation could not create the e↵ect we observe.
The basic SDSS galaxy sample is continuum-flux limited. As is described in Paper
I, we selected emission-line galaxies from that sample using criteria that include a S/N
ratio limit on the H

emission line and a very narrow redshift range, so our AGN

sample is e↵ectively H -luminosity limited. For this reason, it is not surprising that
the observed L(H ) given in Table 3.1 is constant across all of the subsamples. The
strong trend in the dereddened L(H ) therefore comes strictly from the observed trend
in the H↵/H ratio from which we deduce E(B

V ). Can this H↵/H trend be due

to selection e↵ects? For this to be true, we would need to either be systematically
missing the high-ionization galaxies with the strongest H↵ lines or the low-ionization
galaxies with the weakest H↵ lines. The latter seems more probable, but even for no
reddening the H↵ line is three times stronger than H , so we do not see why such
lines would be missed or incorrectly measured.
The actual subsamples used here reflect the AGN locus that was defined by the
MFICA techniques described in Paper I. We are confident that we have correctly
studied the physical properties of that sample.

3.4.2 [O III] and [N II] Temperatures
The major problem with these LOC models is that they do not reproduce the increase
in the temperature-sensitive [O III] 4363/ 5007 and [N II] 5755/ 6584 ratios that
is observed as one moves down the AGN sequence towards lower ionization, nor do the
dust-free models reproduce the very high ratios that are reached. Our large galaxy
sample has allowed this e↵ect to be seen not only for the typical [O III] ratio but
also for the [N II] ratio. Although neither ratio can be measured in the very lowest
ionization case, Figure 3.2 shows that [O III] 4363 is easily measured at all other

137

cases and likewise, [N II] 5755 is easily measured in our high moderate ionization
subsets. In the low-density limit, the observed [O III] ratio for the high-ionization
end of the AGN sequence gives a temperature of ⇠11,000 K and the observed [N II]
ratio gives a temperature of ⇠9,000 K. The lowest-ionization subset for which the [O
III] ratio can be measured (the a11 case) corresponds to a temperature of ⇠20,000 K
in the low density limit. High temperatures have long been known to exist in some
AGN (Shuder & Osterbrock 1981) and have also been discussed by Komossa & Schulz
(1997) and Zhang, Liang & Hammer (2013).
The [O III] 5007 and [N II] 6584 lines become collisionally quenched at high
densities, while [N II] 5755 and [O III] 4363 optimally emit. This could mean that
the so called “temperature problem” is actually a density problem and that the [O
III] and [N II] ratios are giving false temperatures. We used simple power laws to
represent the radial and density distributions in our LOC models, but there could in
principal be some other arrangement of the NLR that emphasizes high density gas.
However, there are several problems with this scenario.
First, the [N II] 5755 and [O III] 4363 optimally emit in di↵erent regions of
the LOC plane. So emphasizing a moderately high density (e.g. log(nH ) ⇠ 6.5)
would solve the [N II] temperature problem but give an even worse fit to the [O III]
temperature. Second, our low-density indicating line ratio, [S II] 6716/ [S II] 6731,
fits the entire AGN sequence to within ⇠10 per cent and our high-density indicating
line ratio, [S II] 4070/ [S II] 6720, fits the AGN sequence to within a factor of 2
or 3 down to all but the lowest ionization observed AGN subset. Finally, additional
high density components would add substantial Balmer emission (since these lines
emit strongly over a wide density range) which would result in underpredicting many
key ratios (e.g. [O III] 5007/H , [S II] 6720/ H↵, etc.). Indeed, Komossa & Schulz
(1997) attempted to solve the “temperature problem” by emphasizing a high density
component, only to find that it severely overestimated [O I] emission. Therefore, it

138

is likely that the low ionization AGN really do contain low density gas with high Te .
Figure 3.7 shows that [O III] 4363 optimally emits at log(r) ⇠ 18.0 while [O
III] 5007 optimally emits at log(r) ⇠ 19. Similarly, [N II] 5755 optimally emits at
log(r) ⇠ 19.5 while [N II] 6584 optimally emits at log(r) ⇠ 21.0. Therefore smaller
values of the radial exponent (which correspond to less extended NLRs) should give
larger [O III] and [N II] ratios and hence larger deduced Te . However, the observations
surprisingly show that the higher ionization cases have smaller [O III] and [N II] ratios
even though the other line ratios indicate that higher ionization cases are the ones
with the more concentrated NLRs (smaller fitted values of ).
One possible explanation for the high temperatures is an increasing contribution of
shock excitation at lower ionization. The shock models presented by Dopita & Sutherland (1995) do reach the high [O III] temperatures found in our subsets. However,
their shock models that produce these high temperatures significantly overestimate
the [O I] emission given in our lower ionization subsets. Most importantly, their shock
models cannot reproduce the emission found in high excitation Seyfert galaxies, so
our hope of unifying the low- and high-ionization AGN would have to be abandoned.
Another factor that could explain the high [O III] ratio in our low-ionization AGN
is a decrease in metallicity. A decrease in metal abundances results in a decrease in
cooling which raises the electron temperature. Komossa & Schulz (1997) studied
this e↵ect in luminous Seyfert 2 galaxies. Their models that resulted from scaling
solar metal abundances by 0.3, which they favored, reproduce the [O III] ratio found
in our low ionization subsets. This e↵ect has an upper limit of log([O III] 4363/
5007) ⇠

1.3, but this is already larger than any [O III] ratio found in our subsets.

However, studies with far larger samples than the one used by Komossa & Schulz
(1997) have shown that Seyferts with subsolar metallicities are rare (Groves et al.
2006), so this explanation may be limited to low ionization AGN.
Photoelectric heating from grains can also increase the electron temperature enough

139

to reproduce the temperature sensitive ratios found in the high-ionization Seyfert
galaxies studied by Komossa & Schulz (1997) and G04b. Our dusty models have
the same e↵ect, producing an [O III] 4363/ 5007 ratio that matches the moderate
ionization subsets, but they do not accurately predict the highest ionization observations.
This problem might be ameliorated in models adding cosmic rays, shock excitation, turbulent heating or an increase in grain abundance in more extended regions
of the NLR. Although a full study must be left for a future paper, we briefly explored
one of these possibilities by running our dust-free and dusty model with a cosmic ray
ionization rate 103 times greater than the Galactic background given by Indriolo et
al. (2007). It is possible that AGN and starburst galaxies see a greater cosmic ray
flux than other galaxies (Papadopoulos 2010). However, the extra heating by cosmic
rays had a negligible e↵ect on the [O III] and [N II] temperatures in both our dust
free and dusty models, leaving the need for a deeper investigation of an additional
component to address the temperature problem in low ionization AGN.

3.4.3 The Possibility of a Dusty Narrow Line Region
We were unable to find a set of free parameters for our dusty models that can account
for systematic variations along the AGN locus. However, our results do not rule out
the dusty models. Describing the NLR as either completely dusty or dust-free is likely
to be a gross simplification for a complicated environment; it is indeed plausible that
both dusty and dust-free regions exist.
A dusty, radiation pressure dominated NLR has previously been proposed by
G04a and provides a di↵erent interpretation of the physical meaning behind the AGN
locus. The basic premise behind that model is that radiation pressure exerted by the
ionizing source upon grains, the photoelectric heating produced by the grains, and
the ionizing photon absorption are non-negligible e↵ects in the NLR. Similar e↵ects

140

occur in individual Galactic H II regions (e.g. the Orion Nebula; Baldwin et al.
1991). Aside from these e↵ects, grains also deplete refractory elements and destroy
resonance lines (e.g. Ly↵ 1216).
In this case, the theoretical diagnostic diagrams computed by G04b indicate that
for Z = 2Z

a decrease in ionization along the AGN sequence represents a decrease

in ionization parameter for a dusty, constant pressure NLR. The G04a models do not
even consider the low ionization end of the AGN locus, nor do they fit the observed
values for those cases. In addition, their interpretation was mostly limited to strong
lines. For the weak line ratios that were included in this analysis, our dust-free models
fit our subsets better than our dusty models and their dusty models clearly do not
reproduce our subsets. In fact, their constant density, dust-free models are in better
agreement with these lines.
Dusty models, as in G04b for example, are only capable of reproducing a narrow
range of high ionization UV lines in NLR Seyferts (Nagao, Maiolino & Marconi 2006).
In particular, dusty models for C IV 1551/ He II 1640 > 1.5 fail to match observations due to stagnation at low metallicity and high ionization parameter, while the
dust-free models from F97 are successful. This suggests that at least part of the high
ionization region of the NLR is dust free. One current physical picture envisions the
NLR as a mixture of a dusty region responsible for lower ionization emission (e.g. [O
III], [S II]) and a dust free component that emits coronal lines (e.g. [S VIII], [Fe X])
(Dopita et al. 2002).
Our analysis results in contradictory conclusions about the dust content. Komossa
& Schulz (1997) favored a dust-free model mainly because their dusty model failed
to predict strong, moderate ionization Fe line intensities (e.g. [Fe III] 4658, [Fe VII]
6087) due to highly depleted Fe abundances. F97 found that lower ionization regions
were dusty and higher ionization regions dust free.
Low ionization lines may suggest that the gas is dusty. The [Fe III] line is not

141

detected, with I([Fe III] 4658)/I(H ) < 0.020.04. The predicted [Fe III]/H ratios
range from 0.05 for our highest ionization dust-free LOC model to 0.01 for our moderate ionization dusty model, which are too close to the observational limits to be
conclusive. But at most, marginal depletion of Fe might be required to fit the observations. The [Ca II]

7291, 7324 doublet is typically strong in dust-free regions

(Kingdon, Ferland & Feibelman 1995). The [Ca II] 7324 line is blended with [O II]
7325, but we do not detect the [Ca II] 7291 line. Our dust-free models predict that
I( 7291)/I(H ) = 0.31 for higher ionization AGN, and they also overpredict the [Ca
II] 7324/H ratio at I( 7324)/I(H ) = 0.21 since it is predicted to be larger than
the [O II] 7325/H↵ emission line ratio. This indicates significant depletion of Ca
onto dust.
On the other hand, high ionization lines suggest that the gas is dust free. The [Fe
VII] 6087 emission line is observed to be moderately strong in our higher-ionization
subsets. The dust free model does correctly predict the [Fe VII] 6087/H↵ ratio over
the range of ionization that it is observed, while the dusty model underpredicts this
line by almost a factor of three.
All of this leaves the overall picture of the NLR’s dust content unclear, and it
likely that both high ionization dust-free and lower ionization dusty clouds exist in
regions of the NLR. An ablating wind, as mentioned above, from inner regions could
be dust-free and lead to a mixture of dust abundances.
F97 came across a similar problem. They found that low-ionization regions appeared to be dusty while high ionization region must be dust free to create strong
high ionization Fe lines. The spectral diagnostics do not clearly indicate the ionization potential where the change occurs, so we cannot say whether any particular line
might form in a dusty or dust free region. The results in Figure 3.13, which show
that the intensity ratios involving low-ionization lines, such as I([N I] 5200)/I([O II]
7325) or I([O I] 6300)/I(H↵), fit the dust-free models better than the dusty models

142

add to the confusion.

3.5 Conclusions
Spectral signatures of the dust content are inconclusive, as found in previous studies.
It is likely that di↵erent regions have a di↵erent dust to gas ratio, perhaps because of
their history. Luckily, the key spectral diagnostics used to reach the main conclusions
of this paper are insensitive to whether or not our LOC models contain dust.
Our LOC models, in which only the radial extent of the gas is varied, successfully reproduce most of the BPT and VO87 diagnostic line ratios observed over the
sequence of high- through moderate-ionization NLRs in narrow-lined AGN. This includes fitting the line ratios which depend most strongly on the abundances, gas
density and spectral energy distribution. Our sequence of models also fit a number
of the observed line ratios even for very low-ionization AGN. The most significant
discrepancies between the observations and our models are that we cannot fit several
important line ratios in the low-ionization objects, and in particular the models do
not reproduce the high Te implied by the large observed [O III] 4363/ 5007 ratio.
We have shown that the high - moderate ionization AGN along the locus generated
by MFICA represent a physically meaningful sequence in NLR properties. The most
viable interpretation describes the sequence from high ionization to low ionization
AGN as a change in the radial distribution of the NLR gas, which is proportional
to changing the flux that is incident on each cloud. This finding states that higher
ionization AGN contain optimally emitting clouds that are more concentrated towards
the central continuum source than in lower ionization AGN, but that the density
distribution is on average the same. Scaling from the H luminosity, we have found
that the central continuum source in low-ionization objects is more luminous than
in high-ionization objects. The ionization sequence might be an age sequence, where
the low-ionization objects are older and have systematically cleared out their central
143

regions by radiation pressure.
The first index of the AGN locus can be understood by the radial distribution of
clouds in the NLR, but the meaning behind the second index is most likely objectto-object scatter. The majority of diagnostic diagrams in Figure 3.8 show that the
observed subsets corresponding to the two “wings” parallel to the AGN sequence
(j = 0 and j = 2) are so similar to the central AGN locus (j = 1) that we cannot
deduce any physical meaning from their di↵erences. This supports describing the
AGN sequence, from high at least through medium ionization, as a single-parameter
sequence in the radial concentration of the NLR.
We arrived at these results by not only matching typical line ratio diagnostics
based on the strongest emission lines but also by taking advantage of weaker lines
for consistency checks. This allowed greater confidence in the temperature, density,
abundances and ionizing continuum constrained from strong line ratios. We will
continue this analysis with the SF sequence in a future paper where the physical
properties of low ionization SF galaxies might overlap with those of low ionization
AGN.

144

4 Theoretical Modeling of the
Ionization Sequence in
Star-Forming Galaxy
Emission-Line Spectra
High ionization star forming galaxies are easily identified with strong emission line
techniques such as the BPT diagram, but for ionization levels below log([O III]/H ) ⇠
0.3 they become confused with low-ionization AGN, making their physical interpretation di cult. Mean field independent component analysis (MFICA) is a novel approach to processing emission line spectra that can disentangle the AGN and starlight
contributions to emission lines allowing the properties of pure AGN and pure starburst galaxies to be interpreted over a wide range of ionization. We applied MFICA
to a large sample of low-z SDSS galaxies and created subsamples of pure star forming
galaxies, resulting in a sequence of varying ionization. We used a locally optimally
emitting cloud (LOC) model to fit emission line ratios that constrain the excitation
mechanism, spectral energy distribution, abundances and physical conditions. Preliminary results in fitting diagrams that constrain the excitation mechanism indicate
that the variation of starburst galaxies is due to a change in the radial distribution
of clouds and confirms that MFICA is a powerful tool to assess di↵erences in emission line properties solely due to starbursts. We briefly discuss future work that will
decipher other properties in star forming galaxies.

145

4.1 Introduction
Starburst galaxies typically exhibit massive star formation, leading to an emission-line
spectrum that is dominated by gas ionized by new O stars formed in molecular clouds.
Early work by Baldwin, Phillips & Terlevich (1981), later expanded upon by Veilleux
& Osterbrock (1987; hereafter VO87), showed that the intensity ratios of strong
emission lines (known as “strong-line ratios”) found in active galactic nuclei (AGN)
and star forming (SF) galaxies can be used to classify the excitation mechanism
in emission line galaxies. In particular, the plot of [O III]

5007/ H

vs. [N II]

6584/ H↵ (hereafter the BPT diagram) has been remarkably successful at separating
galaxies highly ionized by a central AGN from those ionized by intense starbursts. On
the BPT diagram, SF galaxies fall in a well-defined region that spans approximately
an order of magnitude in these ionization-sensitive line ratios.
One viable interpretation for explaining why SF galaxies exhibit such a wide range
of ionization comes from di↵erences in ionization parameter and metallicity (Kewley
et al. 2001). The ionization parameter, U , is a convenient dimensionless number that
describes the degree of ionization in a gaseous cloud and is defined as U = (H)/nH c,
where (H) is the ionizing photon flux and nH is the hydrogen density. Optical line
ratio diagrams found in Veilleux & Osterbrock (1987) were successfully reproduced
over a large range in ionization by Kewley et al (2001). Subsequent work extended
the analysis to include infrared line ratio diagnostics (Snijders, Kewley & van der
Werf 2007) and low metallicity galaxies via the abundance sensitive [N II]/[O II]
ratio (Levesque, Kewley & Larson 2010).
While this interpretation is valid at higher ionization, the fit at lower ionization
becomes more uncertain due to possible entanglement with emission line properties
resulting from an active galactic nucleus. Figure 4.1a shows the BPT diagram for
the sample of galaxies used here (which is the same as was described in Paper I
(Allen et al. 2013) and Paper II (Richardson et al. 2013). The solid and dashed
146

Figure 4.1 (Top panel) The BPT diagram for our sample of galaxies from Paper I.
The red triangles show our sequences of pure SF picked out by MFICA with s4j
representing the highest ionization observations. (Bottom panel) The 3-D MFICA
classification diagram with axes representing the strengths of the first three MFICA
components. In this representation the AGN locus appears as the red line, and the
SF locus is the blue line.

147

lines represent classification curves from Kewley et al. (2001) and Kau↵mann et
al. (2003), and respectively define the upper and lower limits on the number of SF
galaxies. The dotted line represents the Kau↵mann et al. (2003) line that divides
a class of low ionization objects (LINERs) from more classical AGN. The correct
classification of AGN and SF galaxies is readily apparent at higher ionization due
to the large degree of separation in [N II] 6584/ H↵ above log([O III/H ]) ⇠ 0.3.
However below this threshold, the two “branches” of AGN and SF galaxies merge to
form a region containing a mix of star-forming galaxies, AGNs and composite cases.
The actual theoretical boundary between AGN and SF galaxies, and the degree of
mixing between the two, is a matter of debate (Kewley et al 2001; Kau↵mann et al.
2003; Kewley et al. 2006; Stasi´ska et al. 2006). Therefore, it is likely that previous
n
models are limited in their interpretation of galaxies that lie inside the composite
region. Furthermore, optical analysis has largely been limited to classical line ratio
diagrams from Veilleux & Osterbrock (1987). Additional diagrams could provide
additional constraints or consistency checks on models developed from strong line
techniques.
In Allen et al. (2013) we showed that a variant of blind source separation, known
as mean field independent component analysis (MFICA), successfully alleviates the
restriction of tangled emission line components from AGN and starlight excitation.
We direct the reader to Paper I and Paper II for a more detailed description of the
MFICA procedure, but we provide a summary of it here. Using MFICA, we generated
a series of continuum and emission line components using a small subset of our ⇠ 104
low redshift (0.1 < z < 0.12) galaxies sample from the Sloan Digital Sky Survey
(SDSS) to match our full galaxy sample. The unique signatures of AGN and SF
galaxies in our MFICA parameter space allowed the isolation of “pure” SF and AGN
galaxies over a large range of ionization. Our sample size also enabled many weaker
lines to be accurately measured by coadding spectra.

148

Figure 4.1b, taken from Figure 1b in Paper II, again shows the BPT diagram but
in a more abstract sense with the axes representing components from MFICA. The
degree of separation between both sequences is clear. All but the lowest ionization
point of each sequence occupies a unique parameter space. This outlines the key
point that a single parameter is responsible for the variation in the AGN and SF
components in emission line galaxies.
In Paper II we took advantage of this degree of separation by performing plasma
simulations with the spectral synthesis code Cloudy (Ferland et al. 1998) to reproduce
AGN observations over a wide range of ionization with a locally optimally emitting
cloud (LOC) model (Ferguson et al. 1997). This model states the cumulative observed
emission from emission line galaxies is the result of selection e↵ects stemming from
various emission lines optimally emitted from a large number of gas clouds spanning
a vast range in physical conditions. After integrating over clouds with a wide range
of gas density and radial distance from the ionizing continuum source, weighting the
gas clouds by power laws in these two parameters, we found that almost all of our
higher to moderate ionization diagnostic line ratios could be represented by change
in radial integration weighting. In several cases, especially for the classical line ratio
diagrams, we fit the entire sequence of observations from Seyferts all the way down
to low ionization nuclear emitting regions (LINERs). This finding emphasized the
success of MFICA in isolating galaxies with only AGN activity. In addition, we also
used weaker lines in diagnostic diagrams to act as consistency checks on our models.
In the present work, we perform a similar analysis for our subsets of SF galaxies.
LOC models have successfully reproduced observations in the NLR (Ferguson et al.
1997) and BLR (Korista et al. 1997) of AGN and here, for the first time, we extend
the use of LOC models to describe emission from SF galaxies. As was the case in
Paper II for AGN galaxies, finding a physically meaningful parameter to reproduce
variation of SF galaxies over the entire length of the SF locus is the goal of this work.

149

Our overall aim is to test the hypothesis that the range in ionization present in the
SF sequence is either due to the variation in the central concentration of gas clouds
(as in Paper II), starburst age or metallicity.
In §4.2, we describe five roughly equally spaced points along the MFICA SF locus.
Then in §4.3 we describe our current LOC model and the preliminary results for
classical diagrams that constrain excitation mechanisms. Finally, we outline future
work in §4.4.

4.2 Observations
In Paper I, we showed that MFICA allows a degree of separation between SF galaxies
and AGN that we can exploit here in order to define and interpret a pure sample
of low-ionization star forming galaxies over a large range of ionization. Here, we
present five subsamples of SF galaxies identified by MFICA that extend over the
entire length of the SF locus. These subsamples correspond to “pure” SF galaxies
where the emission lines are solely the result of starburst excitation. We also identify
two parallel wing sequences of SF galaxies, on opposite sides of the central sequence.
After the continuum components were fitted to each individual galaxy using MFICA,
a composite emission line spectrum was generated at each locus point by co-adding the
continuum subtracted spectra in each subsample. We adopt the same nomenclature
as in Paper II by naming the points along the sequence in the format sij. The first
index, i, indicates the ionization with i = 4 corresponding to high ionization starburst
galaxies and i = 0 corresponding to low ionization SF galaxies. The second index,
j, describes the position of the sequence orthogonal to the SF locus with j = 1
designating the central sequence, j = 0 the subsets positioned closest to galaxies with

150

Table 4.1: Properties of sij subsets
Subset
E(B

s00
V)

s01

s02

s10

s11

s12

s20

s21

s22

s30

s31

s32

s40

s41

s42

0.69 0.68 0.67 0.37 0.44 0.49 0.28 0.28 0.30 0.15 0.20 0.27 0.14 0.14 0.14

Component 1 0.00 0.02 0.05 0.00 0.01 0.02 0.00 0.04 0.10 0.18 0.24 0.30 0.53 0.52 0.52
Component 2 0.06 0.04 0.02 0.32 0.31 0.29 0.57 0.54 0.50 0.74 0.59 0.43 0.47 0.44 0.40
Component 3 0.93 0.90 0.87 0.68 0.67 0.65 0.42 0.41 0.39 0.05 0.16 0.27 0.00 0.02 0.06

151

AGN and j = 2 the subsets located on the opposite side.
Figure 4.1 displays these sequences on the traditional BPT diagram and our “3D”
BPT diagram. Figure 4.1a shows the highest ionization SF galaxies are isolated, leaving little ambiguity about their excitation mechanism. However, the lowest ionization
subsets descend into regions close to that of composite galaxies, leading to questions
about whether AGN have contributions to emission lines seen in this low ionization
SF subsets. However, as Figure 4.1b shows, MFICA has enabled the isolation of quite
pure subsamples of galaxies. At the low ionization end of the SF sequence, there is
minimal overlap with the AGN sequence.
Table 4.1 lists the MFICA component weights that define each of the subsamples
along the SF sequence. From the large SDSS sample, we identified for each subsample
the 50 galaxies which lie closest to the defining point of that subsample in terms of
their MFICA component weights, and co-added their spectra to make an average
subsample spectrum. Figure 4.2 displays several regions of the co-added spectra for
our central sequence; the wing sequences are not displayed due to their close similarity
to the central sequence.
As was done in Paper II for the AGN in the sample, the emission lines present
in the co-added spectra of these SF galaxies were measured by subtracting a locally
fitted continuum for each spectrum and then integrating the flux within each line
profile. The [S II] doublet was separated by integrating to the lowest point in the
blended region between the two lines. Strong Na D absorption can be seen in our
lowest ionization subsets blending together with He I 5876. We de-blended these
lines by fitting the He I line with the profile from H , yielding an uncertainty of
around 10 per cent. Table 4.2a lists our observed emission line fluxes and Table 4.2b
lists those same lines but dereddened assuming a standard galactic reddening curve
with an E(B

V ) that gives a dereddened I(H↵)/I(H ) = 2.86, appropriate for Case

B recombination lines with an electron temperature Te = 104 K and an electron den-

152

Figure 4.2 Observed, coadded spectra for the five subsamples that fall directly along
central the SF locus. Panel (a) shows the full spectra. The other panels show enlargements that slightly overlap in wavelength. All flux values are shown in units of
the peak H intensity in the particular spectrum.
153

Table 4.2a: The measured emission line strengths for the SF locus. Measurements
are relative to H .
Ion
[O II]

air

s00

s01

s02

s10

s11

s12

s20

s21

s22

s30

s31

s32

s40

s41

s42

3727

0.73

0.74

0.77

1.28

1.34

1.39

2.05

2.01

1.92

2.66

2.36

2.08

2.11

2.10

2.08

HI

3750 <0.03 <0.03 <0.04

0.04

0.03

0.03

0.04

0.03

0.05

0.05

0.04

0.04

0.04

0.04

0.04

HI

3771 <0.01 <0.01 <0.00

0.03

0.04

0.01

0.05

0.04

0.04

0.05

0.05

0.05

0.05

0.05

0.05

HI

3798 <0.00 <0.01 <0.02

0.05

0.06

0.04

0.05

0.06

0.06

0.06

0.06

0.06

0.06

0.06

0.06

HI

3835 <0.02 <0.03 <0.04

0.06

0.05

0.05

0.05

0.05

0.07

0.07

0.07

0.07

0.07

0.07

0.07

[Ne III] 3869

0.05

0.05

0.05

0.07

0.03

0.06

0.05

0.05

0.06

0.15

0.11

0.10

0.22

0.22

0.23

HI

3889

0.08

0.09

0.09

0.14

0.13

0.12

0.16

0.16

0.18

0.19

0.18

0.18

0.20

0.20

0.20

HI

3970

0.11

0.13

0.14

0.11

0.13

0.14

0.13

0.14

0.16

0.19

0.18

0.18

0.22

0.22

0.23

4070 <0.05 <0.04 <0.01 <0.01 <0.02

0.03

0.03

0.03

0.02

0.03

0.02

0.02

0.02

0.02

0.02

4102

0.18

0.18

0.16

0.20

0.21

0.21

0.21

0.21

0.23

0.24

0.23

0.24

0.25

0.25

0.25

[Fe V]? 4229

0.04

0.05

0.06

0.03

0.02

0.34

0.34

0.34

0.40

0.40

[S II]
HI

HI

4340

<0.01 <0.01 <0.01 <0.00 <0.01 <0.00 <0.01 <0.00 <0.00 <0.00
0.40

0.40

0.42

0.42

[O III] 4363 <0.00 <0.00 <0.00 <0.00 <0.00 <0.00 <0.00 <0.00 <0.01
He I

4471 <0.01 <0.00 <0.00 <0.02 <0.00 <0.00

0.03

He II

4686 <0.00 <0.00 <0.01 <0.01 <0.01 <0.00 <0.00 <0.00 <0.01
154

0.02

0.02

0.43

0.43

0.43

0.45

0.45

0.45

0.01

0.01

0.01

0.03

0.03

0.03

0.03

0.03

0.03

0.03

0.03

0.03

0.00

0.01

0.01

0.01

0.01

0.01

Table)4.2a)(cont1d)
Ion

air

s00

s01

s02

s10

s11

s12

[Ar IV] 4711 <0.01 <0.01 <0.01 <0.03 <0.02 <0.02
HI

s20

s21

s22

s30

s31

s32

s40

s41

s42

0.02

0.01

0.01

0.01

0.01

0.01

0.01

0.01

0.01

4861

1.00

1.00

1.00

1.00

1.00

1.00

1.00

1.00

1.00

1.00

1.00

1.00

1.00

1.00

1.00

[O III] 4959

0.07

0.08

0.10

0.09

0.12

0.15

0.18

0.25

0.30

0.68

0.58

0.55

1.05

1.06

1.08

[O III] 5007

0.18

0.20

0.26

0.21

0.33

0.43

0.46

0.68

0.90

1.16

1.76

1.63

2.52

2.54

2.56

0.04

0.05

0.04

0.02

0.03

0.04

0.03

0.02

0.03

0.03

0.03

0.03

0.02

0.02

0.02

[N I]

5200

[N II]

5755 <0.01 <0.02 <0.01 <0.01 <0.02 <0.01 <0.01 <0.00 <0.01

0.00

0.00

0.01

0.00

0.00

0.00

He I

5876

0.09

0.09

0.09

0.09

0.11

0.12

0.12

0.12

0.13

0.12

0.13

0.14

0.13

0.13

0.13

[O I]

6300

0.10

0.10

0.11

0.07

0.08

0.11

0.09

0.08

0.09

0.09

0.09

0.09

0.06

0.06

0.06

[O I]

6363

0.03

0.03

0.05

0.02

0.03

0.03

0.03

0.02

0.03

0.03

0.03

0.03

0.02

0.02

0.02

[N II]

6548

0.70

0.69

0.71

0.43

0.48

0.51

0.32

0.33

0.34

0.14

0.20

0.28

0.13

0.13

0.13

HI

6563

5.86

5.79

5.74

4.21

4.49

4.77

3.82

3.83

3.92

3.34

3.51

3.78

3.31

3.30

3.29

[N II]

6584

2.32

2.30

2.36

1.51

1.61

1.71

1.10

1.11

1.11

0.50

0.66

0.89

0.41

0.41

0.40

He I

6678

0.02

0.02

0.03

0.03

0.03

0.04

0.03

0.04

0.04

0.04

0.04

0.04

0.04

0.04

0.04

[S II]

6716

0.77

0.78

0.79

0.66

0.72

0.78

0.71

0.64

0.59

0.51

0.50

0.49

0.33

0.33

0.32

[S II]

6731

0.59

0.59

0.59

0.50

0.53

0.59

0.52

0.49

0.45

0.38

0.38

0.39

0.26

0.26

0.25

[Ar III] 7135

0.02

0.03

0.04

0.04

0.05

0.05

0.06

0.06

0.07

0.07

0.08

0.09

0.09

0.09

0.09

155

Table)4.2a)(cont1d)
Ion
[O II]

air

s00

s01

s02

s10

s11

s12

s20

s21

s22

s30

s31

s32

s40

s41

s42

7325

0.07

0.07

0.08

0.02

0.06

0.06

0.06

0.05

0.06

0.06

0.06

0.07

0.06

0.06

0.06

156

Table 4.2b: The dereddened emission line strengths for the SF locus. Measurements
are relative to H .
Ion
[O II]

air

s00

s01

s02

s10

s11

s12

s20

s21

s22

s30

s31

s32

s40

s41

s42

3727

1.59

1.59

1.63

1.94

2.17

2.41

2.80

2.75

2.70

3.15

2.94

2.81

2.46

2.45

2.42

HI

3750 <0.06 <0.07 <0.08

0.06

0.05

0.05

0.05

0.04

0.06

0.05

0.05

0.06

0.05

0.05

0.05

HI

3771 <0.03 <0.02 <0.01

0.05

0.06

0.02

0.07

0.05

0.06

0.06

0.06

0.07

0.06

0.06

0.06

HI

3798 <0.01 <0.03 <0.04

0.07

0.10

0.06

0.07

0.07

0.08

0.08

0.08

0.08

0.07

0.07

0.07

HI

3835 <0.04 <0.06 <0.08

0.09

0.07

0.08

0.07

0.07

0.09

0.08

0.08

0.09

0.08

0.08

0.08

[Ne III] 3869

0.11

0.10

0.09

0.10

0.05

0.09

0.07

0.06

0.08

0.17

0.14

0.13

0.25

0.26

0.26

HI

3889

0.17

0.18

0.17

0.20

0.20

0.20

0.21

0.21

0.24

0.22

0.22

0.24

0.23

0.23

0.23

HI

3970

0.20

0.24

0.26

0.15

0.19

0.21

0.17

0.18

0.20

0.22

0.22

0.23

0.25

0.25

0.26

4070 <0.08 <0.06 <0.02 <0.02 <0.02

0.04

0.04

0.04

0.03

0.04

0.03

0.03

0.02

0.02

0.02

4102

0.30

0.30

0.28

0.27

0.29

0.31

0.27

0.26

0.29

0.27

0.27

0.30

0.28

0.28

0.28

[Fe V]? 4229

0.07

0.07

0.10

0.04

0.02

0.01

0.49

0.50

0.49

0.50

0.51

0.52

[S II]
HI

HI

4340

<0.01 <0.01 <0.00 <0.01 <0.00 <0.01 <0.00 <0.00 <0.00
0.46

0.49

0.50

[O III] 4363 <0.00 <0.01 <0.00 <0.00 <0.00 <0.00 <0.00 <0.00 <0.01
He I

4471 <0.01 <0.01 <0.00 <0.02 <0.00 <0.00

0.03

He II

4686 <0.00 <0.01 <0.01 <0.01 <0.01 <0.00 <0.00 <0.00 <0.01

157

0.02

0.03

0.47

0.48

0.50

0.49

0.48

0.48

0.01

0.01

0.01

0.03

0.03

0.03

0.03

0.03

0.03

0.03

0.03

0.03

0.00

0.01

0.01

0.01

0.01

0.01

Table)4.2b)(cont1d)
Ion

air

s00

s01

s02

s10

s11

s12

[Ar IV] 4711 <0.01 <0.01 <0.01 <0.03 <0.02 <0.02
HI

s20

s21

s22

s30

s31

s32

s40

s41

s42

0.02

0.01

0.01

0.01

0.01

0.01

0.01

0.01

0.01

4861

1.00

1.00

1.00

1.00

1.00

1.00

1.00

1.00

1.00

1.00

1.00

1.00

1.00

1.00

1.00

[O III] 4959

0.07

0.08

0.10

0.09

0.11

0.15

0.17

0.24

0.29

0.67

0.57

0.54

1.04

1.05

1.07

[O III] 5007

0.16

0.18

0.24

0.20

0.31

0.41

0.44

0.66

0.87

1.14

1.71

1.58

2.48

2.49

2.52

0.03

0.04

0.04

0.02

0.03

0.03

0.03

0.02

0.03

0.03

0.03

0.03

0.02

0.02

0.02

[N I]

5200

[N II]

5755 <0.01 <0.01 <0.01 <0.01 <0.01 <0.01 <0.00 <0.00 <0.01

0.00

0.00

0.01

0.00

0.00

0.00

He I

5876

0.06

0.05

0.05

0.07

0.08

0.09

0.10

0.10

0.10

0.11

0.11

0.11

0.12

0.12

0.12

[O I]

6300

0.05

0.05

0.06

0.05

0.06

0.07

0.07

0.06

0.07

0.08

0.07

0.07

0.06

0.06

0.06

[O I]

6363

0.01

0.02

0.03

0.01

0.02

0.02

0.02

0.02

0.02

0.03

0.02

0.02

0.02

0.02

0.02

[N II]

6548

0.34

0.34

0.36

0.29

0.31

0.31

0.24

0.25

0.25

0.12

0.16

0.21

0.11

0.11

0.11

HI

6563

2.86

2.86

2.86

2.86

2.86

2.86

2.86

2.86

2.86

2.86

2.86

2.86

2.86

2.86

2.86

[N II]

6584

1.12

1.13

1.17

1.02

1.02

1.02

0.82

0.83

0.81

0.43

0.53

0.67

0.35

0.35

0.35

He I

6678

0.01

0.01

0.02

0.02

0.02

0.02

0.02

0.03

0.03

0.03

0.03

0.03

0.03

0.03

0.03

[S II]

6716

0.36

0.37

0.37

0.44

0.44

0.45

0.52

0.47

0.42

0.44

0.40

0.36

0.28

0.28

0.28

[S II]

6731

0.27

0.27

0.28

0.33

0.33

0.34

0.38

0.35

0.32

0.33

0.31

0.29

0.22

0.22

0.22

[Ar III] 7135

0.01

0.01

0.02

0.03

0.03

0.03

0.04

0.04

0.05

0.06

0.06

0.06

0.07

0.07

0.07

158

Table)4.2b)(cont1d)
Ion
[O II]

air

s00

s01

s02

s10

s11

s12

s20

s21

s22

s30

s31

s32

s40

s41

s42

7325

0.03

0.03

0.03

0.01

0.03

0.03

0.04

0.04

0.04

0.05

0.05

0.05

0.05

0.05

0.05

159

sity ne = 102 cm 3 (Osterbrock & Ferland 2006; hereafter AGN3). We compare our
models below to the dereddened line strengths. Our large sample of pure SF galaxies
has enabled the measurement of several weak lines that are not usually detectable
in lower-quality spectra, including [S II] 4070, [Ar IV] 4711, [N I] 5200, [Ar III]
7135 and [O II] 7325.

4.3 Preliminary Modeling Results
The LOC model is based on a grid of individual models at di↵erent densities and
di↵erent distances from the source. After determining the spectral energy distribution
(SED) and chemical abundances, the main free parameters correspond to the radial
distribution of clouds from the ionizing source and the density distribution of those
clouds.
We used version c10.00 of the plasma simulation code Cloudy to compute our
LOC models. As in Paper II, and Ferguson et al. (1997), we assumed there is no
attenuation of the incident continuum by ISM clouds. Our solar abundance set comes
from Groves et al. (2004a), which were taken from a series of papers (Asplund et al.
2000; Asplund 2000, 2003; Allende Prieto et al. 2001, 2002); otherwise, our heavy
element abundances come from the solar composition in Ferguson et al. (1997).
Typical starburst galaxies have luminosities Lion = 1043.5 ergs s 1 so we used
this value for our central ionizing source. Plasmas with Te < 4000 K or Te > 105 K
do not principally contribute to optimal emission, therefore we use these conditions
as limits on the Cloudy models used in our simulations.
We ran grids of simulations varying the distance from the source r, from 1016
cm to 1023 cm and the total hydrogen density, nH from 100 cm 3 to 1010 cm 3 .
The resolution of our grid resulted in 7171 total clouds evenly distributed in the log
of our parameter space. Each simulation assumed constant total hydrogen density;
however, the molecular and electron density were allowed to change with depth. We
160

assumed that the ionizing source was isotropically emitting so that (H) / r2 . For
our preliminary analysis, all of our simulation grids were dust-free.
The LOC model is based on the assumption that the emission we observe stems
from selection e↵ects; optimal emission line production occurs in di↵erent conditions
for di↵erent emission lines and we only see the cumulative emission from many clouds.
This spectrum is largely determined by radial and density distribution of clouds from
the source. As in Paper II (see also Ferguson et al. 1997; Baldwin et al. 1995), our
model defined the radial and density distributions of gas as f (r) / r and g(n) / nH ,
respectively, where

and

are free parameters. The total line luminosity is then

simply,

Lline /

Z Z

r2 F (r, nH )r nH dnH dr

(4.1)

(Paper II) where F (r, nH ) is the flux of the line and r nH is the spatial distribution
function. For simplicity, we chose the same integration limits as in Paper II, however
we note that this only had small e↵ects on the cumulative observed emission. Other
properties, such as the SED, abundances and gas distributions had far more noticeable
e↵ects. Therefore gas with log(nH ) > 8 was not included due to collisional quenching,
and gas with log(nH ) < 2 was not included due to low optical line intensities.
From here, the predicted emission can constrain the properties of these galaxies
by utilizing line ratio diagnostics in the same manner as VO87, and many other
subsequent papers. Our preliminary diagnostic diagrams will focus on classical line
ratios that have been empirically determined to constrain the excitation mechanism
in emission line galaxies.

161

4.3.1 Spectral Energy Distribution
We used the evolutionary synthesis code Starburst99 to generate synthetic starburst
spectra for our simulations. Several factors influence the radiation field: metallicity,
star formation history (SFH), and stellar population age and evolution. For simplicity,
we chose to adopt the Padova track evolutionary sequence with AGB stars (Bressan
et al. 1993) and Pauldrach / Hillier model atmospheres (Pauldrach et al. 2001; Hillier
& Miller 1998) for all of our SEDs. We also assumed a Kroupa initial mass function
(IMF; Kroupa 2001) with mass intervals of 0.1 M

to 0.5 M

and 0.5 M

to 100

M .
This leaves three additional properties of the SED, SFH, stellar population age
and metallicity. Here we investigate the e↵ects of varying the first two of these. We
used two star formation histories for the incident radiation field of our simulations:
instantaneous and continuous. Instantaneous starbursts assume that a population of
massive stars formed with single burst of star formation while continuous starbursts
continue to evolve until there is balance between newly born stars and the death of
older stars. Our instantaneous starbursts assumed a fixed mass of 106 M , while our
continuous starbursts assumed a star formation rate of 1 M

yr 1 .

Figure 4.3 displays the SEDs for our continuous (Figure 4.3a) and instantaneous
(Figure 4.3b) starbursts as a function of age for solar metallicity. For a stellar population that has only undergone a single burst of star formation at 0 Myr, the hardness
of the SED at high energies begins to decrease as a function of age due to hot OB
stars evolving o↵ the main sequence. Figure 4.3b illustrates this by showing that
after 2 Myr (red line) the intensity of the FUV spectra drops o↵ rapidly. This e↵ect
is less prominent in lower metallicity populations due to less ionizing photons being
absorbed by metals in stellar atmospheres (Snijders et al. 2007). In our preliminary
analysis, we have only considered a solar metallicity cluster of stars for both SFHs.

162

Figure 4.3 The SEDs for continuous (Panel a) and instantaneous (Panel b) star formation histories as a function of age for a cluster of solar metallicity stars. As the
continuous starburst ages, the birth and death of stars eventually comes into equilibrium and the SED ceases to evolve. In contrast, instantaneous starburst undergo
a single period of star formation, as evident by the sharp decrease in intensity past
⇠25 eV for a 10 Myr cluster of stars.

4.3.2 Excitation Mechanism Diagrams
Using the SEDs given above we ran LOC integrations (following the same procedures
as in Paper II) to compare our predicted emission line ratios to the observed values
in our subsets. The VO87 diagnostic diagrams constrain the excitation mechanism
and therefore designate the first order comparison needed to show that MFICA is
correctly isolating SF galaxies. Figure 4.4 shows the predicted emission line ratios
as colored circles connected by a solid line. Each set of circles represents a di↵erent
density distribution in the range -1.8 

 -0.6 in increments of -0.4. The arrow in

the corner indicates the direction of more negative

values. The circles within each

set represent di↵erent radial distributions ranging from -2.0 

 2.0 in increments

of 0.25. The largest circle in each set represents the most negative

value. The white

circle represents our current best fit to the highest ionization subset and is described
below. The observations from the subsets along the central SF locus are also given
163

Figure 4.4 Diagnostic diagrams from VO87 for a continuous SFH as a function of
age. Each quadrant represents a di↵erent starburst age. In each subpanel the squares
connected by the dashed black line are the dereddened measurements from the coadded spectra that constitute the s41-s31-s21-s11-s01 sequence. We do not show the
other SF sequences due to their similarity to the central sequence. The largest shape
in the sequence indicates observations at the top of SF locus (s41). The colored lines
with circles represent LOC integrations. The density-weighting indices in the LOC
integrations, , have more negative values along the direction of the arrow and are
indicated by di↵erent colors. The radial weighting indices, , become more negative
as a function of increasing distance from the largest colored circle for each particular
density weighting. The hollow circle indicates our best fitting set of free parameters
that match the a41 subset (large black square) in the 0.1 Myr starburst. The cross
in the lower right panel of each section indicates the range of acceptable error.

164

Figure)4.4)(cont1d)

165

Figure)4.4)(cont1d)

166

Figure)4.4)(cont1d)

as black squares on these diagrams. We have chosen to only show the central locus
on these diagrams since the wing loci are very similar to the central locus, as can be
seen in Figure 4.2.
Figure 4.4 shows our VO87 diagrams for a starburst with a continuous SFH with
each quadrant showing a di↵erent stellar population age. All emission line ratios
decrease by ⇠0.2 dex between 0.1 Myr and 2 Myr, except for [O II] / [O III], which
increases.
Between these ages, the Balmer emission line equivalent widths, which are relatively flat across the LOC plane, decrease by a smaller amount than the maximum
equivalent widths of the forbidden lines in these diagrams, thus most of the emission
line ratios decrease. However, the [O II] 3727 equivalent widths decrease more slowly
than the [O III] 5007 equivalent widths leading to a larger [O II] / [O III] ratio.
167

Figure 4.5 Diagnostic diagrams from VO87 for a instantaneous SFH as a function of
age, displayed in the same manner as Figure 4.4.

168

Figure)4.5)(cont1d)

169

Figure)4.5)(cont1d)

170

Figure)4.5)(cont1d)

All emission line ratios slightly increase between 2 Myr and 6 Myr except for
[O II] 3727/ [O III] 5007, which slightly decreases, and [N II] 6584/ H↵ which
stays constant. Over this time interval, the Balmer emission line equivalent widths
decrease faster than the maximum equivalent widths of [O III] 5007, [S II] 6720
and [O I] 6300. Here, the [O II] / [O III] ratio decreases due to sharp decrease in
[O II] emission at large radii. After 6 Myr, all emission line ratios stay constant since
the starburst has reached equilibrium (Kewley et al. 2001), as emphasized by our 10
Myr panel in Figure 4.4a.
Figure 4.5 shows the same VO87 diagnostic diagrams as Figure 4.4a but this time
for an instantaneous starburst. Between 0.1 Myr and 2 Myr, all line ratios decrease
by approximately 0.2 dex except for [O III]

5007/ H

and [O II]

3727/ [O III]

5007. In this age range, [O III] 5007/ H , on average, decreases by approximately
171

0.5 dex due to a sharp decrease in [O III] emission. This decrease is larger than that
of [O II] 3727 explaining why the [O II]/[O III] ratio increases between 0.1 Myr and
2 Myr.
Between 2 Myr and 6 Myr, all the emission line ratios in the VO87 diagrams
decrease by ⇠0.2 dex due to a drastic decrease in forbidden line emission over the
entire LOC plane. In contrast, the Balmer line emission also decreases, but by a substantially less amount, and the emission remains flat across the essentially the entire
parameter space. Finally, as Figure 4.3b shows, at 10 Myr the starburst becomes so
soft that only a negligible amount of forbidden line emission is produced. Our current
best fit from the models comes from a continuous starburst (Figure 4.4) with an age
of 0.1 Myr (1st quadrant). The white circle corresponds to

= -0.5 and

= -1.4, and

matches our highest ionization subset to within a factor of 2 for each of the diagrams
shown.

4.4 Future Work
This thesis chapter represents a progress report for an ongoing investigation. The
preliminary conclusions of this work are:
• Currently, our model that best fits the SF subsets results from a 0.1 Myr continuous starburst with a solar metallicity and LOC integration free parameters
= -0.5 and

= -1.4

• This set of free parameters also happens to be very close to our best fit for the
highest ionization AGN subset from Paper II ( = -0.75,

= -1.4)

• We can hypothesize from both instantaneous and continuous starbursts that
the variation in our subsets may not be a function of starburst age, at least for
a dust-free solar metallicity model

172

• We can conclude, at a glance, that our main conclusion from Paper II, the radial
concentration of clouds is the physical parameter responsible for the observed
AGN ionization sequence, also seems to hold for our preliminary results with
SF galaxies
In our preliminary analysis, we assumed that the incident radiation field originated
from a star cluster with either an instantaneous or continuous SFH. We will broaden
this analysis to encompass the SED from 30 Doradus (Selman et al. 1999), which
results from a star cluster that undergoes three discrete periods of star formation and
possesses an IMF slightly flatter than that of the Salpeter IMF (Salpeter 1955).
Here, we have also assumed that the NLR does not increase in metallicity as the
starburst ages, which would be reasonable if the starburst is lighting up pre-existing
gas around it. But an alternate possibility is that the NLR consists of gas produced
by SNe within the starburst, in case its metallicity could evolve as the starburst
ages. Our future work will include assessing the e↵ects of starburst metallicity on
the resulting LOC integrations for both SFHs. Multiple diagnostic diagrams that
constrain the SED will provide insight into the incident continuum and we will include
diagrams with lines as consistency checks on strong line ratios. Starbursts with lower
metallicities have harder ionizing continua, making our observations of He II and Ar
IV emission lines valuable due to their large ionization potential. After constraining
the incident radiation field, we will also constrain the metallicity of the gas clouds
by utilizing metallicity sensitive emission line ratios, in particular the [N II]/ [O II]
ratio. We have standard indicators for electron temperature, [O III] 4363/ [O III]
5007, and low density gas, [S II] 6716/ [S II] 6731, but due to our large galaxy
sample the [S II] 6720/ [S II] 4070 ratio also enables a check on any high density
cloud components.
We will also examine the impact of dust grains by running dusty grids of simulations in the same manner as Paper II, as opposed to the dust-free case studied
173

here. Photoelectric heating from grains will assist in the production of high ionization emission lines. Overall, the incorporation of dust into our models will likely have
around a factor of 2 impact on the equivalent widths (Ferguson et al. 1997) and
may influence the SED optimization (Richardson et al. 2013). Comparisons between
dust-free and dusty models will determine whether a dust-free gas can successfully
recreate the spectra seen in our subsample of starburst galaxies or whether dusty gas,
as in Kewley et al. (2001) and Levesque et al. (2010), is a more accurate description.
We will search for connections between low ionization SF galaxies and low ionization AGN. In Paper II, we fit several low ionization AGN by adjusting the radial
distribution of clouds, however many diagnostic ratios were not adequately fit at low
ionization. Indeed, despite MFICA isolating pure sequences of AGN and SF galaxies, the small degree of cross talk between components could allow enough room for
overlapping AGN and SF properties.
Finally, we will switch to running our simulations with the newest version of
Cloudy (13.02; Ferland et al. 2013). This version features more robust molecular and
chemistry solvers and updated atomic data. In addition, the c10.00 featured an error
that caused the intensities of certain [S II] and [O II] to be incorrectly computed.
This will only have a minor impact on the preliminary results presented here.

174

REFERENCES

175

REFERENCES

Abdo A. A. et al., 2011, Sci, 331, 739
Allen J. T., Hewett P. C., Richardson C. T., Ferland G. J., Baldwin J. A.,
2013, MNRAS, 430, 3510
Allen M. G., Dopita M. A., Tsvetanov Z. I., 1998, ApJ, 493, 571
Allende Prieto C., Lambert D. L., Asplund M., 2001, ApJ, 556, L63
Allende Prieto C., Lambert D. L.; Asplund M., 2002, ApJ, 573, L137
Alonso-Herrero A., Rieke M. J., Rieke G. H., Shields J. C., 2000, ApJ, 532,
845
Asplund M., 2000, A&A, 359, 755
Asplund M., 2003, preprint (astro-ph/0302409)
Asplund M., Nordlund A., Trampedach R., Stein R. F., 2000, A&A, 359,
743
Atoyan A. M., Aharonian F. A., 1996, MNRAS, 278, 525
Baldwin J., Ferland G., Korista K., Verner D., 1995, ApJ, 455, L119
Baldwin J. A., Phillips M. M., Telervich R., 1981, PASP, 93, 5 (BPT)
Baldwin J. A., Ferland G. J., Martin P. G., Corbin M. R., Cota S. A., Peterson B. M., Slettebak A., 1991, ApJ, 374, 580
Baldwin J. A. et al., 1996, ApJ, 468, L115
Baldwin J. A., Hamann F., Korista K. T., Ferland G. J., Dietrich M., Warner C.,
2003, ApJ, 583, 649
Binette L., Wilson A. S., Storchi-Bergmann T., 1996, A&A, 312, 365
Black J. H., van Dishoeck E. F., 1987, ApJ, 322, 412
Bottor↵, M.C., Baldwin, J.A., Ferland, G.J., Ferguson, J.W., Korista, K.T.
2002, ApJ, 581, 932
176

Bradley L. D., Kaiser M. E., Baan W. A., 2004, ApJ, 603, 463
Bressan A., Fagotto F., Bertelli G., Chiosi C., 1993, A&AS, 100, 647
Buehler R. et al., 2012, ApJ, 749, 26
Cazaux S., Tielens, A. G. G. M., 2002, ApJ, 575, 29
Dalgarno A., Yan M., Liu W., 1999, ApJS, 125, 237
Davidson, K., 1978, ApJ, 220, 177
Davidson, K., 1979, ApJ, 228, 179
Davidson K., Fesen R. A., 1985, ARA&A, 23, 119
Dopita M. A., Sutherland R. S., 1995, ApJ, 455, 468
Dopita M. A., Sutherland R. S., 1996, ApJS, 102, 161
Dopita M. A., Groves B. A., Sutherland R. S., Binette L., Cecil G., 2002,
ApJ, 572, 753
Draine, B., 2003, ApJ, 598, 1017
Duyvendak J. J. L., 1942, PASP, 54, 91
Elvis et al., 1994, ApJS, 95, 1
Fabian A. C., Johnstone R. M., Sanders J. S., Conselice C. J., Crawford C.
S., Gallagher J. S. III, Zweibel E., 2008, Nat, 454, 7207
Ferguson J. W., Korista K. T., Baldwin J. A., Ferland G. J., 1997, ApJ,
487, 122 (F97)
Ferland G., 2011, in Colpi M., Gallo L., Grupe D., Komossa S., Leighly K.,
Mathur S., eds, Proc. Sci., Vol. NLS1, Narrow-Line Seyfert 1 Galaxies and Their
Place in the Universe. SISSA, Trieste, p. 013
Ferland, G. J., Osterbrock D. E., 1986, ApJ, 300, 658
Ferland G. J., Fabian A.C., Johnstone R.M., 1994, MNRAS, 266, 399
Ferland G. J., Korista K. T., Verner D. A., Ferguson J. W., Kingdon J. B.,
177

Verner E. M., 1998, PASP, 110, 761
Ferland G. J., Fabian A. C., Hatch N. A., Johnstone R. M., Porter R. L.,
van Hoof P. A. M., Williams R. J. R., 2009, MNRAS, 392,1475
Ferland G. J, et al., 2013, Revista Mexicana de Astronomia y Astrosica, in
press (arXiv:1302.4485)
Fesen R. A., Kirshner R. P., 1982, ApJ, 258, 1
Fesen R., Blair W. P., 1990, ApJ, 351, L45
Fesen R. A., Shull J. M, Hurford, A. P., 1997, AJ, 113, 354
Gallant Y. A., Arons J., 1994, ApJ, 435, 230
Glaccum W., Harper D. A., Loewenstein R. F., Pernic R., Low F. J., 1982,
BAAS, 14, 612
Glassgold A., Langer W., 1974, ApJ, 193, 73
Gomez, H. L. et al., 2012, ApJ, 760, 96
Graham J. R., Wright G. S., Longmore A. J., 1990, ApJ, 352, 172 (G90)
Groves B. A., Dopita, M. A., Sutherland R. S., 2004, ApJS, 153, 9 (G04a)
Groves B. A., Dopita, M. A., Sutherland R. S., 2004, ApJS, 153, 75 (G04b)
Groves B. A., Heckman T. M., Kau↵mann G., 2006, MNRAS, 371, 1559
Hamann F., Ferland G., 1999, ARA&A, 37, 487
Heckmann T. M., 1980, A&A, 87, 142
Hillier D., Miller D. L., 1998, ApJ, 496, 407
Hennessy G. S. et al., 1992, ApJ, 395, L13
Henney W. J., Williams R. J. R., Ferland G. J., Shaw G., O’Dell C. R.,
2007, ApJ, 671, L137
Henry R. B. C., MacAlpine, G. M., 1982, ApJ, 258, 11
Hester J. J., 2008, ARA&A, 46, 127
178

Hester J. J., Graham J. R., Beichman C. A., Gautier T. N. III, 1990, ApJ,
357, 539
Ho L. C., Filippenko A. V., Sargent W. L. W., 1993, ApJ, 417, 63
Holt J., Tadhunter C. N., Morganti R., Emots B. H. C., 2011, MNRAS,
410, 1527
Hoshino M., Arons J., Gallant Y. A., Langdon A. B., 1992, ApJ, 390, 454
Indriolo N., Geballe T. R., Oka T., McCall B. J., 2007, ApJ, 671, 1736
Indriolo N., Fields B. D., McCall B. J., 2009, ApJ, 694, 257
Kau↵mann G., et al., 2003, MNRAS, 346, 1055
Kewley L. J., Dopita M. A., Sutherland R. S., Heisler C. A., Trevena J.,
2001, ApJ, 556, 121
Kewley L. J., Groves B., Kau↵mann G., Heckman T., 2006, MNRAS, 372,
961
Kingdon J., Ferland G. J., Feibelman W. A., 1995, ApJ, 439, 739
Komossa S., Schulz H., 1997, A&A, 323, 31
Korista K., Baldwin J., Ferland G., Verner D., 1997, ApJS, 108, 401
Kroupa, P. 2001, MNRAS, 322, 231
Le Bourlot J., Pineau des Forets G., Roue↵ E., 1995, A&A, 297, 251
Levesque E., Kewley L. J., Larson K., 2010, AJ, 139, 712
Loh E. D., Baldwin J. A., Ferland G. J., 2010, ApJ, 716, L9 (Paper I)
Loh E. D., Baldwin J. A., Curtis Z. K., Ferland G. J., ODell C. R., Fabian
A. C., Salom´ P., 2011, ApJS, 194, 30 (Paper II)
e
Luridiana, V., Sim´n-Diaz, S., Cervi˜o, M., Gonz´lez Delgado, R.M., Porter,
o
n
a
R.L., Ferland, G.J., 2009, ApJ, 691. 1712 Loh E. D., Baldwin J. A., Ferland G. J.,
Curtis Z. K., Richardson C. T., Fabian A. C., Salom´ P., 2012, MNRAS, 421, 789
e
(Paper III)

179

MacAlpine, G. M., Satterfield, T. J., 2008, AJ, 136, 2152
Marsden P. L., Gillett F. C., Jennings R. E., Emerson J. P., de Jong T.,
Olnon F. M., 1984, ApJ, 278, L29
Mathews W. G., Ferland G. J., 1987, ApJ, 323, 456
Mayall N. U., Oort J. H., 1942, PASP, 54, 95
Miller J. S., 1973, ApJ, 180, L83
Mori K., Burrows D.N., Hester J.J., Pavlov G.G., Shibata S., Tsunemi H.
2,004, ApJ, 609, 186
Nagao T., Maiolino R., Marconi A., 2006, A&A, 447, 863
Osterbrock D. E., Ferland G. J., 2006, Astrophysics of Gaseous Nebulae & Active
Galactic Nuclei, 2nd edn. University Science Press, Mill Valley, CA (AGN3)
Papadopoulos, P., 2010, ApJ 720, 226
Pequignot D., Dennefeld M., 1983, A&A, 120, 249
Pellegrini E. W., Baldwin J. A., Brogan C. L. Hanson M. M., Abel N. P.,
Ferland G. J., Nemala H. B., Shaw G., Troland, T. H. 2007, ApJ, 658, 1119
Richardson C. T., Baldwin, J. A. Ferland G. J., Loh E. D., Kuehn C. A.,
Fabian A. C., Salom´ P., 2013, MNRAS, 430, 1257
e
Richardson C. T., Allen J. T., Baldwin J. A., Hewett P. C., Ferland G. J.,
2013, MNRAS (submitted)
Salpeter E. E., 1955, ApJ, 121, 161
Sankrit R. et al., 1998, ApJ, 504, 344
Satterfield, T. J., Katz, A. M., Sibley, A. R., MacAlpine, G. M., Uomoto,
A., 2012, AJ, 144, 27
Selman F., Melnick J., Bosch G., Terlevich R., 1999, A&A, 347, 532
Shaw G., Ferland G. J., Abel N. P., Stancil P. C., van Hoof P. A. M.,
2005, ApJ, 624, 794
Shuder J.M., Osterbrock D.E. 1981, ApJ, 250, 55

180

Snijders L., Kewley L. J., van der Werf P. P., 2007, ApJ, 669, 269
Spinoglio L., Malkan M. A., 1992, ApJ, 399, 504
Spitzer L., 1978, Physical Processes in the Interstellar Medium (New York:
John Wiley & Sons, Inc.)
Stasi´ska G., Cid Fernandes R., Mateus A., Sodr´ L., Asari N. V., 2006,
n
e
MNRAS, 371, 972
Ste↵en A. T., Strateva I., Brandt W. N., Alexander D. M., Koekemoer A.
M., Lehmer B. D., Schneider D. P., Vignali C., 2006, AJ, 131, 2826
Sternberg A., Neufeld, D. A., 1999, ApJ, 516, 371
Storchi-Bergmann T., Pastoriza M. G. 1990, PASP, 102, 1359
Sturm E., Lutz D., Verma A., Netzer H., Sternberg A., Moorwood A. F.
M., Oliva E., Genzel R., 2002, A&A, 393, 821
Tanaka M., 2012, PASJ, 64, 36
Tanaka M., 2012, PASJ, 64, 37
Tavani M. et al., 2011, Sci, 331, 736
Temim T. et al., 2006, AJ, 132, 1610
Temim T., Sonneborn G., Dwek E., Arendt R.G., Gehrz R., Slane P., Roellig T.L., 2012, ApJ, 753, 72
Tielens A. G. G. M., 2005, The Physics and Chemistry of the Interstellar
Medium, (Cambridge, UK: Cambridge University Press)
Tielens A. G. G. M., Hollenbach D., 1985, ApJ, 291, 722
Trimble V., 1968, AJ, 73, 535
Trimble V., 1973, PASP, 85, 579
Veilleux S., Osterbrock D. E., 1987, ApJS, 63, 295 (VO87)
Verner E. M., Verner D. A., Korista K. T., Ferguson J. W., Hamann F.,
Ferland G. J., 1999, ApJS, 120, 101

181

Webber W. R. 1998, ApJ, 506, 329
Wu C. C., 1981, ApJ, 245, 581
Young J. S., Scoville, N.Z., 1991, ARA&A, 29, 581
Zamorani et al., 1981, ApJ, 245, 357
Zhang Z., Liang Y., Hammer F., 2013, preprint (arXiv:1302.3013)

182