Harvard Astronomy 201b

Posts Tagged ‘shock physics’

ARTICLE: The Acceleration of Cosmic Rays in Shock Fronts (Bell 1978)

In Journal Club 2013 on April 9, 2013 at 5:20 am
Cartoon picture of Bell's diffusive shock acceleration, courtesy of Dr. Mark Pulupa's space physics illustration: http://sprg.ssl.berkeley.edu/~pulupa/illustrations/

Cartoon picture of Bell’s diffusive shock acceleration, courtesy of Dr. Mark Pulupa’s space physics illustration: http://sprg.ssl.berkeley.edu/~pulupa/illustrations/

Summary By Pierre Christian


WordPress does not do equations very well!!! For ease of viewing, I attach here the PDF file for this paper review. It is identical in content with the one on WordPress, but since it is generated via Latex, the equations are much easier to read. On the other hand, the WordPress version has prettier pictures – take your pick… Here’s the Latex-ed version: Latex-ed Version

Here is the handout for said Journal Club: Handout

Bell’s paper can be accessed here: The Acceleration of Cosmic Rays in Shock Fronts

Introduction: Cosmic Rays in Context

In order to appreciate Bell’s acceleration mechanism, it is important for us to first learn some background information about the cosmic rays themselves. These are highly energetic particles, mostly protons (with a dash of heavier elements) that are so energetic that at the highest energies one proton has the energy equivalent to the fastest tennis serves. What can produce such energetic particles?

Our first clue comes from the composition of these cosmic rays. First, there is no electrically neutral cosmic ray particle. Electromagnetic forces must therefore be important in their production. Second, there is an overabundance of iron and other heavy elements in the cosmic ray population compared to the typical ISM composition (Drury, 2012).

The favored cosmic ray production sites are supernova remnants (Drury, 2012). Why not extra-galactic sources? Well, it has been shown that cosmic ray intensity is higher in the inner region of the Milky Way and goes down as we move radially to the outer disk. In addition, cosmic ray intensity is also shown to be less in the Magellanic clouds compared to the Milky Way. Additionally, strong shocks amplify magnetic fields to large proportions, and due to stellar nucleosynthesis, there is an overabundance of iron and heavy elements. I should also mention that a recent Fermi result conclusively proves that supernova remnants produce cosmic rays.

The observed cosmic ray spectrum is given in the figure below. We can see that the bulk of cosmic rays follow a broken power law. In particular, there are two breaks in this power law, one at 3 \times 10^{15} eV, called the ‘knee’ and one at 3 \times 10^{18} eV, called the ‘ankle’. Bell’s contribution is the derivation of a cosmic ray power law spectrum via acceleration in the shocks of supernova remnants. Bell’s model did not explain why the ‘ankle’ and the ‘knee’ exist, and to my knowledge the reason for their presence is still an open question. One explanation is that galactic accelerators cannot efficiently produce cosmic rays to arbitrarily high energies. The knee marks the point where galactic accelerators reach their energetic limits. The ankle marks the point where the galactic cosmic ray intensity falls below the intensity of cosmic rays from extragalactic sources, the so called ultra high-energy (UHE) cosmic rays (Swordy, 2001).

Observed cosmic ray spectrum from many experiments. Originally published by Swordy (2001), and modified by Dr. William Hanlon of the University of Utah (http://www.physics.utah.edu/~whanlon/spectrum.html). This image shows the three powerlaw regimes and the corresponding two breaks: the knee at 3 x 10^{15} eV and the ankle at 3 x 10^{18} eV.

Observed cosmic ray spectrum from many experiments. Originally published by Swordy (2001), and modified by Dr. William Hanlon of the University of Utah (http://www.physics.utah.edu/~whanlon/spectrum.html). This image shows the three power law regimes and the corresponding two breaks: the knee at 3 x 10^{15} eV and the ankle at 3 x 10^{18} eV.

Bell’s Big Ideas

Bell saw the potential of using the large bulk kinetic energy produced in objects such as a supernova remnant to power the acceleration of particles to cosmic ray energies. In order to harness the energy in bulk fluid motions, Bell needed a mechanism to transfer this energy to individual particles. Cited in Bell’s paper, Jokipii and Fisk had in the late sixties and early seventies deduced a mechanism using shocks as a method for accelerating particles. Bell modified and perfected Jokipii’s and Fisk’s mechanisms into something that could be applicable in the acceleration of cosmic ray particles in strong magnetized shocks.

General Idea
The general idea of Bell’s mechanism is that a particle with gyroradius, much larger than the shock’s width can move between the upstream (region the shock is moving into) and downstream regions (region that the shock had already passed in). Every crossing increases the particle’s energy slightly, and after many crossings, the particle is accelerated up to cosmic ray energies. This requires the particle to be confined within a certain region around the shock; there must be a mechanism that keeps the particle bouncing around the shock.

As a shock wave wades through the ISM, it produces turbulence in its wake. In the frame of the shock, the bulk kinetic energy of the upstream plasma is processed into smaller scales (turbulence and thermal) in the downstream region. This turbulence cascades down to smaller scales (via a magnetic Kolmogorov spectrum) before being dissipated far downstream. This magnetic turbulence can scatter fast moving charged particles and deflect them. In particular, a particle trying to escape downstream can be deflected by this turbulence back upstream.

If a charged particle saw some magnetic inhomogeneity (a random fluctuation) in the plasma, it is possible for the particle to scatter off these inhomogeneities. It is also known that much like how a particle moving in air produces sound waves, fast moving charged particles in an MHD fluid can produce Alfven waves. In the upstream region, energetic particles will excite Alfven waves. Now, a particle moving much faster than the Alfven speed will essentially see these waves as magnetostatic inhomogeneities. Furthermore, there is a resonant reaction that can scatter a particle with Larmor radius comparable to the wavelength of the inhomogeneities. In the upstream region of the shock, energetic particles will then scatter off Alfven waves that they themselves generate, deflecting them back into the downstream region.

At this point I would like to point out that the Alfven wave scatterings are elastic. Bell used confusing wording that made it seem that Alfven waves decrease the energy of the upstream particles enough so that they get overtaken by the shock. This is not true; magnetostatic interaction necessarily conserve energy because static magnetic fields can do no work (force is always perpendicular to velocity). What happens is that, in the frame of the shock, particles scatter with Alfven waves through a multitude of small angle scatterings, so that the angle between the particle motion and the background magnetic field undergoes a random walk. After several steps in this random walk, some particles will be scattered backwards towards the downstream direction.

Energy Source
Bell emphasizes that there are two different scattering centers in the process. The downstream scattering centers are the turbulence excited by the shock while the upstream scattering centers are the Alfven waves produced by the energetic particles. What is important is that these two waves are moving at different speeds in the frame of the shock. The energy powering the particle acceleration is exactly harnessed from the speed differential between the upstream and downstream scattering centers. Look at the ‘Intuitive Derivation of Equation (4)’ section for more details on this process.

The Power Law Spectra of Cosmic Rays

Bell found the differential energy spectrum of cosmic rays to be:
N(E) dE = \frac{\mu - 1} {E_0} \left( \frac{E} {E_0} \right) ^{-\mu} dE \; ,
where E_0 is the energy of the particle as it is injected in the acceleration process and \mu is defined to be:
\mu = \frac{2 u_2 + u_1} {u_1 - u_2} + O\left(\frac{u_1 - u_2} {c} \right) \; .
Where O\left(\frac{u_1 - u_2} {c} \right) denotes terms of order \left(\frac{u_1 - u_2} {c} \right) and higher.
Now, u_1 and u_2 are respectively the velocity of the scattering centers. Each are given by:
u_1 = v_s - v_A \; ,
u_2 = \frac{v_s} {\chi} + v_W \; ,
where v_A is the Alfven speed, v_W the mean velocity of scattering centers downstream, v_s the shock velocity, and \chi the factor by which the gas is compressed at the shock (\chi = 4 for high Mach number shocks obeying Rankine-Hugoniot jump conditions). Combining these equations we get:
\mu = \frac{ (2 + \chi) + \chi(2 v_W / v_s - 1/M_A)} {(\chi - 1) - \chi (v_W/v_s + 1/M_A)} \; ,
where M_A is the Alfven Mach number of the shock. In particular, for typical shock conditions this gives us a slope of \sim -2 to \sim -2.5, in agreement with the observed power law of the cosmic ray spectrum.

Damping of Alfven Waves

Bell notices that Alfven waves are damped by two important processes: collision of charged particles with neutral particles and the sound cascade. The first effect can be explained by noting that the magnetic field is flux frozen only to the charged particles in the fluid. If there is a significant amount of neutral particles, the charged particles that are ‘waving’ along with the Alfven waves can collide and scatter with these neutral particles (because the neutral particles do not ‘wave’ along with the magnetic field). This transfers energy from the hydromagnetic Alfven waves to heating the fluid, effectively damping the Alfven waves. Bell notes that the spectral index he derived is only good up to an energy cutoff, above which the spectrum becomes more steep. This energy cutoff is:
\frac{E_{crit}}{\rm{GeV}} = 0.07 \left( \frac{100 u_1} {c} \right)^{4/3} \left(\frac{n_e}{cm^{-3}} \right)^{-1/3}\left(\frac{n_H}{cm^{-3}} \right)^{-2/3} \left( \frac{f(0, p) - f_0(p)} {f_{gal}(p)} \right)^{2/3} \; .
Note that f_0(p) , the background particle distribution, is exactly f_{gal}(p) (the Galactic particle distribution) if the shock is moving through undisturbed interstellar gas. This condition is false for objects that generates multiple shock fronts. Also note that due to compression in the shockwave, f(0,p) can be many times the galactic value.

The second damping mechanism is due to the sound cascade. Alfven waves can interact with and lose energy to magnetosonic waves of lower wavelengths. This requires the sound speed to be less than the Alfven speed:
T < 3100 \left( \frac{B} {3 \mu G} \right)^2 \left( \frac{n_e} {cm^{-3}} \right) ^{-1} K \; .
However, this damping mechanism does not completely remove the Alfven waves, but merely limits the wave intensity. If this process is important, it will allow particles upstream to travel further from the shock before crossing back downstream. Bell does not believe that this will hamper the acceleration process in most physical cases.

Injection Problem

Throughout this paper, Bell assumed that the accelerated particles are sufficiently energetic to be able to pass through the shock. In order for this assumption to be valid, the gyroradius of the particle needs to be larger than the shock width. The typical thermal energy of particles in the fluid is much too low to satisfy this condition. Therefore, there is a need for a pre-acceleration phase where thermal particles become energized enough for them to participate in Bell’s acceleration (Bell, 1978).

A lot of work has been put into answering this question and to my knowledge it is still an open problem. One method due to Malkov and Volk is to note that even when the gyroradius of the particle is smaller than the shock width, there is still a leaking of mildly suprathermal particles from downstream to upstream. These particles will then excite Alfven waves and accelerate in much the same way as Bell’s diffusive shock acceleration (Malkov & Voelk, 1995).

What About Electrons?

Why do we observe a deficit of electrons in the cosmic ray composition? There are two primary reasons. The first is that the injection problem is much more severe for an electron. Recall that the gyroradius is:
r_{gyro} = \frac{mv_{\perp}} {|q|B} \; .
As a consequence of electrons being much less massive than protons (protons are heavier by a factor of about 2000), they have a much smaller gyroradius. In particular, this means that an electron needs to be much more energetic than a proton for it to be able to pass through a shock.

The second reason for the electron deficit is efficient radiative cooling processes. Due to its tiny mass, both inverse Compton and synchrotron cooling are extremely efficient is taking kinetic energy away from an electron. In particular, close to a shock where magnetic fields can be amplified to many times their background value, synchrotron radiation becomes extremely effective. As such, not only is it more difficult for electrons to join in the acceleration process, electrons also lose kinetic energy much more efficiently than protons or other heavier nuclei.

However, it should be noted that there are methods for electrons to participate in the acceleration. One of the cutest ones involves electrons hitching a ride on ions. It turns out that the acceleration timescale of ions is comparable to their electron stripping timescale in typical shock conditions. So, ions that still harbor electrons can be accelerated via Bell’s diffusive shock acceleration before the electrons, now also possessing a large kinetic energy due to the acceleration, get stripped from the ions and join in the acceleration party.

Oblique Shocks

Figure 2 of Bell's paper, oblique shock geometry.

Figure 2 of Bell’s paper, showing oblique shock geometry.

What if the shock normal is not parallel to the magnetic field? In general, since the plasma can move both parallel and perpendicular to the magnetic field lines, we have to include the electric field in the description. However, Bell claimed that we can Lorentz transform to a frame where the shock front is stationary and all bulk fluid motion is parallel to the magnetic field as long as the angle between the shock’s normal and the magnetic field is less than cos^{-1} (v_s / c) (better argued in Drury, 1983). In this frame, the electric field would be non-existent and we can use much of our previous calculations. If velocity downstream and upstream are \vec{w_1} and \vec{w_2} respectively (look at Figure 2 of Bell’s paper, reproduced preceding this paragraph), the energy increase due to one crossing is:

E_{k+1} = E_{k} \left( \frac{1 + \vec{v}_{k1} \cdot (\vec{w}_1 - \vec{w}_2)} {1 + \vec{v}_{k2} \cdot (\vec{w}_1 - \vec{w}_2)} \right ) \; ,
Since \vec{w}_1 and \vec{w}_2 are not parallel, the increase per crossing is smaller than that of the parallel shock. However, the particle’s gyration around magnetic field lines would allow the particle to cross the shock multiple times when it drifts close to the shock (remember, the particle’s gyration radius must be large compared to the shock’s width for it to participate in the acceleration process!). In total, Bell claims that these two effects cancel each other out, resulting in the same power law spectrum.


In this paper, Bell described a novel method to use the bulk kinetic energy in shocks around supernova remnants to accelerate particles to cosmic ray energies. The analytically calculated particle energy power law spectrum agrees with the observed spectrum. However, the study of cosmic ray acceleration is still far from over. An issue left open by this paper is the ‘seed’ of the accelerating particles. Because in standard SNR condition thermal particles do not possess enough energy to participate in the acceleration mechanism, Bell’s process requires a pre-acceleration mechanism. This injection mechanism is still unknown, and we would like to direct interested readers to Malkov & Voelk, 1995 for a hypothesis. The subject of ultra high energy cosmic rays (UHECR) is also left undiscussed. In particular, although most observed cosmic rays are accelerated by Milky Way’s SNR’s, these UHECRs are thought to originate from extragalactic sources. Both the site of their acceleration and the mechanism for said acceleration is still a mystery.

Derivation of Equation (4)

Equation (4) in Bell’s paper is a cornerstone equation which Bell’s argument rests upon. For some reason he did not show the derivation of this equation. The easiest way to derive it is to follow Bell’s advice and perform Lorentz transforms of the energy in the rest frame of the scattering center of the new region. This Lorentz boost will be in the direction parallel to the shock’s normal. For a single upstream to downstream crossing, the particle energy is increased by:
E = \frac{E' + vp'} {\sqrt{1 - v^2/c^2}} = \frac{E' + v \gamma m v_{k1}} {\sqrt{1 - v^2/c^2}} = E' \frac{1 + v v_{k1}/c^2} {\sqrt{1 - v^2/c^2}}
where E' is the energy in the original region (upstream) and E' is the energy in the downstream region; v_{k1} is the velocity at which the particle is crossing from upstream to downstream; and v is the difference of the velocity between scattering centers upstream and downstream parallel to the shock’s normal (direction of Lorentz boost) given by: v = (u_1 - u_2) \cos \theta_{k1} , where \theta_{k1} is the angle the motion is making with the shock’s normal. Putting this together gives:
E_{downstream} = E_{upstream} \left( \frac{1 + v_{k1}(u_1 - u_2) \cos \theta_{k1}/c^2} {\sqrt{1 - ((u_1 - u_2) \cos \theta_{k1})^2/c^2}} \right) \; .
Now, equation (4) in Bell’s paper is the total energy change of the particle if it crosses from upstream to downstream and to upstream again. This means we have to once again perform a Lorentz boost, this time from downstream to upstream. This is the exact same procedure, but now we solve for E' instead of E (inverse Lorentz transformation) since we are measuring everything in the frame of the upstream scattering centers. This gives:
E_{final} = E_{initial} \left( \frac{1 + v_{k1}(u_1 - u_2) \cos \theta_{k1}/c^2} {\sqrt{1 - ((u_1 - u_2) \cos \theta_{k1})^2/c^2}} \right) \left(\frac{\sqrt{1 - ((u_1 - u_2) \cos \theta_{k1})^2/c^2}} {1 + v_{k2}(u_1 - u_2) \cos \theta_{k2}/c^2} \right )
= E_{initial} \left( \frac{1 + v_{k1}(u_1 - u_2) \cos \theta_{k1}/c^2} {1 + v_{k2}(u_1 - u_2) \cos \theta_{k2}/c^2} \right ) \; ,
which is exactly equation (4)!

A More Intuitive Derivation of Equation (4)

Although the previous derivation is valid, I think it is not very physically illuminating. Here is a more intuitive derivation of the amount of energy increase a particle receives due to crossing from upstream to downstream. Note that we will not be using Lorentz transforms for this derivation, so it is only accurate to the lowest order in (u_1- u_2)/c .

What is the amount of momentum increase a particle receive due to crossing from upstream to downstream? Suppose p' is the momentum of the particle after it crosses and p is the momentum of the particle before it crosses:
p' = p + \Delta p
= p + \gamma m \Delta v
= p + \gamma m (u_1 - u_2) \cos \theta_{k1} \; .
Now, what is the change of particle energy due to this increase in momentum?
\Delta E = \int \frac{dp}{dt} dl = \int dp v_{k1} \sim v_{k1} \Delta p \;.
The change of energy therefore is:
E' = (E + v_{k1} \Delta p)
= (E + v_{k1} \gamma m (u_1 - u_2) \cos \theta_{k1})
= E ( 1 + v_{k1} m (u_1 - u_2) \cos \theta_{k1}/c^2) \; ,
where we have used E = \gamma m c^2 on the last line. This is the amount of energy the particle possess after it crosses from the upstream to the downstream. The increase in energy is:
\Delta E = E (v_{k1} m (u_1 - u_2) \cos \theta_{k1}/c^2) \; .
Why is this derivation more illuminating than the Lorentz boosts? For one, it is easy to see where the extra energy comes from. It is obvious from our \Delta E equation that it comes from the difference in velocity of the scattering centers upstream and downstream (u_1 - u_2) . Therefore, cosmic ray acceleration is literally a process where energetic particles steal energy from the bulk fluid motion!!!

Another illumination comes from the fact that this equation is exactly the same if we apply it backwards, from downstream to upstream, as long as we change v_{k1} to the particle velocity going back upstream and \theta_{k1} to the angle the motion makes with the shock’s normal (\theta now defined in the downstream region). Now, u_1-u_2 will have to be changed to u_2-u_1 , giving a minus sign. However, the angle must be inverted, \cos\theta \rightarrow - \cos\theta , netting another minus sign. The energy change is always positive no matter if the particle crosses from upstream to downstream or downstream to upstream. Particles always gain and never lose energy when the cross a shock!!! In particular, the total energy gain by a particle after many crossings is just the product of a bunch of these (v_{k} m (u_1 - u_2) \cos \theta_{k}/c^2) terms!

To get to equation (4), we have to perform two crossings, one upstream to downstream and another downstream to upstream.
Some coordinate change trickery is also required. In particular, \theta_{k1} and \theta_{k2} are both measured in the upstream frame. Therefore, the extra minus sign from \cos\theta \rightarrow - \cos\theta is not present. Therefore:
E_{final} = E_{initial} (1 + v_{k1} m (u_1 - u_2) \cos \theta_{k1}/c^2) (1 - v_{k2} m (u_1 - u_2) \cos \theta_{k2}/c^2) \; ,
where the minus sign comes from the u_2 - u_1 term of the second crossing. It is now obvious to see that this is exactly equation (4) if we Taylor expand the denominator to linear order in (u_1 - u_2)/c . Since equation (7) in Bell’s paper only goes to linear order in (u_1 - u_2)/c , this approximation will give us the same differential energy spectrum as equation (9) of Bell’s paper.

Derivation of Equations (1) & (13)

In this section we shall derive Equation (1) of the paper:
\frac{\partial n} {\partial t} + u_2 \frac{\partial n} {\partial x} = \frac{\partial} {\partial x} \left( D(x) \frac{\partial n} {\partial x} \right) \; .
This equation describes the evolution of the particle density, n(x, t) in the downstream region. The two important mechanisms are advection due to the fluid flow (with advection velocity u_2 ) and diffusion.

Start with a flux continuity equation of the particle density n with its flux, n u_2 .
\frac{\partial n} {\partial t} + \nabla \cdot (n \vec{u_2}) = \rm{Source} - \rm{Sink} \; .
Perform a product rule:
\frac{\partial n} {\partial t} + n \nabla \cdot \vec{u_2} + \vec{u_2} \cdot \nabla n = \rm{Source} - \rm{Sink} \; .
Now, n \nabla \cdot \vec{u_2} = 0 since the advection velocity downstream is u_2 regardless of position, so:
\frac{\partial n} {\partial t} + \vec{u_2} \cdot \nabla n = \rm{Source} - \rm{Sink} \; .
This takes care of advection. What about diffusion? If we know how diffusion affects \partial n / \partial t , we can put diffusion into the Source – Sink term.

Let us consider again a flux continuity equation, but only considering diffusion processes (no advection or other source/sink terms) this time:
\frac{\partial n} {\partial t} + \nabla \cdot \vec{J} = 0 \; .
By the phenomenological Fick’s first law:
\vec {J} = -D(\vec{x}) \nabla n \; ,
where D(\vec{x}) is called the \textbf{diffusion coefficient} and can vary with respect to position in the shock. Plugging this to the previous equation gives:
\frac{\partial n} {\partial t} = \nabla \cdot (D(\vec{x}) \nabla n) \; .
We can interpret this equation as an extra contribution to \partial n / \partial t due to diffusion. In particular, this is the Sink-Source due to diffusion. The complete continuity equation with both advection and diffusion is therefore:
\frac{\partial n} {\partial t} + \nabla \cdot (n \vec{u_2}) = \nabla \cdot (D(\vec{x}) \nabla n) \; .
Now, since things do not vary wildly perpendicular to the shock’s normal, the only important terms in these derivatives are the ones parallel to the shock’s normal (the x component in Figure (1) of the paper). Similarly, \vec{u_2} = u_2 \hat{x} and D(\vec{x}) = D(x) , since we assume that the advection velocity and diffusion coefficient do not vary perpendicular to the shock’s normal.
Therefore, the continuity equation becomes:
\frac{\partial n} {\partial t} + u_2 \frac{\partial n} {\partial x} = \frac{\partial} {\partial x} \left( D(x) \frac{\partial n} {\partial x} \right) \; ,
equation (1) of Bell’s paper. Note that I only need to change the advection velocity from u_2 to u_1 if I want to find a similar equation in the upstream region, giving equation (13) of the paper.


Bell, A. R. 1978, MNRAS, 182, 443

Drury, L. O. . 2012, Astroparticle Physics, 39, 52

Drury, L. O. 1983, Reports on Progress in Physics, 46, 973

Malkov, M. A., & Voelk, H. J. 1995, A&A, 300, 605

Swordy, S. P. 2001, Space Sci. Rev., 99, 85

Further Readings

Blandford, R., & Eichler, D. 1987, Phys, Rep., 154, 1

Schure, K. M., Bell, A. R., OC Drury, L., & Bykov, A. M. 2012, Space Sci. Rev., 173, 491

Skilling, J., 1975a. Mon. Not. R. astr. Soc., 172, 557.

Skilling, J., 1975b. Mon. Not. R. astr. Soc., 173, 245.

Skilling, J., 1975c. Mon. Not. R. astr. Soc., 173, 255.

Wentzel, D. G., 1974. A. Rev. Astr. Astrophys., 12, 71.

Cas A, a supernova remnant where shocks are accelerating cosmic rays, courtesy of wikipedia.

Cas A, a supernova remnant where shocks are accelerating cosmic rays, courtesy of wikipedia.

ARTICLE: An L1551 Extravaganza: Three Articles

In Journal Club 2013 on April 1, 2013 at 11:55 am

Wide-Field Near-Infrared Imaging of the L1551 Dark Cloud by Masahiko Hayashi and Tae-Soo Pyo

Observations of CO in L1551 – Evidence for stellar wind driven shocks by Ronald L. Snell, Robert B. Loren & Richard L. Plambeck

Multiple Bipolar Molecular Outflows from the L1551 IRS5 Protostellar System by Po-Feng Wu, Shigehisa Takakuwa, and Jeremy Lim

Summary by Fernando Becerra, Lauren Woolsey, and Walker Lu


Young Stellar Objects and Outflows

In the early stages of star formation, Young Stellar Objects (YSOs) produce outflows that perturb the surrounding medium, including their parental gas cloud. The current picture of star formation indicates that once gravity has overcome pressure support, a central protostar is formed surrounded by an infalling and self-supported gas disk. In this context outflows are powered by the release of gravitational potential energy liberated by matter accreting onto the protostar. Outflows are highly energetic and often spatially extended phenomena, and are observable over a wide range of wavelengths from x-ray to the radio. Early studies of molecular outflows (predominantly traced by CO emission lines, e.g. Snell et al. 1980, see below) have shown that most of their momentum is deposited in the surrounding medium and so provide a mass loss history of the protostar. In contrast, the optical and near-infrared (NIR) emission trace active hot shocked gas in the flow.

Interactions with the surrounding medium: Herbig-Haro objects, bow shocks and knots

When outflows interact with the medium surrounding a protostar, emission can often be produced. One example of this is emission from Herbig-Haro (HH) objects, which can be defined as “small nebulae in star-forming regions as manifestations of outflow activity from newborn stars”. The most common pictures show a HH object as a well-collimated jet ending in a symmetric bow shock. Bow shocks are regions where the jet accelerates the ambient material. The shock strength should be greatest at the apex of the bow, where the shock is normal to the outflow, and should decline in the wings, where the shocks become increasingly oblique. Another interesting feature we can distinguish are knots. Their origin is still unknown but a few theories have been developed over the years. They can formed due to the protostar producing bursts of emission periodically in time, or producing emission of varying intensity. They can also form due to interactions between the jet and the surrounding Interstellar Medium (ISM), or due to different regions of the jet having different velocities.

An exceptional case: The L1551 region

The L1551 system is an example of a region in which multiple protostars exhibiting outflows are seen, along with several HH objects and knots. This system has been catalogued for over fifty years (Lyons 1962), but ongoing studies of the star formation and dynamical processes continue to the present day (e.g. Hayashi and Pyo 2009; Wu et al. 2009). L1551 is a dark cloud with a diameter of ~20′ (~1 pc) located at the south end of the Taurus molecular cloud complex. The dark cloud is associated with many young stellar objects. These YSOs show various outflow activities and characteristics such as optical and radio jets, Herbig-Haro objects, molecular outflows, and infrared reflection nebulae. We will start by giving a broad view of the region based on Hayashi and Pyo 2009, and then we will focus on a subregion called L1551 IRS 5 following Snell et al. 1980 and Wu et al. 2009.

Paper I: An overview of the L1551 Region (Hayashi and Pyo 2009)

The L1551 region is very rich in YSOs, outflows and their interaction with the ISM. The most prominent of the YSOs in this region are HL Tau, XZ Tau, LkHα 358, HH 30, L1551 NE, and L1551 IRS 5 (see Fig. 1), arrayed roughly north to south and concentrated in the densest part (diameter ~10′) of the cloud. The authors based their study on observations using two narrowband filters [Fe II] (\lambda_c = 1.6444 μm, \Delta\lambda = 0.026 μm), H_2 (\lambda_c = 2.116 μm, \Delta\lambda = 0.021 μm) and two broad-band filters: H (\lambda_c = 1.64 μm, \Delta\lambda = 0.28 μm) K_s (\lambda_c = 2.14 μm, \Delta\lambda = 0.31 μm). The choice of [Fe II] and H_2 is motivated by previous studies suggesting that the [Fe II] line has higher velocity than the H_2, and thus arises in jet ejecta directly accelerated near the central object, while H_2 emission may originate in shocked regions. In the particular case of bow shocks, regions of higher excitation near the apex are traced by [Fe II], while H_2 is preferentially found along bow wings. The broadband filters were chosen for comparison with NIR narrowband filters and comparison with previous studies. The total sky coverage was 168 arcmin2, focused on 4 regions of the densest part of the L1551 dark cloud, including HL/XZ Tau, HH30, L1551 IRS5, some HH objects to the west, L1551 NE, and part of HH 262 (see Fig. 1).


Figure 1: An overview of L1551 (Figure 1 of Hayashi and Pyo 2009)

HL/XZ Region

Some of the features the authors identify in this region are:

  • A faint [Fe II] jet emanating from HL Tau to its northeast and southwest. The H_2 emission is hard to identify in the northeast part, but significant H_2 emission blobs are detected in the southwest part (denoted “H2 jet” in Fig. 2)
  • A diffuse feature is also distinguished to the north-northeast of XZ Tau, which may be related to the outflow from one member of the XZ Tau binary.
  • A continuum arc from HL Tau to the north and then bending to the east (“cont arc” in Fig. 2) is also identified. This arc may be a dust density discontinuity where enhanced scattering is observed. Although it is not clear if this arc is related to activities at HL Tau or XZ Tau.
  • Another arc feature to the south from HL Tau curving to the southeast can be identified. Two H_2 features are located in the arc and indicated by arrows in Fig. 2. This may be shocked regions in the density discontinuity.
  • Other H_2 features can be distinguished: “A” (interpreted as a limb-brightened edge of the XZ Tau counter-outflow) and “B”, “C”, “a” (blobs driven by the LkH\alpha 358 outflow and interacting with the southern outflow bubble of XZ Tau).

Figure 2: HL/XZ Region (Figure 2 of Hayashi and Pyo 2009)

HH 30 Region

HH 30 is a Herbig-Haro (HH) object including its central star, which is embedded in an almost edge-on flared disk. Although this object doesn’t have clear signs of large-scale [Fe II] or H_2 emission (see Fig. 3), a spectacular jet was detected in the [S II] emission line in previous studies. Despite that, the authors identify two faint small-scale features based on the [Fe II] frame: one to the northeast (corresponding to the brightest part of the [S II] jet) and one to the south-southeast (corresponding to a reflection nebula)


Figure 3: HH 30 Region (Figure 3 of Hayashi and Pyo 2009)

L1551 NE

L1551 NE is a deeply embedded object associated with a fan-shaped infrared reflection nebula opening toward the west-southwest seen in the broad-band Ks continuum emission. It has an opening angle of 60^o. The most important features in this region are:

  • A needle-like feature connecting L1551 NE and HP2 is distinguished from the continuum-substracted [Fe II] image, associated with an [Fe II] jet emanating from L1551 NE.
  • A diffuse red patch at the southwest end of the nebula (denoted as HP1) is dominated by H_2 emission
  • Five isolated compact features are detected in the far-side reflection nebula: HP3 and HP3E ([Fe II] emission), HP4 (both [Fe II] and H_2 emission) and HP5 and HP6 (H_2 emission). All of them are aligned on a straight line that is extrapolated from the jet connecting NE and HP2, naturally assigned to features on the counter-jet.
  • Comparing this data to previous observations in [S II] and radio we can deduce radial velocities of 160-190 km/s for HP2, and 140-190 km/s for HP4 and HP5. With radial velocities in the range 100-130 km/s for these knots, the inclination of the jet axis is estimated to be 45^o60^o.

L1551 IRS-5

L1551 IRS 5 is a protostellar binary system with a spectacular molecular outflow (Snell et al. 1980; see below) and a pair of jets emanating from each of the binary protostars. A conspicuous fan-shaped infrared reflection nebula is seen in Fig. 4, widening from IRS 5 toward the southwest. At the center of this nebula, the two [Fe II] jets appear as two filaments elongated from IRS 5 to its west-southwest; the northern jet is the brighter of the two. Knots A, B and C located farther west and west-southwest of PHK3 (associated with H_2 line emission) have significant [Fe II] emission.


Figure 4: A close-up of IRS-5 (Figure 5 of Hayashi and Pyo 2009)

A counter-jet only seen in the [Fe II] frame can be distinguished to the northeast of IRS 5. Considering its good alignment with the northern jet, it can be interpreted as the receding part of the jet. Based on brightness comparison between the both jets, and transforming H-band extinction to visual extinction the authors deduce a total visual extinction of Av=20-30 mag. Besides the counter-jet, the authors also detect the northern and southern edge of the reflection nebula that delineate the receding-side outflow cone of IRS5.

A brief summary of the HH objects detected in the IRS5 region:

  • HH29: Consistent with a bow shock, its [Fe II] emission features are compact, while the H_2 emission is diffuse. Both emissions are relatively separate.
  • HH260: Consisted with a bow shock with compact [Fe II] emission knot located at the apex of a parabolic H_2 emission feature.
  • HP7: Its [Fe II] and H_2 emission suggest it is also a bow shock driven by an outflow either from L1551 IRS5 or NE.
  • HH264: It is a prominent H_2 emission loop located in the overlapping molecular outflow lobes of L1551 IRS5 and NE. Its velocity gradients are consistent with the slower material surrounding a high-velocity (~ -200 km/s in radial velocity) wind axis from L1551 IRS 5 (or that from L1551 NE)
  • HH 102: Loop feature dominated by H_2 emission (and no [Fe II] emission) similar to HH264. Considering that the major axes of the two elliptical features are consisted with extrapolated axis of the HL Tau jet, it is suggested that they might be holes with wakes on L1551 IRS5 and/or NE outflow lobe(s) that were bored by the collimated flow from HL Tau.

Comparison of Observations

Near-infrared [Fe II] and H_2 emission show different spatial distributions in most of the objects analyzed here. On one hand the [Fe II] emission is confined in narrow jets or relatively compact knots. On the other hand, the H_2 emission is generally diffuse or extended compared with the [Fe II] emission, with none of the H_2 features showing the well collimated morphology as seen in [Fe II].
These differences can be understood based on the conditions that produce different combinations of [Fe II] and H_2 emission:

  • Case of spatially associated [Fe II] and H_2 emissions: Generally requires fast dissociative J shocks (Hollenbach & McKee 1989; Smith 1994; Reipurth et al. 2000).
  • Case of a strong H_2 emission without detectable [Fe II] emission: Better explained by non-dissociative C shocks

The interpretation of differences in [Fe II] and H_2 emission as a result of distinct types shocks is supported by observational evidence showing that the [Fe II] emission usually has a much higher radial velocity than the H_2 emission. In the case of HH 29, HH 260 and HP 7 the [Fe II] emission arises in the bow tips where the shock velocity is fast (~50 km/s) and dissociative whereas H_2 emission occurs along the trailing edges where the shock is slower (~20 km/s)

Paper II: Landmark Observations of Snell et al. 1980

One of the original papers in the study of L1551 was written by Snell, Loren and Plambeck (1980). In this paper, the authors use 12CO to map what they find to be a double-lobed structure extending from the infrared source IRS-5 (see Figures 1, 4). This system is also associated with several Herbig-Haro objects, which are small dense patches that are created in the first few thousands of years after a star is formed. This star is consistent with a B star reddened by 20 magnitudes of dust extinction along the line of sight, through a distance of 160 pc (Snell 1979). By studying these outflows, we are able to better understand the evolution of YSOs.


Snell et al. (1980) made their observations using the 4.9 meter antenna at Millimeter Wave Observatory in Texas. Specifically, they considered the J = 1-0 and J = 2-1 transitions of ^{12}CO and ^{13}CO. Additionally, they made J = 1-0 observations with the NRAO 11 meter antenna. They found asymmetries in the spectral lines, shown below in Figure 5. To the northeast of IRS-5, the high-velocity side of the line has a broad feature, and the southwest of IRS-5 presents a similar broad feature on the low-velocity side of the spectral line. No such features were found to the NW, SE, or in the central position of IRS-5.

Snell et al. 1980, figure 4

Figure 5: 12CO and 13CO 1-0 transition lines; top is NE of central source, bottom is SW of source (Figure 4 of Snell et al. 1980)

The J = 2-1 ^{12}CO transition is enhanced relative to the J = 1-0 transition of ^{12}CO, suggesting that the ^{12}CO emission is not optically thick. If the emission was optically thick, the J = 1-0 line would be the expected dominant transition as it is a lower level transition. The observations also suggest an excitation temperature for the 2-1 transition of T_{ex} ~ 8-35 K. This would only relate to the gas temperature if the environment is in local thermal equilibrium, but it does set a rough minimum temperature. The ^{13}CO emission for the 1-0 transition is roughly 40 times weaker than the same transition for ^{12}CO, which further suggests both isotopes are optically thin in this region (if the ^{12}CO is already optically thin, the weaker transition means ^{13}CO is even more so). The geometry of the asymmetries in the line profiles seen to the NE and SW combined with the distance to L1551 suggest lobes that extend out 0.5 pc in both directions.


Column density

The authors make a rough estimate of the column density of the gas in these broad velocity features by making the following assumptions:

  • the ^{12}CO emission observed is optically thin
  • the excitation temperature is 15 K
  • the ratio of CO to H2 is a constant 5 \times 10^{-5}

With these assumptions, the authors find a column density of 10^{20} cm^{-2}. This is much lower than the region’s extinction measurement of A_{v} = 20 magnitudes by Snell (1979), as the outflow is sweeping out material around the star(s).

Stellar wind and bow shocks

The model of the wind that Snell et al. (1980) suggest is a bimodal wind that sweeps out material in two “bubble-like” lobes, creating a dense shell and possible ionization front that shocks the gas (More on shocks). The physical proximity of the Herbig-Haro (HH) objects in the southwest lobe coming from IRS-5 suggests a causal relationship. Previous work found that the optical spectra of the HH objects resemble spectra expected of shocked gas (Dopita 1978; Raymond 1979).

There is evidence that the CO lobes are the result of a strong stellar wind, the authors clarify this with the schematic shown 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 shell is moving out at that speed. The Herbig-Haro objects HH29 and HH102 have radial velocities of approximately 50 km/s in the same direction as the SW lobe is expanding (Strom, Grasdalen and Strom 1974). Additionally, Cudworth and Herbig (1979) measured the transverse velocities of HH28 and HH29, and found that the objects were moving at a speed of 150 to 170 km/s away from IRS-5. To have reached these velocities, the HH objects must have been accelerated, most likely by a strong stellar wind at speeds above 200 km/s. The bimodal outflow suggests a thick accretion disk around the young star.

Snell et al. 1980, figure 5

Figure 6: Schematic drawing of stellar outflow (Figure 5 of Snell et al. 1980)

Mass-loss rate

The average density in the region away from the central core is 10^{3} {\rm ~cm}^{-3} (Snell 1979), so the extent and density of the shell implies a swept-up mass of 0.3 to 0.7 solar masses. With the measured velocity of ~15 km/s assumed to be constant during the lifetime of the shell and at the measured distance of 0.5 pc from the star, the shell was created 30,000 years ago. With this age, the authors determined a mass loss using the lower end of the assumed swept-up mass and the observed volume of the shell. They found a mass-loss rate of 8 \times 10^{-7} {\rm ~M}_{Sun}{\rm ~yr}^{-1}, which can be compared to other stars using a chart like that shown in Figure 7. This is not meant to present constraints on the processes that produce the mass loss in the IRS-5 system, but rather to simply provide context for the stellar wind observed. The low-mass main sequence star(s) that will eventually arise from the IRS-5 system will be characterized by much lower mass loss rates, and studies of mass loss rates from other YSOs suggest that this source is at the high end of the range of expected rates.

Cranmer - Ay201a - Stellar Winds

Figure 7: A representative plot of the different types of stellar wind, presented by Steve Cranmer in lectures for Ay201a, Fall 2012

Snell et al. (1980) suggest observational tests of this wind-driven shock model:

  • H2 emission from directly behind the shock
  • FIR emission from dust swept up in the shell; this is a possibly significant source of cooling
  • radio emission from the ionized gas in the wind itself near to IRS-5 with the VLA or similar; an upper limit of 21 mJy at 6 cm for this region was determined by Gilmore (1978), which suggests the wind is completely ionized

The results of some newer observations that support this wind model are presented in the following section.

Paper III: A new look at IRS-5 by Wu et al. 2009

Wu et al. (2009) focus on the outflows in L1551 IRS5, the same region studied by Snell et al. (1980), but at a higher angular resolution (~3 arcsec; <1000 AU) and much smaller field of view (~1 arcmin; ~0.05 pc or 10,000 AU). Using the sub-millimeter array (SMA) right after it was formally dedicated, the authors detected CO(2-1) line and millimeter continuum in this low-mass star formation system. The mm continuum, which comes mostly from thermal dust emission, is used to estimate the dust mass. The CO(2-1) spectral line is used to trace the outflows around the binary or triple protostellar system at the center, revealing complex kinematics that suggest the presence of three possible bipolar outflows. The authors construct a cone-shaped outflow cavity model to explain the X-shaped component, and a precessing outflow or a winding outflow due to orbital motion model to explain the S-shaped component. The third component, the compact central one, is interpreted as newly entrained material by high-velocity jets.

Important concepts

There are several concepts related to radio interferometry that merit some discussion:

1. Extended emission filtered out by the interferometer
This is the well-known ‘missing flux problem’ unique in interferometry. There is a maximum scale over which structure cannot be detected by an interferometer, and this scale is set by the minimum projected baseline (i.e. projected distance between a pair of antennas) in the array. In the channel maps (Fig. 3 of Wu et al. 2009), there is a big gap between 5.8 km/s and 7.1 km/s. It does not indicate that there is no CO gas at these velocities, but rather is the result of very extended and homogenous CO distribution which is unfortunately filtered out. This effect applies to all channels.

2. Visibility vs. image
The data directly obtained from the interferometers are called visibility data, in which the amplitude and phase are stored for each baseline. The amplitude, as it literally means, measures the flux; and the phase derives the relative location with respect to the phase center (a reference position on the antenna). We need to convolve the visibility data with a point spreading function, also called ‘dirty beam’, to get the image we need. Mathematically, visibility and image are related through a Fourier transform. For more information, see this online course.

3. Channel maps and P-V diagram
In radio observations, velocity of spectral lines plays an important role by providing kinematic information inside ISM. Normally, the radio data has (at least) three dimensions, two in spatial (e.g. R.A. and Decl.) and one in frequency or velocity. Velocity itself can be used to identity outflows, turbulence or infall by the analysis of line profile, or it can be combined with spatial distribution of emission, if the spatial resolution allowed as in this paper, in the form of channel maps or P-V diagram, to show the three-dimensional structure. In terms of outflows, we expect to see gas at velocities much different from the systematic velocity, and a symmetric pattern both in red- and blue-shifted sides will be even more persuasive. For the efforts in visualizing the 3-d datacube in a fancier way, see this astrobite.

The classical image of low-mass star formation

Fig. 8 shows a diagram of a simplified four-step model of star formation (Fig. 7 of Shu et al. 1987). First, the dense gas collapses to form a core; second, a disk forms because of conservation of angular momentum; third, a pair of outflows emerge along the rotational axis; finally, a stellar system comes into place. During this process, bipolar outflows are naturally formed when the wind breaks through surrounding gas. Therefore, bipolar outflows are useful tools to indirectly probe the properties of protostellar systems.

Shu et al. 1987, figure 7

Figure 8: The formation of a low-mass star (see Shu et al. 1987)

Protostar candidates in L1551 IRS5

Two protostellar components as well as their own ionized jets and circumstellar disks have been found in this source. In addition, Lim & Takakuwa (2006) found a third protostellar candidate, as seen in Fig. 9. In this paper the authors investigated the possible connection between these protostars and the outflows.

Lim and Takakuwa 2006

Figure 9: Two 7mm continuum peaks in the north and south represent the binary system in IRS 5. Arrows show the direction of jets from each of the protostars. A third protostellar candidate is found close to the northern protostar, marked by a white cross. (see Lim and Takakuwa 2006)

Three outflow components

Based on CO(2-1) emission, the authors found three distinct structures. Here the identification was not only based on morphology, but also on the velocity (see Fig. 5 and Fig. 9 of Wu et al. 2009). In other words, it is based on information in the 3-d datacube, as shown in 3 dimensions by the visualization below.


Figure 10: 3-D Datacube of the observations. Arrows mark the outflows identified in this paper. Red/blue colors indicate the red/blue-shifted components. The solid arrows mark the X-shaped component; the broad arrows are the S-shaped component; the narrow arrows are the compact central component. There is an axes indicator at the lower-left corner (x – R.A., y – Decl., z- velocity). (visualization by Walker Lu)

  • The X-shaped component

The first one is an X-shaped structure, with its morphology and velocity shown in the paper. Four arms comprise an hour-glass like structure, with an opening angle of ~90 degree. The northwest and southwest arms are blue-shifted with respect to the systematic velocity, and the northeast and southeast arms are red-shifted. This velocity trend is the same with the large-scale bipolar outflow (Snell et al. 1980, Fig. 6 above). However the two blue-shifted arms, i.e. the NW and SW arms, are far from perfect: the SW arm is barely seen, while the NW arm consists of two components, and presents a different velocity pattern. This component coincides well the U-shaped infrared emission found to the SW of IRS5 (see Hayashi & Pyo 2009, or Fig. 4 above).

Outflow components

Figure 11: Components of the Outflow. Coordinates are offset from the pointing center. (Figure 7 of Wu et al. 2009)

  • The S-shaped component

The second component is an S-shaped structure. It extends along the symmetry axis of the X-shaped component as well as the large scale outflow, but in an opposite orientation. As it literally means, this component is twisted like an ‘S’, although the western arm is not so prominent.

  • The compact central component

The third component is a compact, high-velocity outflow very close to the protostars. The authors fitted a 2-d Gaussian structure to this component and made a integrated blue/red-shifted intensity map, which shows a likely outflow feature in the same projected orientation with the X-shaped component and the large scale outflow.

Modeling the outflows

  • A Cone-shaped cavity model for the X-shaped component

The authors then move on to construct outflow models for these components. For the X-shaped component, a cone-shaped outflow cavity model is proposed (see Fig. 12 and compare with Fig. 11). By carefully selecting the opening angle of the cone and position angle of the axis, plus assuming a Hubble-like radial expansion, this model can reproduce the X-shaped morphology and the velocity pattern. The origin of this cone is related to a high-velocity and well collimated wind, followed by a low-velocity and wide-angle wind that excavates the cone-shaped cavity. Therefore, what we see as X-shaped structure is actually the inner walls of the cavity. However, this model cannot incorporate the NW arm into the picture.

Outflow Models

Figure 12: Models from the paper (Figure 10 of Wu et al. 2009)

  • Entrained material model of the compact central component

For the compact central component, the authors argue that it is material that has been entrained recently by the jets. After comparing the momentum of this component and that of the optical jets, they found that the jets are able to drive this component. Moreover, the X-ray emission around this component indicates strong shocks, which could be produced by the high-velocity jets as well. Finally, the possibility of infalling motion instead of outflow is excluded, because the velocity gradient is larger than the typical value and rotation should dominate over infall at this scale.

  • Three models of the X-shaped component

For the S-shaped component, the authors present the boldest idea in this paper. They present three possible explanations at one go, then analyze their pros and cons in four pages. But before we proceed, do we need to consider other possibilities, such as that this S-shaped component is not interconnected at all, but instead contains two separate outflows from different objects, or that the western arm of the component is not actually part of the outflow but some background or foreground emission? Although the model can reproduce the velocity pattern, it has so many adjustable parameters that we could use it to reproduce anything we like (which reminds me of the story of ‘drawing an elephant with four parameters‘). Anyway, let’s consider the three explanations individually.

  1. Near and far sides of outflow cavity walls? This possibility is excluded because it cannot be incorporated into the cone-shaped cavity model aforementioned, and cannot explain the S-shaped morphology.
  2. A precessing outflow? The outflow should be driven by a jet. Then if the jet is precessing, the outflow will also be twisted. The authors considered two scenarios: a straight jet and a bent jet, and found with finely tuned precessing angle and period, the bent jet model can best reproduce the velocity pattern along the symmetry axis. Therefore, the orbital motion between the third protostellar component, which is thought to be the driving source of this jet, and northern protostellar component, is proposed to cause the precession of the jet thus the outflow. A schematic image is shown in Fig 12.
  3. A winding outflow due to orbital motion? The difference between this explanation and the previous one is that in a precessing outflow the driving source itself is swung by its companion protostar, so the outflow is point symmetric, while in a winding outflow the driving source is unaffected but the outflow is dragged by the gravity of its companion as it proceeds, so it has mirror symmetry with respect to the circumstellar disk plain. Again, if we fine-tune the parameters in this model, we can reproduce the velocity pattern.
S-shaped component

Figure 13: The best-fit model for the S-shaped component, a bent precessing jet. Note the velocity patterns for red/blue lobes are not symmetric (Figure 11 of Wu et al. 2009)

A problem here, however, is although either the precessing jet or the winding outflow model is assumed to be symmetric, the authors use asymmetric velocity patterns to fit the two arms of the S-shaped component (see Fig. 12 and 13 in the paper). In the winding outflow model for instance, in order to best fit the observed velocities, the authors fit the eastern arm starting at a net velocity of 2 km/s at the center, while they fit the western arm starting at ~1.2 km/s. This means the two arms start at different velocities at the center.


The nature of X-shaped and S-shaped structures interpreted in this paper is based on the analysis of kinematics and comparison with toy models. However, the robustness of their conclusion suffers from several questions: for example, how to explain the uniqueness of the NW arm in the X-shaped structure? Is the X-shaped structure really a bipolar outflow system, or just two crossing outflows? Why is the compact central component filtered out around the systematic velocity? Is the S-shaped structure really a twisted outflow, or it is two outflow lobes from two separated protostars?

All these questions might be caused by the missing flux problem discussed above. Observations from a single-dish telescope could be  combined with the interferometric data to: 1) find the front and back walls of the outflow cavity, given sufficient sensitivity, to confirm that the X-shaped component is interconnected; 2) detect the extended structure around the systematic velocity, thus verify the nature of the compact central component; 3) recover at least part of the flux in the SW arm of the X-shaped component and the west arm of the S-shaped component, and better constrain the models.


Using radio and infrared observations, these three papers together provide a integrated view of jets and outflows around YSOs in L1551. The near infrared observations of Hayashi & Pyo (2009) searched for [Fe II] and H2 features introduced by shocks, and found quite different configurations among the YSOs in this region. Some have complicated IR emission, such as HL/XZ Tau, while others like L1551 NE and IRS5 have well-collimated jets traced by [Fe II]. Among them, L1551 IRS5 is particularly interesting because it shows two parallel jets. The pilot work of Snell et al. (1980) revealed a bipolar molecular outflow traced by 12CO from IRS 5, which is interpreted to be created by a strong stellar wind from the young star. High angular resolution observation by Wu et al. (2009) confirms this outflow component, as well as the presence of another two bipolar outflows originating from the binary, or triple system in IRS 5. All these observations show us that jets and outflows are essential in star formation, no only by transporting the angular momentum so that YSOs can continue accreting, but also by stirring up ambient gas and feeding turbulence into the ISM, which might determine the core mass function as mentioned in Alves et al. 2007.


Cudworth, K. M. and Herbig, G. 1979, AJ, 84, 548.
Dopita, M. 1978, ApJ Supp., 37, 117.
Gilmore, W. S. 1978, Ph.D Thesis, University of Maryland
Hayashi, M. and Pyo, T.-S. 2009, ApJ, 694, 582-592.
Hollenback, D. and McKee, C. F. 1989, ApJ, 342, 306-336.
Lim, J. and Takakuwa, S. 2006, ApJ, 653, 1.
Lyons, B. T. 1962, ApJ Supp., 7, 1.
Raymond, J. C. 1979, ApJ Supp., 39, 1.
Reipurth, B., et al. 2000, AJ, 120, 3.
Pineda, J. et al. 2011, ApJ Lett., 739, 1.
Shu, F. H., Adams, F. C., Lizano, S. 1987, ARA&A, 25, 23-81.
Smith, M. D. 1994, A&A, 289, 1.
Snell, R. L. 1979, Ph.D Thesis, University of Texas.
Strom, S. E., Grasdalen, G. L. and Strom, K. M. 1974, ApJ, 191, 111.
Vrba, F. J., Strom, S. E. and Strom, K. M., 1976, AJ, 81, 958.
Wu, P.-F., Takakuwa, S., and Lim, J. 2009, 698, 184-197. Read the rest of this entry »

ARTICLE: A Theory of the Interstellar Medium: Three Components Regulated by Supernova Explosions in an Inhomogeneous Substrate

In Journal Club 2013 on March 15, 2013 at 11:09 pm

Abstract (the paper’s, not ours)

Supernova explosions in a cloudy interstellar medium produce a three-component medium in which a large fraction of the volume is filled with hot, tenuous gas.  In the disk of the galaxy the evolution of supernova remnants is altered by evaporation of cool clouds embedded in the hot medium.  Radiative losses are enhanced by the resulting increase in density and by radiation from the conductive interfaces between clouds and hot gas.  Mass balance (cloud evaporation rate = dense shell formation rate) and energy balance (supernova shock input = radiation loss) determine the density and temperature of the hot medium with (n,T) = (10^{-2.5}, 10^{5.7}) being representative values.  Very small clouds will be rapidly evaporated or swept up.  The outer edges of “standard” clouds ionized by the diffuse UV and soft X-ray backgrounds provide the warm (~10^{4} K) ionized and neutral components.  A self-consistent model of the interstellar medium developed herein accounts for the observed pressure of interstellar clouds, the galactic soft X-ray background, the O VI absorption line observations, the ionization and heating of much of the interstellar medium, and the motions of the clouds.  In the halo of the galaxy, where the clouds are relatively unimportant, we estimate (n,T) = (10^{-3.3}, 10^{6.0}) below one pressure scale height.  Energy input from halo supernovae is probably adequate to drive a galactic wind.

The gist

The paper’s (McKee and Ostriker 1977) main idea is that supernova remnants (SNRs) play an important role in the regulation of the ISM.  Specifically, they argue that these explosions add enough energy that another phase is warranted: a Hot Ionized Medium (HIM)

A Bit About SNRs…

A basic supernova explosion consists of several phases.  Their characteristic energies are on the order of 10^{51} erg, and indeed this is a widely-used unit.  For a fairly well-characterized SNR, see Cas A which exploded in the late 1600s.

Nearby supernova remnant Cassiopeia A, in X-rays from NuSTAR.

Nearby supernova remnant Cassiopeia A, in X-rays from NuSTAR.

  1. Free expansion
    A supernova explosion begins by ejecting mass with a range of velocities, the rms of which is highly supersonic.  This means that a shock wave propagates into the ISM at nearly constant velocity during the beginning.  Eventually the density decreases and the shocked pressure overpowers the thermal pressure in the ejected material, creating a reverse shock propagating inwards.  This phase lasts for something on the order of several hundred years.  Much of the Cas A ejecta is in the free expansion phase, and the reverse shock is currently located at 60% of the outer shock radius.
  2. Sedov-Taylor phase
    The reverse shock eventually reaches the SNR center, the pressure of which is now extremely high compared to its surroundings.  This is called the “blast wave” portion, in which the shock propagates outwards and sweeps up material into the ISM.  The remnant’s time evolution now follows the Sedov-Taylor solution, which finds R_s \propto t^{2/5}.  This phase ends when the radiative losses (from hot gas interior to the shock front) become important.  We expect this phase to last about 10^3 years.
  3. Snowplow phase
    When the age of the SNR approaches the radiative cooling timescale, cooling causes thermal pressure behind the shock to drop, stalling it.  This phase features a shell of cool gas around a hot volume, the mass of which increases as it sweeps up the surrounding gas like a gasplow.  For typical SNRs, this phase ends at an age of about 10^6 yr, leading into the next phase:
  4. Fadeaway
    Eventually the shock speed approaches the sound speed in the gas, and turns into a sound wave.  The “fadeaway time” is on the order of 10^{6} years.

So why are they important?

To constitute an integral part of a model of the ISM, SNRs must occur fairly often and overlap.  In the Milky Way, observations indicate a supernova every 40 years.  Given the size of the disk, this yields a supernova rate of 10^{-13} pc^{-3} yr^{-1}.

Here we get some justification for an ISM that’s a bit more complicated than the then-standard two-phase model (proposed by Field, Goldsmith, and Habing (1969) consisting mostly of warm HI gas).  Taking into account the typical fadeaway time of a supernova, we can calculate that on average 1 other supernova will explode within a “fadeaway volume” within that original lifetime.  That volume is just the characteristic area swept out by the shock front as it approaches the sound speed in the last phase.  For a fadeaway time of 10^6 yr and a typical sound speed of the ISM, this volume is about 100 pc.  Thus in just a few million years, this warm neutral medium will be completely overrun by supernova remnants!  The resulting medium would consist of low-density hot gas and dense shells of cold gas.  McKee and Ostriker saw a better way…

The Three Phase Model

McKee and Ostriker present their model by following the evolution of a supernova remnant, eventually culminating in a consistent picture of the phases of the ISM. Their model consists of a hot ionized medium with cold dense clouds dispersed throughout. The cold dense clouds have surfaces that are heated by hot stars and supernova remnants, making up the warm ionized and neutral media, leaving the unheated interiors as the cold neutral medium. In this picture, supernova remnants are contained by the pressure of the hot ionized medium, and eventually merge with it. In the early phases of their expansion, supernova remnants evaporate the cold clouds, while in the late stages, the supernova remnant material cools by radiative losses and contributes to the mass of cold clouds.

McKee and Ostriker Three Phase Model
A schematic of the three phase model, showing how supernovae drive the evolution of the interstellar medium.

In the early phases of the supernova remnant, McKee and Ostriker focus on the effects of electron-electron thermal conduction. First, they cite arguments by Chevalier (1975) and Solinger, Rappaport, and Buff (1975) that conduction is efficient enough to make the supernova remnant’s interior almost isothermal. Second, they consider conduction between the supernova remnant and cold clouds that it engulfs. Radiative losses from the supernova remnant are negligible in this stage, so the clouds are evaporated and incorporated into the remnant. Considering this cloud evaporation, McKee and Ostriker modify the Sedov-Taylor solution for this stage of expansion, yielding two substages. In the first substage, the remnant has not swept up much mass from the hot ionized medium, so mass gain from evaporated clouds dominates. They show this mechanism actually modifies the Sedov-Taylor solution to a t^{3/5} dependance. In the second substage, the remnant has cooled somewhat, decreasing the cloud evaporation, making mass sweep-up the dominant effect. The classic t^{2/5} Sedov-Taylor solution is recovered.

The transition to the late stages occurs when the remnant has expanded and cooled enough that radiative cooling becomes important. Here, McKee and Ostriker pause to consider the properties of the remnant at this point (using numbers they calculate in later sections): the remnant has an age of 800 kyr, radius of 180 pc, density of 5 \times 10^{-3} cm^{-3}, and temperature of 400 000 K. Then, they consider effects that affect the remnant’s evolution at this stage:

  • When radiative cooling sets in, a cold, dense shell is formed by runaway cooling: in this regime, radiative losses increase as temperature decreases. This effect is important at a cooling radius where the cooling time equals the age of the remnant.
  • When the remnant’s radius is larger than the scale height of the galaxy, it could contribute matter and energy to the halo.
  • When the remnant’s pressure is comparable to the pressure of the hot ionized medium, the remnant has merged with the ISM.
  • If supernovae happen often enough, two supernova remnants could overlap.
  • After the cold shell has developed, when the remnant collides with a cold cloud, it will lose shell material to the cloud.

Frustratingly, they find that these effects become important at about the same remnant radius. However, they find that radiative cooling sets in slightly before the other effects, and continue to follow the remnant’s evolution.

The mean free path of the remnant’s cold shell against cold clouds is very short, making the last effect important once radiative cooling has set in. The shell condenses mass onto the cloud since the cloud is more dense, creating a hole in the shell. The density left behind in the remnant is insufficient to reform the shell around this hole. The radius at which supernova remnants are expected to overlap is about the same as the radius where the remnant is expected to collide with its first cloud after having formed a shell. Then, McKee and Ostriker state that little energy loss occurs when remnants overlap, and so the remnant must merge with the ISM here.

At this point, McKee and Ostriker consider equilibrium in the ISM as a whole to estimate the properties of the hot ionized medium in their model. First, they state that when remnants overlap, they must also be in pressure equilibrium with the hot ionized medium. Second, the remnants have added mass to the hot ionized medium by evaporating clouds and removed mass from the hot ionized medium by forming shells – but there must be a mass balance. This condition implies that the density of the hot ionized medium must be the same as the density of the interior of the remnants on overlap. Third, they state that the supernova injected energy that must be dissipated in order for equilibrium to hold. This energy is lost by radiative cooling, which is possible as long as cooling occurs before remnant overlap. Using supernovae energy and occurrence rate as well as cold cloud size, filling factor, and evaporation rate, they calculate the equilibrium properties of the hot ionized medium. They then continue to calculate “typical” (median) and “average” (mean) properties, using the argument that the hot ionized medium has some volume in equilibrium, and some volume in expanding remnants. They obtain a typical density of 3.5 \times 10^{-3} cm^{-3}, pressure of 5.0 \times 10^{-13} cm^{-2} dyn, and temperature of 460 000 K.

McKee and Ostriker also use their model to predict different properties in the galactic halo. There are fewer clouds, so a remnant off the plane would not gain as much mass from evaporating clouds. Since the remnant is not as dense, radiative cooling sets in later – and in fact, the remnant comes into pressure equilibrium in the halo before cooling sets in. Supernova thus heat the halo, which they predict would dissipate this energy by radiative cooling and a galactic wind.

Finally, McKee and Ostriker find the properties of the cold clouds in their model, starting from assuming a spectrum of cloud sizes. They use Hobbs’s (1974) observations that the number of clouds with certain column density falls with the column density squared, adding an upper mass limit from when the cloud exceeds the Jeans mass and gravitationally collapses. A lower mass limit is added from considering when a cloud would be optically thin to ionizing radiation. Then, they argue that the majority of the ISM’s mass lies in the cold clouds. Then using the mean density of the ISM and the production rate of ionizing radiation, they can find the number density of clouds and how ionized they are.

Parker Instability

The three-phase model gives little prominence to magnetic fields and giant molecular clouds. As a tangent from McKee and Ostriker’s model, the Parker model (Parker 1966) will be presented briefly to showcase the variety of considerations that can go into modelling the ISM.

The primary motivation for Parker’s model are observations (from Faraday rotation) that the magnetic field of the Galaxy is parallel to the Galactic plane. He also assumes that the intergalactic magnetic field is weak compared to the galactic magnetic field: that is, the galactic magnetic field is confined to the galaxy. Then, Parker suggests what is now known as the Parker instability: that instabilities in the magnetic field cause molecular cloud formation.

Parker’s argument relies on the virial theorem: in particular, that thermal pressure and magnetic pressure must be balanced by gravitational attraction. Put another way, field lines must be “weighed down” by the weight of gas they penetrate: if gravity is too weak, the magnetic fields will expand the gas it penetrates. Then, he rules out topologies where all field lines pass through the center of the galaxy and are weighed down only there: the magnetic field would rise rapidly towards the center, disagreeing with many observations. Thus, if the magnetic field is confined to the galaxy, it must be weighed down by gas throughout the disk.

He then considers a part of the disk, and assumes a uniform magnetic field, and shows that it is unstable to transverse waves in the magnetic field. If the magnetic field is perturbed to rise above the galactic plane, the gas it penetrates will slide down the field line towards the disk because of gravity. Then, the field line has less weight at the location of the perturbation, allowing magnetic pressure to grow the perturbation. Using examples of other magnetic field topologies, he argues that this instability is general as long as gravity is the force balancing magnetic pressure. By this instability, he finds that the end state of the gas is in pockets spaced on the order of the galaxy’s scale height. He suggests that this instability explains giant molecular cloud formation. The spacing between giant molecular clouds is of the right order of magnitude. Also, giant molecular clouds are too diffuse to have formed by gravitational collapse, whereas the Parker instability provides a plausible mode of formation.

In today’s perspective, it is thought that the Parker instability is indeed part of giant molecular cloud formation, but it is unclear how important it is. Kim, Ryu, Hong, Lee, and Franco (2004) collected three arguments against Parker instability being the sole cause:

  • The formation time predicted by the Parker instability is ~10 Myr. However, looking at giant molecular clouds as the product of turbulent flows gives very short lifetimes (Ballesteros-Paredes et al. 1999). Also, ~10 Myr post T Tauri stars are not found in giant molecular clouds, suggesting that they are young (Elmegreen 2000).
  • Adding a random component to the galactic magnetic field can stabilize the Parker instability (Parker & Jokipii 2000, Kim & Ryu 2001).
  • Simulations suggest that the density enhancement from Parker instability is too small to explain GMC formation (Kim et al. 1998, 2001).

Does it hold up to observations?

The paper offers several key observations justifying the model.  First, of course, is the observed supernova rate which argues that a warm intercloud medium would self-destruct in a few Myr.  Other model inputs include the energy per supernova, the mean density of the ISM, and the mean production rate of UV photons.

They also cite O VI absorption lines and soft X-ray emission as evidence of the three-phase model.  The observed oxygen line widths are a factor of 4 smaller than what would be expected if they originated in shocks or the Hot Ionized Medium, and they attribute this to the idea that the lines are generated in the conductive surfaces of clouds — a key finding of their model above.  If one observes soft X-ray emission across the sky, a hot component of T ~ 10^{6.3} K can be seen in data at 0.4-0.85 keV, which cannot be well explained just with SNRs of this temperature (due to their small filling factor).  This is interpreted as evidence for large-scale hot gas.

So can it actually predict anything?

Sure!  Most importantly, with just the above inputs — the supernova rate, the energy per supernova, and the cooling function — they are able to derive the mean pressure of the ISM (which they predict to be 3700 K cm^-3, very close to the observed thermal pressures).

Are there any weaknesses?

The most glaring omission of the three-phase model is that the existence of large amounts of warm HI gas, seen through 21cm emission, is not well explained; they underpredict the fraction of hydrogen in this phase by a factor of 15!  In addition, observed cold clouds are not well accounted for; they should disperse very quickly even at temperatures far below that of the ISM that they predict.

CHAPTER: The Sound Speed

In Book Chapter on February 7, 2013 at 10:00 pm

(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.7 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.