Spatial Structure and Collisionless Electron Heating in Balmer-dominated Shocks

Spatial Structure and Collisionless Electron Heating in Balmer-dominated Shocks

Matthew van Adelsberg11affiliation: JILA, University of Colorado, 440 UCB, Boulder, CO 80309; , Kevin Heng22affiliation: Institute for Advanced Study, Einstein Drive, Princeton, NJ 08540 , Richard McCray11affiliation: JILA, University of Colorado, 440 UCB, Boulder, CO 80309; , John C. Raymond33affiliation: Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138

Balmer-dominated shocks in supernova remnants (SNRs) produce strong hydrogen lines with a two-component profile composed of a narrow contribution from cold upstream hydrogen atoms, and a broad contribution from hydrogen atoms that have undergone charge transfer reactions with hot protons. Observations of emission lines from edge-wise shocks in SNRs can constrain the gas velocity and collisionless electron heating at the shock front. Downstream hydrogen atoms engage in charge transfer, excitation and ionization reactions, defining an interaction region called the shock transition zone. The properties of hot hydrogen atoms produced by charge transfers (called broad neutrals) are critical for accurately calculating the structure and radiation from the shock transition zone. This paper is the third in a series describing the kinetic, fluid and emission properties of Balmer-dominated shocks, and is the first to properly treat the effect of broad neutral kinetics on shock transition zone structure. We use our models to extract shock parameters from observations of Balmer-dominated SNRs. We find that inferred shock velocities and electron temperatures are lower than those of previous calculations by for km s, and by for km s. This effect is primarily due to the fact that excitation by proton collisions and charge transfer to excited levels favor the high speed part of the neutral hydrogen velocity distribution. Our results have a strong dependence on the ratio of electron to proton temperatures, , which allows us to construct a relation between the temperature ratio and shock velocity. We compare our calculations to previous results by Ghavamian et al. (2007b).

shock waves — supernova remnants

1 Introduction

Balmer-dominated shocks in supernova remnants (SNRs) encounter upstream gas containing a substantial fraction of neutral hydrogen atoms. Emission from the shock is characterized by strong Balmer and Lyman lines with two-component profiles, which consist of a narrow contribution from direct excitation of the hydrogen atoms entering the shock, and a broad contribution from excitation of “broad neutrals,” hydrogen atoms that have been produced by charge transfer reactions with protons in the downstream gas (Chevalier & Raymond, 1978; Chevalier et al., 1980). Such shocks are seen in many SNRs, including parts of the Cygnus Loop (Ghavamian et al. 2001, hereafter G01), Tycho and Kepler’s remnants (Kirshner et al. 1987, KCW87; Smith et al. 1991, S91; G01; Fesen et al. 1989, F89), RCW 86 (G01; Ghavamian et al. 2007b, G07b), SN 1006 (KCW87; S91; Ghavamian et al. 2002, G02), and several remnants located in the LMC (Tuohy et al. 1982, T82; S91; Ghavamian et al. 2003, G03; Ghavamian et al. 2007a, G07a).

Profiles of broad emission lines from edge-wise observations of shocks in SNRs can be used to infer shock velocities. These calculations can be compared to analyses of shock front proper motion to compute distances to the objects (Chevalier et al., 1980; Kirshner et al., 1987). Observations which resolve the spatial profile of the combined H emission can constrain the neutral fraction and density of the upstream gas (Raymond et al., 2007).

In addition to diagnosing parameters of SNRs, Balmer-dominated shocks provide an important probe into the plasma physics of collisionless, non-relativistic shocks. Models of the broad component width and integrated broad-to-narrow intensity ratio can be used to derive the ratio of electron to proton temperature as a function of shock velocity. Such a relation can provide insight into the physical mechanisms at work in the collisionless plasma.

This is the third in a series of papers investigating the hydrodynamics, kinetics, and line emission from Balmer-dominated shocks. In Heng & McCray (2007, Paper 1), we calculated velocity distribution functions for the broad neutrals and computed the ratio of broad-to-narrow line emission as a function of shock velocity. In Heng et al. (2007, Paper 2) we calculated the density structure of the shock transition zone, where hydrogen atoms passing through the shock front undergo charge transfer reactions, emit radiation, and become fully ionized.

The approximations employed in Paper 2 limit the validity of the results to shocks entering the upstream gas with velocities  km s, due to our treatment of the broad neutrals as a single fluid with the same bulk velocity as the ions (the “restricted three-component model” of Paper 2, §3). In fact, the broad neutrals in the shock transition zone are not a fluid, but have distinct anisotropic distribution functions depending on the number of charge transfer reactions they engage in (Paper 1). In the present paper, we treat charged species in the shocked gas as fluids, and describe the hydrogen atoms with appropriate kinetic distribution functions. In addition, we provide an improved calculation of the broad line velocity profile for km s. Using this methodology, we calculate the structure of the shock transition zone more accurately than in Paper 2, and characterize hydrogen line emission from shocks with  km s.

We find that we can self-consistently determine shock parameters for most Balmer-dominated SNRs. Our results yield lower values of the inferred shock velocity and proton-electron temperature equilibration than those derived using previous models. We compute the dependence of the electron temperature on shock velocity and compare it to the results of Ghavamian et al. (2007b). We note that our calculations are unable to fit the observations for several SNRs; in these cases, our basic model must be augmented by new physics to account for the data.

In §2 of this paper, we describe our physical model. In §3, we display the equations employed to calculate the structure of the shock transition zone and describe our numerical solution. In §4, we compute the spatial emission profile and hydrogen line spectra from the shock transition zone. In §5, we analyze our results and describe how observations of Balmer-dominated SNRs should be interpreted in light of our new calculations. Finally, in §6, we discuss implications of our results for collisionless electron heating, explore limitations of our model, and identify areas for future research.

2 Physical Model

We consider a shock with velocity km s traveling through the ISM, which consists of cold, partially neutral hydrogen and neutral helium. The pre-shock fraction of helium, relative to hydrogen, is denoted . If the total upstream density of protons and hydrogen atoms is , then the fraction of pre-shock protons is defined to be .

In the frame of the shock, the cold upstream neutrals and charged particles flow downstream with uniform velocity . At the shock, we assume that the protons are heated in a thermal distribution to K for cm s. At the high post-shock temperatures in Balmer-dominated remnants, emission from hydrogen excitation is much stronger than forbidden transitions in metals seen in radiative shocks (e.g., Shull & McKee, 1979). The length scale for thermalization of protons is set by the proton cyclotron gyroradius, cm at G, which is much less than the length scale for the shock transition zone,  cm (set by the mean-free paths for excitation, charge transfer and ionization reactions, and the total upstream density).

Without any energy transfer between shocked protons and electrons, the electron temperature is expected to be a factor of less than that of the protons. Since the collisional timescale for temperature equilibration between electrons and protons is in some cases longer than the age of the remnant (e.g., Spitzer, 1962), we assume that the electrons remain at a fixed temperature with respect to the protons in the downstream material (note, however, that intra-species and inter-ion collisional timescales can be much shorter).111Some collisional heating of electrons by protons will occur. We estimate that such an effect can lower the derived value of by However, several authors have argued that electrostatic instabilities at the shock front can increase the electron temperature, with estimates ranging from to (see Cargill & Papadopoulos, 1988; Ghavamian et al., 2007b, and the references therein). Since the physics of non-relativistic, collisionless shocks is still poorly understood, we parametrize the ratio of electron to proton temperatures using the definition


One of the goals of our work is to use observations of Balmer-dominated remnants to constrain the value of in the shock transition zone as a function of shock velocity.

In the uniform velocity upstream gas, we assume that negligible interactions take place between neutral and charged species. The cold neutral hydrogen atoms passing through the shock are not affected by the discontinuity. In contrast, the bulk velocities of the downstream thermal protons and electrons become approximately , at which point ionization and charge transfer reactions between cold atoms and hot protons occur (see §5.1). However, if a significant amount of the energy dissipated in the shock is used to produce cosmic rays, the ion compression ratio increases, and their bulk velocity decreases. Furthermore, the presence of a cosmic ray precursor can broaden the distribution of upstream neutrals. We discuss the implications of these physical effects in §6.

Downstream charge transfer reactions between cold neutrals and hot protons produce a new population of hot atoms, referred to as broad neutrals. Initially, such an interaction will produce a cold proton, which is re-energized through gyro motion around the magnetic field and intra-species Coloumb collisions with the hot proton population. We assume that such protons rapidly equilibrate back into the thermal pool. However, recent calculations by Raymond et al. (2008) indicate that interactions within the shock front may produce a distinct “pickup” ion population analogous to that in the solar wind.

H emission is produced by excitation reactions and charge transfers to excited states in the shock transition zone. The narrow and broad components of the line are emitted by the cold hydrogen and broad neutrals, respectively. The spatial and velocity distributions of the broad neutrals are critical for calculating the correct structure and emission from the transition zone. In Paper 2, we assumed that the broad neutrals could be characterized as a single fluid, with velocity and temperature equal to that of the ions. This treatment is flawed because, unlike the ionic species, the broad neutrals have negligible interactions among themselves. Therefore, each time a neutral engages in charge exchange, it becomes part of a new distribution with distinct velocity and temperature (see Paper 1) and is not subsumed into the original population. This fact has not, until now, been properly taken into account in the literature.

The detailed kinetic properties of the particle distributions and their reaction rates were calculated in Paper 1. We use these methods, with updated values for the atomic cross sections (see Appendix A), in our fluid calculation of the shock transition zone structure. In this picture, the broad neutrals form an infinite set of atomic populations with distinct reaction rates. In practice, as noted in Paper 1, the velocity distribution function of atoms experiencing many charge transfers rapidly converges to that of the protons, usually after only three or more interactions. Thus, only three distinct broad components interact with the other species. The bulk velocities and temperatures of the broad components are written:


where and are the velocity and temperature of the broad neutral population after charge transfers. The coefficients and are calculated using the methods described in Paper 1 and, in principle, are also functions of and . However, we show in §5.1 that the proton velocity and temperature do not vary significantly enough in the shock transition zone to affect the values of these coefficients. For , and are taken to be unity.

Ionization of helium in the shock transition zone produces singly-ionized He and alpha particles. We assume that the helium ions and protons are coupled through rapid inter-ion Coloumb collisions. In this case, all of the ions share the same bulk velocity () and temperature (). Figure 1 schematically displays the density variation of the different particle species in the shock transition zone.

3 Spatial Structure of the Transition Zone

3.1 Basic Equations

We employ a plane-parallel coordinate system in the frame of the shock, in which the shock front is at . The density structure of the transition zone is determined by conservation of mass flux:


where the rate coefficients for the atomic interactions (in units of cm s) are defined in Appendix A. Charge conservation requires the electron density to obey the relation . The subscripts , , and denote ionization of hydrogen, helium, and singly-ionized helium, while the index runs over the charged particles participating in the reactions, (electrons, protons, singly-ionized helium and alpha particles). The subscript indicates charge transfer of cold hydrogen atoms with protons. As discussed in §2, the subscript denotes the neutral hydrogen population which has engaged in charge transfers, where the index runs from to . While each broad population will have a separate rate coefficient for every atomic interaction, the values are not sensitive to the details of the broad distribution functions. Therefore, we use one rate coefficient to describe the interactions of all the broad populations: the subscripts and denote ionization and charge transfer, respectively. We list the references for our interaction cross sections and rate coefficients in Appendix A.

The system of differential equations is completed by expressions for conservation of momentum and energy flux:


where is the pressure and is the energy density. The ratios and are set by the kinetic calculation of the broad neutral distribution functions through eq. (3), which assumes . Though the ion temperature varies by as much as 10% in the shock transition zone, this produces a negligible effect on the rate coefficients. Thus, the value of is proportional to the local ion temperature. Furthermore, our steady-state fluid calculation neglects kinetic evolution of the broad neutral distributions, which we expect to be a modest effect.

We define a natural length scale for the transition zone, according to the formula:


where is a typical value for the interaction rate coefficents, taken to be cm s. We also define a set of dimensionless variables with , , , and . As discussed above, broad neutrals which have engaged in three or more charge transfers can be described with the thermal proton distribution function. We therefore define . Summing over eqs. (7) with and using these definitions, the formulas for the shock transition zone structure can be written:


where .

3.2 Solution Method

Eqs. (15)–(24) are a set of ten coupled, non-linear ordinary differential equations, which are solved with the constraint of charge conservation. We compute eqs. (15)–(22) directly in terms of the number density flux, denoted . The ion velocity and temperature, and , determine the broad velocity and temperature through the relations (2)-(3). Once and are known, the dimensionless number density can be substituted into the right hand side of the conservation equations.

We solve eqs. (23) and (24) subject to the boundary conditions , , , and . The upstream neutral temperature is taken to be zero, while the protons and electrons are in thermal distributions with temperatures . Integrating these equations, substituting the initial conditions, and using the definition for yields:


Note that the sums over charged species run from protons to alpha particles and do not include electrons, which appear separately in the equations (we neglect terms proportional to ). The system (3.2)–(3.2) can be solved simultaneously to yield a quadratic equation in with coefficients that are functions of . The solution to the quadratic has two positive roots. Only one of the roots is less than unity; it represents the physical solution for the bulk ion velocity, which must be less than the shock velocity. We use a standard Runge-Kutta method to solve eqs. (15)–(24), at each integration step using (3.2)–(3.2) to compute the ion velocity which determines the particle densities.

Figure 2 shows the density structure of the shock transition zone for km s, , and . Using eq. (14), the spatial coordinate has been converted from to physical distance behind the shock front. We use external density cm and helium fraction for all calculations in this paper. The left and right panels of Figure 2 show the dimensionless densities for the neutral and charged species, respectively. The electron and proton densities have been scaled by a factor of . In the left panel, the solid curve shows a monotonic decrease in the density of the cold hydrogen atoms, which are removed by both charge transfer and ionization reactions. Charge transfer produces three populations of broad neutrals. At low velocities, charge transfer dominates over ionization, and many broad neutrals with are produced. The dash-dotted curve shows the density of neutral helium, which is ionized farther downstream than hydrogen due to its smaller ionization rate coefficients. In the right panel, the solid curve shows singly-ionized helium, which is produced downstream the from neutral atoms, and is then ionized to yield a monotonically increasing population of alpha particles, shown by the dotted curve. The proton and electron densities are depicted by the short-dashed and long-dashed curves, which saturate when all of the neutral species have been depleted. The final electron density is slightly higher than that of the protons due to the presence of alpha particles.

4 Line Emission

Once the density structure of the transition zone is determined, we can compute the hydrogen line emission, including the spatial distribution and line profiles for the broad and narrow components. We neglect collisional dexcitation and assume that every atom is excited from the ground state. H photons are produced by transitions from atomic levels and to , as well as from to . In the latter case, the atomic physics is complicated by a possible transition from directly to , which results in a Lyman photon. If the medium is optically thin to Ly photons (Case A conditions), this possibility can be taken into account by proper weighting of the distinct angular momentum states in constructing excitation and change transfer cross sections for H emission:


where the factor is the fraction of transitions from to . If the medium is optically thick to Ly emission (Case B conditions), re-absorption by ground state hydrogen effectively traps Ly photons until they are re-emitted as H photons. In this case, all of the transitions eventually result in H emission, and we set .

Case A and B conditions represent the two extremes of media that are optically thin and thick to Ly scattering. For the stationary atoms (i.e., cold hydrogen) which produce the narrow line, the optical depth to scattering of Ly photons is, at line center, , with cm, cm, and cm (see, e.g., Rybicki & Lightman, 1979; Cox, 2001). The column of upstream neutrals will also line scatter Ly emission, increasing the effective value of . Thus, in conditions appropriate for many Balmer-dominated SNRs, partial scattering of Lyman photons will occur, producing results intermediate between Cases A and B for the narrow line (Ghavamian et al., 2001).

The conversion efficiency of narrow component Ly to H was first computed by Chevalier et al. (1980). Subsequent Monte Carlo calculations were performed by Laming et al. (1996) and Ghavamian et al. (2001). Briefly, models of the neutral hydrogen density and narrow component excitation rate are computed as a function of distance behind the shock. These are used to calculate the profile of excitations to the 3p level. Photons are emitted at frequencies distributed according to the pre-shock temperature in random directions. They are followed as they are absorbed and re-emitted as Ly or H at a different locations until they escape. The conversion fraction depends on shock speed, electron-ion temperature ratio, and pre-shock ionization fraction, but not on total upstream density. We assume that the variation with is modest, and that the temperature ratio decreases from for km s to for km s, and use the calculations of Laming et al. (1996), in which of 3p excitations result in Ly photons. We fit the resulting conversion fraction as a function of shock velocity, yielding:


where is the shock velocity in units of cm s.

4.1 Spatial Emissivity Profile

To determine the emission, we employ excitation rate coefficients calculated using the methods of Paper 1 (see also Appendix A). An emission line photon is produced when a cold hydrogen atom or broad neutral is excited by a charged particle, or undergoes charge transfer to an excited state. In the former case, the rates must be weighted by the probability of repeated excitation. We use the following equations to calculate the spatial emissivity profiles for the narrow () and broad () components:


where the rate coefficients are labeled by the transition , atomic interaction, and particle type. The symbols and denote excitation and charge transfer to an excited state for cold hydrogen, while and denote these interactions for broad neutrals. The probabilities and are calculated using the total reaction rates per atom or broad neutral for ionization, excitation, and charge transfer. These are calculated using the weighted sum of the rate coefficients, for example, .

Figure 3 shows the spatial emissivity profiles for the narrow and broad components as a function of distance behind the shock front. The thin vertical lines indicate the centroids of the emission components, calculated according to the formula:


For the low value of shown in this figure, significant ionization must occur before there are enough charged particles to excite the narrow emission and engage in charge transfers to produce broad neutrals. Therefore, the intensity peaks for both components are shifted downstream from the shock front, and the centroid of the broad line emission is shifted farther downstream than that of the narrow line.

4.2 Line Profiles

The full width half maximum (FWHM) of the broad line can be related to the velocity and temperature equilibration of the shock. The line profile is a convolution of the broad neutral and exciting species distribution functions with the cross sections for excitation and charge transfer to an excited state, projected along the line of sight to the observer. In a cylindrical coordinate system where the axis is along the shock velocity direction, it is straightforward to project the line profile for observers oriented both edge-wise (along the axis) and face-on (along the axis) with respect to the shock front. Due to limb brightening, most observations will be selected for shocks viewed edge-wise or nearly so.

We calculate hydrogen line profiles , and according to the formulas:


where , and the cross sections denote excitation and charge transfer reactions for H emission to appropriately weighted angular momentum states, using eqs. (27) and (28). It should be noted that in eqs. (32)–(35), the quantities denoted by are kinetic distribution functions, and should not be confused with the pre-shock ionization fraction . The edge-wise profile can be calculated as a function of position behind the shock. In practice, however, the observed lines are not spatially resolved. We therefore calculate the profiles as a function of and , spatially averaged over the shock transition zone. Figure 4 shows examples of symmetric, edge-wise broad neutral velocity distributions, which are inputs to eq. (33). Results are depicted in a reference frame where the average ion velocity is zero, for , with km s (left panel), and km s (right panel). For all cases, the pre-shock ionization fraction is set to .

For km s, excitation by electrons dominates over that by protons (including charge transfer into excited states). Even at moderate values of the electron temperature, , the thermal width of the electron distribution function is much larger than that of the broad neutral distribution, and the electrons typically have much higher velocities. Thus, for the majority of the range of integration, , and the electron contribution to eqs. (32)–(33) is approximately equal to the projected velocity distribution of the broad neutrals, multiplied by the rate coefficient for excitation by electrons. This approximation has been used in previous studies of the line profile (e.g., Chevalier et al., 1980; Ghavamian et al., 2001; Heng & McCray, 2007). Nevertheless, for shock velocities km s, excitation rates by electrons and protons are comparable, with the proton contribution increasing for faster shocks. The cross sections for excitation by protons increase with the relative speed of the colliding particles. Since the broad neutrals and ions have comparable speeds, high speed neutrals are more likely to produce H photons. The effect on the line profile is somewhat mitigated by the integration over velocity space, but the observed line width is larger than the velocity width of the broad neutral distribution. This effect is relatively small at low velocity () but is significant for high velocity shocks (see §6).

Figure 5 displays the FWHM of the broad, edge-wise H line profile as a function of , at several values of . The pre-shock ionization is set at . The FWHM increases monotonically with due to the increased temperature of the post-shock ion distributions (and hence broad neutral distributions) behind faster moving shocks. As is increased, energy is transferred from the protons to heat the electrons, leading to lower proton temperatures and a smaller FWHM for the broad component.

For fast shocks, in which proton excitation contributes substantially to eqs. (32)–(33), we expect different transitions (e.g., H and Ly) to have different FWHM relations, since they employ distinct reaction cross sections in the calculation of . This is a unique prediction made by our calculation, and may have consequences for studies of Ly emission from shocks with km s, in which the FWHM of the Ly line increases from 10-60% over that of the H line.

5 Dependence on Shock Parameters

Here we describe how the structure and emission from the transition zone of Balmer-dominated SNRs depend on shock velocities (ranging from km s), pre-shock ionization fractions (ranging from ), and temperature equilibration ratios (ranging from ), as well as Case A and Case B conditions for the narrow line.

5.1 Spatial Structure of the Shock Transition Zone

The basic features of the shock transition zone density structure were shown in Figure 2. The details of broad neutral production are highly sensitive to , , and . We can see the effect of increasing shock velocity by comparing Figure 6 to Figure 2. At high , the ionization rate is greater than the charge transfer rate, so that relatively few broad neutrals are produced. The densities of the broad neutral populations decrease rapidly with each subsequent charge transfer reaction. The helium ionization rates decrease, causing the and alpha particle production to peak further downstream.

In Paper 2, we explored variations in ion velocity throughout the shock transition zone and found that, for km/s, with negligible deviation. We confirm this conclusion with our multi-component models. Figure 7 shows the dimensionless ion velocity and temperature as a function of distance from the shock front. Results are shown for km s and km s, with , . In the left panel, we display the percent deviation of the ion velocity from for the two shock velocities. At  km s, broad neutrals are produced with smaller momentum density than the original ion population, leading to an increase in ion velocity by conservation of momentum (solid curve). For km s, the opposite occurs, and the ion velocity decreases as the broads are produced (dotted curve). The deviation of ion velocity from is less than and has a negligible effect on calculations of the reaction rate coefficients . In the right panel, we plot the ion temperature profile. In this case, the broad neutrals are produced with slightly lower temperature than the ions, leading to slight heating of the ions. At large distances downstream, the presence of alpha particles increases the mean atomic mass in the gas to 1.27 (assuming a helium abundance), leading to a final ion temperature that is slightly higher than . The maximum deviation of from its expected value is of order , which in practice has a negligible effect on the values of the reaction rate coefficients.

5.2 H Emissivity

In Figure 3, we show a typical emissivity profile at relatively low shock velocity and ionization fraction. In the left panel of Figure 8, we show the effect of increasing the initial ionization fraction to . We see that the neutral density decreases rapidly and the narrow line emission peaks at the shock front (left panel). Many charge transfer reactions occur close to the shock, pushing the centroid of the broad emission farther upstream compared to the case. Comparing the right panel of Figure 8 to Figure 3, we can see the effect of increasing the shock velocity. For the high velocity case, both ionization and charge transfer rates are decreased, shifting the centroids of both line components farther downstream in the transition zone.

As discussed by Raymond et al. (2007) and Paper 2, the spatial shift between the broad and narrow line centroids potentially provides a constraint on the pre-shock ionization fraction, , and external density, . In Figure 9, we plot as a function of shock velocity , with cm. Note that the spatial shift scales as . In the left panel, we show how depends on , holding . As the velocity increases, the charge transfer rate decreases relative to the ionization rate, delaying the production of broad neutrals. Consequently, the centroid of shifts downstream and the value of increases. For small , few protons initially exist to engage in charge exchange. Consequently, shifts downstream from the shock front. In the right panel, we show how the results depend on the temperature equilibration parameter, , for fixed . At shock velocities km s, increasing reduces the proton ionization rate for cold hydrogen. This effect causes the peak of broad production to shift downstream from shock front. At high velocities, increasing tends to increase the charge transfer rate for broad neutrals over that for cold hydrogen. This effect causes the peak of broad production to shift closer to the shock front, and decreases the value of . For more discussion on the observational significance of the spatial emissivity profile and the shift, see §6.2.

Another important observational diagnostic for Balmer-dominated SNRs is the ratio of integrated broad-to-narrow line strengths, defined as:


This ratio has a strong dependence on both and . Figure 10 shows versus shock velocity at fixed , for several values of temperature equilibration , using Case A (left panel) and Case B (right panel) conditions for the narrow line. For a given , the variation of the intensity ratio with velocity is the result of competition between charge transfer and ionization, which contribute equally at km s. At very low , the ionization rate begins to decrease precipitously while the charge transfer rate stays roughly constant, leading to the spike in . At a fixed shock velocity, increasing causes a decrease in the proton ionization rate for broad neutrals compared to cold hydrogen, leading to an increase in the intensity ratio. In Case B conditions (right), additional narrow emission due to absorption of trapped Ly photons decreases the values of relative to Case A. As noted in §4, for most Balmer-dominated SNRs, partial line scattering of Ly photons yields emission intermediate between Cases A and B.

According to the picture in Paper 1, the dependence of on is expected to be weak, since the number of broad neutrals produced in each population per cold hydrogen atom is fixed. This intuition is confirmed by our multi-component calculation. When is increased, we expect more charge transfers in the shock transition zone. However, this effect is balanced by higher ionization rates for both broad and cold neutral populations. The net result is a negligible change in the ratio. The presence of neutral helium introduces a weak dependence of on . In Figure 11, we show as a function of for fixed in Case A (left panel) and Case B (right panel) conditions. As decreases, fewer protons engage in charge transfer and ionization reactions close to the shock front, shifting the peak of broad neutral production downstream. In this case, the broad neutrals persist far enough downstream to interact with the charged helium species produced there, altering the ratio . This effect is of order for and remains small for larger helium fractions. However, it should be noted that, when the effects of partial Ly scattering are included, variations in optical depth with can introduce a dependence of on pre-shock ionization fraction (Ghavamian et al., 2001). While we have roughly treated Ly scattering using the prescription of §4, a full calculation of line scattering in the shock transition zone is needed to properly take into account these effects and incorporate the dependence of the integrated line ratio on the pre-shock ionization fraction. This is a source of systematic error in our calculation.

Neglecting the dependence, the broad-to-narrow intensity ratio can be written as a function


where and represent the measured and theoretical values of the intensity ratio, respectively.

5.3 Broad Line Profiles and Observational Interpretations

Eqs. (32) and (33) allow us to model the FWHM of the broad line as a function of and , as shown in Fig. 5. In the extreme scenarios of Case A and Case B conditions, the FWHM has a weak () dependence on the pre-shock ionization after taking a spatial average over the shock transition zone. Figure 12 shows examples of broad neutral velocity distributions [which are inputs to eq. (32)] for the case of face-on orientation of the shock as a function of line-of-sight velocity , in a frame of reference where the ion velocity is zero. Results are shown for fixed at several values of at shock velocities km s (left panel) and km s (right panel). At low shock velocities, charge transfer is extremely efficient, yielding broad neutral distribution functions nearly identical to that of thermal protons, with the broad neutral moving slightly slower than the ions. At high shock velocities, the broad neutral distributions are skewed and offset from the proton distribution, leading to the asymmetric profiles depicted in the right panel, with the broad neutrals moving considerably faster than the ions.

We write the FWHM relation as a function of and :


where and represent the measured and theoretical values of the line FWHM, respectively. For SNRs with measured and , combining eqs. (37) and (38) self-consistently constrains and . This is accomplished by inverting to yield , and calculating the root of the expression


This procedure produces a pair of values for which the theoretical calculations and equal the measured and .

We emphasize that a self-consistent calculation is critical for accurately determining the values of and . In the previous literature, two bracketing values for the temperature equilibration are sometimes chosen (e.g., ), and the FWHM relation is employed to give a range of possible values for the shock velocity (e.g., Heng & McCray, 2007). However, this procedure has a major flaw. In practice, the self-consistent quantity will have a minimum value over the range . If for a particular observed shock is less than this minimum value, then no pair will yield the observed values and . In such a case, the model breaks down and quoting a range of possible shock velocities for two bracketing values of is inappropriate. Additional physics must be invoked to account for the observations. Moreover, when no measurement of exists, the bracketing procedure may or may not yield an accurate estimate of the range of shock velocities.

We use the methodology described above to self-consistently extract shock parameters from observations of Balmer-dominated SNRs. We summarize our results in Table 1, which lists, from left to right, the object, reference, H and , calculated values for and from H, Ly , and calculated value for from Ly. If our models do not yield a fit to the observations, we do not list values for the shock velocity and temperature equilibration ratio. Our calculations show a characteristic range of between and for km s. Therefore, in the case of Ly observations where no measurement of exists, we report derived shock velocities for the range ; further observations are required to confirm the accuracy of these estimates. For SNR 0519-69.0, observations exist in H and Ly. To calculate shock velocities from the Ly observations, we use the derived value of from the H diagnostics.

Our models successfully fit 14 measurements from seven Balmer-dominated SNRs, within the observational uncertainties. In the cases of the Cygnus Loop (Ghavamian et al., 2001) and one measurement from 0519-69.0 by Tuohy et al. (1982), the observed ratio is too low to be accounted for by our calculations. In addition, we are unable to fit the majority of measurements from the object DEM L71/0505-67.9, which we have omitted from the table. Below, we discuss a possible explanation for these discrepancies in our model.

We note that our inferred shock velocities and temperature equilibration ratios are systematically lower than in previous studies, and show an increased sensitivity to compared to Papers 1 and 2. The ability to sensitively probe both shock velocity and temperature equilibration is of interest for studies of collisionless electron heating in shocks, as described below. Typically, inferred values for are smaller than those quoted in Paper 1, primarily because of the contribution of broad neutral velocities to relative speeds in fast neutral-ion interactions. The shock speeds in Tycho’s SNR and SN1006 are smaller by about 15% and 27%, respectively, than those reported previously. One implication is that the distances to these SNRs derived from the shock speeds and proper motions are correspondingly smaller. The inferred distance to SN1006 is reduced from 2.18 kpc (Winkler et al., 2003), to 1.6 kpc. The corresponding brightness of the SN is squarely in the middle of the Type I SN distribution at 2.18 kpc, but is 0.7 magnitudes fainter at 1.6 kpc. The smaller distance is more comfortable in comparison with that of the S-M star which lies between 1.05 kpc and 2.1 kpc and whose spectrum shows absorption by SN1006 ejecta (Burleigh et al., 2000). On the other hand, the ejecta are observed to expand at 7026 km s (Hamilton et al., 2007), and the requirement that this material lies within the remnant places a conservative lower limit to the distance of 1.6 kpc, which is just consistent with that derived here.

6 Discussion

6.1 Collisionless Shock Heating

One of the motivations for studying Balmer-dominated SNRs is to probe the physics of collisionless shocks. While various mechanisms for transferring energy from proton to electron populations have been proposed in the literature, there is no consensus on how to predict the value of for a given set of shock parameters in astrophysical contexts. In a recent paper, Ghavamian et al. (2007b) attempted to address this problem by deriving a relation for temperature equilibration ratio versus shock velocity from observations of shocks in SNRs. They reported that their results could be fit by a curve , with for km s. Given that the proton temperature scales roughly as , this conclusion implies that the electrons are heated to a constant temperature, independently of the shock velocity. While Ghavamian et al. (2007b) discussed a possible mechanism for this dependence, it has not yet been established that theoretical models can produce such an effect.

Using our new model for emission from Balmer-dominated shocks, we have calculated an updated relation, which is displayed in Figure 13. The solid curve depicts the proposed dependence, with for km s. For all the points shown in the plot, we have self-consistently fit measurements of both and . We fix the pre-shock ionization at , and exclude the majority of measurements from DEM L71/0505-67.9, as well as one measurement each from Cygnus and 0519-69.0 which we are unable to account for with our models. We find from our calculations that the data are not well fit by a power-law relation ; if we set , the fit to the curve yields a reduced chi-square of .222Note that such a large value of should be interpreted in the context of uncertainties associated with the atomic cross sections used. For shock velocities km s, our results are fairly close to those of previous models. At higher velocities km s, the deviations are more significant. The minimum value of the temperature equilibration ratio is at velocities  km s. This value is greater by several orders of magnitude than the theoretical minimum , but is smaller than the value predicted by some collisionless heating models (e.g., Cargill & Papadopoulos, 1988).

6.2 Limitations to our Model and Future Work

As demonstrated by SNR 0505-67.9, the Cygnus Loop, and 0519-69.0, an additional physical mechanism is needed to account for the observed H emission seen in some SNRs. One possibility is that a significant amount of the dissipated energy in the shock is transferred to cosmic rays, producing a precursor which can heat and accelerate the upstream gas, altering the shock jump conditions (Smith et al., 1991; Hester et al., 1994; Sollerman et al., 2003). In addition, the precursor can “push” on the upstream protons, leading to a velocity differential between neutrals and charged species (e.g., Berezhko & Ellison, 1999). Furthermore, all previous models have assumed that the kinetic distribution functions for protons and electrons in the shock transition zone are Maxwellian. However, recent calculations have shown that the proton distribution can significantly depart from Maxwellian behavior, forming a distinct “pickup” ion population (Raymond et al., 2008). This effect will change the structure of the transition zone, broad neutral distributions, and kinetic reaction rates.

The inclusion of a cosmic ray precursor will affect the predicted broad-to-narrow line strength in several ways. Broadening of the cold neutral distribution function effectively reduces the optical depth to Ly scattering, which will decrease the conversion efficiency of narrow Ly to H and increase . In contrast, additional excitation of cold neutral H atoms in the precursor will increase narrow line emission and decrease . While it is not obvious which effect is dominant, evidence for such precursor effects exists in the anomalously large widths of H lines in most Balmer-dominated SNRs (Sollerman et al., 2003). In addition, H emission from a spatially resolved precursor in Tycho’s SNR has recently been reported by Lee et al. (2007).

The spatial structure of the shock transition zone provides a way to infer the external density and pre-shock ionization fraction . As noted in §5.2 the spatial emissivity profile and centroid shift between the narrow and broad components have a strong dependence on and . Raymond et al. (2007) were able to spatially resolve the total H emission from SN1006 using the ACS camera on the Hubble Space Telescope. Raymond et al. (2007) adopted the values of and from those inferred by Ghavamian et al. (2002) using observations of the FWHM and , and then attempted to fit for the spatial structure by varying and . Future instruments similar to ACS with high spatial resolution (e.g., WFC3 onboard Hubble with the narrow H component isolated by a narrow band filter) will allow for more accurate measurements of the transition zone structure and provide additional constraints on shock parameters.

The models presented in our series of papers generically describe the physics of non-radiative shocks interacting with cold, pre-shock gas. In addition to using them to understand nearby SNRs, we can also calculate hydrogen line emission from SNRs in young, distant galaxies (with redshifts ), a subject first explored by Heng & Sunyaev (2008). We improve on their estimates of the luminosity ratios of Ly and Ly to H, denoted and , respectively (see Figure 14; c.f. Figure 1 of Heng & Sunyaev 2008). Heng & Sunyaev (2008) underestimated the broad Ly emission at high shock velocities due to their neglect of the multi-component shock transition zone. Nevertheless, their conclusions remain intact: the sensitivity of and to by a factor over the velocity range km s is a direct and unique way to measure the temperature ratio. A valuable extension to their work will incorporate models of Ly and Ly scattering, taking into account geometric effects.

The biggest systematic uncertainties in our current model are: i) a proper treatment of Lyman line scattering in the shock transition zone; ii) the effect of a cosmic ray precursor on ; and iii) inclusion of a non-thermal population of protons in calculating the reaction kinetics. The incorporation of these physical effects is important for future models. Our new results yield substantial differences in derived shock speed (and hence inferred distances) for measurements of Balmer-dominated SNRs compared to previous models. While our results are not consistent with a power-law relation between shock velocity and electron temperature, future work is needed to resolve the remaining approximations in our model and make a more definitive statement about shock heating of electrons.

We thank Carles Badenes, Anatoly Spitkovsky and Eliot Quataert for useful conversations, and Parviz Ghavamian for helping clarify several points in the observational literature. K.H. thanks the Institute for Advanced Study for their generous support.

Appendix A Shock Kinetics

We summarize important details of our calculations of rate coefficients, broad neutral velocity distributions and broad line profiles. For a more detailed description of our methods, see Paper 1. In this work, we treat the atomic sub-levels separately instead of considering a single, summed level as in Paper 1. We treat charge transfers, excitation and ionization events between electrons, protons, alpha particles and hydrogen atoms using the cross sections of Barnett et al. (1990), Belkić et al. (1992), Janev & Smith (1993), Balança et al. (1998), and Harel et al. (1998). Fitting functions for some of these cross sections are provided in Heng & Sunyaev (2008).

We also consider the ionization of helium atoms and singly-ionized helium by electrons and protons, using the cross sections of Peart et al. (1969), Angel et al. (1978), Rudd et al. (1983), Shah & Gilbody (1985), Rinn et al. (1986) and Shah et al. (1988). Following Barnett et al. (1990), we fit the cross sections to the function


where the components of are the fitting parameters, and the quantities are the Chebyshev polynomials:


We define the fitting variable as


where is the relative energy between the interacting particles and and are the minimum and maximum energies for which data are available. We assume a fiducial error of 10 for the data. The cross sections and corresponding fits are displayed in Figure 15 and the fitting parameters are presented in Table 2.

We neglect charge transfer reactions of neutral and singly-ionized helium with protons due to scarcity of cross sections for these processes in our velocity range of interest. While such interactions may affect helium line emission, we do not expect them to strongly couple to the Balmer radiation. In addition, for ease of computation, we approximate excitation and ionization of hydrogen by He using relevant proton rate coefficients.333Our assumption is based on the Weizsacker-Williams approximation in which scattering dynamics are dominated by the charge of the impacting particle as opposed to its mass (e.g., Jackson, 1998). This approximation is valid when the relative velocity of the collision is greater than that of an electron orbiting the hydrogen atom. At low velocities km s, this assumption breaks down. However, cross sections for this process do exist in the literature and should be used in future calculations (Barnett et al., 1990). Nevertheless, we expect corrections to these cross sections to have a small quantitative effect on our results.

We make several approximations to speed up our computations. For temperature ratios , the velocity width of the electron distribution is generally broader than that of the broad neutrals (characteristic electron velocities are greater than proton velocities by a factor of ), implying that rate coefficients for interactions involving electrons are insensitive to changes in the broad neutral velocity distribution. We therefore approximate electron rate coefficients to be the same for reactions involving cold atoms and broad neutrals. Calculations of the FWHM for line profiles should include excitations by electrons, protons, singly-ionized helium and alpha particles. However, alpha particles contribute due to the relatively smaller density , and can typically be omitted.

The reaction rate coefficients used in eqs. (4)–(11) are defined as:


where denotes the interaction (ionization or charge transfer), the atomic species (broad or cold neutrals), and the interacting charged particle (electrons, protons, singly-ionized helium, or alpha particles). In contrast to eqs. (33)–(35), the interaction cross sections used here represent the sum of the reactions to all levels for which atomic data is available [c.f., eq. (27)]. Using the definition of eq. (A6), the quantity gives the number of interactions , between the species and , per unit volume.

Appendix B Typographical Errors in Papers 1 & 2

We point out several minor typographical errors in Heng & McCray (2007) and Heng et al. (2007). In §5.1 of Heng & McCray (2007), should be . In §5 of Heng et al. (2007), should be , should be and should be . In §6.2, paragraph 2 of Heng et al. (2007), the sentence that begins “Before CKR80 and HM07, …” should be changed to “Before (CKR80 and HM07), …”.


  • Angel et al. (1978) Angel, G. C., Dunn, K. F., Sewell, E. C., & Gilbody, H. B. 1978, Journal of Physics B Atomic Molecular Physics, 11, L49
  • Balança et al. (1998) Balança, C., Lin, C. D., & Feautrier, N. 1998, Journal of Physics B Atomic Molecular Physics, 31, 2321
  • Barnett et al. (1990) Barnett, C. F., Hunter, H. T., Fitzpatrick, M. I., Alvarez, I., Cisneros, C., & Phaneuf, R. A. 1990, NASA STI/Recon Technical Report N, 91, 13238
  • Belkić et al. (1992) Belkić, D., Gayet, R., & Salin, A. 1992, Atomic Data and Nuclear Data Tables, 51, 59
  • Berezhko & Ellison (1999) Berezhko, E. G. & Ellison, D. C. 1999, ApJ, 526, 385
  • Burleigh et al. (2000) Burleigh, M. R., Heber, U., O’Donoghue, D., & Barstow, M. A. 2000, A&A, 356, 585
  • Cargill & Papadopoulos (1988) Cargill, P. J. & Papadopoulos, K. 1988, ApJ, 329, L29
  • Chevalier et al. (1980) Chevalier, R. A., Kirshner, R. P., & Raymond, J. C. 1980, ApJ, 235, 186
  • Chevalier & Raymond (1978) Chevalier, R. A. & Raymond, J. C. 1978, ApJ, 225, L27
  • Cox (2001) Cox, A. 2001, Allen’s Astrophysical Quantities (Springer Verlag)
  • Fesen et al. (1989) Fesen, R. A., Becker, R. H., Blair, W. P., & Long, K. S. 1989, ApJ, 338, L13
  • Ghavamian et al. (2007a) Ghavamian, P., Blair, W. P., Sankrit, R., Raymond, J. C., & Hughes, J. P. 2007a, ApJ, 664, 304
  • Ghavamian et al. (2007b) Ghavamian, P., Laming, J. M., & Rakowski, C. E. 2007b, ApJ, 654, L69
  • Ghavamian et al. (2003) Ghavamian, P., Rakowski, C. E., Hughes, J. P., & Williams, T. B. 2003, ApJ, 590, 833
  • Ghavamian et al. (2001) Ghavamian, P., Raymond, J., Smith, R. C., & Hartigan, P. 2001, ApJ, 547, 995
  • Ghavamian et al. (2002) Ghavamian, P., Winkler, P. F., Raymond, J. C., & Long, K. S. 2002, ApJ, 572, 888
  • Hamilton et al. (2007) Hamilton, A. J. S., Fesen, R. A., & Blair, W. P. 2007, MNRAS, 381, 771
  • Harel et al. (1998) Harel, C., Jouin, H., & Pons, B. 1998, Atomic Data and Nuclear Data Tables, 68, 279
  • Heng & McCray (2007) Heng, K. & McCray, R. 2007, ApJ, 654, 923
  • Heng & Sunyaev (2008) Heng, K. & Sunyaev, R. A. 2008, A & A, in press, 481, 117
  • Heng et al. (2007) Heng, K., van Adelsberg, M., McCray, R., & Raymond, J. C. 2007, ApJ, 668, 275
  • Hester et al. (1994) Hester, J. J., Raymond, J. C., & Blair, W. P. 1994, ApJ, 420, 721
  • Jackson (1998) Jackson, J. D. 1998, Classical Electrodynamics, 3rd Edition (New York: Wiley & Sons)
  • Janev & Smith (1993) Janev, R. & Smith, J. 1993, Cross Sections for Collision Processes of Hydrogen Atoms with Electrons, Protons, and Multiply Charged Ions (Vienna: International Atomic Energy Agency)
  • Kirshner et al. (1987) Kirshner, R., Winkler, P. F., & Chevalier, R. A. 1987, ApJ, 315, L135
  • Laming et al. (1996) Laming, J. M., Raymond, J. C., McLaughlin, B. M., & Blair, W. P. 1996, ApJ, 472, 267
  • Lee et al. (2007) Lee, J. J., Koo, B. C., Raymond, J. C., Ghavamian, P., Pyo, T. S., Tajitsu, A., & Hayashi, M. 2007, A&A, 569, L133
  • Peart et al. (1969) Peart, B., Walton, D. S., & Dolder, K. T. 1969, Journal of Physics B Atomic Molecular Physics, 2, 1347
  • Raymond et al. (2008) Raymond, J. C., Isenberg, P. A., & Laming, J. M. 2008, astro-ph/0804.3808
  • Raymond et al. (2007) Raymond, J. C., Korreck, K. E., Sedlacek, Q. C., Blair, W. P., Ghavamian, P., & Sankrit, R. 2007, ApJ, 659, 1257
  • Rinn et al. (1986) Rinn, K., Melchert, F., Rink, K., & Salzborn, E. 1986, Journal of Physics B Atomic Molecular Physics, 19, 3717
  • Rudd et al. (1983) Rudd, M. E., Goffe, T. V., Dubois, R. D., Toburen, L. H., & Ratcliffe, C. A. 1983, Phys. Rev. A, 28, 3244
  • Rybicki & Lightman (1979) Rybicki, G. B. & Lightman, A. P. 1979, Radiative Processes in Astrophysics (New York: Wiley-Interscience)
  • Shah et al. (1988) Shah, M. B., Elliott, D. S., McCallion, P., & Gilbody, H. B. 1988, Journal of Physics B Atomic Molecular Physics, 21, 2751
  • Shah & Gilbody (1985) Shah, M. B. & Gilbody, H. B. 1985, Journal of Physics B Atomic Molecular Physics, 18, 899
  • Shull & McKee (1979) Shull, J. M. & McKee, C. F. 1979, ApJ, 227, 131
  • Smith et al. (1991) Smith, R. C., Kirshner, R. P., Blair, W. P., & Winkler, P. F. 1991, ApJ, 375, 652
  • Sollerman et al. (2003) Sollerman, J., Ghavamian, P., Lundqvist, P., & Smith, R. C. 2003, A&A, 407, 249
  • Spitzer (1962) Spitzer, L. 1962, Physics of Fully Ionized Gases (Physics of Fully Ionized Gases, New York: Interscience (2nd edition), 1962)
  • Tuohy et al. (1982) Tuohy, I. R., Dopita, M. A., Mathewson, D. S., Long, K. S., & Helfand, D. J. 1982, ApJ, 261, 473
  • Winkler et al. (2003) Winkler, P. F., Gupta, G., & Long, K. S. 2003, ApJ, 585, 324
Object Reference H FWHM H H Ly FWHM Ly
(km s) (km s) (km s) (km s)
Cygnus G01
RCW 86 G07b