# Harvard Astronomy 201b

## Bubbles, Bubbles: Visualization of Milky Way Project and GBT HII Region Survey

In Special Topics Modules on May 8, 2013 at 5:01 pm

If we are lucky to be able to look up on a clear summer night at a place with little light pollution, we would see Milky Way with its faint stripe of stars. If we are a little bit luckier that we can see infrared light with our naked eyes, we would see Milky Way in its hidden and more beautiful form. In infrared, our mother galaxy would reveal itself as a long and broader stripe full of gas and dust, with sometimes complicated structures like twisted filaments. Likely your eyes would be caught by the many spherical shapes with reddish hues surrounded by a round blue circle (here the “red” and “blue” are the two new colors received by your far/near-infrared color receptive cones). These are what astronomers (loosely) call “bubbles” or “shells,” which are indeed almost everywhere across our Milky Way. If you ever participate in the citizen science project, the “Milky Way Project,” you are probably familiar with these bubbles on the infrared maps that pop up on your screen.

## Milky Way Project

Citizen science projects took off with the fast development of internet speed and home computers in 2000s and have since evolved from the rotating numbers and curves in SETI@home to the pretty maps which you can actually work on in Galaxy Zoo. The Milky Way Project is part of the Zooniverse program which first started its crowd sourcing project with Galaxy Zoo in 2007. It now includes two sub-projects: BUBBLES and CLOUDS, where different maps are used to identify, unsurprisingly, bubbles and the so-called infrared dark clouds. In both projects, images taken by NASA Spitzer Space Telescope are used as the maps where citizen scientists look for these structures just as they do when identifying, say, a park on Google Maps.

M20 is a typical bubble in our own galaxy, with a huge HII region related to it. This picture is taken by Spitzer Space Telescope. Image credit: NASA, JPL-Caltech, J. Rho (SSC/Caltech)

BUBBLES, serving as the trailblazer for the Milky Way Project, ask participants to describe the shape and the size of a bubble on Spitzer maps covering the inner part of our galaxy. On the interface, citizen scientists are able to draw a circle on what they think is a bubble and to subsequently change the thickness of the circle, the shape (or, the ellipticity) of the circle and, when needed, to report an incompleteness in the circle, to match the geometry of the bubble. Like most other citizen science projects, BUBBLES propel the participants by giving scores and letting them compete with each other on the quality and the number of bubbles they identify. Like non-citizen science, citizen scientists are encouraged to review their results and discuss with their fellow scientists. If you want to have a taste of being an astronomer (or you ARE an astronomer), BUBBLES also provide advanced tools for citizen scientists to search in astronomy catalogues such as SIMBAD to look for previously known astronomical objects on the maps.

A screen shot of the Milky Way Project: BUBBLES interface. On the left are tools to adjust the shape of the circle.

### For Astronomers

In the Milky Way Project, both Spitzer and Herschel data are being shown to the citizen scientists. In order to achieve the best result from crowd sourcing, different combinations of wavelengths are used respectively in BUBBLES and CLOUDS. While CLOUDS set their goal on identifying infrared dark clouds and use longer-wavelength maps from Herschel, BUBBLES use 4.5/8.0/24-µm images from Spitzer on blue/green/red color maps. Data in BUBBLES are taken from the Galactic Legacy Infrared Mid-Plane Survey Extraordinaire (GLIMPSE) and the Multiband Imaging Photometer for Spizter Galactic Plane Survey (MIPSGAL), from Spitzer Space Telescope, with a coverage of ±65˚ in galactic longitudes and ±1˚ in galactic latitudes.

The first result of the Milky Way Project was released in 2012, with a catalogue of 5,106 bubbles after combing through a total of 520,120 circles drawn by participants. The combination process pushes the user-drawn circles through two sets of 2˚ by 2˚ grids, with one offset by 1˚ from the other. An automatic clustering algorithm would start to group user-drawn circles if there are more than five circles in a 2˚ by 2˚ box, and a “hit rate” is calculated from the number of user-drawn circles over the number of times the map is shown to participants. The offset grids prevent circles on the border of the first set of grids from being left out. The whole process automatically leaves out “lonely circles” which are less likely to be real bubbles.

An example of user drawings and reduced, cleaned circles. The opacities of the cleaned circles are 2 times their hit rates. Simpson et al. 2012

## Bubbles

So what are these bubbles? To answer this question, we have to know what are giving light in these infrared wavelengths first. Now, imagine that you are holding a dust grain in your hand and shrink it tenfold. That tiny little “dust” is what floats between the cold  and vacuous spaces between stars. Just as we see in infrared cameras things that are too cold to emit light in optical wavelengths, we “see” these cold dust grains in infrared. Bubbles we identify on infrared images such as those taken by Spitzer are composed of these dust grains. The shapes of bubbles are the result of “stellar winds,” which are often related to young massive stars that are blowing “winds” and shining strong radiation. The stellar winds push these dust grains outward until there are too many that the force of winds is balanced by the interstellar pressure. What is left behind would be a more vacuous interior with heated up dust and gas, surrounded by a thick cold “ring” or, here we are, “bubbles.”

The reason why we study these bubbles is that they are often related to star formation, of which we still do not have a full picture although a simple scheme can be derived from a few decades of observations. We know that the star formation tends to happen in clusters and that it is affected by forces and radiation from stars. As massive stars live shorter lives, those we see on the sky are young and thus close to where they are born. The stellar winds from these massive stars thus cast strong impacts on the star forming sites and can either trigger or disturb the next generation of star formation (or more often, trigger and disturb at the same time). This makes bubbles interesting places that astronomers want to study, and identifying them would be the first step.

### For Astronomers

The massive stars that push the expansion of bubbles do not only drive momentum-loaded materials but strong radiation as well. All O-type stars and the most massive B-type stars emit photons strong enough to ionize the hydrogen atoms, dissociated by cosmic rays and photons from nearby stars at the surface of denser and thus heavily self-shielded molecular regions. Due to the complicated geometry of stellar clusters, shapes of HII regions are often irregular and can sometimes be filamentary. Bubbles in the Milky Way Projects are not necessarily related to HII regions, and less massive stars and young stellar objects (YSOs) can create bubble-like cavities as well.

When stellar winds from massive stars blow outward in the dust and gas medium, they often move in supersonic speed and thus create discontinuities of physical properties (temperatures, densities, pressures, etc.) at the front of expanding bubbles. These discontinuities are called shocks and provide a genuine environment for characteristic emissions and chemical reactions. Masers and molecular line transitions of certain species like SiO are observed in this environment. People use these emissions to trace strong interaction between stars/YSOs and the interstellar medium to find out what effects this feedback process can have on the subsequent star formation.

## “Bubbles, Bubbles” in World Wide Telescope

This module allows you to play around bubbles identified in the Milky Way Project in World Wide Telescope. By fading between different all-sky maps, you will be able to explore the correlation of emissions at different wavelengths and the bubbles. The complete GLIMPSE+MIPSGAL survey map is readily imported in World Wide Telescope in the same set of color combination as that in the Milky Way Project, so you can get a sense of how people identify bubbles on the exact map shown to them. Other infrared surveys taken by various instruments such as Wide-Field Infrared Survey Explorer (WISE) and Infrared Astronomical Satellite (IRAS) that can be interesting to overlay on are also available in World Wide Telescope. You can also play with the dust map and the H-Alpha emission map, which are also closely related to star formation.

In this module, catalogues taken from HII Region Discovery Survey (HRDS) is used to show positions of galactic HII regions. Although the HRDS catalogue does not intend to be complete, covering only a sub-region of -17˚ to 67˚ in galactic longitudes within ±1˚ in galactic latitudes, it does give a sense of where typical HII regions are distributed. By exploring this catalogue in the module together with the data from Milky Way Project, you can see how these two types of objects are (not) correlated with each other.

Screen shot of the interactive map shown on the HRDS website. When this WWT module is finished, the WWT html5 implementation could provide a better visualization of the HRDS data.

### How to Use Data Tables in WWT

Here are the WWT tour and the excel files that are shown in the tour. To view the tour, simply download the file in the link below. You can open the tour file in both web client and Windows client by going to the pull-down menu under “Explore,” finding “Open…” and then clicking on “Tours.” You are also welcome to download the data tables and explore them on your own in WWT. To view these data tables in WWT, first you have to install Windows client of the latest WWT version on your machine, and go here for WWT Excel Add-in. When you finish the installation, you will see a WWT tab in Excel when you open these tables. Select the columns you want to import into WWT and click on “Visualize Selection” button. After setting the correct labels for each column, you can click on “View in WWT” and start your exploration!

WWT saved view of imported data: red circles indicate bubbles identified in the Milky Way Project, and purple dots indicate positions of HII region observed in HRDS.

### For Astronomers

The HII Region Discovery Survey released the result in a series of five papers. The targets were selected based mainly on previous results from 21-cm HI emission and continuum observation at Very Large Array (VLA), including a piloting effort on radio survey of the galactic plane: the Multi-Array Galactic Plane Imaging Survey (MAGPIS). The 24-µm image from Spitzer MIPSGAL is used as a reference in selecting targets as well. Seven radio recombination lines (RRLs, from H87α to H93α) with frequencies falling in the GBT X-band receiver are observed toward each target. The velocity resolution of such observations can reach ~ 0.4 km/s per channel and thus allows us to derive kinematic distances for these objects.

Distribution of HII RRLs FWHM widths. Anderson et al. 2011

## What You Can Do!

Pick your favorite project from the zoo of citizen science projects and take part! Zooniverse would be a good place to start, with projects in fields ranging from astronomy to meteorology to biology. One advantage of Zooniverse projects over other citizen science projects (besides the pretty pictures) is their well-documented background explanation. By reading through the project introduction, you can pick up what is really intriguing to scientists and understand why people want to study these objects. As the Galaxy Zoo co-founder Kevin Schawinski said, “we prefer to call this [Galaxy Zoo] citizen science because it’s a better description of what you’re doing; you’re a regular citizen but you’re doing science. …You’re pro-actively involved in the process of science by participating.” If you are interested in science but was daunted by equations, give these projects a try. Science is actually a few clicks away!

### For Astronomers

Data sharing is the future of science, and luckily for those who are interested in bubbles, the data presented in this module are all publicly available. The Milky Way Project released its first result of BUBBLES in 2012, which includes one catalogue of “large bubbles” and one of “small bubbles.” The “large” and “small” here indicate whether or not the bubbles/”knots” are large enough to be simulated by at least the smallest possible circular shapes available on the identification interface. In the catalogues, columns of galactic coordinates and geometric parameters are included together with “hit rate” and “dispersion” to indicate the reliabilities of identified bubbles.

The HII Region Discovery Survey data are available via the project website. The Arecibo and GBT parts of the survey are separated in two independent files. Correlation can be found using the source id in the catalogues. Here in this module only the GBT part is presented as the derived kinematic distances are available in the GBT catalogue but not the Arecibo one, due to the fact that Arecibo did not bear the capability to resolve the confusion of overlapped velocity components toward inner Milky Way with the settings of this survey. Besides coordinates, velocities and velocity-based distances are available in the catalogues, and those who are interested can refer to their papers for details about the line surveys.

## ARTICLE: Interpreting Spectral Energy Distributions from Young Stellar Objects

In Journal Club, Journal Club 2013 on May 3, 2013 at 8:02 pm

*** Still a draft! ***

Posted by: Meredith MacGregor (4/18/13)

#### 1. Introduction

When discussing young stellar objects and protoplanetary disks, astronomers continually throw out the term ‘SED.’ For the past week, I have been at a conference titled ‘Transformational Science with ALMA: From Dust to Rocks to Planets– Formation and Evolution of Planetary Systems.’ And, I can attest to the fact that the term ‘SED’ has come up in a very significant fraction of the contributed talks. For those not intimately familiar with this sub-field, this rampant use of abbreviations can be confusing, making it difficult to glean any useful take-aways from a talk. So, in addition to summarizing the article by Robitaille et al. (2007), the goal of this post is to give a bit of an introduction to the terminology and motivation for the growing field of star and planetary system formation.

SEDs: What They Are and Why We Care So Much

The abbreviation ‘SED’ stands for ‘Spectral Energy Distribution.’ If you want to sound like an expert, this should be pronounced exactly as it appears (analogous to the non-acronym ‘said’). A SED is essentially a graph of how the flux of emission from an object depends on the wavelength the object is observed at– flux plotted versus wavelength. In the context of the Robitaille et al. article, we are most interested in the SEDs for young stars and the envelopes and disks surrounding them. So, why exactly does the flux from a young stellar object (YSO) vary with wavelength? As it turns out, different regions of the YSO emit at different wavelengths. This means that when we observe at different wavelengths, we are actually probing distinct regions of the star, envelope, disk, etc. By tracing out the entire SED for a YSO, we can actually determine what the geometry, structure, and constituents are for that object.

SEDs for pre-main sequence stars are often seen to have a bump or excess in the infrared. This bump is generally interpreted as being due to thermal emission from warm dust and gas surrounding the central star. As an illustrative example, let’s consider a protoplanetary disk around a young star. The image below is taken from a review article by Dullemond & Monnier (2010) and shows a graphical representation of the emission spectrum from such a disk.

In this case, emission in the Near-IR traces the warmer inner regions of the disk. As you move into submillimeter wavelengths, you start to probe the outer, cooler regions of the disk. By modeling the observed SED of a pre-main sequence star, you can determine what components of the source are contributing to the observed flux. The second figure below is taken from a paper by Guarcello et al. (2010).   The left panel of this figure shows the observed SED of a BWE star (in line with other ridiculous astronomy acronyms, this stands from ‘blue stars with excesses’), while the right panel shows a model fit to this data. The authors describe the observed SED using a model with four distinct components: (1) emission from the reddened photosphere of the central star, (2) radiation scattered into the line of sight from dust grains in the disk, (3) emission from a collapsing outer envelope, and (4) thermal emission from a circumstellar disk. The summation of these four component makes up the complete SED for this BWE.

Describing and Classifying the Evolution of a Protostar

As a protostar evolves towards the Zero Age Main Sequence (ZAMS), the system geometry and, as a result, the SED will evolve as well. Thus, the stage of evolution of such a protostar is often classified according to features seen in the SED. A graphical overview of the four stages of protostar evolution are shown below (Andrea Isella’s thesis, 2006). Class 0 objects are characterized by a very embedded central core in a much larger accreting envelope. The mass of the central core grows in Class I objects and a flattened circumstellar accretion disk develops. For Class II objects, the majority of circumstellar material is now found in a disk of gas and dust. Finally, for Class III objects, the emission from the disk becomes negligible and the SED resembles a pure stellar photosphere.

#### 2. The Article Itself

One common method of fitting SEDs is to assume a given gas and circumstellar dust geometry and set of dust parameters and then use radiative transfer models to find a set of parameters that best reproduce the observations.  However, fitting SEDs by trial and error is a time consuming way to explore a large parameter space.  The problem is even worse if you want to consider thousands of sources.  So, what’s to be done? Enter Robitaille et al.  In order to attempt to make SED fitting more efficient, they have pre-calculated a large number of radiative transfer models that cover a reasonable amount of parameter space.  Then, for any given source, they can quickly find all of the combinations of parameters that could explain the observations.

Let’s Get Technical

The online version of this fitting tool draws from 20,000 sets of physical parameters and 10 viewing angles (if you are particularly curious, the online tool is available here).  A brief overview of the parameter space covered is as follows:

• Stellar mass between 0.1 and 50 solar masses
• Stellar ages between 10^3 and 10^7 years
• Stellar radii and temperatures from evolutionary tracks given mass
• Disk and envelope parameters sampled given age of source

However, there are truly a vast number of parameters that could be varied in models of YSOs.  Thus, for simplicity, the authors are forced to make a number of assumptions.  What are some of the biggest assumptions involved?

1. All stars form via accretion through a disk and an envelope.
2. The gas-to-dust ratio in the disk is 100.
3. The apparent size of the source is not larger than the given aperture.

The last constraint is not required, but it allows a number of model SEDs to be cut out and thus speeds up the process.  Furthermore, the authors make a point of saying that the results can always be scaled to account for varying gas-to-dust ratios, since only the dust is taken into account in the actual radiative transfer calculations.

Does This Method Really Work?

If this tool works well, it should be able to correctly reproduce previous results.  In order to test this out, the authors turn to the Taurus-Auriga star forming region.  They select a sample of 30 sources from Kenyon & Hartman (1995) that are spatially resolved, meaning that there is prior knowledge of their evolutionary stage from direct observations (i.e. there is a known result that astronomers are fairly certain of to compare the model fits against).  In the context of this particular star forming region, the authors throw in a few additional assumptions:

1. All sources are within a distance range of 120 – 160 AU.
2. The foreground interstellar extinction is no more that 20.

The authors then assign an arbitrary cut-off in chi-squared for the best-fit model.  They even say themselves that this is arbitrary: ‘Athough this cut-off is arbitrary, it provides a range of acceptable fits to the eye.’  After taking Jim Moran’s Noise and Data Analysis class, I for one would like to see the authors try a Monte Carlo Markov Chain (MCMC) analysis of their 14-dimensional space.  That might make the analysis a bit less ‘by eye’ and more ‘by statistics.’

The upshot of this study is that for the vast majority of the sources considered, the best-fit values obtained by this new SED fitting tool are close to the previously known values. Check.

It is also worth mentioning here, that there are many other sets of SED modeling codes.  One set of codes of particular note are those written by Paola D’Alessio (D’Alessio et al., 1998; D’Alessio et al., 1999; D’Alessio et al, 2001).  These codes were the most frequently used in the results presented at the ALMA conference I attended.  The distinct change in the D’Alessio models is that they solve for the detailed hydrostatic vertical disk structure in order to account for observations of ‘flared’ disks around T Tauri stars (flaring refers to an increase in disk thickness at larger radii).

But, Wait! There are Caveats!

Although the overall conclusion is that this method fits SEDs with reasonable accuracy, there are a number of caveats that are raised.  First of all, the models tend to overestimate the mid-IR fluxes for DM Tau and GM Aur (two sources known to have inner regions cleared of dust).  The authors explain that this is most likely due to the fact that their models for central holes assume that there is no dust remaining in the hole.  In reality, there is most likely a small amount of dust that remains.  Second, the models do not currently account for the possibility of young binary systems and circumbinary disks (relevant for CoKu Tau 1).

The paper also addresses estimating parameters such as stellar temperature, disk mass, and accretion rate from SED fits.  And, yes, you guessed it, these calculations raise several more issues.  For very young objects, it is difficult to disentangle the envelope and disk, making it very challenging to estimate a total disk mass.  To make these complications clearer, the set of two plots below from the paper show calculated values for the disk mass plotted against the accepted values from the literature.  It is easily seen that the disk masses for the embedded sources are the most dissimilar from the literature values.

Furthermore, even if the disk can be isolated, the dust mass in the disk is affected by the choice of dust opacity.  That’s a pretty big caveat!  A whole debate was started at the ALMA conference over exactly this issue and the authors have swept it under the rug in just one sentence.  In addition, the calculated accretion rates from SED fitting are systematically larger than what is presented in the literature.  The authors conclude that future models should include disk emission inside the dust destruction radius, the radius inside which it is too hot for dust to survive.  A great example of the complications that arise from a disk with a central hole can be seen in LkCa 15 (Espaillat et al., 2010Andrews et al., 2011).   The figure below shows the observed and simulated SEDs for the source (left) as well as the millimeter image (right).  The double ansae seen in the millimeter contours are indicative of a disk with a central cavity.

In this case, a population of sub-micron sized dust within the hole is needed in order to produce the observed silicate feature at 10 microns.  Furthermore, an inner ring is required to produce the strong near-IR excess shortward of 10 microns.  A cartoon image of the predicted disk geometry is shown below.  To make things even more complicated, the flux at shorter wavelengths appears to vary inversely with the flux at longer wavelengths over time (Espaillat et al., 2010).  This phenomenon is explained by changing the height of the inner disk wall over time.

Finally, the authors discuss how well parameters are constrained given different combinations of data points for two example sources: AA Tau and IRAS 04361+2547.  In both sources, if only IRAC (Infrared Array Camera on the Spitzer Space Telescope) fluxes are used, the stellar mass, stellar temperature, disk mass, disk accretion rate, and envelope accretion rate are all poorly constrained.  Things are particularly bad for AA Tau in this scenario, where only using IRAC data results in ~5% of all SED models meeting the imposed goodness of fit criterion (yikes!).  Adding in optical data to the mix helps to rule out models that have low central source temperatures and disk accretion rates.  Adding data at wavelengths longer than ~ 20 microns helps to constrain the evolutionary stage of the YSO, because that is where any infrared excess is most apparent.  And, adding submillimeter data helps to pin down the disk mass, since the emission at these longer wavelengths is dominated by the dust.  This just goes to show how necessary it is to obtain multi-wavelength data is we really want to understand YSOs, disks, and the like.

#### 3. Sources:

Where to look if you want to read more about anything mentioned here…

## ARTICLE: Stellar feedback in galaxies and the origin of galaxy-scale winds (Hopkins et al. 2012)

In Journal Club 2013 on April 27, 2013 at 4:21 pm

Summary by Kate Alexander

Link to paper: Hopkins et al. (2012)

## Abstract

Feedback from massive stars is believed to play a critical role in driving galactic super-winds that enrich the intergalactic medium and shape the galaxy mass function, mass–metallicity relation and other global galaxy properties. In previous papers, we have introduced new numerical methods for implementing stellar feedback on sub-giant molecular cloud (sub-GMC) through galactic scales in numerical simulations of galaxies; the key physical processes include radiation pressure in the ultraviolet through infrared, supernovae (Type I and Type II), stellar winds (‘fast’ O star through ‘slow’ asymptotic giant branch winds), and HII photoionization. Here, we showthat these feedback mechanisms drive galactic winds with outflowrates as high as ∼10–20 times the galaxy star formation rate. The mass-loading efficiency (wind mass-loss rate divided by the star formation rate) scales roughly as $\dot{M}_{wind}/\dot{M}_* \propto V_c^{-1}$ (where $V_c$ is the galaxy circular velocity), consistent with simple momentum-conservation expectations. We use our suite of simulations to study the relative contribution of each feedback mechanism to the generation of galactic winds in a range of galaxy models, from Small Magellanic Cloud-like dwarfs and Milky Way (MW) analogues to z ∼ 2 clumpy discs. In massive, gas-rich systems (local starbursts and high-z galaxies), radiation pressure dominates the wind generation. By contrast, for MW-like spirals and dwarf galaxies the gas densities are much lower and sources of shock-heated gas such as supernovae and stellar winds dominate the production of large-scale outflows. In all of our models, however, the winds have a complex multiphase structure that depends on the interaction between multiple feedback mechanisms operating on different spatial scales and time-scales: any single feedback mechanism fails to reproduce the winds observed.We use our simulations to provide fitting functions to the wind mass loading and velocities as a function of galaxy properties, for use in cosmological simulations and semi- analytic models. These differ from typically adopted formulae with an explicit dependence on the gas surface density that can be very important in both low-density dwarf galaxies and high-density gas-rich galaxies.

## Introduction

Galaxy evolution cannot be properly understood without accounting for strong feedback from massive stars. Specifically, in cosmological models that don’t include feedback processes, the star formation rates in simulated galaxies are much too high, as gas quickly cools and collapses. Additionally, these simulations find that the total amount of gas present in galactic disks is too high. Both of these problems can be solved by including local outflows and galactic superwinds that remove baryons from the disks, slowing star formation and bringing simulations in line with observations. Such winds are difficult to include in simulations, however, because they have their origins in stellar feedback processes, which occur on small scales. Most simulations are either too low resolution to properly model these processes, or they make simplifying assumptions about the physics that prevent accurate modeling of winds. Thus, although we have seen observational evidence of such winds in real galaxies (for example, Coil et al. 2011; Hall et al. 2012), until recently simulations have not been able to generate galactic winds from first principles and have instead added them in manually. Hopkins, Quataert, and Murray for the first time present a series of numerical simulations that successfully reproduce galactic winds that are consistent with observations for a wide range of galaxy types. Unlike previous work, their simulations have both sufficiently high resolution to focus on small-scale processes in giant molecular clouds (GMCs) and star forming regions and the physics to account for multiple types of stellar feedback, not just thermal heating from supernovae. These simulations are also discussed in two companion papers (Hopkins et al. 2011 and Hopkins et al. 2012), which focus on the star formation histories and properties of the galactic ISM of simulated galaxies and outline rigorous numerical tests of the models. The 2011 paper was discussed in the 2011 Ay 201b journal club and is summarized nicely here.

## Key Points

1. Simulations designed to study stellar feedback processes have, for the first time, succeeded in reproducing galactic winds capable of removing material from galaxies at several times the star formation rate when multiple feedback mechanisms are included. They also reproduce the observed inverse scaling of wind mass loading with galactic circular velocity, $\dot{M}_{wind}/\dot{M}_* \propto V_c^{-1}$.
2. Radiation pressure is the primary mechanism for the generation of winds in massive, gas-rich galaxies like local starburst galaxies and high redshift galaxies, while supernovae and stellar winds that shock-heat gas are more important in less gas-rich Milky Way-like galaxies and dwarf galaxies.
3. The wind mass loading and velocity are shown to depend on the gas surface density, an effect which has not previously been quantified.

## Models and Methods

The authors used the parallel TreeSPH code GADGET-3 (Springel 2005) to perform their simulations. The simulations include stars, gas, and dark matter and accounts for cooling, star formation, and stellar feedback. The types of stellar feedback mechanisms they include are local deposition of momentum from radiation pressure, supernovae, and stellar winds; long-range radiation pressure from photons that escape star forming regions; shock heating from supernovae and stellar winds; gas recycling; and photoheating of HII regions. The populations of young, massive stars responsible for most of these feedback mechanisms are evolved using standard stellar population models.

These feedback mechanisms are considered for four different standard galaxy models, each containing a bulge, a disk consisting of stars and gas, and a dark matter halo. These four models are:

1. HiZ: a massive, starburst galaxy at a redshift of 2, with properties chosen to resemble those of non-merging submillimeter galaxies.
2. Sbc: a gas-rich spiral galaxy, with properties chosen to resemble those of luminous infrared galaxies (LIRGs).
3. MW: a Milky Way-like spiral galaxy.
4. SMC: a dwarf galaxy, with properties similar to those of the Small Magellanic Cloud.

Simulations were run for each of these models at a range of resolutions (ranging from 10 pc to sub-pc smoothing lengths) to ensure numerical convergence before settling on a standard resolution. The standard simulations include about $10^8$ particles with masses of $500M_{\odot}$ and have smoothing lengths of about 1-5 pc. (For more details, see the companion papers Hopkins et al. 2011 and Hopkins et al. 2012 or the appendix of this paper). The authors then ran a series of simulations with one or more feedback mechanisms turned off, to assess the relative importance of each mechanism to the properties of the winds generated in the standard model containing all of the feedback mechanisms.

## Results

When all feedback mechanisms are included, the simulations produce the galaxy morphologies seen below. The paper also considers what each galaxy would look like in the X-ray, which traces the thermal emission from the outflows.

The range of morphologies produced with all feedback mechanisms active for the four different model galaxies studied (from left to right: HiZ, Sbc, MW, and SMC). The top two rows show mock images of the galaxies in visible light and the bottom two rows show the distribution of gas at different temperatures (blue=cold molecular gas, pink=warm ionized gas, yellow=hot X-ray emitting gas). Taken from figure 1 of the paper.

## Wind properties and dependence on different feedback mechanisms

As shown above, all four galaxy models have clear outflows when all of these feedback mechanisms are included. When individual feedback mechanisms are turned off to study the relative importance of each mechanism in the different models, the strength of the outflows diminishes. For the HiZ and Sbc models, radiation pressure is shown to be the most important contributing process, while for the MW and SMC models gas heating (from supernovae, stellar wind shock heating, and HII photoionization heating) is more important. The winds are found to consist of mostly (mostly ionized) warm (2000 < T < 400000) and (diffuse) hot (T > 400000) gas, with small amounts of (mostly molecular) colder gas (T < 2000K). Particles in the wind have a range of velocities, differing from simple simulations that often assume a wind with a single, constant speed.

For the purpose of studying galaxy formation, the most important property of the wind is the total mass outflow rate, $\dot{M}$. This is often expressed in terms of the galactic wind mass-loading efficiency, defined as $M_{wind}/M_{new} = \int{\dot{M}_{wind}}/\int{\dot{M}_*}$, where $\dot{M}_{wind}$ is the wind outflow rate and $\dot{M}_*$ is the star formation rate. The galactic mass-loading efficiency for each galaxy model is shown below. By comparing the mass-loading efficiency produced by simulations with all feedback mechanisms turned on (the “standard model”) to simulations with some feedback mechanisms turned off, the importance of each mechanism becomes clear. While the standard model cannot be replicated without all of the feedback mechanisms turned on, radiation pressure is clearly much more important than heating for the HiZ case and less important than heating for the MW and SMC cases. The Sbc case is intermediate, with radiation pressure and heating being of comparable importance.

Galactic wind-mass loading efficiency for each of the four galaxy models studied, taken from figure 8 of the paper.

Derivation of a new model for predicting the wind efficiency of a galaxy

After doing some basic plotting of the wind mass-loading efficiency versus global properties of the galaxy models studied (such as star formation rate and total galaxy mass), the authors explore whether there exists a better model for predicting what the wind mass-loading efficiency should be for a given galaxy. After studying the relations between the wind mass loss rate $\dot{M}_{wind}$ and a range of galaxy properties as a function of radius R, time t, and model type, they conclude that the mass loss rate is most directly dependent on the star formation rate $\dot{M}_*$, the circular velocity of the galaxy $V_c(R)$, and the gas surface density $\Sigma_{gas}$. They find that the mass-loading efficiency can be described by:

$\left<\frac{\dot{M}_{wind}}{\dot{M}_*}\right>_R \approx 10\eta_1\left(\frac{V_c(R)}{100 \text{km/s}}\right)^{-(1+\eta_2)}\left(\frac{\Sigma_{gas}(R)}{10 M_{\odot}\text{pc}^{-2}}\right)^{-(0.5+\eta_3)}$

where $\eta_1\sim 0.7-1.5, \eta_2\sim\pm0.3$, and $\eta_3\sim\pm0.15$ are the uncertainties from the fits of individual simulated galaxies to the model. This relationship is plotted below along with instantaneous and time-averaged properties of simulated galaxies. The dependence of the wind mass loss rate on the star formation rate and the circular velocity of the galaxy match previous results and are easily understandable in terms of conservation of momentum, but the dependence on the surface density of the gas initially seems more surprising. Hopkins et al. posit that this is due to the effects of density on supernovae remnants: for low-density galaxies the expanding hot gas from the supernova will sweep up material with little resistance, increasing its momentum over time, while for high-density galaxies radiative cooling of this gas becomes more important, so it will impart less momentum to swept up material. Therefore supernovae in denser environments contribute less to the wind, all other factors being equal, introducing a dependence of the wind mass loss rate on gas surface density.

## Discussion and Caveats

The results from this paper are fairly robust, as the detailed treatment of multiple feedback mechanisms allows the authors to avoid making some of the simplifying assumptions that are often necessary in galaxy simulations (artificially turning off cooling to slow star formation rates, etc). The combination of high resolution simulations and more realistic physics does a good job of confirming previous numerical work and observational results. Needless to say, however, there is still room for improvement.

One major caveat of these results is that all of the model galaxies are assumed to exist in isolation, with no surrounding intergalactic medium (IGM). In reality, galactic outflows will interact with the IGM and hot coronal gas already present in a halo around the galaxy, which affects the structure of the wind. Additionally, feedback effects from black holes and AGN are not discussed, nor are galactic inflows. Comparisons between observations and simulations have shown that AGN-driven winds cannot alone explain the observed star formation rates in real galaxies (Keres et al. 2009), but they may still be an important contributing factor.

Furthermore, the authors note that their method of modeling radiation pressure is “quite approximate” and could be improved. Cosmic rays are not included and the scattering and absorption of UV and IR photons has been simplified. Computational limits (unresolvable processes) also place constraints on the robustness of the results.

Many of the quantities discussed here are easily derivable from high-resolution simulations, but harder to estimate from observations of real galaxies or simulations that have lower resolution. A good discussion of how simulations compare to the observed galaxy population can be found in Keres et al. 2009. Measurements of hydrogen-alpha emission in galaxies can be used to infer their star formation rate and measurements of their X-ray halos can be used to infer the mass-loss rate from galactic winds, but this requires high-quality observational data that becomes increasingly difficult to capture for galaxies at non-zero redshift (Martin 1998). Depending on the resolution of a simulation or telescope, determining quantities like the galactic rotation curve and gas surface density may not be directly possible. When seeking to apply these results to understand the formation history of galaxies in observational data, these limitations should be taken into account.

## References and Further Reading

Hopkins, P. F., Quataert, E., & Murray, N. 2012, MNRAS, 421, 3522

Keres, D. et al., 2009, MNRAS, 396, 2332

Martin, C. L., 1998, ApJ, 513, 156

## Details on radio observations

In Uncategorized on April 26, 2013 at 9:18 pm

How were the observations done?

If you are a theorist, you probably don’t like words like “bandpass”, and hate words like “calibration.”  If you are an observer, you probably know and love these words.  But even theorists have to understand how observations are done: they’re what keep science science.  Here I want to briefly explain some aspects of the paper’s discussion of how the data were taken that were on first read opaque.

First of all, given that the $NH_3$ rest frequency is roughly 24 GHz, we can find the energy of the transition: $E=h\nu = 10^{-4}\;\rm{eV}$.  This corresponds to an excitation temperature of 1.2 K, which we can also translate into a velocity via $v\sim (E/m)^{1/2} \sim 24\; {\rm m/s}$.  A typical molecular cloud temperature is on the order of 20 K, high enough that there should be many molecules in the upper energy state (think of it as a toy-model 2 level system) and hence many upper to lower transitions occurring (leading to emission at 24 GHz).

A 1 MHz window means the observations were done between 24 GHz – 500 MHz and 24 GHz + 500 MHz, and dividing this interval by 256 gives the 3.9 kHz per channel quoted.   This can be translated into a velocity by thinking of it using the small-z redshift formula $z\approx v/c$. Here, we have a fractional frequency shift $\Delta \nu / \nu \simeq v/c$, and setting $\Delta \nu / \nu$ =  3.9  kHZ/24 GHz leads to the velocity quoted (.049 km/s).

The baselines given (35 meters up to 1 km) determine the spatial resolution: just good old Carroll and Ostlie diffraction limit, $\theta = 1.22 \lambda/D$, where $\theta$ is the angular resolution, $\lambda$ the wavelength of light, and $D$ the diameter of the aperture (here, the baseline).  24 GHz has $\lambda = 1.2$ cm, so with a 1 km baseline the angular resolution is about .05’’ (for comparison, this is 1.3 times the angular resolution of the Hubble Space Telescope at 500 nm).

Nepers is a unit of optical depth; .0689 is optically thin (and measures the optical depth due to the Earth’s atmosphere, since the authors note this corresponds to fair weather; there is negligible optical depth from sources in space at this frequency).  It is interesting to note that the optical depth is given at 22 GHz, but the observations are made at 24 GHz; this is because 22 GHz is a standard fiducial frequency to allow comparison with the weather conditions under which other radio observations might have been made.

Now, calibration—theorists’ bane.  The paper’s discussion is compact and therefore somewhat confusing to those uninitiated into the secrets of radio astronomy, but understanding it is worthwhile for the insight it can yield into what actually goes on at these remote “large arrays” in the desert (alien storage?)

Possible very large array user? Image credit: http://www.123rf.com

In a radio interferometric observation, what you actually measure is the amplitude and phase of the signal in each channel (here, there are 256): i.e., at each different frequency corresponding to the channels in the window about 24 GHz.  (Note: interferometric means that phases are being compared between signal received at spatially separated points on Earth: this allows the difference in path length from the source to two antennae to be computed, and hence the direction to the source). The amplitude and phase at these different frequencies is desirable because it gives (eventually) intensity as a function of frequency, or equivalently a spectral energy distribution (SED).  This in turn allows inferences about the properties of the source (e.g. dust grain size from SEDs of debris disks, cf. problem set 2, problem 4 ($\beta$ describes the shape of the SED)).There will be some intrinsic noise in both amplitude and phase, which you want to eliminate.  For a source that has a flat spectrum, you know any bumps you see at a particular frequency are due to noise in the channel at that frequency.  This is the purpose of having an amplitude calibrator: typically quasars are used because they have flat spectra.

Now, the amplitude is measured by the amount of electrical current coming from a given channel: the more photons at that frequency hit the antenna, the more current.  But astronomers care about flux, not current.  The absolute flux calibrator is a source with a known flux. Measuring the voltage this source produces allows the conversion factor between voltage and flux to be deduced.  Using a source with known flux also allows one to calibrate the bandpass: the flux should be zero outside the window, but to know this is because the window is working, you have to be sure the flux is not just coincidentally zero at frequencies outside the window anyway.

Finally, the astute reader may notice that noise is reported as mJy/beam.  The noise should be linear in the area of the beam, so a larger beam would have more noise.  Thus two sets of observations done with different beam sizes would have different noise levels due to this; reporting noise/beam allows direct comparison of how good observations done with different beam sizes are.

## ARTICLE: On the Density of Neutral Hydrogen in Intergalactic Space

In Journal Club, Journal Club 2013 on April 20, 2013 at 12:25 am

Author: Yuan-Sen Ting, April 2013

In the remembrance of Boston bombing victims: M. Richard, S. Collier, L. Lu, K. Campell

“Weeds do not easily grow in a field planted with vegetable.
Evil does not easily arise in a heart filled with goodness.”

Dharma Master Cheng Yen.

1. For more information and interactive demonstrations of the concepts discussed here, see the interactive online software module I developed for the Harvard AY201b course.

### Introduction

Although this seminal paper by Gunn & Peterson (1965) comprises only four pages (not to mention, it is in single column and has a large margin!), the authors suggested three ideas that are still being actively researched by astronomers today, namely

1. Lyman alpha (Lyα) forest
2. Gunn-Peterson trough
3. Cosmic Microwave Background (CMB) polarization

Admittedly, they put them in a slightly different context, and the current methods on these topics are much more advanced than what they studied, but one could not ask for more in a four pages note! In the following, we will give a brief overview on these topics, the initial discussion in Gunn & Peterson (1965) and the current thinking/results on these topics. Most of the results are collated from the few review articles/papers in the reference list and we will refer them herein.

### Lyman alpha forest

Gunn and Peterson propose using Lyα absorption in the distant quasar spectrum to study the neutral hydrogen in the intergalactic medium (IGM). The quasar acts as the light source, like shining a flashlight across the Universe, and provides a characteristic, relatively smooth spectrum. During the transmission of the light from the quasar to us, intervening neutral hydrogen clouds, along the line of sight to the quasar, cause absorptions on the quasar continuum at the neutral hydrogen HI transition 1215.67Å in the rest frame of the absorbing clouds. Because the characteristic quasar spectrum is well-understood, it’s relatively easy to isolate and study this absorption.

More importantly, in our expanding Universe, the spectrum is continuously redshifted. Therefore the intervening neutral hydrogen at different redshifts will produce absorption at many different wavelengths blueward of the quasars’ Lyα emission. After all these absorptions, the quasar spectrum will have a “forest” of lines (see figure below), hence the name Lyα forest.

Figure 1. Animation showing a quasar spectrum being continuously redshifted due to cosmic expansion, with various Lyman absorbers acting at different redshifts. From the interactive module.

For a quasar at high redshift, along the line of sight, we are sampling a significant fraction of the Hubble time. Therefore, by studying the Lyα forests, we are essentially reconstructing the history of structure formation. To put our discussion in context, there are two types of Lyα absorbers:

1. Low column density absorbers: These absorbers are believed to be the sheetlike-filamentary neutral hydrogen distributed in the web of dark matter formed in the ΛCDM. In this cosmic web, the dark matter filaments are sufficiently shallow to avoid gravitational collapse of the gas but deep enough to prevent these clouds of the inter- and proto-galactic medium from dispersing.
2. High column density absorbers: these are mainly due to intervening galaxies/proto-galaxies such as the damped Lyα systems. As opposed to the low column density absorbers, since galaxies are generally more metal rich than the filamentary gas; the absorption of Lyα due to these systems will usually be accompanied by some characteristic metal absorption lines. There will also be an obvious damping wing visible in the Voigt profile of the Lyα line, hence the name damped Lyα system.

Since our discussion here focuses on the IGM, although the high column density absorbers are extremely important to detect high redshift galaxies, for the rest of the discussion in this post, we will only discuss the former.

### Low column density absorbers

So, what can we learn about neutral hydrogen in the Universe by looking at Lyα forests? With a large sample of high redshift quasar spectra, we can bin up the absorption lines due to intervening clouds of the same redshift. This allows us to study the properties of gas at that distance (at that time in the Universe’s history), and learn about the evolution of the IGM. In this section, we summarize some of the crucial properties that are relevant in our later discussion:

Figure 2. Animation illustrating the effect of the redshift index on the quasar absorption spectrum. A higher index means more closely spaced absorption near the emission peak. From the interactive module.

It turns out that the properties of the IGM can be rather well described by power laws. If we define n(z) as the number density of clouds as a function of redshift z, and N(ρ) as the number density of clouds as a function of column density ρ. It has been observed that [see Rauch (1998) for details],



Note that, in order to study these distributions, ideally, high resolution spectroscopy is needed (FWHM < 25 kms-1) such that each line is resolved.

There is an upturn in the redshift power law index at high redshift. This turns out to be relevant in the discussion of the Gunn-Peterson trough. We will discuss this later. It has also been shown, via cosmological hydrodynamic simulations, that the ΛCDM model is able to reproduce to the observed power laws.

But we can learn about more than just gas density. The Doppler widths of Lyα line are consistent with the photoionization temperature (about 104 K). This suggests that the neutral hydrogen gas in the Lyα clouds is in photoionization equilibrium. In other words, the photoionization (i.e. heating) due to the UV background photons is balanced by cooling processes via thermal bremsstrahlung, Compton cooling and the usual recombination. Photoionization equilibrium allows us to calculate how many photons have to be produced, either from galaxies or quasars, in order to maintain the gas at the observed temperature. Based on this number of photons, we can deduce how many stars were formed at each epoch in the history of the Universe.

We’ve been relying on the assumption that these absorbers are intergalactic gas clouds. How can we be sure they’re not actually just outflows of gas from galaxies? This question has been tackled by studying the two point correlation function of the absorption lines. In a nutshell, two point correlation function measures the clustering of the absorbers. Since galaxies are highly clustered, at least at low redshift, if the low density absorbers are related to the outflow of galaxies, one should expect clustering of these absorption signals. The studies of Lyα forest suggest that the low column density absorbers are not clustered as strongly as galaxies, favoring the interpretation that they are filamentary structures or intergalactic gas [see Rauch (1998) for details].

However this discussion is only valid in the relatively recent Universe. When we look further, and therefore study the earliest times, there are less photons and correspondingly more neutral hydrogen atoms. With more and more absorption, eventually the lines in the forest become so compact and deep that the “forest” becomes a “trough”.

### Gunn-Peterson trough

In their original paper in 1965, Gunn and Peterson estimated how much absorption there should be in a quasar spectrum due to the amount of hydrogen in the Universe. The amount of absorption can be quantified by optical depth which is a measure of transparency. A higher optical depth indicates more absorption due to the neutral hydrogen. Their cosmological model during that time is obsolete now, but their reasoning and derivation are still valid, one can still show that the optical depth



With a modern twist on the cosmological model, the optical depth of neutral hydrogen at each redshift can be estimated to be



where z is the redshift, λ is the transition wavelength, e is the electric charge, me is the electron mass, f is the oscillatory strength of the electronic state transition from N=1 to 2, H(z) is the Hubble constant, h is the dimensionless Hubble’s parameter, Ωb is the baryon density of the Universe, Ωm is the dark matter density, nHI is the number density of neutral hydrogen and nH is the total number density of both the neutral and ionized hydrogen.

Since optical depth is proportional to the density of neutral hydrogen, and the transmitted flux ratio decreases exponentially with respect to the optical depth, even a tiny neutral density fraction, nHI/nH=10-4 (i.e. one neutral hydrogen atom in ten thousands hydrogen atoms) will give rise to complete Lyα absorption. In other words, the transmitted flux should be zero. This complete absorption is known as the Gunn-Peterson trough.

In 1965, Gunn and Peterson performed their calculation assuming that most hydrogen atoms in the Universe are neutral. From their calculation, they expected to see a trough in the quasar (redshift z∼3) data they were studying. This was not the case.

In order to explain the observation, they found that the neutral hydrogen density has to be five orders of magnitude smaller than expected. They came to the striking conclusion that most of the hydrogen in the IGM exists in another form — ionized hydrogen. Lyman series absorptions occur when a ground state electron in a hydrogen atom jumps to a higher state. Ionized hydrogen has no electrons, so it does not absorb photons at Lyman wavelengths.

Figure 3. Animation illustrating Gunn-Peterson trough formation. The continuum blueward of Lyα is drastically absorbed as the neutral hydrogen density increases. From the interactive module.

Gunn and Peterson also showed that since the IGM density is low, radiative ionization (instead of collisional) is the dominant source to maintain such a level of ionization. What an insight!

### The history of reionization

In modern astronomy language, we refer to Gunn and Peterson’s observation by saying the Universe was “reionized”. Combined with our other knowledge of the Universe, we now know the baryonic medium in the Universe went through three phases.

1. In the beginning of the Universe (z>1100), all matter is hot, ionized, optically thick to Thomson scattering and strongly coupled to photons. In other words, the photons are being scattered all the time. The Universe is opaque.
2. As the Universe expands and cools adiabatically, the baryonic medium eventually recombines, in other words, electrons become bound to protons and the Universe becomes neutral. With matter and radiation decoupled, the CMB radiation is free to propagate.
3. At some point, the Universe gets reionized again and most neutral hydrogen becomes ionized. It is believed that the first phase of reionization is confined to the Stromgren spheres of the photoionization sources. Eventually, these spheres grow and overlap, completing the reionization. Sources of photoionizing photons slowly etch away at the remaining dense, neutral, filamentary hydrogen.

When and how reionization happened in detail are still active research areas that we know very little about. Nevertheless, studying the Gunn-Peterson trough can tell us when reionization completed. The idea is simple. As we look at quasars further and further away, there should be a redshift where the quasar sources were emitting before reionization completed. These quasar spectra should show Gunn-Peterson troughs.

### Pinning down the end of reionization

The first confirmed Gunn-Peterson trough was discovered by Becker et al. (2001). As shown in the figure below, the average transmitted flux in some parts of the quasar spectrum is zero.

Figure 4. The first confirmed detection of the Gunn-Peterson trough (top panel, where the continuum goes to zero flux due to intervening absorbers), from a quasar beyond z = 6. The bottom panel shows that the quasar after z = 6 has non-zero flux at all parts. Adapted from Becker et al. (2001).

Figure 5. The close up spectra of the redshift z>6 quasar. In the regions corresponding to the absorbers redshift 5.95 <z< 6.16 (i.e. the middle regions), the transmitted flux is consistent with zero flux. This indicates a strong increment of neutral fraction at z>6. Note that, near the quasar (the upper end), there is transmitted flux. This is mostly owing to the fact that the neutral gas near the quasar is ionized by the quasar emission and therefore the Gunn-Peterson trough is suppressed. This is also known as the proximity effect. As we go to the regions corresponding to lower redshift, there is also apparent transmitted flux, indicating the neutral fraction decreases significantly.

Could this be due to the dense interstellar medium of a galaxy along the line of sight, rather than IGM absorbers? A few lines of evidence refute that idea. First, metal lines typically seen in galaxies (even in the early Universe) are not observed at the same redshift as the Lyman absorptions. Second, astronomers have corroborated this finding using several other quasars. We see these troughs not just a few particular wavelengths in some of the quasar spectra (as you might expect for scattered galaxies), but rather the average optical depth steadily grows, versus a simple extrapolation of how it grows at low redshifts, starting at redshift of z∼6 (see figure below). The rising optical depth is due to the neutral hydrogen fraction rising dramatically beyond z>6.

Figure 6. Quasar optical depth as a function of redshift. The observed optical depth exceeds the predicted trend (solid line), suggesting reionization. Adapted from Becker et al. (2001).

Similar studies using high redshift z>5 quasars confirm a dramatic increases in the IGM neutral fraction from nHI/nH=10-4 at z<6, to nHI/nH>10-3–0.1 at z∼6. At z>6, complete absorption troughs begin to appear.

Figure 7. A more recent compilation of Figure 6 from Fan et al. (2006). Note that the sample variance increases rapidly with redshift. This could be explained by the fact that high redshift star-forming galaxies, which likely provide most of the UV photons for reionization, are highly clustered, therefore reionization is clumpy in nature.

Figure 8. A similar plot, but translating the optical depth into volume-averaged neutral hydrogen fraction of the IGM. The solid points and error bars are measurements based on 19 high redshift quasars. The solid line is inferred from the reionization simulation from Gnedin (2004) which includes the overlapping stage to post-overlapping stage of the reionization spheres.

But it’s important to note that, because even a very small amount of neutral hydrogen in the IGM can produce optical depth τ≫1 and cause the Gunn-Peterson trough, this observational feature saturates very quickly and becomes insensitive to higher densities. That means that the existence of a Gunn-Peterson trough by itself does not prove that the object is observed prior to the start of reionization. Instead, this technique mostly probes the later stages of cosmic reionization. So the question now is: when does reionization get started?

### The start of reionization

Gunn and Peterson told us not only how to probe the end of reionization through the Gunn-Peterson trough, but also how to probe the early stages of reionization using the polarization of the CMB. They discussed the following scenario: if the Universe starts to be reionized at certain redshift, then beyond this, the background distant light sources should experience Thomson scattering.

Thomson scattering is the scattering between free electrons and photons. Recall the fact that due to Thomson scattering by the surrounding nebulae, the starlight undergoes this scattering will be polarized (due to the angle dependence of Thomson scattering). By the same analogy, the CMB, after Thomson scattering off the ionized IGM, is polarized.

As opposed to the case of the Gunn-Peterson troughs, the Thomson scattering results from light interacting with the ionized fraction of hydrogen, not the neutral fraction. Therefore, the large scale polarization of the CMB can be regarded as a complimentary probe of reionization. The CMB polarization detection by Wilkinson Microwave Anisotropy Probe (WMAP) space telescope suggests a significant ionization fraction extending to much earlier in the cosmic history, with z∼11±3 (see figure below). This implies that reionization is a process which starts at some point between z∼8–14 and ends at z∼6–8. This means that cosmic reionization appears to be quite extended in time. This is in contrast to recombination, when electrons became bound to protons, which occurred over a very narrow range in time (z=1089±1).

In summary, a very crude chart on reionization using Gunn-Peterson trough and CMB polarization is shown below.

Figure 9. The volume-averaged neutral fraction of the IGM versus redshift measured using various techniques, including CMB polarization and Gunn-Peterson troughs.

So far, we have left out the discussion of Lyβ, Lyγ absorptions as well as the cosmic Stromgren sphere technique as shown in the figure above. We will now fill in the details before proceeding to the suspects of reionization.

### How about Lyman beta and Lyman gamma

When people talk about the Lyα forest, sometimes it could be an over-simplification. Besides Lyα, other neutral hydrogen transitions such as Lyβ and Lyγ, play a crucial role. Sometimes we can even study the neutral helium Lyα (note that Lyα only designates the transition from electronic state N=1 to 2). We now first look at the advantages and also impediments to include Lyβ and Lyγ absorption lines in the study of the IGM.

1. Recall the formula for the optical depth: τ∝f λ. Due to the decrease in oscillator strength and also the wavelength, for the same neutral fraction density, the Gunn-Peterson optical depths of Lyβ and Lyγ are 6.2 and 17.9 times smaller than Lyα. Having smaller optical depth means that the absorption lines will only saturate at higher neutral fraction/density. Note that once a line is saturated, a large change in column density will only induce a small change in apparent line width. In other words, we are in the flat part of the curve of growth. In this case, the column densities are relatively difficult to measure. Since Lyβ and Lyγ have a smaller optical depths, they could provide more stringent constraints on the IGM ionization state and density when Lyα is saturated. Especially in the case of Gunn-Peterson trough, as we have shown in the figures above, the ability to detect Lyβ and Lyγ is crucial to study reionization.

2. Figure 10. The left-hand panel shows the Voigt profile of Lyα absorption line. The animation shows that when the line is saturated, a large change in column density will only induce a small change in the equivalent width (right-hand panel). This degeneracy makes the measurement of column density very difficult when the absorption line is saturated. From the interactive module.

3. Note that, given the same IGM clouds, the absorption of Lyα will also be accompanied by Lyβ and Lyγ absorptions. Since Lyα has a stronger optical depth, it has broader wings in the Voigt profile and has a higher chance to blend with the other lines. Therefore, simultaneous fits to higher order Lyman lines could help to deblend these lines.
4. As the redshift of quasar increases, there is a substantial overlapping region between Lyβ with Lyα, and so forth. Disentangling which absorption line is due to Lyα and which is due to higher order absorptions is difficult within these regions. Due to this limitation, the study of Lyγ, and higher order lines, is rare.

5. Figure 11. Animation showing the Lyα forest including Lyα and Lyβ absorption lines. One can see that there is a substantial overlapping of Lyα absorption lines onto the the Lyβ region. Therefore, the study of higher order absorptions is very difficult in practice. From the interactive module.

### Probing even higher neutral fraction

In some of the plots above, we have seen that as the neutral fraction increases, Lyα starts to saturate and can only give a moderate lower bound on the optical depth and the neutral fraction. For higher neutral fraction, we have to use Lyβ and Lyγ. It is also clear that even Lyγ will get saturated very quickly. Are there other higher density probes?

Yes, there are other probes! At high neutral fraction, the quasar in this kind of environment resembles a scaled-up version of a O/B type star in a dense neutral hydrogen region (although the latter medium is mostly in molecules, and the former is in atoms). Luminous quasars will produce a highly ionized HII region and create the so called cosmic Stromgren spheres. These cosmic Stromgren spheres could have physical size of about 5 Mpc owing to the extreme luminosity of the quasars.

Measurements of cosmic Stromgren spheres are sensitive to a much larger neutral fraction than the one that Lyγ can probe. However, on the flip side of the coin, they are subjected to many systematics that one could easily imagine, including the three-dimension geometry of the quasar emission, age of the quasar, the clumping factor of the surrounding materials, etc. This explains the large error bar in the reionization chart above.

### How about helium absorption

Although hydrogen is the most abundant element in the Universe, we also know that helium constitutes a large fraction of the Universe. A natural question to ask: how about helium?

As neutral helium has a similar ionization potential and recombination rate compared to neutral hydrogen, it is believed that neutral helium reionization is likely to happen at a similar time together with neutral hydrogen. On the other hand, the ionized HeII has a higher ionization potential, and its reionization is observed to happen at lower redshift, about z∼3. The fact that ionized helium has a higher ionization potential is crucial in disentangling the UV background sources. We will discuss this in detail when we talk about the suspects. In this section, we will discuss some potentials and impediments of the usage of helium Lyα.

1. Helium and hydrogen have different atomic masses. Note that the bulk motion and thermal motion in gas both contribute to the line broadening. However bulk motions are insensitive to the atomic masses whereas the thermal motion velocities go as 1/√μ, where μ is the mean atomic mass. By comparing the broadening of the HeII 304Å and HI 1215Å lines, it is in principle possible to use the difference in atomic masses to measure the contribution from the bulk and thermal motions to the line broadening, separately. Therefore, the theories of the IGM kinematics can be, in principle, be tested. The usefulness of this approach however is remained to be seen.
2. Since ionized helium has a higher ionization potential (Lyman Limit at 228Å) than the neutral hydrogen (Lyman Limit at 912Å), the relative ionization fraction is a sensitive probe of the spectral shape of the UV background. For instance, the contribution from quasars will provide harder spectra than the soft stellar sources and is able to doubly ionize helium.
3. Although helium Lyα forest is awesome, only a small fraction of all quasars are suitable for the search of HeII absorption. The short wavelength of 304Å of Lyman alpha line for ionized helium requires the quasar to be redshifted significantly in order to enter the optical band such that we have enough high resolution spectroscopic instrument to detect it. Furthermore, it is important to note that the helium Lyα line has smaller wavelength than the neutral hydrogen Lyman Limit. As shown in the interactive module, the bound-free absorption beyond the Lyman Limit, as especially at the presence of high column density absorber, will create a large blanketing of lines. This renders the large majority of quasars useless for a HeII search even they might be redshifted enough into the optical band.

### The suspects — who ionized the gas?

We have discussed enough on the technologies to catch the suspects. So we will now go back to the main question on reionization. Who did it? As we have discussed earlier in this post, the IGM is believed to have been reionized early in the Universe’s history via photoionization. Who produced all the high energy, ionizing photons? There are two natural suspects:

1. Star-forming galaxies
2. Quasars

It has been argued that soft sources like star-forming galaxies have to be the dominant sources owing to the fact that there are not many quasars in the early Universe. The Sloan Digital Sky Survey (SDSS) has shown a steep decline in the number of quasars after its peak at z=2.5. If most quasars formed after this epoch, could they play a significant role in the reionization of the Universe, which seems to have completed as early as z∼6.

Taking a closer look at this problem, Faucher-Giguère et al. (2008) made an independent estimate of the photoionization contribution from quasars. In their study, they consider Lyα forests optical depth at redshift 2<z<4.2 from 86 high resolution quasar spectra.

The idea to estimate the fractional contribution from quasars has been described earlier in this post, we will iterate with more details here. Only quasars can produce the very high energy (∼50 eV) photons necessary to doubly ionize helium. Therefore, the luminosity of sources at ∼50 eV can directly tell us the contribution from quasars. By assuming a spectral index for the spectral energy distribution, one can then extrapolate to the lower energies and infer the photoionization rate at the energy range that is relevant to the HI ionization. They show that at most only about 20 per cent of these photons could come from quasars.

Beside showing quasars cannot be the main suspect, Faucher-Giguère et al. also use the Lyα forests to derive the star formation rate indirectly. We have discussed that the absorbers are in photoionzation equilibrium. With this assumption, the authors use the Lyα opacity depth to derive the photoionization rate at different redshifts. With the photoionization rate, they turn that into UV emissivity from the star-forming galaxies and then infer the star formation rate in galaxies at different times in the history of the Universe.

The figure below shows that the derived hydrogen photoionization rate is remarkably flat over the history of the Universe. This suggests that there should be a roughly constant photoionizing source over a large range of cosmic history. This is only possible if the star formation rate in galaxies continues to be high in the early Universe, at high redshift (z∼2.5–4) since the contribution from quasars only begins at this redshift. As the figure shows, such a star formation rate is consistent with the simulation of Hernquist & Springel (2003).

Figure 12. Observations of the hydrogen photoionization rate compared to the best-fitting model from Hernquist & Springel (2003). This suggests that the combination of quasars and stars alone may account for the photoionization. Adapted from Faucher-Giguère et al. (2008).

But there is another way to trace the history of star formation — by directly counting up the number and size of galaxies at different redshifts, using photometric surveys. It turns out that the results from these direct surveys, as performed by Hopkins & Beacom (2006) for example, are in tension with indirect approach. The photometric surveys suggest that the star formation rate decreases after z>2.5 [Prof Hans-Walter Rix once wittily quoted this as the Universe gets tenured at redshift z∼2.5], like in the figure below. The results says that both the major photoionization machines (quasars and stars in galaxies) would be in decline after z>2.5.

Figure 13. Like Figure 12 above, but now using the models of Hopkins & Beacom (2006) where the star formation history is inferred from photometric surveys.

If the photometric surveys are right, Faucher-Giguère et al. strongly suggest that the Universe could not provide enough photons to reionize the intergalactic medium. So, why are there two different observational results?

To sum up this section, by studying Lyα forests, Faucher-Giguère et al. argue that

1. Star-forming galaxies are the major sources of photoionization, although quasars can contribute as much as 20 per cent.
2. The star formation rate has to be higher than the one estimated observationally from the photometric surveys to provide the required photons.

### How the suspect did it

Although the Faucher-Giguère et al. results suggest that star-forming galaxies are the major sources of photoionization, exactly how these suspects actually did it is not clear. The discrepancy between direct and indirect tracers of the star formation history of the Universe might be reconcilable. Observations of star-forming galaxies are plagued with uncertainties, many of which are still areas of active research. Among the major uncertainties are:

1. The star formation rate of galaxies at high redshift. We cannot observe galaxies that are faint and far away. How many photons could they provide? Submillimetre observations in near future could shed more light on this due to their amazing negative K-correction and the ability to observe gas and dust in high redshift galaxies.
2. Starburst galaxies create dust, and dust can obscure observation. Infrared and submillimetre observations are crucial to detect these faint galaxies. It is believed that perhaps most of the ionizing photons come from these unassuming, hard-to-detect galaxies, because they are likely individually small and dim but collectively numerous.
3. IGM clumping factor. In other words, how concentrated is the material surrounding the galaxies? This factor could affect the fraction of ionizing photons escaping into the IGM.

In short, it is possible that the star formation did indeed begin way back in the cosmic history, at z>2.5, but we underestimated it in the photometric surveys due to these uncertainties. In this case, all the results from simulations, Lyα forests, and photometric surveys could fit together.

### Any other suspects on loose

We have now two suspects in custody. The remaining question: are there any suspects on loose? To answer this question, some other interesting proposals were raised, among them:

1. High-energy X-ray photons from supernova remnants or microquasars at high redshift. This is possible, but its contribution might overestimate the soft X-ray background that we observe today.
2. Neutrino is believed to gain their mass through the seesaw mechanism. One of the potential ingredient of this mechanism is the sterile neutrino. Reionization by decaying sterile neutrinos cannot be completely ruled out yet. That said, since there is no confirmed detection of sterile neutrino in particle physics experiments, one should not take this possibility too seriously and blow one own mind.

In summary, we cannot be 100 per cent sure, but other suspect other that quasars and stellar sources is unlikely.

### So what have we learned?

To summarize what we have discussed in this post:

1. Non-detection can be a good thing. The non-detection of the Gunn-Peterson trough in the early days demonstrated that the IGM is mostly in an ionized state.
2. The Lyα forest gives us a 3D picture of the IGM.
3. The Gunn-Peterson trough provides a direct probe on the end stage of the reionization.
4. CMB radiation polarization suggests when the process of reionization began, and shows that reionization is a process extended in time.
5. Star-forming galaxies are likely to be responsible for the reionization, but how they exactly did it is a question still under investigation.

## ARTICLE: Interpreting Spectral Energy Distributions from Young Stellar Objects

In Journal Club, Journal Club 2013 on April 15, 2013 at 8:02 pm

Posted by: Meredith MacGregor

#### 1. Introduction

When discussing young stellar objects and protoplanetary disks (as well as many other topics in astronomy), astronomers continually throw out the term ‘SED.  In early April, I attended a conference titled ‘Transformational Science with ALMA: From Dust to Rocks to Planets– Formation and Evolution of Planetary Systems.’  And, I can attest to the fact that the term ‘SED’ has come up in a very significant fraction of the contributed talks.  For those not intimately familiar with this sub-field, this rampant use of abbreviations can be confusing, making it difficult to glean any useful take-aways from a talk.  So, in addition to summarizing the article by Robitaille et al. (2007), the goal of this post is to give a bit of an introduction to the terminology and motivation for the growing field of star and planetary system formation.

SEDs: What They Are and Why We Care So Much

The abbreviation ‘SED’ stands for ‘Spectral Energy Distribution.’  If you want to sound like an expert, this should be pronounced exactly as it appears (analogous to the non-acronym ‘said’).  A SED is essentially a graph of flux versus wavelength.  In the context of the Robitaille et al. article, we are most interested in the SEDs for young stars and the envelopes and disks surrounding them.  So, why exactly does the flux from a young stellar object (YSO) vary with wavelength?  As it turns out, different regions of the YSO emit at different wavelengths.  This means that when we observe at different wavelengths, we are actually probing distinct regions of the star and its surrounding disk and envelope. By tracing out the entire SED for a YSO, we can determine what the geometry, structure, and constituents are for that object.

Before we go too far, it is worth taking a moment to clarify the terminology used to classify young stars.  A young stellar object or YSO is typically defined as any star in the earliest stages of development.  YSOs are almost always found in or near clouds of interstellar dust and gas.  The broad YSO class of objects is then divided into two sub-classes: protostars and pre-main sequence stars.  Protostars are heavily embedded in dust and gas and are thus invisible at optical wavelengths.  Astronomers typically use infrared, sub-millimeter, and millimeter telescopes to explore this stage of stellar evolution, where a star acquires the bulk of its material via infall and accretion of surrounding material.  Pre-main sequence (PMS) stars are defined to be low mass stars in the post-protostellar phase of evolution.  These stars have yet to enter the hydrogen-burning phase of their life and are sill surrounded by remnant accretion disks called ‘protoplanetary’ disks.  A more detailed a summary of this classification scheme can be found here.  However, because the evolution of young stars is far from well-understand, astronomers often typically these words interchangeably.

SEDs for pre-main sequence stars are often seen to have a bump or excess in the infrared. This bump is generally interpreted as being due to thermal emission from warm dust and gas surrounding the central star. As an illustrative example, let’s consider a protoplanetary disk around a young star. The image below is taken from a review article by Dullemond & Monnier (2010) and shows a graphical representation of the emission spectrum from such a disk.

A graphical representation of the emission spectrum from a protoplanetary disk and the telescopes that can be used to probe different regions of the disk. Taken from Dullemond & Monnier (2010).

In this case, emission in the near-infrared traces the warmer inner regions of the disk.  As you move into submillimeter wavelengths, you start to probe the outer, cooler regions of the disk.  By modeling the observed SED of a pre-main sequence star, you can determine what components of the source are contributing to the observed flux.  The second figure below is taken from a paper by Guarcello et al. (2010).  The left panel of this figure shows the observed SED (solid line) of a BWE star (in line with other ridiculous astronomy acronyms, this stands from ‘blue stars with excesses’) and the unreddened photospheric flux we would expect to see if the star did not have a disk (dashed line).  The right panel shows a model fit to this data.  The authors describe the observed SED using a model with four distinct components, each represented by a colored line in the figure: (1) emission from the reddened photosphere of the central star, (2) radiation scattered into the line of sight from dust grains in the disk, (3) emission from a collapsing outer envelope, and (4) thermal emission from a circumstellar disk.  The summation of these four component makes up the complete SED for this BWE star.

Left panel: The observed SED (solid line) of a BWE star and the unreddened photospheric flux we would expect to see if the star did not have a disk (dashed line). Right panel: A four component model fit to this data. Taken from Guarcello et al. (2010).

Describing and Classifying the Evolution of a Protostar

As a protostar evolves towards the Zero Age Main Sequence (ZAMS), the system geometry (and thus the SED) will evolve as well. Therefore, the stage of evolution of a protostar is often classified according to both the general shape and the features of the SED.  A graphical overview of the four stages of protostellar evolution are shown below (Andrea Isella’s thesis, 2006).  Class 0 objects are characterized by a very embedded central core in a much larger accreting envelope.  The mass of the central core grows in Class I objects and a flattened circumstellar accretion disk develops.  For Class II objects, the majority of circumstellar material is now found in a disk of gas and dust.  Finally, for Class III objects, the emission from the disk becomes negligible and the SED resembles a pure stellar photosphere.  The distinction between these different classes was initially defined by the slope of the SED (termed the ‘spectral index’) at infrared wavelengths.  Class I sources typically have SEDs that rise in the far- and mid-infrared, while Class II sources have flat or falling SEDs in the mid-infrared.  However, this quantitative distinction is not always clear and is not an effective way to unambiguously distinguish between the different object classes (Protostars and Planets V, page 127-128).

A graphical overview of the four stages of protostar evolution taken from Andrea Isella’s thesis (2006). A typical SED of each class is shown in the left column and a cartoon of the corresponding geometry is shown in the right column.

#### 2. The Article Itself

One common method of fitting SEDs is to assume a given gas and circumstellar dust geometry and set of dust properties (grain size distribution, composition, and opacity), and then use radiative transfer models to predict the resulting SED and find a set of parameters that best reproduce the observations.  However, fitting SEDs by trial and error is a time consuming way to explore a large parameter space.  The problem is even worse if you want to consider thousands of sources.  So, what’s to be done?  Enter Robitaille et al.  In order to attempt to make SED fitting more efficient, they have pre-calculated a large number of radiative transfer models that cover a reasonable amount of parameter space.  Then, for any given source, one can compare the observed SED to this set of models to quickly find the set of parameters that best explains the observations.

Let’s Get Technical

The online version of this fitting tool draws from 20,000 combinations of physical parameters and 10 viewing angles (if you are particularly curious, the online tool is available here).  A brief overview of the parameter space covered is as follows:

• Stellar mass between 0.1 and 50 solar masses
• Stellar ages between $10^3$ and $10^7$ years
• Stellar radii and temperatures (derived directly from stellar mass using evolutionary tracks)
• Disk parameters (disk mass, accretion rate, outer radius, inner radius, flaring power, and scale height) and envelope parameters (the envelope accretion rate, outer radius, inner radius, cavity opening angle, and cavity density) sampled randomly within ranges dictated by the age of the source

However, there are truly a vast number of parameters that could be varied in models of YSOs.  Thus, for simplicity, the authors are forced to make a number of assumptions.  Here are some of the biggest assumptions involved:

1. All stars form via accretion through a disk and an envelope.
2. The gas-to-dust ratio in the disk is 100.
3. The apparent size of the source is not larger than the given aperture.

The last constraint is not required, but it allows a number of model SEDs to be cut out and thus speeds up the process.  Furthermore, the authors make a point of saying that the results can always be scaled to account for varying gas-to-dust ratios, since only the dust is taken into account in the actual radiative transfer calculations.

Does This Method Really Work?

If this tool works well, it should be able to correctly reproduce previous results.  In order to test this out, the authors turn to the Taurus-Auriga star forming region.  They select a sample of 30 sources from Kenyon & Hartman (1995) that are spatially resolved, meaning that there is prior knowledge of their evolutionary stage from direct observations (i.e. there is a known result that astronomers are fairly certain of to compare the model fits against).  When fitting their model SEDs with the observed SEDs for this particular star forming region, the authors throw in a few additional assumptions:

1. All sources are within a distance range of 120 – 160 AU (helps to rule out models that are too faint or too luminous).
2. The foreground interstellar extinction is no more than $A_v$ = 20.
3. None of the sources appeared larger than the apertures used to measure fluxes.

The authors then assign an arbitrary cut-off in chi-squared for acceptable model fits: $\chi^2 - \chi^2_\text{best} < 3$.  Here, $\chi^2_\text{best}$ is the $\chi^2$ of the best-fit model for each source.  Robitaille et al. acknowledge that this cut-off has no statistical justification: ‘Athough this cut-off is arbitrary, it provides a range of acceptable fits to the eye.’  After taking Jim Moran’s Noise and Data Analysis class, I for one would like to see the authors try a Monte Carlo Markov Chain (MCMC) analysis of their 14-dimensional space (for more detail on MCMC methods see this review by Persi Diaconis).  That might make the analysis a bit less ‘by eye’ and more ‘by statistics.’

The upshot of this study is that for the vast majority of the sources considered, the best-fit values obtained by this new SED fitting tool are close to the previously known values. Check.

It is also worth mentioning here, that there are many other sets of SED modeling codes.  One set of codes of particular note are those written by Paola D’Alessio (D’Alessio et al., 1998; D’Alessio et al., 1999; D’Alessio et al, 2001).  These codes were the most frequently used in the results presented at the ALMA conference I attended.  The distinct change in the D’Alessio models is that they solve for the detailed hydrostatic vertical disk structure in order to account for observations of ‘flared’ disks around T Tauri stars (flaring refers to an increase in disk thickness at larger radii).

But, Wait! There are Caveats!

Although the overall conclusion is that this method fits SEDs with reasonable accuracy, there are a number of caveats that are raised.  First of all, the models tend to overestimate the mid-IR fluxes for DM Tau and GM Aur (two sources known to have inner regions cleared of dust).  The authors explain that this is most likely due to the fact that their models for central holes assume that there is no dust remaining in the hole.  In reality, there is most likely a small amount of dust that remains.  Second, the models do not currently account for the possibility of young binary systems and circumbinary disks (relevant for CoKu Tau 1).

The paper also addresses estimating parameters such as stellar temperature, disk mass, and accretion rate from SED fits.  And, yes, you guessed it, these calculations raise several more issues.  For very young objects, it is difficult to disentangle the envelope and disk, making it very challenging to estimate a total disk mass.  To make these complications clearer, the set of two plots below from the paper show calculated values for the disk mass plotted against the accepted values from the literature.  It is easily seen that the disk masses for the embedded sources are the most dissimilar from the literature values.

Two plots from Robitaille et al. (2007) that show calculated values for the disk mass plotted against the accepted values from the literature. It is easily seen that the disk masses for the embedded sources (right) are the most dissimilar from the literature values.

Furthermore, even if the disk can be isolated, the dust mass in the disk is affected by the choice of dust opacity.  That’s a pretty big caveat!  A whole debate was started at the ALMA conference over exactly this issue and the authors have simply stated the problem and swept it under the rug in just one sentence.  In 2009, David Hogg conducted a survey of the rho Ophiucus region and used the models of Robitaille et al. (2007) to determine the best-fit dust opacity index, $\beta$ for this group of sources.   Hogg found that $\beta$ actually decreases for Class II protostars, a possible indication of the presence of larger grains in the disk.  Robitaille et al. also mention that the calculated accretion rates from SED fitting are systematically larger than what is presented in the literature.  The authors conclude that future models should include disk emission inside the dust destruction radius, the radius inside which it is too hot for dust to survive.  A great example of the complications that arise from a disk with a central hole can be seen in LkCa 15 (Espaillat et al., 2010Andrews et al., 2011).   The figure below shows the observed and simulated SEDs for the source (left) as well as the millimeter image (right).  The double ansae (bright peaks or ‘handles‘ apparent on either side of the disk) seen in the millimeter contours are indicative of a disk with a central cavity.

Left: The observed and simulated SEDs for Lk Ca 15. The sharp peak seen at 10 microns is due to silicate grains within the inner hole of the disk. Right: The millimeter image of the disk. (Espaillat et al., 2010; Andrews et al., 2011)

In this case, a population of sub-micron sized dust within the hole is needed in order to produce the observed silicate feature at 10 microns.  Furthermore, an inner ring is required to produce the strong near-IR excess shortward of 10 microns.  A cartoon image of the predicted disk geometry is shown below.  To make things even more complicated, the flux at shorter wavelengths appears to vary inversely with the flux at longer wavelengths over time (Espaillat et al., 2010).  This phenomenon is explained by changing the height of the inner disk wall over time.

A cartoon image of the predicted disk geometry for Lk Ca 15 showing the outer ring, silicate grains within the hole, and the inner ring. (Espaillat et al., 2010)

Finally, Robitaille et al. discuss how well parameters are constrained given different combinations of data points for two example sources: AA Tau and IRAS 04361+2547.  In both sources, if only IRAC (Infrared Array Camera on the Spitzer Space Telescope) fluxes obtained between 3.6 and 8 microns are used, the stellar mass, stellar temperature, disk mass, disk accretion rate, and envelope accretion rate are all poorly constrained.  Things are particularly bad for AA Tau in this scenario, where only using IRAC data results in ~5% of all SED models meeting the imposed goodness of fit criterion (yikes!).  Adding in optical data to the mix helps to rule out models that have low central source temperatures and disk accretion rates.  Adding data at wavelengths longer than ~ 20 microns helps to constrain the evolutionary stage of the YSO, because that is where any infrared excess is most apparent.  And, adding submillimeter data helps to pin down the disk mass, since the emission at these longer wavelengths is dominated by the dust.  This just goes to show how necessary it is to obtain multi-wavelength data if we really want to understand YSOs, disks, and the like.

#### 3. Sources:

Where to look if you want to read more about anything mentioned here…

## A model of the emergence of coherence

In Uncategorized on April 14, 2013 at 4:50 pm

## Ostriker’s 1964 Isothermal Cylinder model

In Uncategorized on April 14, 2013 at 4:43 pm

Particularly impressive is that, if you read the fine print, Jerry Ostriker was an NSF graduate fellow (meaning he was a grad student!) when he wrote this paper, which solves differential equations analytically at a level of technical virtuosity undoubtedly beyond anyone reared in the age of Mathematica.  Ostriker obtains solutions for polytropic cylinders, of which an isothermal cylinder is the case where $n=\infty$.  One has

$P=K_n \rho ^{1+1/n}$,

which leads to the fundamental equation

$K_n (n+1) \nabla ^2 \rho ^{1/n} = -4\pi G\rho$.

Using the transformation

$r \equiv \bigg[ \frac{(n+1)K_n}{4\pi G \lambda ^{1-1/n} } \bigg]^{1/2} \xi$,

$\rho \equiv \lambda \theta ^n$,

one has a version of the Lane-Emden equation

$\frac {1}{\xi}\frac{d}{d\xi} \big(\xi \frac{d\theta}{d\xi} \big) =\theta''+\frac{1}{\xi}\theta'=-\theta ^n.$

For n=0, 1, and infinity there is a closed form solution, for other n a power series solution.  In the particular case of an isothermal cylinder, the EOS is that for an ideal gas, and one has

$r\equiv \bigg[ \frac{K_I}{4\pi G \lambda} \big]^{1/2}\xi$

and

$\rho=\lambda e^{-\psi}$.

Using the fundamental equation noted above and manipulating yields

$\psi''+\frac{1}{\xi}\psi'=e^{-\psi}$.

Letting $z=-\psi +2\ln \xi$ and $t=\sqrt{2}\ln\xi$, we find

$2\frac {d^2z}{dt^2}+e^z=0$.

Letting $z=\ln y$ the equation may be integrated to give

$\psi(\xi)=2\ln\big(1+\frac{1}{8} \xi^2\big)$.

Using these results in expressions Ostriker provides in the paper that give general forms for density and mass, one has

$\rho= \frac{\rho_0}{\big(1+\xi^2/8 \big)^2}$

and $M(\xi)=\frac{2 k_B T}{\mu m_0 G} \frac{1}{1+8/\xi^2}$.

The expression for density should look familiar from the Pineda paper: it gets integrated to yield his eqn. (1) for the surface density of the cylinder.

## 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/

Summary By Pierre Christian

## Disclaimer

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

Confining
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, 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.

## Conclusion

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.

## References

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

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.

## Module Prototype: Director’s Cut of a WorldWide Telescope Tour

In Special Topics Modules on April 2, 2013 at 4:46 pm

#### W5: Multigenerational Star Formation

This post, prepared by Alyssa Goodman and used in class on 4/2/13, is intended to give AY201b students an idea of  how ready their interactive module  should be when it is presented, and the level of detail to offer in a presentation.
Click here to view original WWT Tour, or here for Tour description.  Click here to download “Director’s Cut.” In the (Windows or) web version of WWT, go to “Explore, Open…, Tour” and select the file that you’ve downloaded in order to view the Tour. Or (warning, beta!) watch the Director’s Cut tour in WWT/HTML5.

Background: Who are the original Tour’s authors?

• Xavier Koenig: A finishing graduate student at the Harvard-Smithsonian Center for Astrophysics when this Tour was created.  His PhD thesis concerned analysis of Spitzer Space Telescope observations of the star-forming region W5.  As of 2013, Dr. Koenig is a postdoctoral fellow at Yale University.

• Lori Allen: Thesis advisor at the Center for Astrophysics to Xavier Koenig when this Tour was made.  Today, Dr. Allen, is Deputy Director of the Kitt Peak National Observatory.

• Sanjana Sharma: High-School student at the Winsor School in Boston when this Tour was made.  Sharma was an intern with the WorldWide Telescope Ambassadors group at the Center for Astrophysics before moving to New Haven, where she is in the class of 2014 at Yale.

What points are raised in the WWT Tour narration that could be explored more deeply by an interested viewer? (text in purple, concerning how we know stars ages, are now clickable links within the Tour, as a prototype)

• small and faint (how small (angle), how faint, #’s)

• faint diffuse glow, hot gas  (how hot, what does this region look like at other wavelengths)

• red glow(=warm dust, how warm, and how do we know?)

• one burst of star formation can cause another (triggered star formation)

• may have been 3 successive generations of star formation (how do we know which generation is which?)

• stars…dispersed over time (how fast, how much time? how do we know?)

• large clusters (what’s the definition of a “cluster” and is it different for young stars?)

• changing light they emit to infer the presence of multiple objects (how?)

• disk…that could maybe form planets (discuss how & whether this “hostile” environment matters)

• comet-shaped tail that glows in the infrared (how does that happen?)

• pillars compressed from outside…squeezed on inside by internal gravity (how, which forces do what on what time scales?)

• brand-new stars are emerging (how do we know they are new?)

• comparison of pillars/mountains W5/Eagle nebula same scale (angular/linear?…turns out both, as these sources are coincidentally at similar distances from us!)

• Video: Spitzer “Hidden Universe” interviews with Allen & Koenig about W5
• PhD Thesis: Xavier Koenig’s thesis (PDF)

• Journal Article: Koenig et al. 2008, Clustered and Triggered Star formation in W5: Observations with Spitzer (ADS link)  [Abstract: We present images and initial results from our extensive Spitzer Space Telescope imaging survey of the W5 H II region with the Infrared Array Camera (IRAC) and Multiband Imaging Photometer for Spitzer (MIPS). We detect dense clusters of stars, centered on the O stars HD 18326, BD +60 586, HD 17505, and HD 17520. At 24 μm, substantial extended emission is visible, presumably from heated dust grains that survive in the strongly ionizing environment of the H II region. With photometry of more than 18,000 point sources, we analyze the clustering properties of objects classified as young stars by their IR spectral energy distributions (a total of 2064 sources) across the region using a minimal-spanning-tree algorithm. We find ~40%-70% of infrared excess sources belong to clusters with >=10 members. We find that within the evacuated cavities of the H II regions that make up W5, the ratio of Class II to Class I sources is ~7 times higher than for objects coincident with molecular gas as traced by 12CO emission and near-IR extinction maps. We attribute this contrast to an age difference between the two locations and postulate that at least two distinct generations of star formation are visible across W5. Our preliminary analysis shows that triggering is a plausible mechanism to explain the multiple generations of star formation in W5 and merits further investigation.]

## 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

## Introduction

### 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^o$$60^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.

### Observations

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.

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.

### Interpretations

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.

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.

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.

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.

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

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.

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.

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.

### Discussion

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.

## Conclusions

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.

## References

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 Filament Runs Through It: EVLA observations of B5

In Journal Club, Journal Club 2013 on April 1, 2013 at 4:37 am

“EVLA OBSERVATIONS OF THE BARNARD 5 STAR-FORMING CORE: EMBEDDED FILAMENTS REVEALED”

by JAIME E. PINEDA, ALYSSA A. GOODMAN, HÉCTOR G. ARCE, PAOLA CASELLI, STEVEN LONGMORE, and STUARTT CORDER

Post by Zachary Slepian

0. Main points

Like a good poem, this paper is short but thought-provoking. For a cosmologist, perhaps the most interesting take-aways relate to star formation. First, the paper helps with a long-standing problem (trace it back to Larson’s Laws!) in how stars form—lots of molecular clouds have velocity dispersions too high to allow star formation (or to be explained by thermal motions alone).  But Pineda et al. find that in at least one dense core (within a molecular cloud), velocity dispersions are small enough to make it likely that a Jeans-like instability could form stars.  As evidence for this, the authors find that the Jeans length is on the order of the separation between a young stellar object in the core and a starless condensation to its north, the latter perhaps being a precursor for another star.  The fact that these two objects are separated by a Jeans length suggests they both stem from Jeans-like collapse and fragmentation.

The paper’s second major conclusion is that there is a 5000 AU long filament at the center of the core.  This is interesting on its own, because the filament turns out to be well fit by Ostriker’s 1964 model of an isothermal cylinder in hydrostatic equilibrium, in contrast to a different filament observed in 2011 in another star-forming region, which is not fit by this model and is also much smaller than Pineda et al.’s.   More significantly, though, the filament Pineda et al. find is an additional piece of evidence for the efficacy of a Jeans-like instability, as it could easily have resulted from Jeans collapse.

Ultimately, the low velocity dispersion is the common factor between these two pieces of evidence for Jeans collapse: it is only in the absence of turbulent support that Jeans collapse occurs in cores, and the low velocity dispersion implies this absence.

In what follows, I 1) offer some questions, with linked answers on separate pages, to help you think through the paper’s text and figures, and 2) take-aways summarizing the main point of each figure. Details on how the observations were done, which aren’t essential to understanding the physics,  are here, while I go through Ostriker’s 1964 isothermal filament model here.  I close with some recommended further reading and a linked summary of Alyssa and collaborators’ model of how coherence emerges in cores.

1. Introduction

As noted earlier, molecular clouds previously have been found to have “supersonic” velocity dispersions—what does this mean? (1)

These motions must be dissipated to allow collapse and star formation.  Why?  After all, such motions will increase the Jeans length.   But can’t we just treat them as producing an extra effective pressure, and presume that above the Jeans length these regions still collapse? (2)

The technique used to study these dense cores is $NH_3$ mapping, specifically the (1,1) transition.  Where is another place you can find $NH_3$? Why use $NH_3$ rather than our more popular friend $CO$?  What the heck is the (1,1) transition? (3)

The authors note that previous work with the Green Bank Telescope (GBT) has already mapped the ammonia in B5—why do it again? Hint: they used another array, the Expanded Very Large Array (EVLA), which has better resolution.  Why would they want better resolution if they are interested in studying velocity dispersions? (4)

2. Results

When in a rush, it is apparently standard astronomy practice to scan an article’s figures, introduction, and conclusion—after all, astronomers, unlike particle physicists, are visual people.  So I’ll focus my discussion on guiding you through the figures.

Figure 1

Figure 1 from Pineda et al. 2011.

. . . explains why this paper was done even though observations had already been made of B5 with the GBT.  Compare left to right panel: what is the difference? (5)

The figures say they are integrated intensity, but the right-hand panel is in “mJy beam inverse km/s”.  Help!  That does not look like units of intensity.  Try resolving this yourself before clicking the answer! (6)

• Take-aways: The regions with greater integrated intensity are those with higher density, so this figure is showing us where the structures are likely to form or have formed (e.g. the filament), since the densest regions will collapse first ($\tau_J\propto 1/\sqrt{\rho}$).  The right panel has zoomed in on the regions with subsonic velocity dispersions, which are bounded in the left panel by the orange contour.

Figure 2

Figure 2 from Pineda et al. 2011.

What is the most important part of this figure? I’ll give you a hint: it is black and white, and not so exciting looking! (7)

What are the two panels actually showing?  This is not so obvious! (8)

• Take-aways: the left panel shows that the filament is on order the Jeans length, and also that the dense starless condensation and young stellar object (YSO) are separated by on order a Jeans length.  Hence, Jeans collapse is the likely culprit.  The right panel shows the velocity dispersion: regions with darker blue have smaller $\sigma_v$ and hence more coherent velocities; they become less coherent near the YSO, possibly due to feedback.

Figure 3

Figure 3 from Pineda et al. 2011.

Is Figure 3 redundant? (9)

Why do the authors do two different histograms in Figure 3?  Why divide into bits of the region near the YSO and not near it? Why might the red histogram (bits of the region near the YSO) have a higher centroid value of $\sigma_v$ and width than that for bits of the region not near the YSO? Extra credit: can you use Figure 3 to estimate how many other stars we should eventually expect to form near the YSO already observed (assume there are no stellar outflows or winds)? (10)

And why do they use the criterion of two beam widths as the cut between objects close to the YSO and far from the YSO? (11)

Finally, in the caption for Figure 3, they note $\mu=2.33$?  Why?  There’s a very simple answer! (12)

• Take-aways: This Figure is a different way of seeing the same information as the right panel of Figure 2.  It shows that the velocity dispersion is by and large subsonic, emphasizing that the velocity is fairly coherent, especially in regions away from the YSO.  The red histogram in the Figure emphasizes that near the YSO the velocity dispersion is higher, though still susbsonic, as noted earlier likely due to feedback, e.g. radiation from the YSO or interaction between an outflow or stellar wind and the dense surrounding gas.

Figure 4

Figure 4 from Pineda et al. 2011.

Is perhaps the most difficult figure at first glance.  Radius of zero on the horizontal axis is the center of the filament, + and – move right and left in the yellow box around the filament in Figure 1. What is the key point of panel a? (13)

Panel b is perhaps more interesting, because it shows the filament is isothermal.  The model with p=4 is a better fit to the filament, which is Ostriker’s predicted value for p in eqn. (1) if the filament is isothermal.  What is the role of the blue curve, the beam response—why should we care about that? (14)

• Take-aways: Imagine the filament as a cylinder. The top panel shows that as one goes out radial to concentric shells of the cylinder, the velocity dispersions remain sub-sonic  and roughly constant until one is well away from the filament, which is isothermal and resolved.

Pineda et al. note that their filament contrasts with the ones detected by Herschel in other star-forming regions (Arzoumanian et al. 2011), which are not isothermal, and are fit with p=2 instead (see below).

Showing how Herschel filament is better fit with p=2, non-isothermal profile. From Arzoumanian et al. 2011 Fig. 4.

3. Handouts from JC discussion and further reading

Alyssa and collaborators proposed a model of why coherence emerges in “COHERENCE IN DENSE CORES. II. THE TRANSITION TO COHERENCE”, which I go through here.

How does a turbulently-supported core turn into, for instance, a filament? Stella Offner knows!  See “Observing Turbulent Fragmentation in Simulations: Predictions for CARMA and ALMA” by Stella S.R. Offner, John Capodilupo, Scott Schnee, and Alyssa A. Goodman.

For the discovery of the sharp transition to a dense, coherent (velocity) core, alluded to in the paper, see “Direct observation of a sharp transition to coherence in Dense Cores” by Jaime E. Pineda, Alyssa A. Goodman, Héctor G. Arce, Paola Caselli, Jonathan B. Foster, Philip C. Myers, and Erik W. Rosolowsky.

For extremely recent (March 2013) discussion of core formation and filaments in another region, see “Cores, filaments, and bundles: hierarchical core formation in the L1495/B213 Taurus region” by A. Hacar, M. Tafalla, J. Kauffmann, and A. Kovacs.

For the Herschel core observations alluded to in the paper, see “Characterizing interstellar filaments with Herschel in IC 5146″ by D. Arzoumanian et al., 2011.

For context to  “COHERENCE IN DENSE CORES. II. THE TRANSITION TO COHERENCE”,  you may want also to see the companion paper.

“COHERENT DENSE CORES. I. NH3 OBSERVATIONS” by J.A. Barranco and A.A. Goodman.

## Answer to question (14)

In Uncategorized on March 28, 2013 at 7:19 am

The blue curve, beam response, shows that the beam has the resolution to resolve the two different models.

## Answer to question (13)

In Uncategorized on March 28, 2013 at 7:18 am

That even near the filament the velocity is pretty much sub-sonic—the upper dotted line is the sound speed. This further supports the claim that Jeans collapse may be the mechanism behind the filament, although there may have been post-processing or feedback after the filament formed, so we can take this with a grain of salt.

## Answer to question (12)

In Uncategorized on March 28, 2013 at 7:17 am

Molecular cloud—particles are mostly molecular not atomic hydrogen, meaning $\mu\sim 2$.  The .33 accounts for bigger things, such as Helium.  $\mu=2.33$ corresponds to assuming a typical, Milky Way elemental abundance.

## Answer to question (11)

In Uncategorized on March 28, 2013 at 7:16 am

At first glance the cut seems arbitrary. If the thought is that the YSO is heating a region around it, leading to bigger velocity dispersions, the size of that region should be dictated by something like the physics of a Stromgren sphere, and have nothing to do with the beam size.  A similar comment applies to the case where the region is heated by a stellar wind or outflow.

However, if the beam has width 6”, and we assume that represents the full-width at half maximum, then 1 sigma is 3”, and 12” represents 2-sigma on either side of the Gaussian beam profile.  So, while technically structure >6” would be resolved, anything within 12” of the beam will be smoothed by a Gaussian profile.  Hence there is potential for contamination of the black histogram by the YSO if one is closer to it than 12”.

## Answer to question (10)

In Uncategorized on March 28, 2013 at 7:14 am

Radiation from the YSO, interaction between outflow/stellar wind and gas.

## Answer to question (9)

In Uncategorized on March 28, 2013 at 7:13 am

To some extent, yes—it shows the velocity dispersions are mostly subsonic, which is also evident from Figure 2 b if you calculate $c_{s,ave}=.2\;\rm{km/s}$ and see that most of the region has lower values than this. But Figure 3 shows you just how much of the region doesn’t (i.e. is supersonic), and that most of the bits that are supersonic are near the young stellar object (YSO).

## Answer to question (8)

In Uncategorized on March 28, 2013 at 7:12 am

The left-hand panel shows the centroid velocity: in each little beam-size cell, one has a Gaussian-esque line profile, with some centroid whose velocity can be calculated using the line’s red or blue shift.  The right-hand panel shows, in each little beam size cell, what the width of this profile is.  Interestingly, one can compare the sigma of the centroid velocities in the boxed region of the left panel occupied by the filament with the sigma in each beam-sized cell and use this ratio to test the geometry of the structure in the region.  This is essentially because of Larson’s Laws: Philip Mocz’s post on Larson’s laws explains how this ratio depends on the geometry, and calculates it for several idealized cases.  The most relevant one for us is the long sheet, which in 2-d projection would look similar to a filament.  The ratio predicted for such a geometry is 2.67.  Estimating by eye the needed quantities from Figure 2 (try this yourself), we find $\sigma_v/\Delta \sigma_v\sim 3$—offering somewhat independent confirmation that we have a filament.  Note this calculation will be somewhat sensitive to how you choose the region over which to calculate the ratio—but we already know what we are looking for (the filament), so can choose the region accordingly. However, this is why I qualified this with “somewhat” above.

## Answer to question (7)

In Uncategorized on March 28, 2013 at 7:11 am

The Jeans length drawn on, between the YSO and the starless condensation—it is one of the key pieces of evidence that Jeans collapse may be happening.

## Answer to question (6)

In Uncategorized on March 28, 2013 at 7:10 am

mJy are a unit of specific intensity, or flux at a particular frequency.  Integrated intensity is therefore mJy times Hz (or 1/s).  One can convert Hz to km/s by determining what velocity is implied by a given frequency shift $\Delta\nu/\nu$, and that is what has been done here, with $\nu \approx 23.$ GHz.

## Answer to question (5)

In Uncategorized on March 28, 2013 at 7:07 am

Resolution!  And thus probing details of velocity dispersion. The GBT image does not really resolve the filament, and only marginally shows the starless condensation—which is separated from the young stellar object by a Jeans length, a key piece of evidence for the efficiency of Jeans collapse in the region.

## Answer to question (4)

In Uncategorized on March 28, 2013 at 7:05 am

Because the GBT couldn’t resolve spatial variations in the velocity dispersion or column density.  Resolving the former offers more information on the  problem of supersonic vs. subsonic dispersions, already discussed as key for star formation, and column density spatial variations probe whether there are structures along the line of sight—as indeed the authors find (the filament!)

## Answer to question (3)

In Uncategorized on March 28, 2013 at 7:04 am

In your bathroom, if you clean it—$NH_3$ is actually just ammonia!  Ammonia is common in regions near the galactic center (Kaifu et al. 1975). And CO  becomes optically thick before ammonia—and optically thick is optically useless when you want to find densities!  Ammonia allows study of denser regions—exactly what is needed to probe star formation.

The (1,1) transition is actually quite exotic: no mere rotational line here.  Ammonia is shaped like a triangular pyramid, with the three H’s at the base and the N on top.  The N can quantum mechanically tunnel through the potential barrier of the base, inverting the pyramid.  So the (1,1) transition is also known as an “inversion” transition.

You might wonder why there is a potential barrier at all—after all, each H atom is neutral, and so is the N on top.  But, if you were an electron on one of the H’s, were would you want to be? Certainly far from the other 2 Hs’ electrons!  Hence each H’s electron will spend most of its time outside the base, meaning the triangle formed by the H’s will be slightly positively charged on the inside. Similarly, the electron on the N will want to be as far from the other three electrons as it can be, so it will hover above the N, meaning the bit of the N facing the pyramid’s triangular base will be slightly positively charged.  Ergo, potential barrier.

Making simple assumptions, we can estimate the energy of the inversion.  Assuming the distance to the base’s center for each H’s electron is (a+l), a the Bohr radius and l the ammonia bond length, 1 angstrom, and that the N’s electrons are (a+l) above this center, we calculate the potential where the 7 protons in N are.  Converting to energy and thence frequency, we obtain $\nu\sim 1 \;\rm{GHz}$. Incidentally, to go further one could use the WKB approximation to estimate the tunneling probability.  Given a minimum flux per beam width to which the telescope is sensitive (14 mJy here is the noise), one could even then place a lower bound on the column density observable with this transition by a particular instrument, and assuming isotropy, one could get a density from the column density.

## Answer to question (2)

In Uncategorized on March 28, 2013 at 7:02 am

“Inevitably when large enough” is the key.  One can indeed treat the supersonic motions as contributing an additional pressure, and raising the Jeans length—the problem is, they raise it above the typical size of a dense core!  Hence for dense cores to collapse, the turbulence must be dissipated so that the Jeans length goes down below the size of the core.

## Answer to question (1)

In Uncategorized on March 28, 2013 at 7:00 am

Thermal velocity dispersions mean you have a spectral line with some width, and the width is given by thermal broadening, so that $\sigma_T=\sqrt{k_B T/\mu m_H}$ from the Equipartition Theorem. This also happens to be the sound speed!  Is it mere coincidence that thermal velocities are on order the sound speed? No!  Thermal motions are (no surprise) set by the temperature.  The sound speed is set by pressure, since sound waves are just pressure-density waves, and the pressure is also ultimately set by the temperature. So it’s not coincidence that thermal motions are sub-sonic, and supersonic motions cannot be explained by thermal broadening.

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

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

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.

## ARTICLE: The mass function of dense molecular cores and the origin of the IMF (2007)

In Journal Club, Journal Club 2013, Uncategorized on March 12, 2013 at 9:30 pm

## The mass function of dense molecular cores and the origin of the IMF

by  Joao Alves Marco Lombardi & Charles Lada
Summary by Doug Ferrer

## Introduction

One of the main goals of researching the ISM is understanding the connection between the number and properties of stars and the properties of the surrounding galaxy. We want to be able to look at a stellar population and deduce what sort of material it came from, and the reverse–predict what sort of stars we should expect to form in some region given what we know about its properties. The basics of this connection have been known for a while (eg. Bok 1977). Stars form from the gravitational collapse of dense cores of molecular clouds. Thus the properties of stars are the properties of these dense cores modulated by the physical processes that happen during this collapse.

One of the key items we would like to be able to derive from this understanding of star formation is the stellar initial mass function (IMF)–the number of stars of a particular mass as a function of mass. Understanding how the IMF varies from region to region and across time would be extremely useful in many areas of astrophysics, from cosmology to star clusters. In this paper, Alves et al. attempt to explain the origin of the IMF and provide evidence for this explanation by examining the molecular cloud complex in the Pipe nebula. We will look at some background on the IMF, then review the methods used by Alves et al. and asses the implications for star formation.

## The IMF and its Origins

The IMF describes the probability density for any given star to have a mass M, or equivalently the number of stars in a given region with a mass of M. Early work done by Salpeter (1955) showed that the IMF for relatively high mass stars ( $M > 1 M_{sun}$) follows a power law with the index $\alpha \approx -1.3 - 1.4$. For lower masses, the current consensus is for a break below $1 M_{sun}$, and another break below $.1 M_{sun}$ , with a peak around $.3 M_{sun}$. A cartoon of this sort of IMF is shown in Fig. 1.  The actual underlying distribution may in fact be log-normal. This is consistent with stars being formed from collapsing over-densities caused by supersonic turbulence within molecular clouds (Krumholz 2011). This is not particularly strong evidence, however, as the log-normal distribution can result from any process that depends on the product of many random variables.

Figure 1: An illustration of the accepted form of the IMF. There is a power law region for less than $1 M_{sun}$, and a broad peak for lower masses. Image from Krumholz (2011) produced using Chabrier (2003)

The important implication of these results is that there is a characteristic mass scale for star formation of around $.2-.3 M_{sun}$ . There are two (obvious) ways to explain this:

1. The efficiency of star formation peaks at certain scales.
2. There is a characteristic scale for the clouds that stars form from

There has been quite a lot of theoretical work examining option 1 (see Krumhoz 2011 for a relatively recent, accessible review). There are many different physical processes at play in star formation–turbulent MHD, chemistry, thermodynamics, radiative transfer, and gravitational collapse.  Many of these processes are not separately well understood, and each is occurring in its complex, and highly non-linear regime.  We are not even close to a complete description of the full problem. Thus it would not be at all surprising if there were some mass scale, or even several mass scales that are singled out by the star formation process, even if none is presently known. There is some approximate analytical work showing that feedback from stellar outflows may provide such a scale (eg. Shu et al. 1987) . More recent work (e.g. Price and Bate 2008) has shown that magnetic fields cause significant effects on the structure of collapsing cloud cores in numerical simulations, and may reduce or enhance fragmentation depending on magnetic field strength and mass scale.

Nevertheless, the authors are skeptical of the idea that star-formation is not a scale-free process. Per Larson (2005), they do not believe that turbulence or magnetic fields are likely to be very important for the smallest, densest scale clouds that starts ultimately form from. Supersonic turbulence is quickly dissipated on these scales, and magnetic fields are dissipated through ambipolar diffusion–the decoupling of neutral molecules from the ionic plasma. Thus Larson argues that thermal support is the most important process in small cores, and the Jeans analysis will be approximately correct.

The authors thus turn to option 2. It is clear that if the dense cores of star forming clouds already follow a distribution like the IMF, then there will be no need for option 1 as an explanation. Unfortunately though, the molecular cloud mass function (Fig. 2) does not at first glance show any breaks at low mass and has too shallow of a power law index .  But what if we look at only the smallest, densest clumps?

Figure 2: The cumulative cloud mass function (Ak is proortional to mass) for several different cloud complexes from Lombardi et al. (2008). While this is not directly comparable to the IMF, the important take away is that there are no breaks at low mass.

## Observations

Observations using dense gas emission tracers like $C^{18} O$ and $HC^{13}O$ produces mass distributions more like the stellar IMF (eg. Tachihara et al. 2002, Onishi 2002) . However, there are many systematic uncertainties in emission based analysis. To deal with these issues, the authors instead probed dense cloud mass using wide field extinction mapping (this work).   An extinction map was constructed of the nearby Pipe nebula using the NICER method of Lombardi & Alves (2001), which we have discussed in class. This produced the extinction map shown in Fig. 3 below.

Figure 3: NICER Extinction map of the Pipe nebula. A few dense regions are visible, but the noisey, variable background makes it difficult to seperate out seperate cores in a consistent way.

## Data Analysis

The NICER map of the pipe nebula reveals a complex filamentary structure with very high dynamic range in column density ($~ 10^{21} - 10^{22.5} cm^{-2}$). It is difficult to assign cores to regions in such a data set in a coherent way (Alves et al. 2007) using standard techniques–how do we determine what is a core and what is substructure to a core? Since it is precisely the number and distribution of cores that we are interested in, we cannot use a biased method to identify the cores. To avoid this, the authors used a technique called  multiscale wavelet decomposition. Since the authors do not give much information on how this works, we will give a brief overview following a description of a similar technique from Portier-Fozzani et al. (2001).

### Wavelet Decomposition

Wavelet analysis is a hybrid of Fourier and coordinate space techniques. A wavelet is a function that  has a characteristic frequency, position and spatial extent, like the one in Fig.  4. Thus if we convolve a wavelet of a given frequency and length with a signal, it will tell us how the spatial frequency of the signal varies with position. This is the type of analysis used to produce a spectrogram in audio visualization.

Figure 4: A wavelet. This is one is the product of a sinusoid with a gausian envelope. This choice is a perfect tradeoff between spatial and frequency resolution, but other wavelets may be ideal for some other resolution goal. Note that $\Delta x \Delta k \ge 1/2$

The authors used this technique to separate out dense structures from the background. They performed the analysis using wavelets with spatial scales close to three different length scales separated by factors of 2 (.08 pc, .15 pc, .3 pc). Then they combined the structures identified at each length scale into trees based on spatial overlap, and called only the top level of these trees separate cores. The resulting identified cores are shown in fig. 5.  While the cores are shown as circles in fig. 5, this method does not assume spherical symmetry, and the actual identified core shapes are filamentary, reflecting the qualitative appearance of the image

Figure 5: Dense cores identified in the Pipe nebula (circles) from Alves et al. (2007). The radii of the circles is proportional to the core radius determined by the wavelet decomposition analysis

## Results

After obtaining what they claim to be a mostly complete sample of cores, the authors calculate the mass distribution for them. This is done by converting dust extinction to a dust column density. This gives a dust mass for each clump, which can then be extrapolated to a total mass by assuming a set dust fraction.  The result of this is shown in fig. 6 below. The core mass function they obtain is qualitatively similar in shape to the IMF, but scaled by a factor of ~4 in mass. The analysis is only qualitative, and no statistics are done or functions fit to the data. The authors claim that this result evinces a universal star formation efficiency of ~30%, and that this is good agreement with that calculated analytically by (Shu 2004) and numerical simulations. This is again only a qualitative similarity, however.

We should also note that the IMF is hypothesized to be a log-normal distribution. This sort of distribution can arise out of any process that depends multiplicatively on many independent random factors. Thus the fact that dense cores have a mass function that is a scaled version of the IMF is not necessarily good evidence that they share a simple causal link, in the same way that two variables both being normally distributed does not mean that they are any way related.

Figure 6: The mass function of dense molecular cores (points) and the IMF (solid grey line). The dotted gray line is the IMF with mass argument scaled by a factor of 4. The authors note the qualitative agreement, but do not perform any detailed analysis.

## Conclusions

Based on this, the authors conclude that there is no present need for a favored mass scale in star formation as an explanation of the IMF.  Everything can be explained by the distribution produced by dense cloud cores as they are collapsing.  There are a few points of caution however. This is a study of only one cluster, and the data are analyzed using an opaque algorithm that is not publicly available. Additionally, the distributions are not compared statistically (such as with KS), so we have only qualitative similarity.  It would be interesting to see these results replicated for a different cluster using a more transparent statistical method

## References

Bok, B. J., Sim, M. E., & Hawarden, T. G. 1977, Nature, 266, 145

Krumholz, M. R. 2011, American Institute of Physics Conference Series, 1386, 9

Lombardi, M., Lada, C. J., & Alves, J. 2008, A&A, 489, 143

Krumholz, M. R., & Tan, J. C. 2007, ApJ, 654, 304

Alves, J., Lombardi, M., & Lada, C. J. 2007, A&A, 462, L17

Lombardi, M., Alves, J., & Lada, C. J. 2006, A&A, 454, 781

Larson, R. B. 2005, MNRAS, 359, 211

Shu, F. H., Li, Z.-Y., & Allen, A. 2004, Star Formation in the Interstellar Medium: In Honor of
David Hollenbach, 323, 37

Onishi, T., Mizuno, A., Kawamura, A., Tachihara, K., & Fukui, Y. 2002, ApJ, 575, 950

Tachihara, K., Onishi, T., Mizuno, A., & Fukui, Y. 2002, A&A, 385, 909

Portier-Fozzani, F., Vandame, B., Bijaoui, A., Maucherat, A. J., & EIT Team 2001, Sol. Phys.,
201, 271

Shu, F. H., Adams, F. C., & Lizano, S. 1987, ARA&A, 25, 23

Salpeter, E. E. 1955, ApJ, 121, 161

## CHAPTER: The Virial Theorem

In Book Chapter on March 12, 2013 at 3:21 pm

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

## Detecting the phases of the ISM: wild ideas

In Journal Club, Journal Club 2013 on March 7, 2013 at 7:45 pm

These are the ideas reported by the student groups in our 7 March discussion of McKee & Ostriker (1977).

## Group 1

• Use line ratios to discriminate between CNM/WNM
• Use SNR light echoes to determine dust 3D structure

## Group 2

• The key is to disentangle v_thermal from v_bulk
• If the velocity dispersion is dominated by thermal motion, then v ~ m^1/2, and therefore different species should have different line widths. If instead you’re looking at a cold cloud with a velocity distribution dominated by bulk motion, the line width will be independent of the mass of the species.
• Use different distributions along a line of sight. If the implied temperature variance changes as a function of distance along the line of sight, then it may be bulk velocity. if it doesn’t, you’re looking at thermal velocity.
• Use face on spiral
• Use dust content as proxy, since it will only sublimate at high temperatures
• Use hyperfine lines kin a transition which have different optical depth to look into different shells of a cloud

## Group 3

• conduct a large survey of 21 cm emission
• goal: see whether WNM is a distinct phase, or just the product of several CNM clouds moving quickly
• method: try to distinguish between line widths and line shifts
• need very high-resolution spectrograph!
• look for dust towards a region that might be WNM
• dust sublimates at ~2000 K
• so, dust would exist if the region is actually CNM, but not if it’s WNM
• probe a cloud in a “layered” approach
• test the proposed structure of CNM core with WNM, WIM in successive layers
• observe at a range of wavelengths that probe different layers (perhaps at a range of wavelengths that become optically thick at different layers)
• collaborate with other alien astronomers
• get different lines of sight!

## Group 4

• 21-cm observations at different galactic latitudes
• count clouds in a volume to get a filling factor (use CO or IR emission from dust)
• look at other galaxies
• look for bubbles along the line of sight

## CHAPTER: Ion-Neutral Reactions

In Book Chapter on March 7, 2013 at 3:20 pm

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

## CHAPTER: Neutral-Neutral Interactions

In Book Chapter on March 7, 2013 at 3:19 pm

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

## CHAPTER: Excitation Processes: Collisions

In Book Chapter on March 7, 2013 at 3:18 pm

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

## CHAPTER: Definitions of Temperature

In Book Chapter on March 7, 2013 at 3:27 am

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

## CHAPTER: Important Properties of Local Thermodynamic Equilibrium

In Book Chapter on March 7, 2013 at 3:25 am

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

## The “True” Column Density Distribution in Star-Forming Molecular Clouds

In Journal Club 2013 on March 5, 2013 at 8:25 pm

## CHAPTER: The Saha Equation

In Book Chapter on March 5, 2013 at 3:21 am

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

## CHAPTER: Spitzer Notation

In Book Chapter on March 5, 2013 at 3:19 am

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

## CHAPTER: Thermodynamic Equilibrium

In Book Chapter on February 28, 2013 at 3:13 am

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

## CHAPTER: Introductory Remarks on Radiative Processes

In Book Chapter on February 28, 2013 at 3:10 am

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

## CHAPTER: Relevant Velocities in the ISM

In Book Chapter on February 28, 2013 at 3:06 am

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

## CHAPTER: Energy Density Comparison

In Book Chapter on February 26, 2013 at 3:04 am

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

## CHAPTER: Measuring States in the ISM

In Book Chapter on February 26, 2013 at 3:00 am

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

#### SEDs

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)

#### Spectra

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

Key

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)

## ARTICLE: Dark Nebulae, Globules, and Protostars

In Journal Club, Journal Club 2013 on February 19, 2013 at 10:45 pm

### Dark Nebulae, Globules, and Protostars by Bart Bok (1977)

Summary by George Miller

### Introduction

In Bart Bok’s 1977 paper Dark nebulae, globules, and protostars (Bok 1977), largely based on a lecture given upon acceptance of the Astroomical Society of the Pacific’s Bruce Medal, he presents two fundamentally different pictures of star formation. The first, constituting the majority of the paper’s discussion, occurs in large Bok globules which are compact, rounded and remarkably well-defined regions of high-extinction ranging from 3′ to 20′.  The globules show a strong molecular hydrogen and dust component and relatively little signs of higher neutral HI concentrations than its surroundings. In contrast, Bok briefly examines star formation in the Magellanic Clouds which show a vast amount of neutral atomic hydrogen and a comparatively small amount of cosmic dust. In this review, I will summarize a number of key points made by Bok, as well as provide additional information and modern developments since the paper’s original publishing.

### Large Bok Globules

#### A history of observations

In 1908, Barnard drew attention to “a number of very black, small, sharply defined spots or holes” in observations of the emission nebula Messier 8 (Barnard 1908).  39 years later Bok published extensive observations of 16 “globules” present in M8 as well others in $\eta$ Carinae, Sagittarius, Ophiuchus and elsewhere, making initial estimates of their distance, diameter and extinction (Bok & Reilly 1947). He further claimed that these newly coined “globules” were gravitationally contracting clouds present just prior to star formation, comparing them to an “insect’s cocoon” (Bok 1948). As we will see, this bold prediction was confirmed over 40 years later to be correct. Today there over 250 globules known within roughly 500 pc of our sun and, as Bok claims in his 1977 paper, identifying more distant sources is difficult due to their small angular diameter and large number of foreground stars.  There are currently four chief methods of measuring the column density within Bok Globules: extinction mappings of background stars, mm/sub-mm dust continuum emission, absorption measurements of galactic mid-IR background emission, and mapping molecular tracers.  See Figure 1 for a depiction of the first three of these methods.  At the time Bok published his paper in 1977, only extinction mapping and molecular tracer methods were readily available, thus I will primarily discuss these two.  For a more in depth discussion, see Goodman et. al. 2009 and the subsequent AST201b Journal Club review.

Figure 1.  Three methods of determining column density of starless molecular cores or Bok globules. (a) K-band image of Barnard 68 and plot of the $A_K$ as a function of radius from the core.  This method measures the H–K excess, uses the extinction law to convert into $A_V$, and then correlated to the $H_2$ column density from UV line measurements, parameterized by f. (b) 1.2-mm dust continuum emission map and ﬂux versus radius for L1544.  $\kappa_{\nu}$ is the dust opacity per unit gas mass, ρ is the dust density, and m the hydrogen mass (corrected for He). (c) 7-μm ISOCAM image and opacity versus radius for ρ Oph D.  In this method the absorbing opacity is related to the hydrogen column via the dust absorption cross section, $\sigma_{\lambda}$.  Figure taken from Bergin & Tafalla 2007.

#### Measuring photometric extinction

Measuring the photometric absorption, and thus yielding a minimum dust mass, for these globules is itself an arduous process. For globules with $A_v<10$   mag, optical observations with large telescopes can be used to penetrate through the globules and observe the background stars.  Here $A_{\lambda} \equiv m_{\lambda}-m_{\lambda, 0} = 2.5 \, log(\frac{F_{\lambda,0}}{F_{\lambda}})$.  Thus an extinction value of $A_v=10$ mag means the flux is decreased by a factor of $10^4$.  By using proper statistics of the typical star counts and magnitudes seen within a nearby unobstructed field of view, extinction measurements can be made for various regions.  It is important to note that the smaller an area one tries to measure an extinction of, the greater the statistical error (due to a smaller number of background stars).  This is one of the key limitations of extinction mappings.  For the denser cores or more opaque globules with $10 < A_V < 20$ mag, observations in the near infrared are needed (which is relatively simple by today’s standards but not so during Bok’s time). This is further complicated due to imprecisely defined BVRI photometric standard sequences for fainter stars, a problem still present today with various highly-sensitive space telescopes such as the HST. Bok mentions two methods. In the past a Racine (or Pickering) prism was used to produce fainter companion images of known standards, yet as discussed by Christian & Racine 1983 this method can produce important systematic errors. The second, and more widely used, method is to pick an easily accessible progression of faint stars and calibrate all subsequent photographic plates (or ccd images) from this. See Saha et. al. 2005 for a discussion of this problem in regards to the Hubble Space Telescope.

Obtaining an accurate photometric extinction for various regions within the globule is valuable as it leads an estimate of the dust density. Bok reports here from his previous Nature paper (Bok et. al. 1977) that the extinction $A_v$ within the Coalsack Globule 2 varies inversely as the square of distance, thus also implying the dust density varies inversely as the cube of distance from the core.  Modern extinction mappings, as seen in Figure 1(a) of Barnard 68,  show that at as one approaches the central core the extinction vs. distance relation actually flattens out nearly to $r^{-1}$.  This result was a key discovery, for the Bonnor-Ebert (BE) isothermal sphere model predicts a softer power law at small radii.  In his paper, Bok remarks “The smooth density gradient seems to show that Globule 2 is […] an object that reminds one of the polytropic models of stars studied at the turn of the century by Lane (1870) and Emden (1907)”.  It is truly incredible how accurate this assessment was.  The Bonnor-Ebert sphere is a model derived from the Lane-Emden equation for an isothermal, self-gravitating sphere which remains in hydrostatic equilibrium.  Figure 2 displays a modern extinction mapping of Barnard 68 along with the corresponding BE sphere model, showing that the two agree remarkably well.  There are, however, a number of detractors from the BE model applied to Bok globules.  The most obvious is that globules are rarely spherical, implying that some other non-symmetric pressure must be present.  Furthermore, the density gradient between a globule’s core and outer regions often exceeds 14 ($\xi_{max} > 6.5$) as required for a stable BE sphere (Alves, Lada & Lada 2001).

Figure 2.  Dust extinction mapping for Barnard 68 plotted against an isothermal Bonnor-Ebert sphere model.  Figure taken from Alves, Lada & Lada 2001.

#### Using CO as a tracer

Important tracer molecules, such as CO, are used to study the abundance of $H_2$, temperatures and kinematics of these globules. Because the more common $^{12}CO$ isotope tends to become optically thick and saturate in regions of higher column density such as globules, the strength of $^{13}CO$ emission is usually used to indicate the density of H2.  The conversion factor of $N_{H_2} = 5.0 \pm 2.5 \times 10^5 \times N_{13},$ from Dickman 1978 has changed little in over three decades. The column density of $H_2$, combined with its known mass and radius of the globule, can then be used to estimate the globule’s total mass. Furthermore, the correlation of $^{13}CO$ density with photometric extinction, $A_v = 3.7 \times 10^{-16} \times N_{13},$ is another strong indication that $^{13}CO$ emission is an accurate tracer for H$_2$ and dust. Further studies using $C^{17}O$ and $C^{18}O$ have also been used to trace even higher densities when even $^{13}CO$ can become optically thick(Frerking et. al. 1982).  As an example, Figure 3 shows molecular lines from the central region of the high-mass star forming region G24.78+0.08.  In the upper panel we can see the difference between the optically thick $^{12}CO$ and thin $C^{18}O$.  The $^{12}CO$ line shows obvious self-absorption peaks associated with an optically thick regime, and one clearly can not make a Gaussian fit to determine the line intensity.   $^{12}CO$, due to the small dipole moment of its $J=1 \rightarrow 0$ transition and thus ability to thermalize at relatively low densities, is also used to measure the gas temperature within globules. These temperatures usually range from 7K to 15K. Finally, the width of CO lines are used to measure the velocity dispersion within the globule. As Bok states, most velocities range from 0.8 to 1.2 km/s. This motion is often complex and measured excess line-widths beyond their thermal values are usually attributed to turbulence (Bergin & Tafalla 2007). Importantly, the line-width vs. size relationship within molecular clouds first discovered by Barnard 1981 does not extend to their denser cores (which have similar velocity motions as Bok globules).  Instead, a “coherence” radius is seen where the non-thermal component of a linewidth is approximately constant (Goodman et. al. 1998).  In the end, as Bok surmises, the subsonic nature of this turbulence implies it plays a small role compared to thermal motions.

Figure 3.  Spectra taken from the core of the high-mass star forming region G24.78+0.08.  The solid line corresponds to $^{12}CO (1\rightarrow 0)$, $^{12}CO (2\rightarrow 1)$, and $C^{32}S (3\rightarrow 2)$, the dashed line to $^{13}CO (1\rightarrow 0)$$^{13}CO (2\rightarrow 1)$, and $C^{34}S (3\rightarrow 2)$ and the dotted line to $C^{18}O (1\rightarrow 0)$.  From the top panel, one can clearly see the difference between the optically thick, saturated $^{12}CO (1\rightarrow 0)$ line and the optically thin $C^{18}O (1\rightarrow 0)$ transition.  Figure taken from Cesaroni et. al. 2003.

#### The current status of Bok globules

Today, the majority of stars are thought to originate within giant molecular clouds or larger dark cloud complexes, with only a few percent coming from Bok globules. However, the relative simplicity of these globules still make them important objects for studying star formation. While an intense debate rages today regarding the influence of turbulence, magnetic fields, and other factors on star formation in GMCs, these factors are far less important than simple gravitational contraction within Bok globules. The first list of candidate protostars within Bok globules, obtained by co-adding IRAS images, was published in 1990 with the apropos title “Star formation in small globules – Bart Bok was correct” (Yun & Clemens 1990).  To conduct the search, Yun & Clemens first fit a single-temperature modified blackbody model the the IRAS 60 and 100 μm images (after filtering out uncorrelated background emission) to obtain dust temperature and optical depth values.  This result was then used as a map to search for spatially correlated 12 and 25 μm point sources (see Figure 4.).  More evidence of protostar outflows (Yun & Clemens 1992), Herbig-Haro objects due to young-star jets (Reipurth et al. 1992) and the initial stages of protostar collapse (Zhou et. al. 1993) have also been detected within Bok Globules. Over 60 years after Bok’s pronouncement that these globules were “insect cocoons” encompassing the final stages of protostar formation, his hypothesis remains remarkably accurate and validated. It is truly “pleasant indeed that globules are there for all to observe!”

Figure 4.    (a) Contour map of the dust temperature $T_{60/100}$ of the Bok Globule CB60 derived from 60 and 100 μm IRAS images.  (b) 12 μm IRAS image of CB60 after subtracting background emission using median-filtering.  This source is thought to be a young stellar object or protostar located within the globule.  The other 12 μm field sources seen in (b) are thought not to be associated with the globule. Figure taken from Yun & Clemens 1990.

### Magellanic Cloud Star Formation

At the end of his paper, Bok makes a 180 degree turn and discusses the presence of young stars and blue globulars within the Magellanic Clouds. These star formation regions stand in stark contrast to the previously discussed Bok globules; they contain a rich amount of HI and comparatively small traces of dust, they are far larger and more massive, and they form large clusters of stars as opposed to more isolated systems. Much more is known of the star-formation history in the MCs since Bok published this 1977 paper. The youngest star populations in the MCs are found in giant and supergiant shell structures which form filamentary structures throughout the cloud. These shells are thought to form from supernova, ionizing radiation and stellar wind from massive stars which is then swept into the cool, ambient molecular clouds. Further gravitational, thermal and fluid instabilities fragment and coalesce these shells into denser star-forming regions and lead to shell-shell interactions (Dawson et. al. 2013). The initial onset of this new ($\sim$ 125 Myr) star formation is thought to be due to close encounters between the MCs, and is confirmed by large-scale kinematic models (Glatt et al. 2010).

### References

Alves, J. F., Lada, C. J., & Lada, E. A. 2001, Nature, 409, 159
Barnard, E. E. 1908, Astronomische Nachrichten, 177, 231
Bergin, E. A., & Tafalla, M. 2007, ARA&A, 45, 339
Bok, B. J. 1948, Harvard Observatory Monographs, 7, 53
—. 1977, PASP, 89, 597
Bok, B. J., & Reilly, E. F. 1947, ApJ, 105, 255
Bok, B. J., Sim, M. E., & Hawarden, T. G. 1977, Nature, 266, 145
Cesaroni, R., Codella, C., Furuya, R. S., & Testi, L. 2003, A&A, 401, 227
Christian, C. A., & Racine, R. 1983, PASP, 95, 457
Dawson, J. R., McClure-Griffiths, N. M., Wong, T., et al. 2013, ApJ, 763, 56
Dickman, R. L. 1978, ApJS, 37, 407
Frerking, M. A., Langer, W. D., & Wilson, R. W. 1982, ApJ, 262, 590
Glatt, K., Grebel, E. K., & Koch, A. 2010, A&A, 517, A50
Goodman, A. A., Barranco, J. A., Wilner, D. J., & Heyer, M. H. 1998, ApJ, 504, 223
Goodman, A. A., Pineda, J. E., & Schnee, S. L. 2009, ApJ, 692, 91
Reipurth, B., Heathcote, S., & Vrba, F. 1992, A&A, 256, 225
Saha, A., Dolphin, A. E., Thim, F., & Whitmore, B. 2005, PASP, 117, 37
Yun, J. L., & Clemens, D. P. 1990, ApJ, 365, L73
—. 1992, ApJ, 385, L21
Zhou, S., Evans, II, N. J., Koempe, C., & Walmsley, C. M. 1993, ApJ, 404, 232

## CHAPTER: Chemistry

In Book Chapter on February 13, 2013 at 10:04 pm

See Draine Table 1.4 for elemental abundances for the proto-solar environment. H:He:C = $1:0.1:3 \times 10^{-4}$ by number. $1:0.4:3.5 \times 10^{-3}$ by mass 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 > Carbon) is twice as low at the sun’s position than in the Galactic center. Even though metals account for 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. Generally, it is easier (i.e. requires less energy) to dissociate a molecule than to ionize something. The lower the electronic state you are trying to ionize, the more energy is needed. The Lyman Limit is the minimum photon energy needed to ionize Hydrogen from the ground state (13.6 eV, 912 Angstrom).

## CHAPTER: Hydrogen Slang

In Book Chapter on February 12, 2013 at 10:02 pm

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 \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 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)$.

Note that the Lyman (or Balmer, Paschen, etc.) limit can be computed by inserting $n_i=\infty$.

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.

## ARTICLE: The Physical State of Interstellar Hydrogen

In Journal Club 2013 on February 12, 2013 at 9:57 pm

The Physical State of Interstellar Hydrogen by Bengt Strömgren (1939)

Summary by Anjali Tripathi

Abstract

In 1939, Bengt Strömgren published an analytic formulation for the spatial extent of ionization around early type stars.  Motivated by new H-alpha observations of sharply bound “diffuse nebulosities,” Strömgren was able to characterize these ionized regions and their thin boundaries in terms of the ionizing star’s properties and abundances of interstellar gas.  Strömgren’s work on these regions, which have come to be eponymously known as Strömgren spheres, has found longstanding use in the study of HII regions, as it provides a simple analytic approach to recover the idealized properties of such systems.

Background: Atomic Physics in Astronomy & New Observations

Danish astronomer Bengt Strömgren (1908-87) was born into a family of astronomers and educated during a period of rapid development in our understanding of the atom and modern physics.  These developments were felt strongly in Copenhagen where Strömgren studied and worked for much of his life.  At the invitation of Otto Struve, then director of Yerkes Observatory, Strömgren visited the University of Chicago from 1936 to 1938, where he encountered luminaries from across astrophysics, including Chandrasekhar and Kuiper.  With Struve and Kuiper, Strömgren worked to understand how photoionization could explain observations of a shell of gas around an F star, part of the eclipsing binary $\epsilon$ Aurigae (Kuiper, Struve and Strömgren, 1937).  This work laid out the analytic framework for a bounded region of ionized gas around a star, which provided the theoretical foundation for Strömgren’s later work on HII regions.

The observational basis for Strömgren’s 1939 paper came from new spectroscopic measurements taken by Otto Struve.  Using the new 150-Foot Nebular Spectrograph (Struve et al, 1938) perched on a slope at McDonald Observatory, pictured below, Struve and collaborators were able to resolve sharply bound extended regions “enveloped in diffuse nebulosities” in the Balmer H-alpha emission line (Struve and Elvey, 1938).  This emission line results from recombination when electrons transition from the n = 3 to n = 2 energy level of hydrogen, after the gas was initially ionized by UV radiation from O and B stars.  Comparing these observations to those of the central parts of the Orion Nebula led the authors to estimate that the number density of hydrogen with electrons in the n=3 state is $N_3 = 3 \times 10^{-21} cm^{-3}$, assuming a uniform concentration of stars and neglecting self-absorption (Struve and Elvey, 1938).  From his earlier work on $\epsilon$ Aurigae, Strömgren had an analytic framework with which to understand these observations.

Instrument used to resolve HII Regions in H-alpha (Struve et al, 1938)

Putting it together – Strömgren’s analysis

To understand the new observations quantitatively, Strömgren worked out the size of these emission nebulae by finding the extent of the ionized gas around the central star.  As in his paper with Kuiper and Struve, Strömgren considered only neutral and ionized hydrogen, assumed charge neutrality, and used the Saha equation with additional terms:

${N'' N_e \over N'} = \underbrace{{(2 \pi m_e)^{3/2} \over h^3} {2q'' \over q'} (kT)^{3/2}e^{-I/kT}}_\text{Saha} \cdot \underbrace{\sqrt{T_{el} \over T}}_\text{Temperature correction} \cdot \underbrace{R^2 \over 4 s^2}_\text{Geometrical Dilution}\cdot \underbrace{e^{-\tau_u}}_\text{Absorption}\\ N': \text{Neutral hydrogen (HI) number density}\\ N'':\text{Ionized hydrogen (HII) number density}\\ N_e:\text{Electron number density, }N_e = N''\text{ by charge neutrality}\\ x: \text{Ionization fraction}, x = N''/(N'+N'')$

Here, the multiplicative factor of $\sqrt{T_{el} \over T}$ corrects for the difference between the stellar temperature($T$) and the electron temperature($T_{el}$) at a distance $s$ away from the star.  The dilution factor ${R^2 \over 4 s^2}$, where $R$ is the stellar radius and $s$ is the distance from the star, accounts for the decrease in stellar flux with increasing distance.  The factor of $e^{-\tau_u}$, where $\tau_u$ is the optical depth, accounts for the reduction in the ionizing radiation due to absorption.  Taken together, this equation encapsulates the physics of a star whose photons ionize surrounding gas.  This ionization rate is balanced by the rate of recombination of ions and electrons to reform neutral hydrogen.  As a result, close to the star where there is abundant energetic flux, the gas is fully ionized, but further from the star, the gas is primarily neutral.  Strömgren’s formulation allowed him to calculate the location of the transition from ionized to neutral gas and to find the striking result that the transition region between the two is incredibly sharp, as plotted below.

Plot of ionization fraction vs. distance for an HII Region (Values from Table 2 of Strömgren, 1939)

Strömgren found that the gas remains almost completely ionized until a critical distance $s_0$, where the ionization fraction sharply drops and the gas becomes neutral due to absorption.  This critical distance has become known as the Strömgren radius, considered to be the radius of an idealized, spherical HII region.  The distance over which the ionization fraction drops from 1 to 0 is small (~0.01 pc), corresponding to one mean free path of an ionizing photon, compared to the Strömgren radius(~100pc).  Thus Strömgren’s analytic work provided an explanation for sharply bound ionized regions with thin transition zones separating the ionized gas from the exterior neutral gas.

Strömgren also demonstrated how the critical distance depends on the total number density $N$, the stellar effective temperature $T$, and the stellar radius $R$:

$\log{s_0} = -6.17 + {1 \over 3} \log{ \left( {2q'' \over q'} \sqrt{T_{el} \over T} \right)} - {1 \over 3} \log{a_u} - {1 \over 3} \frac{5040K}{T} I + {1 \over 2} \log{T} + {2 \over 3} \log{R} - {2 \over 3} \log{N},$

where $a_u$ is the absorption coefficient for the ionizing radiation per hydrogen atom (here assumed to be frequency independent) and $s_0$ is given in parsecs.  From this relation, we can see that for a given stellar radius and a fixed number density, $s_0 \propto T^{1/2}$, so that hotter, earlier type stars have larger ionized regions.  Plugging in numbers, Strömgren found that for a total number density of $3~cm^{-3}$, a cluster of 10 O7 stars would have a critical radius of 100-150 parsecs, in agreement with estimates made by the Struve and Elvey observations.

To estimate the hydrogen number density from the H-alpha observations, Strömgren also considered the excitation of the n=3 energy levels of hydrogen.  Weighing the relative importance of various mechanisms for excitation – free electron capture, Lyman-line absorption, Balmer-line absorption, and collisions – Strömgren found that their effects on the number densities of the excited states and electron number densities were comparable.  As a result, he estimated from Struve’s and Elvey’s $N_3$ that the number density of hydrogen is 2-3 $cm^{-3}$.

Strömgren’s analysis of ionized regions around stars and neutral hydrogen in “normal regions” matched earlier theoretical work by Eddington into the nature of the ISM (Eddington, 1937).  “With great diffidence, having not yet rid myself of the tradition that ‘atoms are physics, but molecules are chemistry’,” Eddington wrote that “presumably a considerable part” of the ISM is molecular.  As a result, Strömgren outlined how his analysis for ionization regions could be modified to consider regions of molecular hydrogen dissociating, presciently leaving room for the later discovery of an abundance of molecular hydrogen.  Instead of the ionization of atomic hydrogen, Strömgren worked with the dissociation of molecular hydrogen in this analysis.   Given that the energy required to dissociate the bond of molecular hydrogen is less than that required to ionize atomic hydrogen, Stromgren’s analysis gives a model of a star surrounded by ionized atoms, which is surrounded by a sharp, thin transition region of atomic hydrogen, around which molecular hydrogen remains.

In addition to HI and HII, Strömgren also considered the ionization of other atoms and transitions.  For example, Strömgren noted that if the helium abundance was smaller than that of hydrogen, then most all of the helium will be ionized out to the boundary of the hydrogen ionization region.  From similar calculations and considering the observations of Struve and Elvey, Strömgren was able to provide an estimate of the abundance of OII, a ratio of $10^{-2}-10^{-3}$ oxygen atoms to each hydrogen atom.

Strömgren Spheres Today

Strömgren’s idealized formulation for ionized regions around early type stars was well received initially and  has continued to influence thinking about HII regions in the decades since.  The simplicity of Strömgren’s model and its assumptions, however, have been recognized and addressed over time.  Amongst these are concerns about the assumption of a uniformly dense medium around the star.  Optical and radio observations, however, have revealed that the surrounding nebula can have clumps and voids – far from being uniformly dense (Osterbrock and Flather, 1959).  To address this, calculations of the nebula’s density can include a ‘filling factor’ term.  Studies of the Orion Nebula (M42), pictured below, have provided examples of just such clumpiness.  M42 has also been used to study another related limitation of Strömgren’s model – the assumption of a central star surrounded by spherical symmetry.

Orion Nebula, infrared image from WISE. Credit: NASA/JPL/Caltech

Consideration of the geometry of Strömgren spheres has been augmented by blister models of the 1970s whereby a star ionizes surrounding gas but the star is at the surface or edge of a giant molecular cloud (GMC), rather than at the center of it.  As a result, ionized gas breaks out of the GMC, like a popping blister, which in turn can prompt “champagne flows” of ionized gas leaching into the surrounding medium.  In a review article of Strömgren’s work, Odell (1999) states that due to observational selection effects, many HII regions observed in the optical actually are more akin to blister regions, rather than Strömgren spheres, since Strömgren spheres formed at the center or heart of a GMC may be obscured so much that they are observable only at radio wavelengths.

In spite of its simplifying assumptions, Strömgren’s work remains relevant today.   Given its abundance, hydrogen dominates the physical processes of emission nebulae and, thus, Strömgren’s idealized model provides a good first approximation for the ionization structure, even though more species are involved than just atomic hydrogen.  Today we can enhance our understanding of these HII regions using computer codes, such as CLOUDY, to calculate the ionization states of various atoms and molecules.  We can also computationally model the  hydrodynamics  of shocks radiating outwards from the star and use spectral synthesis codes to produce mock spectra.  From these models and the accumulated wealth of observations over time, we have come to accept that dense clouds of molecular gas, dominated with molecular hydrogen, are the sites of star formation.  Young O and B-type stars form out of clumps in these clouds and their ionizing radiation will develop into an emission nebula with ionized atomic hydrogen, sharply bound from the surrounding neutral cloud.  As the stars age and the shocks race onwards, the HII regions will evolve.  What remains, however, is Strömgren’s work which provides a simple analytic basis for understanding the complex physics of HII regions.

Selected references (in order of appearance in this article)

Strömgren, “The Physical State of Interstellar Hydrogen”, ApJ (1939)

Kuiper et al, “The Interpretation of $\epsilon$  Aurigae”, ApJ (1937)

Struve et al, “The 150-Foot Nebular Spectrograph of the McDonald Observatory”, ApJ (1938)

Eddington, “Interstellar Matter”, The Observatory (1937)

Osterbrock and Flather, “Electron Densities in the Orion Nebula. II”, ApJ (1959)

O’Dell, “Strömgren Spheres”, ApJ (1999)

## CHAPTER: Stromgren Sphere: An example “chalkboard derivation”

In Book Chapter on February 12, 2013 at 8:31 pm

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

## ARTICLE: On the Dark Markings in the Sky

In Journal Club, Journal Club 2013 on February 8, 2013 at 2:46 pm

On the Dark Markings in the Sky by Edward E. Barnard (1919)

Summary by Hope Chen

#### Abstract

By examining photographic plates of various regions on the sky, Edward E. Barnard concluded in this paper that what he called “dark markings” were in fact due to the obscuration of nearby nebulae in most cases. This result had a significant impact on the debate regarding the size and the dimension of the Milky Way and also the research of the interstellar medium, particularly work by Vesto Slipher, Heber Curtis and Robert Trumpler. The publication of  Photographic Atlas of Selected Regions of the Milky Way after Barnard’s death, which included many of the regions mentioned in the paper, further provided a new method of doing astronomy research. In this paper and the Atlas, we are also able to see a paradigm very different from that of today.

It is now well-known that the interstellar medium causes extinction of light from background stars. However, think of a time when the infrared imaging was impossible, and the word “photon” meant nothing but a suspicious idea. Back in such a time in the second decade of the twentieth century, Edward Edison Barnard, by looking at hundreds of photographic plates, proposed an insightful idea that “starless” patches of the sky were dark because they are obscured by nearby nebulae. This idea not only built the foundation of the modern concept of the interstellar medium, but also helped astronomers figure out that the Universe extended so much farther beyond the Milky Way.

#### Young Astronomer and His Obsession of the Sky

In 1919, E. E. Barnard published this paper and raised the idea that what he called “dark markings” are mostly obscuration from nebulae close to us. The journey, however, started long before the publication of this paper. Born in Nashville, Tennessee in 1857, Barnard was not able to receive much formal education owing to poverty. His first interest, which became important for his later career, was in photography. He started working as a photographer’s assistant at the age of nine, and the work continued throughout most of his teenage years. He then developed an interest in astronomy, or rather, “star-gazing,” and would go watch the sky almost every night with his own telescope. He took courses in natural sciences at Vanderbilt University and started his professional career as an astronomer at the Lick Observatory in 1888. He helped build the Bruce Photographic Telescope at the Lick Observatory and there he started taking pictures of the sky on photographic plates. He then moved on to his career at the Yerkes Observatory at Chicago University and worked there until his death in 1922. (Introduction of the Atlas, Ref. 2)

Fig. 1 One of the many plates in the Atlas including the region around Rho Ophiuchii, which was constantly mentioned in many of Barnard’s works. (Ref. 2)

Fig. 1 is one of the many plates taken at the Yerkes Observatory. It shows the region near Rho Ophiuchii, which was a region constantly and repetitively visited by Barnard and his telescope. Barnard noted in his description of this plate, “the [luminous] nebula itself is a beautiful object. With its outlying connections and the dark spot in which it is placed and the vacant lanes running to the East from it, … it gives every evidence that it obscures the stars beyond it.” Numerous similar comments spread throughout his descriptions of various regions covered in A Photographic Atlas of Selected Regions of the Milky Way (hereafter, the Atlas). Then finally in his 1919 paper, he concluded, “To me these are all conclusive evidence that masses of obscuring matter exist in space and are readily shown on photographs with the ordinary portrait lenses,” although “what the nature of this matter may be is quite another thing.” The publication of these plates in the Atlas (unfortunately after his death, put together by Miss Mary R. Calvert, who was Barnard’s assistant at the Yerkes Observatory and helped publish many of Barnard’s works after his death) also provided a new way of conducting astronomical research just as the World Wide Telescope does today. The Atlas for the first time allowed researchers to examine the image and the astronomical coordinates along with lists of dominant objects at the same time.

Except quoting Vesto Slipher’s work on spectrometry measurements of these nebulae, most of the evidences in Barnard’s paper seemed rather qualitative than quantitative. So, as of today’s standard, was the “evidence” really conclusive? Again, the question cannot be answered without knowing the limits of astronomical research at the time. Besides an immature understanding of the underlying physics, astronomers in the beginning of the twentieth century were limited by the lack of tools on both the observation and analysis fronts. Photographic plates as those in the Atlas were pretty much the most advanced imaging technique at the time, on which even a quantitative description of “brightness” was not easy, not to mention an estimation of the extinction of these “dark markings.” However, this being said, a very meaningful and somewhat “quantitative” assumption was drawn in Barnard’s paper: the field stars were more or less uniformly distributed. Barnard came to this assumption by looking at many different places, both in the galactic plane and off the plane, and observing the densities of field stars in these regions. Although numbers were not given in the paper, this was inherently similar to a star count study. Eventually, this assumption lead to what Barnard thought as the conclusive evidence of these dark markings being obscuring nebulae instead of “vacancies.” Considering the many technical limits at the time, while the paper might not seem to be scientific in today’s standard, this paper did pose a “conclusion” which was strong enough to sustain many of the more quantitative following examinations.

#### The “Great Debate”

Almost at the same time, perviously mentioned Vesto Slipher (working at the Lowell Observatory) began taking spectroscopic measurements of various clouds and tried to understand the constituents of these clouds. Although limited by the wavelength range and the knowledge of different radiative processes (the molecular transition line emission used largely in the research of the interstellar medium today was not observed until half a century later in 1970, by Robert Wilson, who, on a side note, also discovered the Cosmic Microwave Background), Slipher was able to determine the velocities of clusters by measuring the Doppler shifts and concluded that many of these clusters move at a faster rate than the escape velocity of the Milky Way (Fig. 2). This result, coupled with Barnard’s view of intervening nebulae, revolutionized the notion of the Universe in the 1920s.

Fig. 2 The velocity measurements from spectroscopic observations done by Vesto Slipher. (Ref. 3)

On April 26, 1920 (and in much of the 1920s), the “Great Debate” took place between Harlow Shapley (the Director of Harvard College Observatory at the time) and Curtis Heber (the Lick Observatory, 1902 to 1920). The general debate concerned the dimension of the Universe and the Milky Way, but the basic issue was simply whether distant “spiral nebulae” were small and lay within the Milky Way or whether they were large and independent galaxies. Besides the distance and the velocity measurements, which suffered from large uncertainties due to the technique available at the time, Curtis Heber was able to “win” the debate by claiming that dark lanes in the “Great Andromeda Nebula” resemble local dark clouds as those observed by Barnard (Fig. 3, taken in 1899). The result of the debate then sparked a large amount of work on “extragalactic astronomy” in the next two decades and was treated as the beginning of this particular research field.

Fig. 3 The photographic plate of the “Great Andromeda Nebula” taken in 1988 by Isaac Roberts.

#### The Paper Finally Has a Plot

Then after the first three decades of the twentieth century, astronomers were finally equipped with a relatively more correct view of the Universe, the idea of photons and quantum theory. In 1930, Robert J. Trumpler (the namesake of the Trumpler Award) published his paper about reddening and reconfirmed the existence of local “dark nebulae.” Fig. 4 shows the famous plot in his paper which showed discrepancies between diameter distances and photometric distances of clusters. In the same paper, Trumpler also tried to categorize effects of the ISM on light from background stars, including what he called “selective absorption” or reddening as it is known today. This paper, together with many of Trumpler’s other papers, is one of the first systematic research results in understanding the properties of Barnard’s dark nebulae, which are now known under various names such as clouds, clumps, and filaments, in the interstellar medium.

Fig. 4 Trumpler’s measurements of diameter distances v. photometric distances for various clusters.

#### Moral of the Story

As Alyssa said in class, it is often more beneficial than we thought to understand what astronomers knew and didn’t know at different periods of time and how we came to know what we see as common sense today, not only in the historically interesting sense but also in the sense of better understanding of various ideas. In this paper, Barnard demonstrated a paradigm which we may call unscientific today but made a huge leap into what later became the modern research field of the interstellar medium.

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