Harvard Astronomy 201b

Course Notes

In Uncategorized on April 12, 2011 at 11:26 pm

Stromgren Sphere: An example “chalkboard derivation”

(updated for 2013)

The Stromgren sphere is a simplified analysis of the size of HII regions. Massive O and B stars emit many high-energy photons, which will ionize their surroundings and create HII regions. We assume that such a star is embedded in a uniform medium of neutral hydrogen. A sphere of radius r around this star will become ionized; is called the “Stromgren radius”. The volume of the ionized region will be such that the rate at which ionized hydrogen recombines equals the rate at which the star emits ionizing photons (i.e. all of the ionizing photons are “used up” re-ionizing hydrogen as it recombines)

The recombination rate density is \alpha n^2, where \alpha is the recombination coefficient (in \mathrm{cm}^3~\mathrm{s}^{-1}) and n=n_e=n_\mathrm{H} is the number density (assuming fully ionized gas and only hydrogen, the electron and proton densities are equal). The total rate of ionizing photons (in photons per second) in the volume of the sphere is N^*. Setting the rates of ionization and recombination equal to one another, we get

\frac43 \pi r^3 \alpha n^2 = N^*, and solving for r,

r = ( \frac {3N^*} {4\pi\alpha n^2})^{\frac13}

Typical values for the above variables are N^* \sim 10^{49}~\mathrm{photons~s}^{-1}, \alpha \sim 3\times 10^{-13}\; \mathrm{cm}^3 \; \mathrm s^{-1} and n \sim 10\; \mathrm {cm}^{-3}, implying Stromgren radii of 10 to 100 pc. See the journal club (2013) article for discussion of Stromgren’s seminal 1939 paper.

How do we know there is an ISM?

(updated for 2013)

Early astronomers pointed to 3 lines of evidence for the ISM:

  • Extinction. The ISM obscures the light from background stars. In 1919, Barnard (JC 2011, 2013) called attention to these “dark markings” on the sky, and put forward the (correct) hypothesis that these were the silhouettes of dark clouds. A good rule of thumb for the amount of extinction present is 1 magnitude of extinction per kpc (for typical, mostly unobscured lines-of-sight).
  • Reddening. Even when the ISM doesn’t completely block background starlight, it scatters it. Shorter-wavelength light is preferentially scattered, so stars behind obscuring material appear redder than normal. If a star’s true color is known, its observed color can be used to infer the column density of the ISM between us and the star. Robert Trumpler first used measurements of the apparent “cuspiness” and the brighnesses of star clusters in 1930 to argue for the presence of this effect. Reddening of stars of “known” color is the basis of NICER and related techniques used to map extinction today.
  • Stationary Lines. Spectral observations of binary stars show doppler-shifted lines corresponding to the radial velocity of each star. In addition, some of these spectra exhibit stationary (i.e. not doppler-shifted) absorption lines due to stationary material between us and the binary system. Johannes Hartmann first noticed this in 1904 when investigating the spectrum of \delta Orionis: “The calcium line at \lambda 3934 [angstroms] exhibits a very peculiar behavior. It is distinguished from all the other lines in this spectrum, first by the fact that it always appears extraordinarily week, but almost perfictly sharp… Closer study on this point now led me to the quite surprising result that the calcium line… does not share in the periodic displacements of the lines caused by the orbital motion of the star”

Helpful References: Good discussion of the history of extinction and reddening, from Michael Richmond.

A Sense of Scale

(updated for 2013)

How dense (or not) is the ISM?

  • Dense cores: n \sim 10^5 ~{\rm cm}^{-3}
  • Typical ISM: n \sim 1 ~{\rm cm}^{-3}
  • This room: 1 mol / 22.4L \sim 3 \times 10^{19}~ {\rm cm}^{-3}
  • XVH (eXtremely High Vacuum) — best human-made vacuum: n \sim 3 \times 10^{4}~ {\rm cm}^{-3}
  • Density of stars in the Milky Way: 2.8~{\rm stars/pc}^3 \approx 0.125~M_\odot/{\rm pc}^3 = 8.5 \times 10^{-24} ~{\rm g / cm}^3 \sim 5~{\rm cm}^{-3}

In other words, most of the ISM is at a density far below the densities and pressures we can reproduce in the lab. Thus, the details of most of the microphysics in the ISM are still poorly understood. We also see that the density of stars in the Galaxy is quite small – only a few times the average particle density of the ISM.

See also the interstellar cloud properties table and conversions between angular and linear scale.

Density of the Milky Way’s ISM

(updated for 2013)

How do we know that n \sim 1 ~{\rm cm}^{-3} in the ISM? From the rotation curve of the Milky Way (and some assumptions about the mass ratio of gas to gas+stars+dark matter), we can infer

M_{\rm gas} = 6.7 \times 10^{9} M_\odot

Maps of HI and CO reveal the extent of our galaxy to be

D = 40 kpc

h = 140 pc (scale height of HI)

This applies an approximate volume of

V = \pi D^2 h / 4 = 5 \times 10^{66} ~{\rm cm}^{3}

Which, yields a density of

\rho = 2.5 \times 10^{-24} ~{\rm g cm}^{-3}

Density of the Intergalactic Medium

(updated for 2013)

From cosmology observations, we know the universe to be very nearly flat (\Omega = 1). This implies that the mean density of the universe is \rho = \rho_{\rm crit} = \frac{3 H_0^2}{8 \pi G} = 7 \times 10^{-30} ~{\rm g~ cm}^{-3} \Rightarrow n<4.3 \times 10^{-6}~{\rm cm}^{-3}.

This places an upper limit on the density of the Intergalactic Medium.

Composition of the ISM

(updated for 2013)

  • Gas: by mass, gas is 60% Hydrogen, 30% Helium. By number, gas is 88% H, 10% He, and 2% heavier elements
  • Dust: The term “dust” applies roughly to any molecule too big to name. The size distribution is biased towards small (0.2 \mum) particles, with an approximate distribution N(a) \propto a^{-3.5}. The density of dust in the galaxy is \rho_{\rm dust} \sim .002 M_\odot ~{\rm pc}^{-3} \sim 0.1 \rho_{\rm gas}
  • Cosmic Rays: Charged, high-energy (anti)protons, nuclei, electrons, and positrons. Cosmic rays have an energy density of 0.5 ~{\rm eV ~ cm}^{-3}. The equivalent mass density (using E = mc^2) is 9 \times 10^{-34}~{\rm g cm}^{-3}
  • Magnetic Fields: Typical field strengths in the MW are 1 \mu G \sim 0.2 ~{eV ~cm}^{-3}. This is strong enough to confine cosmic rays.

Bruce Draine’s List of constituents in the ISM:

(updated for 2013)

  1. Gas
  2. Dust
  3. Cosmic Rays*
  4. Photons**
  5. B-Field
  6. Gravitational Field
  7. Dark Matter

*cosmic rays are highly relativistic, super-energetic ions and electrons

**photons include:

  • The Cosmic Microwave Background (2.7 K)
  • starlight from stellar photospheres (UV, optical, NIR,…)
  • h\nu from transitions in atoms, ions, and molecules
  • “thermal emission” from dust (heated by starlight, AGN)
  • free-free emission (bremsstrahlung) in plasma
  • synchrotron radiation from relativistic electrons
  • \gamma-rays from nuclear transitions

His list of “phases” from Table 1.3:

  1. Coronal gas (Hot Ionized Medium, or “HIM”): T> 10^{5.5}~{\rm K}. Shock-heated from supernovae. Fills half the volume of the galaxy, and cools in about 1 Myr.
  2. HII gas: Ionized mostly by O and early B stars. Called an “HII region” when confined by a molecular cloud, otherwise called “diffuse HII”.
  3. Warm HI (Warm Neutral Medium, or “WNM”): atomic, T \sim 10^{3.7}~{\rm K}. n\sim 0.6 ~{\rm cm}^{-3}. Heated by starlight, photoelectric effect, and cosmic rays. Fills ~40% of the volume.
  4. Cool HI (Cold Neutral Medium, or “CNM”). T \sim 100~{\rm K}, n \sim 30 ~{\rm cm}^{-3}. Fills ~1% of the volume.
  5. Diffuse molecular gas. Where HI self-shields from UV radiation to allow H_2 formation on the surfaces of dust grains in cloud interiors. This occurs at 10~{\rm to}~50~{\rm cm}^{-3}.
  6. Dense Molecular gas. “Bound” according to Draine (though maybe not). n >\sim 10^3 ~{\rm cm}^{-3}. Sites of star formation.  See also Bok Globules (JC 2013).
  7. Stellar Outflows. T=50-1000 {\rm K}, n \sim 1-10^6 ~{\rm cm}^{-3}. Winds from cool stars.

These phases are fluid and dynamic, and change on a variety of time and spatial scales. Examples include growth of an HII region, evaporation of molecular clouds, the interface between the ISM and IGM, cooling of supernova remnants, mixing, recombination, etc.

Topology of the ISM

(updated for 2013)

A grab-bag of properties of the Milky Way

  • HII scale height: 1 kpc
  • CO scale height: 50-75 pc
  • HI scale height: 130-400 pc
  • Stellar scale height: 100 pc in spiral arm, 500 pc in disk
  • Stellar mass: 5 \times 10^{10} M_\odot
  • Dark matter mass: 5 \times 10^{10} M_\odot
  • HI mass: 2.9 \times 10^9 M_\odot
  • H2 mass (inferred from CO): 0.84 \times 10^9 M_\odot
  • HII mass: 1.12 \times 10^9~M_\odot
  • -> total gas mass = 6.7 \times 10^9~M_\odot (including He).
  • Total MW mass within 15 kpc: 10^{11} M_\odot (using the Galaxy’s rotation curve). About 50% dark matter.

So the ISM is a relatively small constituent of the Galaxy (by mass).

The Sound Speed

(updated for 2013)

The speed of sound is the speed at which pressure disturbances travel in a medium. It is defined as

c_s \equiv \frac{\partial P}{\partial \rho} ,

where P and \rho are pressure and mass density, respectively. For a polytropic gas, i.e. one defined by the equation of state P \propto \rho^\gamma, this becomes c_s=\sqrt{\gamma P/\rho}. \gamma is the adiabatic index (ratio of specific heats), and \gamma=5/3 describes a monatomic gas.

For an isothermal gas where the ideal gas equation of state P=\rho k_B T / (\mu m_{\rm H}) holds, c_s=\sqrt{k_B T/\mu}. Here, \mu is the mean molecular weight (a factor that accounts for the chemical composition of the gas), and m_{\rm H} is the hydrogen atomic mass. Note that for pure molecular hydrogen \mu=2. For molecular gas with ~10% He by mass and trace metals, \mu \approx 2.3 is often used.

A gas can be approximated to be isothermal if the sound wave period is much higher than the (radiative) cooling time of the gas, as any increase in temperature due to compression by the wave will be immediately followed by radiative cooling to the original equilibrium temperature well before the next compression occurs. Many astrophysical situations in the ISM are close to being isothermal, thus the isothermal sound speed is often used. For example, in conditions where temperature and density are independent such as H II regions (where the gas temperature is set by the ionizing star’s spectrum), the gas is very close to isothermal.

Hydrogen “Slang”

(updated for 2013)

Lyman limit: the minimum energy needed to remove an electron from a Hydrogen atom. A “Lyman limit photon” is a photon with at least this energy.

E = 13.6~{\rm eV} = 1~{\rm `Rydberg'} = hcR_{\rm H} ,

where R_{\rm H}=1.097 \times 10^{7} {\rm m}^{-1} is the Rydberg constant, which has units of 1/\lambda. This energy corresponds to the Lyman limit wavelength as follows:

E = h\nu = hc/\lambda \Rightarrow \lambda=912~{\rm \AA} .

Lyman series: transitions to and from the n=1 energy level of the Bohr atom. The first line in this series was discovered in 1906 using UV studies of electrically excited hydrogen gas.

Balmer series: transitions to and from the n=2 energy level. Discovered in 1885; since these are optical transitions, they were more easily observed than the UV Lyman series transitions.

There are also other named series corresponding to higher n. Examples include Paschen (n=3), Brackett (n=4), and Pfund (n=5). The lowest energy (longest wavelength) transition of a series is designated \alpha, the next lowest energy is \beta, and so on. For example, the transition from n=2 to 1 is Lyman alpha, or {\rm Ly}\alpha, while the transition from n=7 to 4 is Brackett gamma, or {\rm Br}\gamma. The wavelength of a given transition can be computed via the Rydberg equation:

\frac{1}{\lambda}=R_{\rm H} \big|\frac{1}{n_f^2}-\frac{1}{n_i^2}\big| ,

where n_i and n_f are the initial and final energy levels of the electron, respectively. See this handout for a pictorial representation of the low n transitions in hydrogen. Note that the Lyman (or Balmer, Paschen, etc.) limit can be computed by inserting n_i=\infty in the above equation.

The Lyman continuum corresponds to the region of the spectrum near the Lyman limit, where the spacing between energy levels becomes comparable to spectral line widths and so individual lines are no longer distinguishable. Such continua exist for each series of lines.


(updated for 2013)

See Draine Table 1.4 for elemental abundances for the Sun (and thus presumably for the ISM near the Sun).

By number: {\rm H:He:C} = 1:0.1:3 \times 10^{-4} ;

by mass: {\rm H:He:C} = 1:0.4:3.5 \times 10^{-3} .

However, these ratios vary by position in the galaxy, especially for heavier elements (which depend on stellar processing). For example, the abundance of heavy elements (Z ≥ 6, i.e. carbon and heavier) is twice as low at the sun’s position than in the Galactic center. Even though metals account for only 1% of the mass, they dominate most of the important chemistry, ionization, and heating/cooling processes. They are essential for star formation, as they allow molecular clouds to cool and collapse.

Dissociating molecules takes less energy than ionizing atoms, in general. For example:

E_{I,{\rm H}}=13.6~{\rm eV}

E_{D,{\rm H}_2}=4.52~{\rm eV} \Rightarrow \lambda=2743~{\rm \AA} (UV transition)

E_{D,{\rm CO}}=11.2~{\rm eV},

where E_I and E_D are the ionization and dissociation energies, respectively. We can see that it is much easier to dissociate molecular hydrogen than to ionize atomic hydrogen; in other words, atomic H will survive a harsher radiation field than molecular H. The above numbers thus set the structure of molecular clouds in the interstellar radiation field; a large amount of molecular gas needs to gather together in order to allow it to survive via the process of self-shielding, in which a thick enough column of gas exists such that at some distance below the surface of the cloud all of the energetic photons have already been absorbed. Note that the high dissociation energy of CO is a result of the triple bond between the carbon and oxygen atoms. CO is a very important coolant in molecular clouds.

Measuring States in the ISM

(updated for 2013)

There are two primary observational diagnostics of the thermal, chemical, and ionization states in the ISM:

  1. Spectral Energy Distribution (SED; broadband low-resolution)
  2. Spectrum (narrowband, high-resolution)


Very generally, if a source’s SED is blackbody-like, one can fit a Planck function to the SED and derive the temperature and column density (if one can assume LTE). If an SED is not blackbody-like, the emission is the sum of various processes, including:

  • thermal emission (e.g. dust, CMB)
  • synchrotron emission (power law spectrum)
  • free-free emission (thermal for a thermal electron distribution)


Quantum mechanics combined with chemistry can predict line strengths. Ratios of lines can be used to model “excitation”, i.e. what physical conditions (density, temperature, radiation field, ionization fraction, etc.) lead to the observed distribution of line strengths. Excitation is controlled by

  • collisions between particles (LTE often assumed, but not always true)
  • photons from the interstellar radiation field, nearby stars, shocks, CMB, chemistry, cosmic rays
  • recombination/ionization/dissociation

Which of these processes matter where? In class (2011), we drew the following schematic.

A schematic of several structures in the ISM


A: Dense molecular cloud with stars forming within

  • T=10-50~{\rm K};~n>10^3~{\rm cm}^{-3} (measured, e.g., from line ratios)
  • gas is mostly molecular (low T, high n, self-shielding from UV photons, few shocks)
  • not much photoionization due to high extinction (but could be complicated ionization structure due to patchy extinction)
  • cosmic rays can penetrate, leading to fractional ionization: X_I=n_i/(n_H+n_i) \approx n_i/n_H \propto n_H^{-1/2}, where n_i is the ion density (see Draine 16.5 for details). Measured values for X_e (the electron-to-neutral ratio, which is presumed equal to the ionization fraction) are about X_e \sim 10^{-6}~{\rm to}~10^{-7}.
  • possible shocks due to impinging HII region – could raise T, n, ionization, and change chemistry globally
  • shocks due to embedded young stars w/ outflows and winds -> local changes in Tn, ionization, chemistry
  • time evolution? feedback from stars formed within?

B: Cluster of OB stars (an HII region ionized by their integrated radiation)

  • 7000 < T < 10,000 K (from line ratios)
  • gas primarily ionized due to photons beyond Lyman limit (E > 13.6 eV) produced by O stars
  • elements other than H have different ionization energy, so will ionize more or less easily
  • HII regions are often clumpy; this is observed as a deficit in the average value of n_e from continuum radiation over the entire region as compared to the value of ne derived from line ratios. In other words, certain regions are denser (in ionized gas) than others.
  • The above introduces the idea of a filling factor, defined as the ratio of filled volume to total volume (in this case the filled volume is that of ionized gas)
  • dust is present in HII regions (as evidenced by observations of scattered light), though the smaller grains may be destroyed
  • significant radio emission: free-free (bremsstrahlung), synchrotron, and recombination line (e.g. H76a)
  • chemistry is highly dependent on nT, flux, and time

C: Supernova remnant

  • gas can be ionized in shocks by collisions (high velocities required to produce high energy collisions, high T)
  • e.g. if v > 1000 km/s, T > 106 K
  • atom-electron collisions will ionize H, He; produce x-rays; produce highly ionized heavy elements
  • gas can also be excited (e.g. vibrational H2 emission) and dissociated by shocks

D: General diffuse ISM

  • UV radiation from the interstellar radiation field produces ionization
  • ne best measured from pulsar dispersion measure (DM), an observable. {\rm DM} \propto \int n_e dl
  • role of magnetic fields depends critically on XI(B-fields do not directly affect neutrals, though their effects can be felt through ion-neutral collisions)

Energy Density Comparison

(updated for 2013)

See Draine table 1.5. The primary sources of energy present in the ISM are:

      1. The CMB (T_{\rm CMB}=2.725~{\rm K}
      2. Thermal IR from dust
      3. Starlight (h\nu < 13.6 {\rm eV}
      4. Thermal kinetic energy (3/2 nkT)
      5. Turbulent kinetic energy (1/2 \rho \sigma_v^2)
      6. Magnetic fields (B^2 / 8 \pi )
      7. Cosmic rays

All of these terms have energy densities within an order of magnitude of 1 ~{\rm eV ~ cm}^{-3}. With the exception of the CMB, this is not a coincidence: because of the dynamic nature of the ISM, these processes are coupled together and thus exchange energy with one another.

Relevant Velocities in the ISM

(updated for 2013)

Note: it’s handy to remember that 1 km/s ~ 1 pc / Myr.

  • Galactic rotation: 18 km/s/kpc (e.g. 180 km/s at 10 kpc)
  • Isothermal sound speed: c_s =\sqrt{\frac{kT}{\mu}}
    • For H, this speed is 0.3, 1, and 3 km/s at 10 K, 100 K, and 1000 K, respectively.
  • Alfvén speed: The speed at which magnetic fluctuations propagate. v_A = B / \sqrt{4 \pi \rho} Alfvén waves are transverse waves along the direction of the magnetic field.
    • Note that v_A = {\rm const} if B \propto \rho^{1/2}, which is observed to be true over a large portion of the ISM.
    • Interstellar B-fields can be measured using the Zeeman effect. Observed values range from 5~\mu {\rm G} in the diffuse ISM to 1 mG in dense clouds. For specific conditions:
      • B = 1~\mu{\rm G}, n = 1 ~{\rm cm}^{-3} \Rightarrow v_A = 2~{\rm km~s}^{-1}
      • B = 30~\mu {\rm G}, n = 10^4~{\rm cm}^{-3} \Rightarrow v_A = 0.4~{\rm km~s}^{-1}
      • B = 1~{\rm mG}, n = 10^7 {\rm cm}^{-3} \Rightarrow v_A = 0.5~{\rm km~s}^{-1}
    • Compare to the isothermal sound speed, which is 0.3 km/s in dense gas at 20 K.
      • c_s \approx v_A in dense gas
      • c_s < v_A in diffuse gas
  • Observed velocity dispersion in molecular gas is typically about 1 km/s, and is thus supersonic. This is a signature of the presence of turbulence. (see the summary of Larson’s seminal 1981 paper)

Introductory remarks on Radiative Processes and Equilibrium

(updated for 2013)

The goal of the next several sections is to build an understanding of how photons are produced by, are absorbed by, and interact with the ISM. We consider a system in which one or more constituents are excited under certain physical conditions to produce photons, then the photons pass through other constituents under other conditions, before finally being observed (and thus affected by the limitations and biases of the observational conditions and instruments) on Earth. Local thermodynamic equilibrium is often used to describe the conditions, but this does not always hold. Remember that our overall goal is to turn observations of the ISM into physics, and vice-versa.

The following contribute to an observed Spectral Energy Distribution:

      • gas: spontaneous emission, stimulated emission (e.g. masers), absorption, scattering processes involving photons + electrons or bound atoms/molecules
      • dust: absorption; scattering (the sum of these two -> extinction); emission (blackbody modified by wavelength-dependent emissivity)
      • other: synchrotron, brehmsstrahlung, etc.

The processes taking place in our “system” depend sensitively on the specific conditions of the ISM in question, but the following “rules of thumb” are worth remembering:

      1. Very rarely is a system actually in a true equilibrium state.
      2. Except in HII regions, transitions in the ISM are usually not electronic.
      3. The terms Upper Level and Lower Level refer to any two quantum mechanical states of an atom or molecule where E_{\rm upper}>E_{\rm lower}. We will use k to index the upper state, and j for the lower state.
      4. Transitions can be induced by photons, cosmic rays, collisions with atoms and molecules, and interactions with free electrons.
      5. Levels can refer to electronic, rotational, vibrational, spin, and magnetic states.
      6. To understand radiative processes in the ISM, we will generally need to know the chemical composition, ambient radiation field, and velocity distribution of each ISM component. We will almost always have to make simplifying assumptions about these conditions.

Thermodynamic Equilibrium

(updated for 2013)

Collisions and radiation generally compete to establish the relative populations of different energy states. Randomized collisional processes push the distribution of energy states to the Boltzmann distribution, n_j \propto e^{-E_j / kT}. When collisions dominate over competing processes and establish the Boltzmann distribution, we say the ISM is in Thermodynamic Equilibrium.

Often this only holds locally, hence the term Local Thermodynamic Equilibrium or LTE. For example, the fact that we can observe stars implies that energy (via photons) is escaping the system. While this cannot be considered a state of global thermodynamic equilibrium, localized regions in stellar interiors are in near-equilibrium with their surroundings.

But the ISM is not like stars. In stars, most emission, absorption, scattering, and collision processes occur on timescales very short compared with dynamical or evolutionary timescales. Due to the low density of the ISM, interactions are much more rare. This makes it difficult to establish equilibrium. Furthermore, many additional processes disrupt equilibrium (such as energy input from hot stars, cosmic rays, X-ray background, shocks).

As a consequence, in the ISM the level populations in atoms and molecules are not always in their equilibrium distribution. Because of the low density, most photons are created from (rare) collisional processes (except in locations like HII regions where ionization and recombination become dominant).

Spitzer Notation

(updated for 2013)

We will use the notation from Spitzer (1978). See also Draine, Ch. 3. We represent the density of a state j as

n_j(X^{(r)}), where

      • n: particle density
      • j: quantum state
      • X: element
      • (r): ionization state
      • For example, HI = H^{(0)}

In his book, Spitzer defines something called “Equivalent Thermodynamic Equilibrium” or “ETE”. In ETE, n_j^* gives the “equivalent” density in state j. The true (observed) value is n_j. He then defines the ratio of the true density to the ETE density to be

b_j = n_j / n_j^*.

This quantity approaches 1 when collisions dominate over ionization and recombination. For LTE, b_j = 1 for all levels. The level population is then given by the Boltzmann equation:

\frac{n_j^\star(X^{(r)})}{n_k^\star(X^{(r)})} = (\frac{g_{rj}}{g_{rk}})~e^{ -(E_{rj} - E_{rk}) / kT },

where E_{rj} and g_{rj} are the energy and statistical weight (degeneracy) of level j, ionization state r. The exponential term is called the “Boltzmann factor”‘ and determines the relative probability for a state.

The term “Maxwellian” describes the velocity distribution of a 3-D gas. “Maxwell-Boltzmann” is a special case of the Boltzmann distribution for velocities.

Using our definition of b and dropping the “r” designation,

\frac{n_k}{n_j} = \frac{b_k}{b_j} (\frac{g_k}{g_j})~e^{-h \nu_{jk} / kT }

Where \nu_{jk} is the frequency of the radiative transition from k to j. We will use the convention that E_k > E_j, such that E_{jk}=h\nu_{jk} > 0.

To find the fraction of atoms of species X^{(r)} excited to level j, define:

\sum_k n_k^\star (X^{(r)}) = n^\star(X^{(r)})

as the particle density of X^{(r)} in all states. Then

\frac{ n_j^* (X^{(r)}) } { n^* (X^{(r)})} = \frac{ g_{rj} e^{-E_{rj} / kT} } {\sum_k g_{rk} e^{ -E_{rk} / kT} }

Define f_r, the “partition function” for species X^{(r)}, to be the denominator of the RHS of the above equation. Then we can write, more simply:

\frac{n_j^\star}{n^\star} = \frac{g_{rj}}{f_r} e^{-E_{rj}/kT}

to be the fraction of particles that are in state j. By computing this for all j we now know the distribution of level populations for ETE.

The Saha Equation

(updated for 2013)

How do we deal with the distribution over different states of ionization r? In thermodynamic equilibrium, the Saha equation gives:

\frac{ n^\star(X^{(r+1)}) n_e } { n^\star (X^{(r)}) } = \frac{ f_{r+1} f_e}{f_r},

where f_r and f_{r+1} are the partition functions as discussed in the previous section. The partition function for electrons is given by

f_e = 2\big( \frac{2 \pi m_e k T} {h^2} \big) ^{3/2} = 4.829 \times 10^{15} (\frac{T}{K})^{3/2}

For a derivation of this, see pages 103-104 of this handout from Bowers and Deeming.

If f_r and f_{r+1} are approximated by the first terms in their sums (i.e. if the ground state dominates their level populations), then

\frac{ n^\star ( X^{ (r+1) } ) n_e } {n^\star ( X^{ (r) } ) } = 2 \big(\frac{ g_{r+1,1} }{g_{ r,1}}\big) \big( \frac{ 2 \pi m_e k T} {h^2} \big)^{3/2} e^{-\Phi_r / kT},

where \Phi_r=E_{r+1,1}-E_{r,1} is the energy required to ionize X^{(r)} from the ground (j = 1)  level. Ultimately, this is just a function of n_e and T. This assumes that the only relevant ionization process is via thermal collision (i.e. shocks, strong ionizing sources, etc. are ignored).

Important Properties of Local Thermodynamic Equilibrium

(updated for 2013)

For actual local thermodynamic equilbrium (not ETE), the following are important to keep in mind:

      • Detailed balance: transition rate from j to k = rate from k to j (i.e. no net change in particle distribution)
      • LTE is equivalent to ETE when b_j = 1 or \frac{b_j}{b_k} = 1
      • LTE is only an approximation, good under specific conditions.
      • Radiation intensity produced is not blackbody illumination as you’d want for true thermodynamic equilibrium.
      • Radiation is usually much weaker than the Planck function, which means not all levels are populated.
      • LTE assumption does not mean the Saha equation is applicable since radiative processes (not collisions) dominate in many ISM cases where LTE is applicable.

Definitions of Temperature

(updated for 2013)

The term “temperature” describes several different quantities in the ISM, and in observational astronomy. Only under idealized conditions (i.e. thermodynamic equilibrium, the Rayleigh Jeans regime, etc.) are (some of) these temperatures equivalent. For example, in stellar interiors, where the plasma is very well-coupled, a single “temperature” defines each of the following: the velocity distribution, the ionization distribution, the spectrum, and the level populations. In the ISM each of these can be characterized by a different “temperature!”

Brightness Temperature

T_B = the temperature of a blackbody that reproduces a given flux density at a specific frequency, such that

B_\nu(T_B) = \frac{2 h \nu^3}{c^2} \frac{1}{{\rm exp}(h \nu / kT_B) - 1}

Note: units for B_{\nu} are {\rm erg~cm^{-2}~s^{-1}~Hz^{-1}~ster^{-1}}.

This is a fundamental concept in radio astronomy. Note that the above definition assumes that the index of refraction in the medium is exactly 1.

Effective Temperature

T_{\rm eff} (also called T_{\rm rad}, the radiation temperature) is defined by

\int_\nu B_\nu d\nu = \sigma T_{{\rm eff}}^4 ,

which is the integrated intensity of a blackbody of temperature T_{\rm eff}. \sigma = (2 \pi^5 k^4)/(15 c^2 h^3)=5.669 \times 10^{-5} {\rm erg~cm^{-2}~s^{-1}~K^{-4}} is the Stefan-Boltzmann constant.

Color Temperature

T_c is defined by the slope (in log-log space) of an SED. Thus T_c is the temperature of a blackbody that has the same ratio of fluxes at two wavelengths as a given measurement. Note that T_c = T_b = T_{\rm eff} for a perfect blackbody.

Kinetic Temperature

T_k is the temperature that a particle of gas would have if its Maxwell-Boltzmann velocity distribution reproduced the width of a given line profile. It characterizes the random velocity of particles. For a purely thermal gas, the line profile is given by

I(\nu) = I_0~e^{\frac{-(\nu-\nu_{jk})^2}{2\sigma^2}},

where \sigma_{\nu}=\frac{\nu_{jk}}{c}\sqrt{\frac{kT_k}{\mu}} in frequency units, or

\sigma_v=\sqrt{\frac{kT_k}{\mu}} in velocity units.

In the “hot” ISM T_k is characteristic, but when \Delta v_{\rm non-thermal} > \Delta v_{\rm thermal} (where \Delta v are the Doppler full widths at half-maxima [FWHM]) then T_k does not represent the random velocity distribution. Examples include regions dominated by turbulence.

T_k can be different for neutrals, ions, and electrons because each can have a different Maxwellian distribution. For electrons, T_k = T_e, the electron temperature.

Ionization Temperature

T_I is the temperature which, when plugged into the Saha equation, gives the observed ratio of ionization states.

Excitation Temperature

T_{\rm ex} is the temperature which, when plugged into the Boltzmann distribution, gives the observed ratio of two energy states. Thus it is defined by

\frac{n_k}{n_j}=\frac{g_k}{g_j}~e^{-h\nu_{jk}/kT_{\rm ex}}.

Note that in stellar interiors, T_k = T_I = T_{\rm ex} = T_c. In this room, T_k = T_I = T_{\rm ex} \sim 300K, but T_c \sim 6000K.

Spin Temperature

T_s is a special case of T_{\rm ex} for spin-flip transitions. We’ll return to this when we discuss the important 21-cm line of neutral hydrogen.

Bolometric temperature

T_{\rm bol} is the temperature of a blackbody having the same mean frequency as the observed continuum spectrum. For a blackbody, T_{\rm bol} = T_{\rm eff}. This is a useful quantity for young stellar objects (YSOs), which are often heavily obscured in the optical and have infrared excesses due to the presence of a circumstellar disk.

Antenna temperature

T_A is a directly measured quantity (commonly used in radio astronomy) that incorporates radiative transfer and possible losses between the source emitting the radiation and the detector. In the simplest case,

T_A = \eta T_B( 1 - e^{-\tau}),

where \eta is the telescope efficiency (a numerical factor from 0 to 1) and \tau is the optical depth.

Excitation Processes: Collisions

(updated for 2013)

Collisional coupling means that the gas can be treated in the fluid approximation, i.e. we can treat the system on a macrophysical level.

Collisions are of key importance in the ISM:

      • cause most of the excitation
      • can cause recombinations (electron + ion)
      • lead to chemical reactions

Three types of collisions

      1. Coulomb force-dominated (r^{-1} potential): electron-ion, electron-electron, ion-ion
      2. Ion-neutral: induced dipole in neutral atom leads to r^{-4} potential; e.g. electron-neutral scattering
      3. neutral-neutral: van der Waals forces -> r^{-6} potential; very low cross-section

We will discuss (3) and (2) below; for ion-electron and ion-ion collisions, see Draine Ch. 2.

In general, we will parametrize the interaction rate between two bodies A and B as follows:

{\frac{\rm{reaction~rate}}{\rm{volume}}} = <\sigma v>_{AB} n_a n_B

In this equation, <\sigma v>_{AB} is the collision rate coefficient in \rm{cm}^3 \rm{s}^{-1}. <\sigma v>_{AB}= \int_0^\infty \sigma_{AB}(v) f_v~dv, where \sigma_{AB} (v) is the velocity-dependent cross section and f_v~dv is the particle velocity distribution, i.e. the probability that the relative speed between A and B is v. For the Maxwellian velocity distribution,

f_v~dv = 4 \pi \left(\frac{\mu'}{2\pi k T}\right)^{3/2} e^{-\mu' v^2/2kT} v^2~dv,

where \mu'=m_A m_B/(m_A+m_B) is the reduced mass. The center of mass energy is E=1/2 \mu' v^2, and the distribution can just as well be written in terms of the energy distribution of particles, f_E dE. Since f_E dE = f_v dv, we can rewrite the collision rate coefficient in terms of energy as

\sigma_{AB}=\left(\frac{8kT}{\pi\mu'}\right)^{1/2} \int_0^\infty \sigma_{AB}(E) \left(\frac{E}{kT}\right) e^{-E/kT} \frac{dE}{kT}.

These collision coefficients can occasionally be calculated analytically (via classical or quantum mechanics), and can in other situations be measured in the lab. The collision coefficients often depend on temperature. For practical purposes, many databases tabulate collision rates for different molecules and temperatures (e.g., the LAMBDA databsase).

For more details, see Draine, Chapter 2. In particular, he discusses 3-body collisions relevant at high densities.

Neutral-Neutral Interactions

(updated for 2013)

Short range forces involving “neutral” particles (neutral-ion, neutral-neutral) are inherently quantum-mechanical. Neutral-neutral interactions are very weak until electron clouds overlap (\sim 1 \AA\sim 10^{-8}cm). We can therefore treat these particles as hard spheres. The collisional cross section for two species is a circle of radius r1 + r2, since that is the closest two particles can get without touching.

\sigma_{nn} \sim \pi (r_1 + r_2)^2 \sim 10^{-15}~{\rm cm}^2

What does that collision rate imply? Consider the mean free path:

mfp = \ell_c \approx (n_n \sigma_{nn})^{-1} = \frac{10^{15}} {n_H}~{\rm cm}

This is about 100 AU in typical ISM conditions (n_H = 1 {\rm cm^{-3}})

In gas at temperature T, the mean particle velocity is given by the 3-d kinetic energy: 3/2 m_n v^2 = kT, or

v = \sqrt{\frac{2}{3} \frac{kT}{m_n}}, where m_n is the mass of the neutral particle. The mean free path and velocity allows us to define a collision timescale:

\tau_{nn} \sim \frac{l_c}{v} \sim (\frac{2}{3} \frac{kT}{m_n})^{-1/2} (n_n \sigma_{nn})^{-1} = 4.5 \times 10^3~n_n^{-1}~T^{-1/2}~{\rm years}.

      • For (n,T) = (1~{\rm cm^{-3}, 80~K}), the collision time is 500 years
      • For (n,T) = (10^4~{\rm cm^{-3}, 10~K}), the collision time is 1.7 months
      • For (n,T) = (1~{\rm cm^{-3}, 10^4~K}), the collision time is 45 years

So we see that density matters much more than temperature in determining the frequency of neutral-neutral collisions.

Ion-Neutral Reactions

(updated for 2013)

In Ion-Neutral reactions, the neutral atom is polarized by the electric field of the ion, so that interaction potential is

U(r) \approx \vec{E} \cdot \vec{p} = \frac{Z e} {r^2} ( \alpha \frac{Z e}{r^2} ) = \alpha \frac{Z^2 e^2}{r^4},

where \vec{E} is the electric field due to the charged particle, \vec{p} is the induced dipole moment in the neutral particle (determined by quantum mechanics), and \alpha is the polarizability, which defines \vec{p}=\alpha \vec{E} for a neutral atom in a uniform static electric field. See Draine, section 2.4 for more details.

This interaction can take strong or weak forms. We distinguish between the two cases by considering b, the impact parameter. Recall that the reduced mass of a 2-body system is \mu' = m_1 m_2 / (m_1 + m_2) In the weak regime, the interaction energy is much smaller than the kinetic energy of the reduced mass:

\frac{\alpha Z^2 e^2}{b^4} \ll\frac{\mu' v^2}{2} .

In the strong regime, the opposite holds:

\frac{\alpha Z^2 e^2}{b^4} \gg\frac{\mu' v^2}{2}.

The spatial scale which separates these two regimes corresponds to b_{\rm crit}, the critical impact parameter. Setting the two sides equal, we see that b_{\rm crit} = \big(\frac{2 \alpha Z^2 e^2}{\mu' v^2}\big)^{1/4}

The effective cross section for ion-neutral interactions is

\sigma_{ni} \approx \pi b_{\rm crit}^2 = \pi Z e (\frac{2 \alpha}{\mu'})^{1/2} (\frac{1}{v})

Deriving an interaction rate is tricker than for neutral-neutral collisions because n_i \ne n_n in general. So, let’s leave out an explicit n and calculate a rate coefficient instead, in {\rm cm}^3 {\rm s}^{-1}.

k = <\sigma_{ni} v> (although really \sigma_{ni} \propto 1/v, so k is largely independent of v). Combining with the equation above, we get the ion-neutral scattering rate coefficient

k = \pi Z e (\frac{2 \alpha}{\mu'})^{1/2}

As an example, for C^+ - H interactions we get k \approx 2 \times 10^{-9} {\rm cm^{3} s^{-1}}. This is about the rate for most ion-neutral exothermic reactions. This gives us

\frac{{\rm rate}}{{\rm volume}} = n_i n_n k.

So, if n_i = n_n = 1, the average time \tau between collisions is 16 years. Recall that, for neutral-neutral collisions in the diffuse ISM, we had \tau \sim 500 years. Ion-neutral collisions are much more frequent in most parts of the ISM due to the larger interaction cross section.

The Virial Theorem

(Transcribed by Bence Beky). See also these slides from lecture

See Draine pp 395-396 and appendix J for more details.

The Virial Theorem provides insight about how a volume of gas subject to many forces will evolve. Lets start with virial equilibrium. For a surface S,

0 = \frac12 \frac{\mathrm D^2I}{\mathrm Dt^2} = 2\Gamma + 3\Pi + \mathscr M + W + \frac1{4\pi}\int_S(\mathbf r \cdot \mathbf B_ \mathbf B \cdot \mathrm d \mathbf s - \int_S \left(p+\frac{B^2}{8\pi}\right)\mathbf r \cdot \mathrm d\mathbf s,

see Spitzer pp.~217–218. Here I is the moment of inertia:

I = \int \varrho r^2 \mathrm dV

\Gamma is the bulk kinetic energy of the fluid (macroscopic kinetic energy):

\Gamma = \frac12 \int \varrho v^2 \mathrm dV,

\Pi is \frac23 of the random kinetic energy of thermal particles (molecular motion), or \frac13 of random kinetic energy of relativistic particles (microscopic kinetic energy):

\Pi = \int p \mathrm dV,

\mathscr M is the magnetic energy within S:

\mathscr M = \frac1{8\pi} \int B^2 \mathrm dV

and W is the total gravitational energy of the system if masses outside S don’t contribute to the potential:

W = - \int \varrho \mathbf r \cdot \nabla \Phi \mathrm dV.

Among all these terms, the most used ones are \Gamma, \mathscr M and W. But most often the equation is just quoted as 2\Gamma+W=0. Note that the virial theorem always holds, inapplicability is only a problem when important terms are omitted.

This kind of simple analysis is often used to determine how bound a system is, and predict its future, e.g. collapse, expansion or evaporation. Specific examples will show up later in the course, including instability analyses.

The virial theorem as Chandrasekhar and Fermi formulated it in 1953 is the following:

\underbrace {2T_m} _{2\Gamma} + \underbrace {2T_k} _{3\Pi} + \underbrace {\Omega} _{W} + \mathscr M = \underbrace {0} _{\frac {\mathrm D^2 I} {\mathrm D t^2}}.

This uses a different notation but expresses the same idea, which is very useful in terms of the ISM.

Radiative Transfer

The specific intensity of a radiation field is defined as the energy rate density with respect to frequency, cross sectional area, solid angle and time:

\mathrm dI_\nu = \frac {\mathrm dE} {\mathrm d\nu \mathrm dA \mathrm d\Omega \mathrm dt}

[\mathrm dI_\nu] = 1 \;\mathrm{erg} \; \mathrm {Hz}^{-1} \; \mathrm {cm}^{-2} \; \mathrm {sr}^{-1} \; \mathrm s ^{-1}

in cgs units. It is important to note that specific intensity does not change during the propagation of a ray, no matter what its geometry is, unless there is extinction or emission in the path.

Specific intensity is a function of position, frequency, direction and time. Integrating over all directions, we get the specific flux density, which is a function of position, frequency and time:

F_\nu = \int_{4\pi} I_\nu \cos \theta \mathrm d\Omega

[F_\nu] = 1 \;\mathrm{erg} \; \mathrm {Hz}^{-1} \; \mathrm {cm}^{-2} \; \mathrm s ^{-1}

where \theta is usually assumed to be zero.

A conventional unit of specific flux density, especially in radioastronomy, is jansky, named after the American radio astronomer Karl Guthe Jansky:

1 \; \mathrm {Jy} = 10^{-23} \;\mathrm{erg} \; \mathrm {Hz}^{-1} \; \mathrm {cm}^{-2} \; \mathrm s ^{-1} = 10^{-26} \;\mathrm{W} \; \mathrm {Hz}^{-1} \; \mathrm {m}^{-2}

The specific energy density of a radiation field is

u_\nu = \frac1c \int_{4\pi} I_\nu \mathrm d\Omega

[u_\nu] = 1 \;\mathrm{erg} \; \mathrm {Hz}^{-1} \; \mathrm {cm}^{-3}

The mean specific intensity is the specific intensity at a given position and time averaged over all directions:

J_\nu = \frac1{4\pi} \int_{4\pi} I_\nu \mathrm d\Omega = \frac {cu_\nu} {4\pi} \left( = \frac {F_\nu}{4\pi} \textrm{ if } \theta=0 \right)

[J_\nu] = 1 \;\mathrm{erg} \; \mathrm {Hz}^{-1} \; \mathrm {cm}^{-2} \; \mathrm {s}^{-1}

The above “specific” quantities all have their frequency integrated counterparts: intensity, flux density, energy density and mean intensity.

The emission and absorption coefficients j_\nu and \alpha_\nu are defined as the coefficients in the differential equations governing the change of intensity in a ray transversing some medium:

\mathrm dI_\nu = ( j_\nu - \alpha_\nu I_\nu) \mathrm ds.

The emissivity \varepsilon_\nu and opacity \kappa_\nu are defined by

j_\nu = \frac {\varepsilon_\nu \varrho} {4\pi}, \alpha_\nu = \varrho \kappa_\nu.

To find the integral equation determining the specific intensity as a function of path, we define the source function S_\nu and the differential optical depth \mathrm d\tau_\nu by

S_\nu = \frac {j_\nu} {\alpha_\nu}; \mathrm d\tau_\nu = \alpha_\nu \mathrm ds.

It is left as an exercise to the reader that these lead to the Radiative Transfer Equation

I_\nu (\tau_\nu) = I_\nu (0) e^{-\tau_nu} + \int_0^{\tau_\nu} e^{-(\tau_\nu-\tau'_\nu)} S_\nu(\tau'_\nu) \mathrm d \tau'_\nu

If the source function is constant in a medium, then the two limiting cases are

      • optically thin \tau_\nu \to 0: I_\nu (\tau_\nu) \sim I_\nu(0) + S_\nu \tau_\nu
      • optically thick \tau_\nu \to \infty : I_\nu (\tau_\nu) \to S_\nu

The Planck function, named after the German physicist Max Karl Ernst Ludwig Planck, is

B_\nu (T) = \frac {2h\nu^3}{c^2} \cdot \frac1 {e^{\frac{h\nu}{kT}}-1}

Planck’s law says that a black body, that is, an optically thick object in thermal equilibrium, will emit radiation of specific intensity given by the Planck function, also known as the blackbody function. Any radiation with this specific intensity is called blackbody radiation. A similar concept is thermal emission, which in fact is defined as S_\nu=B_\nu, therefore only implies blackbody radiation if the emitting medium is optically thick. Note that in case of thermal emission, one can substitute the definition of source function to write j_\nu=\alpha B_\nu, which is known as Kirchoff’s law.

The brightness temperature T_b of a radiation of specific intensity I_\nu at a given frequency \nu is defined as the temperature that a black body would have to have in order to have the same specific intensity:

I_\nu = B_\nu (T_\nu)

There are two asymptotic approximations of the blackbody radiation: if h\nu \ll kT, that is, small frequency or high temperature, we have the Rayleigh–Jeans approximation for the specific intensity:

I_\nu^{\mathrm {RJ}} (T) = \frac {2\nu^2}{c^2} \cdot kT

This is valid in most areas of radio astronomy, except for example some cases of thermal dust emission.

In the other limit, where h\nu \gg kT, that is, large frequency or low temperature, we have the Wien approximation, named after the German physicist Wilhelm Carl Werner Otto Fritz Franz Wien who derived it in 1983:

I_\nu^{\mathrm {W}} (T) = \frac {2h\nu^3}{c^2} \cdot e^{-\frac{h\nu}{kT}}

The brightness temperature T_\mathrm b of a radiation of specific intensity I_\nu at a given frequency \nu is defined as the temperature that a black body would have to have in order to have the same specific intensity:

I_\nu = B_\nu (T_\mathrm b)

Now assume the background radiation of brightness temperature T_\mathrm{bg} traverses a medium of temperature T, source function S_\nu=B_\nu(T), and optical depth \tau_\nu. In the Rayleigh–Jeans regime B\propto T, therefore the differential and integral forms of the radiative transfer equation become

\frac {\mathrm dT_\mathrm b}{\mathrm d\tau_\nu} = T - T_\mathrm b

T_\mathrm b = T_\mathrm{bg} e^{-\tau_\nu} + T ( 1-e^{-\tau_\nu})

where T_\mathrm b is the brightness temperature of the radiation leaving the medium.

This can be applied to dust emission observations assuming an “emissivity-modified” blackbody radiation. Then the contribution of particles of a given linear size a is

F_\lambda = N_a \frac {\pi a^2} {D^2} Q_\lambda B_\lambda (T)

where N_a is the number of such particles, D is the distace to the observer, and Q_\lambda is the emissivity. Recall that the blackbody intensity with respect to wavelength is

B_\lambda (T) = \frac {2hc^2} {\lambda ^5} \cdot \frac1 {e^{\frac{hc}{\lambda kT}} - 1 }.

We refer the reader to Hilebrand 1983.

It turns out that the emissivity in far infrared follows a power law:

Q_\mathrm{FIR} \propto \lambda ^{-\beta}

      • \beta=0 for blackbody
      • \beta=1 for amorphous lattice-layer materials
      • \beta=2 for metals and crystalline dielectrics}

The observed flux will be

F_\textrm{observed} = \sum_a N_a \frac {\pi a^2} {D^2} Q_\lambda B_\lambda (T).

The same emission can be a result of different values of T, Q, and N.

To determine the total mass of dust, we write

M_\mathrm{dust} = \frac {4\varrho_\mathrm{dust} F_\lambda D^2} {3B_\lambda(T_\mathrm{dust})} \cdot \left\langle \frac a {Q_\lambda} \right\rangle,

where \langle\cdot\rangle is an average weighted appropriately.

If \tau\ll1, then emission essentially depends on grain surface area, and we can write

F_\nu = \kappa_\nu B_\nu \frac {M_\mathrm{dust}} {D^2}.

Hildebrand 1983 gives us an empirical formula for dust opacity, implicitly assuming a gas-to-dust ratio of 100, which can be extended with the power law described above to arrive at

\kappa_nu = 0.1\;\frac{\mathrm{cm}^2}{\mathrm g} \cdot \left( \frac \nu {1200\;\mathrm{GHz}} \right)^\beta.

A typical modern value for interstellar dust is \beta=1.7, but it is highly disputed. It can be determined based on the SED slope in the Rayleigh–Jeans regime, where

F_\nu \propto \nu^{\beta+2}.

ISM of the Milky Way

The ISM in the Milky way can be divided into cold stuff and hot stuff. Cold stuff are dust and gas. Katherine will talk briefly about the importance of cooling via CO emission. Hot stuff will be discussed by Ragnhild and Vicente (also see Spitzer 1958), and Tanmoy already talked about supernovae. Suggested reading is Chapters 5 and 19 of Draine. In particular, if you ever need a reference of term symbols (uppercase Greek letters with subscripts and superscripts), see page 39.

Cold ISM

The molecular gas is mostly composed of H_2 molecules. This, however, rarely can be observed directly: it has no dipole moment, therefore no rovibrational energy levels. (As exceptions, UV absorption lines can be observed in hot H_2 gas, vibrational lines can be observed in shocked ISM, and NIR lines can be observed in “cold” hot ISM.) Instead of observing hydrogen lines, trace species are used as proxies to infer hydrogen density. The choice of preferred trace species depends on their abundance. To quantify this, we define the critical density as

n_\mathrm{critical} = \frac {A_\mathrm{ul}} {\gamma_\mathrm{ul}},

where A is the Einstein coefficient and \gamma is the collision rate for a given transition. This of course will depend on temperature, as the collision rate is the product of the temperature-independent collisional cross section \sigma and the temperature-dependent particle velocity v. For a detailed explanation of critical density, see pp.~81–87 of Spitzer, and Chapter 19 of Draine.

Typically assumed values are T\sim100\;\mathrm K, v\sim10^5\;\frac{\mathrm{cm}}{\mathrm s}, \sigma\sim10^{-15}\;\mathrm{cm}^{2}, \gamma_\mathrm{ul}\sim10^{-10}\;\frac{\mathrm{cm}^3}{\mathrm s}. The Einstein coefficient for the J=1\to0 transition of CO is A_{10}\sim6\cdot10^{-8}\;\mathrm s^{-1}, yielding a critical density of n_\mathrm{critical}\sim6\cdot10^2\;\mathrm {cm}^{-3}. Spitzer gives us fiducial values for critical densities around T=100\;\mathrm K:

      • CO, J = 1 \to 0, \lambda = 2.6 mm, A = 6 \cdot 10^{-8} \rm s^{-1}, n_{\rm crit} = 4 \cdot 10^3 {~\rm cm^{-3}}
      • NH_3, J = 1 \to 1 \lambda = 12.65 mm, A = 1.7\cdot10^{-7}, n_{\rm crit} = 1.1\cdot10^4 {~\rm cm^{-3}}
      • CS, J=1\to 0, \lambda = 6.12 mm, A = 6.12 {\rm s}^{-1}, n_{\rm crit} = 1.1\cdot10^5 {\rm ~cm^{-3}}
      • HCN, J=1\to 0, \lambda = 3.38 mm, A = 2.5 \cdot10^{-5} {\rm s^{-1}}, n_{\rm crit} = 1.6\cdot10^6 {\rm~ cm^{-3}}

In practice, the following tracers are used for cold clouds in the range of 3\;\mathrm K \lessapprox T \lessapprox 100 \; \mathrm K:

      • Low density (n=10 cm^-3): 12CO
      • Dark Cloud (n=300): 13CO, OH
      • Dense Core (n=10^3): C18O, CS
      • Dense Core (n=5 * 10^3): NH3, N2H+, CS
      • Very Dense (n=10^8): OH Masers
      • Very Very Dense (n=10^10): H20 masers

See handout on dense core multi-line half power contours from Myers 1991, and note that sizes do not proceed as critical densities would suggest.

As a reminder, electronic energy levels are much further apart, resulting in higher frequency transition lines. Vibrational levels are a few orders of magnitude more tightly spaced, then rotational energy levels are even closer. Here we give a brief overview of rotational line structure, see Chapter 5.1.5 of \cite{draine:2010} for more details.

Quantum mechanically, a diatomic molecule has energy levels identified by J=0,1,2,\ldots with energy E_{\rm rot} = \frac{J(J+1) \hbar^2}{2I_v} = B_v J (J+1) (which would have J^2 instead in the quasiclassical theory). Here I_v=r_vm_\mathrm r is the moment of inertia, r_v is the distance between the atoms, and m_\mathrm r is their reduced mass. B_v is the rotation constant. The v index signifies that these values are valid for a given vibrational state only. Therefore the transition energy between the states J-1 and J is

\Delta E_\mathrm{rot} = 2J B_v.

In case of ^{12}C^{16}O, B_v=5.75\cdot10^{16} in some mysterious units. For the J=1\to0 transition,

\Delta E_\mathrm{rot}=4.7\cdot10^{-4}\;\mathrm{eV}, corresponding to \nu=115.271\;\mathrm{GHz}, or \frac{h\nu}k=5.5\;\mathrm K. The J=2\to1 and J=3\to2 transitions have twice and three times higher energy differences, respectively.

We can estimate which transition gives maximum emission (in terms of number of photons) at a given temperature:

E_\mathrm{rot} = B_v J (J+1)

k T_\mathrm{rot} = B_v J (J+1)

J_\mathrm{max} \approx \sqrt {\frac {kT_\mathrm{rot}}{B_v}}

      • The SMA and ALMA are sensitive to J=4-7
      • KAO and Sophia are sensitive to J=20-40

Note that 1\to0 at 2.6 mm is more visible from Earth than 2\to1 due to atmospheric extinction.

Collisional Excitation

In LTE, the transition rates for a collisional u \rightleftarrows l transition are

C_{ul} = n\gamma_{ul}

C_{lu} = C_{ul} \frac{g_u}{g_l} e^{-\frac{h\nu}{T_K}}

where \gamma_{ul}=\langle \sigma_{ul} v \rangle is the rate coefficent for the transition, \sigma is the collisional cross section and v is the particle velocity.

In case of Maxwellian velocity distribution, we have

\gamma_{ul} = \frac 4 {\sqrt \pi} \left( \frac \mu {2kT_K} \right)^{\frac32} \int_0^\infty \sigma_{ul} v^3 e^{-\frac{\frac12\mu v^2}{kT_K}} \mathrm dv

where \mu is the reduced mass. For neutral-neutral state transitions, \gamma is typically 10^{-11}\sim10^{-10}\;\mathrm{cm}^3\;\mathrm s^{-1}, and for an ionizing transition, \gamma\sim10^{-9}\;\mathrm{cm}^3\;\mathrm s^{-1}.

In equilibrium, we have

\dot n_u = n_l C_{lu} - n_u C_{ul} - n_u A_{ul}

0 = n_l C_{ul} \frac{g_u}{g_l} e^{-\frac{h\nu}{T_K}} - n_u C_{ul} - n_u A_{ul}

0 = (n-n_u) \frac{g_u}{g_l} e^{-\frac{h\nu}{T_K}} - n_u - n_u \frac {A_{ul}}{n C_{ul}}

n \frac{g_u}{g_l} e^{-\frac{h\nu}{T_K}} = n_u \frac{g_u}{g_l} e^{-\frac{h\nu}{T_K}} + n_u + n_u \frac {n_\mathrm{crit}}{n}

\frac {n_u} n = \frac {\frac{g_u}{g_l} e^{-\frac{h\nu}{T_K}}} {\frac{g_u}{g_l} e^{-\frac{h\nu}{T_K}} + 1 + \frac {n_\mathrm{crit}}{n}}

where n_\mathrm{crit}=\frac{A_{ul}}{C_{ul}} is the critical density. When n\gg n_\mathrm{crit}, collisions dominate over spontaneous emission, resulting in

\frac {n_u} {n_l} = \frac {g_u}{g_l} e^{-\frac{h\nu}{T_K}}.

However, in case of n\ll n_\mathrm{crit}, spontaneous emission dominates decay, so each collisionally excitation results in an emission:

\frac {n_u} {n_l} = \frac n {n_\mathrm{crit}} \frac {g_u}{g_l} e^{-\frac{h\nu}{T_K}}.

In case of optically thin emission,

F_\mathrm{line} = \frac {h\nu}{4\pi} A_{ul} \Omega \int n_u \mathrm ds,

where \Omega is the apparent solid angle of the source. Then

\textrm {if } n\ll n_\mathrm {crit}: \quad F_\mathrm{line} = \frac {h\nu}{4\pi} \gamma_{ul} \Omega \frac{g_u}{g_l} e^{-\frac{h\nu}{kT_K}} \int n^2 \mathrm ds \propto n^2

\textrm {if } n\gg n_\mathrm {crit}: \quad F_\mathrm{line} = \frac {h\nu}{4\pi} A_{ul} \Omega \frac{g_u}{g_l} e^{-\frac{h\nu}{kT_K}} \int n \mathrm ds \propto n

Recombination of Ions with Electrons

For more details, see Draine Chapter 14.

In HII regions, most recombination happens radiatively: X^+ + e \to X + h \nu.

An electron with kinetic energy “E” can recombine to any level of hydrogen. The energy of the photon emitted is then given by h\nu = E + I_{nl}, where I_{nl} = 13.6 ev / n^2 is the binding energy of quantum state nl.

There are two extreme cases in recombination (see Baker and Menzel 1938):

      • Case A: The medium is opcially thin to ionizing radiation. Appropriate in shock-heated regions (T ) where density is very low.
      • Case B: Optically thick to ionizing radiation. When a atom recombines to the n=1 state, the emitted lyman photon is immediately reabsorbed, so that recombinations to n=1 does not change the ionization state. This is called the “on the spot” approximation

Optical tracers of recombination: H-\alpha and other lines

Radio: Radio recombination lines (see Draine 10.7). Rydberg States are recombinations to very high (n>100) hydrogen energy levels. Spontaneous decay gives

\nu_{n\alpha} = \frac{2n+1}{[n(n+1)]^2} \frac{I_H}{h} \sim 6.479 \big(\frac{100.5}{n+0.5}\big)^3 {\rm GHz}

A popular line is \nu_{166\alpha} = 1425 {\rm MHz}. This is often observed because its proximity to the 1420 MHz line of HI.

Radio recombination lines often involve masing.

Star Formation in Molecular Clouds

Topics to be covered include:

      • The Jeans Mass
      • Free Fall Time
      • Virial Theorem
      • Instabilities
      • Magnetic Fields
      • Non-Isolated Systems
      • “Turbulence”

An overview of the steps of star formation.

Basic properties of a GMC

Mass: 10^{5-6} M_\odot

Lifetime: Uncertain. Probably 10 Myr (maybe as long as 100 Myr). “Lifetime” is not easily defined or inferred, since clouds are constantly fragmenting, exchanging mass with their surroundings, etc.

We roughly describe the hierarchy and fragmentation of clouds via the following terms:

      • Clump. 10-100 M_\odot. 1pc. The progenitor of stellar clusters
      • Core. 1-10 M_\odot. 0.1 pc. The progenitor of individual stars, and small stellar systems.
      • Star. The end-product of fragmentation and collapse. 1 M_\odot.

Note that, across these scales, density increases by tens of orders of magnitude:

\frac{\rho_{\rm star}}{\rho_{\rm core}} \propto \bigg(\frac{R_{\rm star}}{R_{\rm core}} \bigg)^3 = \bigg(\frac{.1 \times 3 \times 10^{18} cm}{7 \times 10^{10} cm} \bigg)^3 \sim 10^{20}



\caption{Schematic of how a GMC fragments into clumps, cores, and stars}


The Jeans Mass

For further details, see Draine ch 41.

Let’s analyze the stability of a sphere of gas, where thermal pressure is balanced by self-gravity.

Start with the basic hydro equations (conservation of mass, momentum, and Poisson’s equation for the gravitational potential):

\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \vec{v}) = 0

\frac{\partial v}{\partial t} + (\vec{v} \cdot \nabla) \vec{v} = -\frac{1}{\rho}\nabla P - \nabla \phi

\nabla^2 \phi = 4 \pi G \rho

Consider an equilibrium solution rho_0(\vec{r}), P_0(\vec{r}), etc such that time derivatives are zero. Let’s perturb that solution slightly, and analyze when that perturbation grows unstably.

\vec{v} = \vec{v_0} + \vec{v_1}, \rho = \rho_0 + \rho_1, P = P_0 + P_1, \phi = \phi_0 + \phi_1 latex.

The linear hydro equations, to first order in the perturbations, are

\frac{\partial \rho_1}{\partial t} + \vec{v_0} + \nabla \rho_1 + \vec{v_1} \cdot \nabla \rho_0 = -\rho_1 \nabla \cdot \vec{v_0} - \rho_0 \nabla \cdot \vec{v_1}

\frac{\partial v_1}{\partial t} + (\vec{v_0} + \cdot \nabla) \vec{v_1} + (\vec{v_1} \cdot \nabla) \vec{v_0} = \frac{\rho_1}{\rho_0^2} \nabla P_0 - \frac{1}{\rho_0} \nabla P_1 -\nabla \phi_1

\nabla^2 \phi_1 = 4 \pi G \rho_1

Lets restrict our attention to an isothermal gas, so the equation of state is P = \rho c_s^2, where c_s is the isothermal sound speed. Then, the momentum becomes

\frac{\partial \vec{v_1}}{\partial t} + (\vec{v_0} + \nabla) \vec{v_1} + (\vec{v_1} \cdot \nabla)\vec{v_0} = -c_s^2 \nabla(\frac{\rho_1}{\rho_0}) - \nabla \phi_1

Jeans took these equations and added:

      • Uniform density to start with (\nabla \rho_0 = 0)
      • Stationary gas (v_0 = 0)
      • Gradient-free equilibrium potential \nabla \phi_0 = 0)

Then, the solution becomes (after taking \nabla \cdot the momentum equation)

\frac{\partial^2 \rho_1}{\partial t^2} = c_s^2 \nabla^2 \rho_1 + (4 \pi G \rho_0) \rho_1

Now consider plane wave perturbations

\rho_1 \propto \times exp(i (\vec{k} \cdot \vec{r}) - \omega t)

\omega^2 = k^2 c_s^2 - 4 \pi G \rho_0

Define k_j^2 = 4 \pi G \rho_0 / c_s^2, so

\omega^2 = (k^2 - k_J^2) c_s^2

\omega is reall IFF k . Otherwise, \omega is imaginary, and there is exponential growth of the instability. This then leads to a Jeans Length:

\lambda_J = 2 \pi / k_J = \bigg(\frac{\pi c_s^2}{G \rho_0} \bigg)^{1/2}

Gas is “Jeans Unstable” when \lambda . Exponentially growing perturbations will cause the gas to fragment into parcels of size \lambda \sim \lambda_J.

Converting the Jeans length into a radius (assuming a sphere) yields

M_J = 0.32 M_\odot \big(\frac{T}{10K}\big)^{3/2} \big(\frac{m_H}{\mu}\big)^{3/2} \big(\frac{10^6 cm^{-3}}{n_H}\big)^{1/2}

This defines a “preferred mass” for substructures within a cloud.

Let’s plug in values for a dense core: T=10~K,~\mu = 2.33~{\rm amu},~n_H = 2 \times 10^5~{\rm cm}^{-3}. This yields M_J = 0.2~M_\odot. If we instead plug in numbers appropriate for the mean conditions in a GMC, T=50~{\rm K},~\mu=2.33~{\rm amu},~n_{\rm H}=200~{\rm cm}^{-3}, we get M_J=70~M_\odot.

Note that, once gravitaional collapse and heating set in, our isothermal sphere assumptions are no longer valid.

Collapse Timescale

For large scales, the growth time for the Jeans instability is

\tau_J \sim \frac{1}{k_J c_s} \sim \frac{1}{\sqrt{4 \pi G \rho_0}} = 2.3 \times 10^4 yr \big(\frac{10^6 cm^{-3}}{n_H} \big)^{1/2}

For n_H = 1000, this is about 0.7 Myr. Compared to a “free fall time” (collapse timescale for a pressure-less gas),

\tau_{ff} = \big( \frac{3 \pi} {32 G \rho_0} \big)^{1/2} = 4.4 \times 10^4 yr \big( \frac{10^6 cm^{-3}}{n_H} \big)^{1/2}

For n_H=1000 this is 1.4 Myr — slightly longer than growth time.

The Jeans Swindle

There is a sinister flaw in Jean’s analysis. Assuming \nabla \phi = 0 implies that \nabla^2 \phi = 0. The only way to satisfy this everywhere is for \rho_0 = 0. However, more rigorous analysis verifies that Jean’s approach still yields approximately correct results.

More Realistic Treatment

The Bonnor Ebert Sphere is the “largest mass an isothermal gas sphere can have in a \textbf{pressurized medium} while staying in hydrostatic equilibrium.

Numerical Star Formation Simulations

(Notes from Guest Lecture by Stella Offner)

We start by writing down the conservation laws for mass, momentum and energy, each in the form of “rate of change of conserved quantity + flux = source density”:

\frac {\partial \varrho}{\partial t} + \nabla (\varrho v) = 0

\frac {\partial(\varrho t)}{\partial t} + \nabla (\varrho v^2 + p) = - \varrho \nabla \phi

\frac {\partial(\varrho E)}{\partial t} + \nabla (\varrho v^3 + pv) = - \varrho v \nabla \phi

Also, the Poission equation for gravity is

\nabla^2 \phi = - 4 \pi G \varrho

The initial conditions are determined by T, \varrho(\mathbf r), \mu and bulk \mathbf v(\mathbf r), which further determine the total mass, total angular momentum, turbulence and other global properties.

Grid-based codes

One type of star formation simulations is grid-based. Many of these feature adaptive mash refinement, increasing spatial resolution in areas where parameters vary on smaller scales. Examples of such codes are Orion, Ramses, Athena, Zeus and Enzo.

This algorithm stores \varrho and \mathbf v values at nodes indexed by j. The values at time step n+1 are calculated from those at time step n by discretized versions of the conservation equations. For instance, a discretized mass equation can be written as

\frac {\varrho_j^{n+1}-\varrho_j^n}{\Delta t} + \frac{\varrho_{j+1}^n v_{j+1}^n - \varrho_{j-1}^n v_{j-1}^n} {2\Delta x} = 0

First, a homogeneous grid is created, with resolution satisfying the Truelove condition, also called as the Jeans condition:

\Delta x \leqslant J \lambda_\mathrm J

where the empirical value for the Jeans number J is \frac14. Then the simulation is carried on, repartitioning a square by adding nodes if deemed necessary based on the spatial variation of \varrho or \mathbf v. The timestep is global, but it can also be adaptively changed, as long as it satisfies the Courant condition

\Delta t \leqslant C \cdot \min \left( \frac {\Delta x}{v+c_\mathrm s} \right)

ensuring that gas particles do not move more than half a cell in one time step.

Particle-Based Codes

The other family of codes is particle-based, also referred to as smooth particle hydrodynamics codes. Here individual pointlike particles are traced. Examples codes are Gadget, Gasoline and Hydra.

To solve the equations, density has to be smoothed with a smoothing length h that has to be larger than the characteristic distance between the particles. Usually we want approximately twenty particles within the smoothing radius. The formula for smoothing is

\langle \varrho (\mathbf r) \rangle = \sum_\mathrm{particles} m \omega (\mathbf r-\mathbf r',h)

where an example of the smoothing kernel is

\omega (\mathbf r- \mathbf r', h) = \frac1{Th^3}

      • 1-\frac32u^2+\frac34u^3 \textrm{if } 0 \leqslant u \leqslant 1
      • \frac14(2-u)^2 \textrm{if } 1
      • 0 ~\textrm{if } 2

where u = \frac{|\mathbf r-\mathbf r'|}h.

The Lagrangian of the system is

L = \sum_i \frac12 m_i v_i^2 - \sum_i m_i \phi(\mathbf r_i).

The Euler–Lagrange equations governing the motion of the particles are

\frac {\partial L}{\partial \mathbf r_i} - \frac {\mathrm d}{\mathrm dt} \frac {\partial L}{\partial \dot {\mathbf r}_i} = 0.

It turns out that solving these equations leads to

\frac {\partial v_i}{\partial t} = -\sum_j m_j \left( \frac {p_i}{\varrho_i^2} + \frac {p_j}{\varrho_j^2} + Q_{ij} \right) \nabla_i \omega(\mathbf r_i-\mathbf r_j, h),

where \mathbf Q is an viscosity term added artificially for more accurate modeling of transient behavior. The nature of the simulation requires that particles are stored in a data structure that facilitates easy search based on their position, for example a tree. As particles move, this structure needs to be updated.

The Jeans condition for these kind of simulations is that

M_\mathrm {min} \leqslant M_\mathrm J

where M_\mathrm {min} is the minimum resolved fragmentation mass, equal to the typical total mass within smoothing length from any particle.


The advantages of the grid-based codes is that they are

      • accurate for simulating shocks,
      • more stable at instabilities.

Their disadvantages are that

      • grid orientation may cause directional artifacts,
      • grid imprinting might appear at the edges.

Also, they are more likely to break if there’s a bug in the code. On the other hand, particle-based codes have the advantages of being

      • inherently Lagrangian,
      • inherently Galilei invariant, that is, describing convection and flows well,
      • accurate with gravity,
      • good with general geometries.

At the same time, they suffer from

      • resolution problems in low density regions,
      • the need for an artificial viscosity term,
      • the need for post processing to extract information on density and momentum distribution,
      • statistical noise.

See slides for the effect of too large J for a grid-based simulation, too small h for a particle-based simulation, for examples involving shock fronts, Kelvin–Helmholtz instability and Rayleigh–Taylor instability, and star formation. Note that the simulated IMF matches theory and observations within an order of magnitude, but this is very sensitive to the initial temperature of the molecular cloud.

Some thoughts on Jeans scales

Recall the equation for the Jeans Mass:

M_J = \frac1{8} \big( \frac{\pi k T}{G \mu} \big)^{3/2} \frac1{\rho^{1/2}} \propto \frac{T^{3/2}}{\rho^{1/2}}

In common units, this is

M_J = 0.32 M_\odot \big(\frac{T}{10K} \big)^{3/2} \big(\frac{m_H}{\mu}\big)^{3/2} \big( \frac{10^6 {\rm cm^{-3}}}{n_H} \big)^{1/2}

      • 10K, ~2.33 {\rm amu}, ~ n_H = 2 \times 10^5 {\rm cm^{-3}} \to 0.2 M_\odot
      • 50K, ~2.33 {\rm amu}, ~ n_H = 200 {\rm cm^{-3}} \to 70 M_\odot

This raises a question — what jeans mass do we choose if a gas cloud is hierarchical, with many different densities and temperatures? This motivates the idea of turbulent fragmentation, wherein a collapsing cloud fragments at several scales as it collapses. We will return to this

Also recall that the Jeans growth timescale for structures much larger than the Jeans size is

\tau_J = \frac1{k_J c_s} = \frac1{\sqrt{4 \pi G \rho_0}} = \frac{2.3 \times 10^4 {\rm yr}}{\sqrt{n_H / 10^6 {\rm cm^{-3}}}}

for n_H = 1000, ~ \tau_J \sim 0.7 {\rm Myr}

Compare this to the pressureless freefall time:

\tau_{ff} = \big(\frac{3 \pi}{32 G \rho_0} \big)^{1/2} = \frac{4.4 \times 10^4 {\rm yr}}{\sqrt{n_H / 10^6 {\rm cm^{-3}}}}

From which we see

\frac{\tau_{ff}}{\tau_J} = \pi \sqrt{3/8} = 1.92

The Jeans growth time is about 1/2 the free-fall time

Finally, compare this to the crossing time in a region with n = 1000 {\rm cm^{-3}}. The sound speed is c_s = \sqrt{kT / \mu}, which is 0.27 km/s for molecular hydrogen, and .08 km/s for 13CO.

However, note that the observed linewidth in 13CO for such regions is of the order 1 km/s. Clearly, non-thermal energies dominate the motion of gas in the ISM.

Recall that the jeans length is \lambda_J = \big( \frac{\pi c_s^2}{G \rho} \big)^{1/2}

Plugging in the thermal H2 sound speed yields \lambda_J = 1 pc. The crossing time for CO gas moving at 1 km/s is thus 1Myr in this region.

So, in a cloud with n \sim 10^3 and T \sim 20K, the Jeans growth time, free fall time, and crossing time are all comparable.

There is a problem with our derivation, though! We assumed the relevant sound speed that sets the Jeans length is the thermal hydrogen sound speed. However, these clouds are not supported thermally, as the non-thermal, turbulent velocity dispersion dominates the observed linewidths. This suggests we use and equivalent Jeans length where we replace the sound speed by the turbulent linewidth.

This is often done in the literature, but its sketchy. Using turbulent linewidths in a virial analysis assumes that turbulence acts like thermal motion — i.e., it provides an isotropic pressure. This may not be the case, as turbulent motions can be partially ordered in a way that provides little support against gravity.

The question of how a cloud behaves when it is not thermally supported leads to a discussion of Larsons Laws (JC 2013)

Larsons Legacy

see these slides

Stellar Winds in Star Forming Regions

Guest Lecture by Hector Arce. See this post for notes. See also this movie

Introduction to shocks

Shocks occur when density perturbations are driven through a medium at a speed greater than the sound speed in that medium. When that happens, sound waves cannot propagate and dissipate these perturbations. Overdense material then piles up at a shock front

Collisions in the shock front will alter the density, temeprature, and pressure of the gas. The relevant size scale over which this happens (i.e., the size of the shock front) is the size scale over which particles communicate with each other via collisions:

l = (\sigma n)^{-1}. If n=10^2 {\rm cm^{-3}}, ~\sigma \sim (10^{-9} {\rm cm})^2 \to l = 0.01 {\rm pc}

A shock

An excellent reference offering a three-page summary of important shock basics (including jump conditions) is available in this copy of a handout from Jonathan Williams ISM class at the IfA in Hawaii.

Shock De-Jargonification

See this handout and this list of examples.

Shock(ing facts about) viscosity

A perfect shock is a discontinuity, but real shocks are smoothed out by viscosity. Viscosiity is the resistance of a fluid to deformation (it is the equivalent of friction in a fluid)

Inviscid or ideal fluids have zero viscosity. Colder gasses are less viscious, because collisions are less frequent.

Simulations include numerical viscosity by accident — these are the result of computational approximations / discretization that violate Eulers equations and lead to a momentum flow that acts like true viscosity. However, the behavior of numerical viscosity can be unrealistic and of inappropriate magnitude

True shocks are a few mean free paths thick, which is typically less than the resolution of simulations. Thus, simulations also add artificial viscosity to smooth out and resolve shock fronts.

Rankine Hugoniot Jump Conditions

This section is based on the 2013 Wikipedia entry for “Rankine-Hugoniot Conditions,” archived as a PDF here.

The Jump conditions describe how properties of a gas change on either side of a shock front. They are obtained by integrating the Euler Equations over the shock front. As a reminder, the mass, momentum, and energy equations in 1 dimension are

      • \frac{d \rho}{d t} = - \frac{d}{dx} (\rho u)
      • \frac{d \rho u} {d t} = - \frac{d}{dx} (\rho u^2 + P)
      • \frac{d \rho E}{d t} = - \frac{d}{dx} \big[ \rho u (e + 1/2 u^2 + P / \rho) \big]

Where E = e + 1/2 u^2 is the fluid specific energy, and e is the internal (non-kinetic) specific energy.

We supplement these equations with an equation of state. For a adiabatic shocks (a bad term, but describes shocks where radiative losses are small), the equation of state is P = (\gamma - 1) \rho e

Generally, what a jump condition is a condition that holds at a near-discontinuity. It ignores the details about how a quantity changes across the discontinuity, and instead describes the quantity on either side.

The general conservation law for some quantity w is

\frac{d w}{d t} + \frac{d}{dx} f(w) = 0

A shock is a jump in w at some point x = x_s(t). We call x1 the point just upstream of the shock, and x2 the point just downstream

A schematic defining the notation for the shock jump conditions

Integrating the conservation equation over the shock,

\frac{d}{dt} \big ( \int_{x_1}^{x_s(t)} w dx + \int_{x_s(t)}^{x_2} w dx \big ) = - \int_{x_1}^{x_2} \frac{d}{dx} f(w) dx

w1 \frac{dx_s}{dt} - w2 \frac{dx_s}{dt} + \int_{x_1}^{x_s(t)} w_t dx + \int_{x_s(t)}^{x_2} w_t dx = -f(w) |_{x_1}^{x_2}

This holds because \frac{ dx_1}{dt} = 0, \frac{dx_2}{dt} = 0 in the frame moving with the shock

Let x_1 \to x_s(t), ~ x_2 \to x_s(t) when

\int_{x_1}^{x_s(t)} w_t dx \to 0, ~ \int_{x_s(t)}^{x_2} w_t dx \to 0

and in the limit S(w_1 - w_2) = f(w_1) - f(w_2)

Where S = \frac{d_{x_s(t)}}{dt} = \text{characteristic shock speed}

S = \frac{f(w_1) - f(w_2)}{w_1 - w_2}

The upstream and downstream speeds are constrained via


Moving away from the generic w formalism and to the Euler equations, we find

      • S(\rho_2 - \rho_1) = \rho_2 u_2 - \rho_1 u_1
      • S(\rho_2 u_2 - \rho_1 u_1) = (\rho_2 u_2^2 + P_2) - (\rho_1 u_1^2 + P_1)
      • S(\rho_2 E_2 - \rho_1 E_1) = \big[ \rho_2 u_2 (e_2 + 1/2 u_2^2 + P_2 / \rho_2) \big] - \big[ \rho_1 u_1 (e_1 + 1/2 u_1^2 + P_1 / \rho_1) \big]

These are the Rankine Hugoniot conditions for the Euler Equations.

In going to a frame co-moving with the shock (i.e. v = S – u), one can show that

S = u_1 + c_1 \sqrt{1 + \frac{\gamma + 1}{2 \gamma} (\frac{P_2}{P_1 - 1})}

where c1 is the upstream sound speed c_1 = \sqrt{\gamma P_1 / \rho_1}

For a stationary shock, S = 0 and the 1D Euler equations become

\rho_1 u_1 = \rho_2 u_2

\rho_1 u_1^2 + P_1 = \rho_2 u_2^2 + P_2

\rho_1 u_1( e_1 + 1/2 u_1^2 + P_1 / \rho_1) = \rho_2 u_2 (e_2 + 1/2 u_2^2 + P_2 / \rho_2)

After more substitution and defining h = P / \rho + e \to 2(h_2 - h_1) = (P_2 - P_1) (\frac1{\rho_1} + \frac1{\rho_2}),

      • \frac{\rho_2}{\rho_1} = \frac{P_2/P_1 (\gamma + 1) + (\gamma - 1)}{(\gamma + 1) + P_2 / P_1 (\gamma - 1)} = \frac{u_1}{u_2}
      • \frac{P_2}{P_1} = \frac{ \rho_2 / \rho_1 (\gamma + 1) + (\gamma - 1)}{ (\gamma + 1) + \rho_2 / \rho_1 (\gamma - 1)}
      • \frac{\rho_2}{\rho_1}

For strong (adiabatic) shocks and an monatomic gas (\gamma = 5/3),

\frac{\rho_2}{\rho_1} \to 4

Introduction to Instabilities

See this handout

Rayleigh-Taylor Instability

The Rayleigh-Taylor instability occurs when a more dense fluid rests on top of a less dense fluid, in the presence of a gravitational potential \phi = g z (though note that accelerations mimic gravity, so many driven jets, etc exhibit the RT instability without gravity)

Notation for the RT instability

To derive the instability, we will solve the momentum and continuity equations while satisfying boundary conditions at the a/b interface

The momentum equation is

\rho \frac{DV}{Dt} = \rho [ \frac{dv}{dt} + v \cdot \nabla v] = - \nabla P - \frac1{8 \pi} \nabla B^2 + \frac1{4 \pi} B \cdot \nabla B - \rho \nabla \phi

The continuity equation is

\frac{D \rho}{D t} = \frac{d \rho}{d t} + v \cdot \nabla \rho = -\rho \nabla \cdot v

To satisfy the continuity equation, we introduce a velocity potential

v = - \nabla \psi

The density is constant (in time) above and below the interface, so the continuity equation implies \nabla^2 \psi = 0

A solution satisfying \psi = 0 for very large Z is

      • \psi = \psi_a = k_a e^{i \omega t - kz} sin(kx) above
      • \psi = \psi_b = k_b e^{i \omega t + kz} sin(kx) below

Where k = 2 \pi / \lambda is the wave number

Take the real part of \psi as the physical answer. Let B = 0, and ignore v \cdot \nabla v in the momentum equation. By integrating the momentum equation over space, we find

\rho \frac{d\psi}{d t} = P + \rho g z

which holds separately above and below the interface

The boundary conditions are that v_z,~ P need to be continuous across the interface

We can find the z of the interface, z_i, by integrating v_z over time, so

      • z_i = - \frac1{i \omega} \frac{d \psi}{dz}
      • z_i = \frac{k}{i \omega} \psi above
      • z_i = - \frac{k}{i \omega} \psi below

To first order in \psi

If v_z is continuous across the boundary then z_i must be the same above and below so

k \psi_a = - k \psi_b at z = z_i

This implies that k_a = -k_b

Continuity of pressure across the intervace gives

P = i \omega \rho \psi - \rho g z

Using the equation for z_i and noting P_a = P_b at z=0 gives

\omega^2 = g K \frac{rho_b - \rho_a}{\rho_b + \rho_a}

If \rho_b then \omega is imaginary, and the instability grows with time

Ionization fronts, HII regions, Stellar Winds, and Disk Structure

Although these 4 objects exist at different scales, their evolutions are coupled. Lets explore this connection

Ionization fronts are caused by stars with large fluxes of ionizing photons. We see them around HII regions. We observe them emission from several ionizaed specias, especially as they recombine and cool

We would like to know how ionization fronts form, expand, and emit. We would also like to know their lifetime, and how different ionization fronts map to different stars in a crowded environment.

The expansion of HII regions is influenced by the time dependent evolution of both the ionization front and the stellar wind with the surrounding medium. Draine 37.2 calculates the expansion of an HII region into a uniform medium

Stellar Winds

We talked about bi-polar flows, but these are a special type of stellar wind. Other winds are less collimated. Different kinds of winds are associated with stars at different evolutionary states:

Pre main sequence stars: T-tauri winds

Main Sequence winds (like the solar wind)

Post-Main Sequence winds / planetary nebulae

Disk Structure

Disks around stars are accretion disks. The structure of a gravitationally bound disk is governed by Keplerian rotation. However, accretion disks have a temperature structure giving rise to thermal forces (plus magnetic/turbulent forces) that can change their internal structure and dynamics

Ionization Fonts

See chapter 20 of Shus Physics of Astrophysics

Note that the mean free path of photons in an ionized gas is much larger than the mean free path in neutral gas:

      • Neutral Gas: Cross section for interaction is size of atom, or \sigma \sim 10^{-17} {\rm cm^{-2}}.
      • Ionized Gas: Interaction is via Thompson scattering, where the cross section is the size of a free electron. \sigma \sim 10^{-24} {\rm cm^{-2}}

We already discussed the Stromgren radius, which sets the size of an idealized HII region by balancing the stellar ionizing photon flux with the rate at which the ionized volume recombines and thus quenches this flux:

R_s = (\frac{3}{4 \pi} \frac{N_\star}{n_0^2 \alpha})^{1/3}

However, this ionized region will have a pressure much higher than the cool, neutral exterior. Thus, HII regions expand over time, and both pressure and ionization fronts are driven into the neutral medium.

Expansion of HII regions

Stage 1


      • Stage begins when a star turns on in a neutral medium
      • Rapid expansion of ionization front into static HI cloud
      • Very little fluid motion (photons travel much faster than shock waves!)
      • This stage ends when the radius approaches the Stromgren radius, and the ionizing photons are quenched by recombinations

Stage 2

      • The hot ionized HII region is over-pressured compared to the ambient exterior. This drives a radiative shock
      • Sweeps up HI gas into a shell. The ionization front is interior to this shell
      • Ends either the star explodes as a supernova, or the radius expands to R_final, where recombinations balance ionizing photons and there is pressure balance between the interior and exterior
      • Realistically, inhomogeneities in the surrounding medium mean that the HII region will evolve in more complicated ways than this

Stellar Winds

See Lamers and Cassivelli 1999

Wind driving mechanisms:

For cool stars (T = 5000K), winds are driven by the pressure expansion of the hot corona. This is like our own sun

For hot stars (T = 10,000 – 100,000K at the surface), the lack of strong convection inhibits the creation of a hot corona. The temperature at these stars surfaces is not hot enough for particles to escape the gravitational potential. Instead, these stars drive winds directly via radiation pressure. These winds travel at v , and have mass loss rates of \sim 10^{-4} M_\odot {\rm yr^{-1}}

The wind mechanism for pre main sequence stars is debated, but probably involves some interaction between an accretion disk and magnetic field

Disk Structure

The initial evidence for accretion disks around young stars was the SED, which shows an excess of IR emission due to accretion-powered luminosity.

The pioneering paper in this field is Adams, Lada and Shu 1987. See also the review paper by Lada 1999

Young stars were initially divided into 3 evolutionary classes, via the shape of there IR SEDS:

Class 0: SED looks like a single, cold blackbody. Emission is due entirely to a cold core, with no central source

Class 0 YSO

Class I: A large infrared excess, with a rising SED in the infrared. Strong IR excess (compared to the blackbody of a hot, compact central source) is due to a massive disk and outflow

Class I YSO

Class II: A flat IR SED, as the flux of the central source rises, and the disk dissipates

Class III: A faling IR SED. The accretion disk is mostly depleted, and adds only a small infrared excess on top of the SED of the central object.


Disks in the context of Star and Planet Formation

For excellent references on star and planet formation, see “Protostars and Planets V”, as well as Stahler and Palla’s “The Formation of Stars”

What is a protostar?

I like the wikipedia entry:

“A contracting mass of gass that represents an early stage in the formation of a star before nucleosynthesis has begun

Since fusion is a negligible energy source in a protostar, its luminosity comes from gravitational contraction. From the virial theorem, half of the gravitational potential energy gain from contraction is converted into kinetic energy, while the other half is radiated away. In other words, for a homogeneous sphere, E_{\rm rad} = \frac{3}{10} \frac{GM^2}{R}.

For a 1 solar mass star contracted to 500 R_\odot,~ E = 2 \times 10^{45} erg.

What is a planet?

Wikipedia says: “A celestial body moving in an elliptic orbit around a star.” This isn’t quite precise enough — could be an asetroid, KBO, binary star, brown dwarf, etc.

How do planets form?

See this post

Key and Open Questions Regarding the Intergalactic Medium

See also this post

The Kenticutt-Schmidt relation empirically links the surface brightness and star forming rate of galaxies. What is the physical explanation for this? In particular, is there a more prescriptive form or interpretation of the KS relation, that is more astrophysical (and perhaps more accurate at predicting star formation rate)?

Not all galaxies are created equal(ly)

What is the origin of Spiral galaxies? Is the spiral density wave theory of Lin and Shu corre t, or might the GMC perturbation theory of D’Onghia and Hernquist be more apropos? Is the formation of grand design spirals different from that of flocculent spirals?

Where to ellipticals come from? They most likely are the remnants of mergers:

      • They they have no ordered rotation. (suggests random mixing)
      • Old stars, little or no dust/gas
      • Often at the center of clusters
      • Always have central Black Holes

“Young” merger galaxies (i.e. before a galaxy has undergone several mergers and formed an elliptical) have extreme star formation rates. SFR can range from 100s to 1000s of stars/year (by comparison, the SFR in the Milky Way is 1 M_\odot/yr

Significant Variations in the Gas/Dust Ratio and Metallicity

This is true both within galaxies and among galaxies.

How does this affect the K-S relationship?

What are the Kennicutt-Schmidt relations?

See also this post

The “Schmidt Law” is due to Marten Schmidt (1959): “The mean density of a galaxy may determine, as a first aproximation, its present gas content and thus its evolution stage”

\dot{M} \propto n_{\rm gas}^2

Here, n is the volume density.

The “Kennicut Law” is similar to the Schmidt law, but uses the surface density:

\Sigma_{\rm SFR} = a \Sigma_{\rm gas}^q

Where, typically, q~1.4 (The Schmidt Law above would give q=1.5). Notably q is not 1, which would imply a constant efficiency function.

Some tracers of star formation rate:

      • UV flux (from unobscured massive young stars)
      • H-alpha flux (also from young stars)
      • FIR flux – reprocessed starlight (but– that’s also dependent on gas surface density, so its sketchy to use in the KS law. That doesn’t mean people don’t do it!)

Some tracers of surface density:

      • HI
      • CO
      • HI+CO
      • HCN and other tracers that preferentially probe high surface densities. Gao and Solomon 2004 quote Log_{IR} = 1 \pm 0.05 Log_{HCN} + 2.9

A big issue in using the KS law is how “independent” the tracers of surface density and star formation rate are. Also, the KS relationship extends over many orders of magnitude. Because of this, its unclear whether KS says something important of the conversion of gas into stars, or rather something trivial like “bigger galaxies have more gas, and form more stars”.

The ISM at very high redshift

See this post

Special Topics

See this post

  1. Should we add links to the (scanned) handwritten notes for each “anchored” section here?

  2. […] in Figure 6. They suggest that the wind is creating a bow shock and a shell of swept-up material (More on outflows in star-forming regions). The broad velocity features on the CO emission line wings reach up to 15 km/s, suggesting the […]

  3. […] in Figure 6. They suggest that the wind is creating a bow shock and a shell of swept-up material (More on outflows in star-forming regions). The broad velocity features on the CO emission line wings reach up to 15 km/s, suggesting the […]

  4. […] derivation from treating small perturbations to the equations of hydrostatic equilibrium (found here) produces a Jeans Length of times the Jeans radius I derived […]

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: