© ESO, 2013

1. Introduction

Optical spectra of a few young supernova remnant (SNR) shocks essentially consist of Balmer emission. Observations show that the Hα line produced in these so-called Balmer dominated shocks is made of two components (see Heng 2010, for a review): a narrow line, with a full width at half maximum (FWHM) of few tens of km s^-1 produced by excitation of cold atoms entering the shock region, and a broad line with a FWHM resembles the shock speed produced by atoms that undergo a charge exchange (CE) process before being excited (Chevalier & Raymond 1978; Chevalier et al. 1980). Measuring the shape of this complex Balmer line allows one to infer information on the conditions around the shock, which in turn are affected by the presence of accelerated particles. Balmer emission from young SNR shocks can be used to measure the cosmic ray (CR) acceleration efficiency at such shocks and even to gather information on a CR induced precursor upstream of the shock as predicted by nonlinear theories of diffusive shock acceleration.

Efficient CR acceleration at SNR shocks changes the Balmer line emission in (at least) two ways. For CR modified shocks, the temperature of the ions downstream of the shock is expected to be lower than in the absence of particle acceleration. This leads to a narrower width of the broad Balmer line. On the other hand, the pressure carried by accelerated particles also produces a CR-induced precursor in which ions can be heated by adiabatic heating (the plasma slows down) and interaction with magnetic turbulence, possibly induced by the same accelerated particles. This second effect leads to a broadening of the narrow component of the Balmer line.

An anomalous width of the narrow Balmer line was first reported as early as 1994 (Smith et al. 1994; Hester et al. 1994): a FWHM ranging from 30 to 50 km s^-1 was detected in several SNRs, implying a preshock temperature around 25 000–50 000 K, which is clearly in excess of the temperature for which hydrogen would stay neutral, in the presence of collisional ionization. This suggests that the line is broadened in some sort of preshock region small enough to make collisional ionization unimportant. This is exactly the effect expected from a CR-induced precursor, as already speculated by Smith et al. (1994). From the theoretical point of view, the first attempt at calculating the broadening of the narrow Balmer lines due to a CR-induced precursor was done by Wagner et al. (2009), who used a two-fluid model to treat ions and CRs but neglected the dynamical role of neutrals. A different model was proposed by Raymond et al. (2011) where momentum and energy transfer between ions and neutrals are included, but both the CR spectrum and the profile of the CR-precursor are assumed rather than calculated. Both works conclude that the observed width of 30–50 km s^-1 can be explained by postulating some level of CR acceleration at the shock. A more detailed calculation has been developed in Morlino et al. (2013) to take the mutual interplay between CRs, neutrals, ionized plasma and magnetic turbulence into account. The conclusion of this latter work is that, as far as the broadening of the narrow line is concerned, the main physical effect to take into account is the damping of magnetic turbulence in the CR precursor, while adiabatic compression alone is ineffective.

A smaller FWHM of the broad Balmer line, such as the one expected from the presence of CRs (van Adelsberg et al. 2008), has been inferred in two SNRs: RCW 86 inside our own Galaxy (Helder et al. 2009), and SNR 0509–67.5 in the LMC (Helder et al. 2010). In the first case, the authors combined proper motion measurements of nonthermal X-ray filaments observed by Chandra and optical measurements by the VLT of the width of the broad Hα component from the same region to show that the temperature behind the shock is lower than expected based on Rankine-Hugoniot relations. The conclusion was that a sizable fraction of the energy (≥50%) is being channeled into CRs. In fact this conclusion has somehow been retracted by the same authors in a recent paper (Helder et al. 2013) where they use new measurements of the proper motion of optical filaments, which imply a lower shock velocity than inferred from the X-ray proper motion. Using the new measurements, the CR acceleration efficiency that the authors estimate is lower and depends on the assumed distance of the remnant, which is compatible with zero for d ≲ 2.5 kpc. For the second remnant, SNR 0509–67.5, Helder et al. (2010) estimate that the CR pressure behind the SW region of the shock is >15% of the total post-shock pressure.

In both these papers, the CR acceleration efficiency was estimated through a relation between the FWHM of the broad line and the shock velocity found by van Adelsberg et al. (2008, hereafter vAHMR) using the formalism developed in Heng & McCray (2007). This relation was based on the assumption that after three CE reactions, the distributions of atoms and ions reach thermal equilibrium (namely the same temperature and bulk velocity). Here we use a kinetic approach to show that this condition is not fulfilled and this fact substantially changes the expected width of the broad Balmer line.

An accurate calculation of the Balmer emission requires a proper description of the distribution function of the neutral particles in the shock region. Neutrals may undergo a number of CE events before being ionized, but in general they do not thermalize with ions, so that the neutral distribution function can only be determined using the full Boltzmann equation. Moreover, when CRs are also present, the shock structure is determined by the mutual interaction between ions, neutrals, CRs and magnetic field, which are all coupled together through different processes, resulting in a highly nonlinear system. In a series of three papers (Blasi et al. 2012; Morlino et al. 2012, 2013, hereafter Papers I, II, and III), we developed a semi-analytical method for describing the diffusive shock acceleration process in the presence of neutral particles. Our method also estimates the Balmer emission produced by excited hydrogen atoms in the shock region.

In the present paper we use the previously developed formalism to assess the relation between broad Balmer line width, shock speed, and CR acceleration efficiency. Compared with our previous work, we improved our calculations by better approximating of the ion-neutral scattering rate. We show that, for shock speed ≳2500 km s^-1, the neutrals experience only a few CE events before ionization. As a consequence, the assumption that the neutral distribution function converges to the ion distribution function (quasi-Maxwellian assumption) is a bad approximation that leads to overestimating the width of the broad line by a factor up to ~20%. When measuring the Balmer line width in real sources, one can then easily obtain unrealistically high values for the CR acceleration efficiency, simply as a result of an erroneous prediction of the line width for any given shock velocity. We find that for an acceleration efficiency of ~20%, the FWHM of the broad Balmer line only decreases by 10% compared to the case in which CRs are absent. Because the effect is so modest, it is clear that if one erroneously expects a very broad line for a given shock velocity, based on the assumption of equilibration between ions and neutrals, a very strong acceleration efficiency (>50%) could easily be inferred.

In Sect. 2 we summarize the calculations of the neutral distribution function and discuss the improvements we implemented with respect to our previous work. In Sect. 3 we discuss our results for the distribution of neutrals and the FWHM of the broad Balmer line and compare them with those in the literature. We conclude in Sect. 4.

2. Improved calculation of atomic rates

We consider a plane-parallel shock wave propagating in a partially ionized proton-electron plasma with velocity V_sh along the z direction. The ionization fraction is fixed at upstream infinity (ISM) where ions and neutrals are assumed to be in thermal equilibrium with each other. The shock structure is determined by the interaction of CRs and neutrals with the background plasma, as discussed in detail in Paper III. In this section we limit ourselves to summarizing the computational method we used to determine the distribution function of neutral atoms in phase space and illustrate an improved approximation for the ion-neutral scattering rates.

The basic equations that describe the behavior of neutral hydrogen have been discussed in our Paper I: neutral hydrogen atoms interact with protons through CE and ionization and with electrons through ionization alone. The hydrogen distribution function, f_N(v,z) satisfies the stationary Boltzmann equation ${\mathit{v}}_{\mathit{z}}\frac{\mathit{\partial }{\mathit{f}}_{\mathit{N}}\mathrm{\left(}{v}\mathit{,z}\mathrm{\right)}}{\mathit{\partial z}}\mathrm{=}{\mathit{\beta }}_{\mathit{N}}{\mathit{f}}_{\mathit{i}}\mathrm{\left(}{v}\mathit{,z}\mathrm{\right)}\mathrm{-}\left[{\mathit{\beta }}_{\mathit{i}}\mathrm{+}{\mathit{\beta }}_{\mathrm{e}}\right]{\mathit{f}}_{\mathit{N}}\mathrm{\left(}{v}\mathit{,z}\mathrm{\right)}\mathit{,}$(1)where z is the distance from the shock (which is located at the origin, z = 0), v_z is the velocity component along the z-axis, and where the electron and proton distribution functions, f_e(v,z) and f_i(v,z), are assumed to be Maxwellian at each position.

The collision terms, β_kf_l, describe the interaction (due to CE and/or ionization) between the species k and l. The interaction rate β_k can be formally written as ${\mathit{\beta }}_{\mathit{k}}\mathrm{\left(}{v}\mathit{,z}\mathrm{\right)}\mathrm{=}\mathrm{\int }\mathrm{d}{w}\hspace{0.17em}{\mathit{v}}_{\mathrm{rel}}\hspace{0.17em}\mathit{\sigma }\mathrm{\left(}{{v}}_{\mathrm{rel}}\mathrm{\right)}{\mathit{f}}_{\mathit{k}}\mathrm{\left(}{w}\mathit{,z}\mathrm{\right)}\mathit{,}$(2)where v_rel =  | v − w | is the relative velocity, and σ is the cross section for the relevant interaction process. More precisely, β_N is the rate of CE of an ion that becomes a neutral, β_i is the rate of CE plus ionization of a neutral due to collisions with protons, while β_e is the ionization rate of neutrals due to collisions with electrons. A full description of the cross sections used in the calculations can be found in our Paper II.

Since the dynamics of electrons are not well known, we parametrize the electron temperature as a given fraction of the proton temperature and investigate the dependence of the results on such a parameter. In particular, we define β_down ≡ T_e/T_p in the region downstream of the shock and a similar parameter upstream of the shock.

A semi-analytical method for solving Eq. (1) was discussed in Paper I, where we showed how to calculate f_N starting from the distribution of protons and electrons. The method is based on decomposing the distribution function f_N at each location, upstream and downstream, as the sum of the neutrals that have undergone 0, 1, 2, ..., k CE processes. Each distribution function is named ${\mathit{f}}_{\mathit{N}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}$, and clearly ${\mathit{f}}_{\mathit{N}}\mathrm{\left(}\mathit{z,}{v}\mathrm{\right)}\mathrm{=}{\sum }_{\mathit{k}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}{\mathit{f}}_{\mathit{N}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}\mathrm{\left(}\mathit{z,}{v}\mathrm{\right)}$. Using this formalism, Eq. (1) can be rewritten as a set of k + 1 equations, one for each component: $\begin{array}{ccc}& & {\mathit{v}}_{\mathit{z}}{\mathit{\partial }}_{\mathit{z}}{\mathit{f}}_{\mathit{N}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{=}\mathrm{-}\left[{\mathit{\beta }}_{\mathit{i}}\mathrm{+}{\mathit{\beta }}_{\mathrm{e}}\right]{\mathit{f}}_{\mathit{N}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}} \mathrm{for}\mathit{k}\mathrm{=}\mathrm{0}\mathit{,} \mathrm{and}\\ & & {\mathit{v}}_{\mathit{z}}{\mathit{\partial }}_{\mathit{z}}{\mathit{f}}_{\mathit{N}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}\mathrm{=}{\mathit{\beta }}_{\mathit{N}}^{\mathrm{(}\mathit{k}\mathrm{-}\mathrm{1}\mathrm{)}}{\mathit{f}}_{\mathit{i}}\mathrm{-}\left[{\mathit{\beta }}_{\mathit{i}}\mathrm{+}{\mathit{\beta }}_{\mathrm{e}}\right]{\mathit{f}}_{\mathit{N}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}} \mathrm{for}\mathit{k}\mathit{>}\mathrm{0.}\end{array}$Each one of these equations is linear, and the solution for ${\mathit{f}}_{\mathit{N}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}$ can be obtained from ${\mathit{f}}_{\mathit{N}}^{\mathrm{(}\mathit{k}\mathrm{-}\mathrm{1}\mathrm{)}}$. The boundary conditions are imposed at upstream infinity where ions and neutrals are assumed to start with the same bulk velocity and temperature, so that CE occurs at equilibrium and the distributions do not change. The resulting boundary conditions are ${\mathit{f}}_{\mathit{N}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{(}\mathrm{-}\mathrm{\infty }\mathit{,}{v}\mathrm{)}\mathrm{=}\mathrm{(}{\mathit{n}}_{\mathit{N,}\mathrm{0}}\mathit{/}{\mathit{n}}_{\mathit{i,}\mathrm{0}}\mathrm{)}{\mathit{f}}_{\mathit{i}}\mathrm{(}\mathrm{-}\mathrm{\infty }\mathit{,}{v}\mathrm{)}$ and ${\mathit{f}}_{\mathit{N}}^{\mathrm{(}\mathit{k}\mathit{>}\mathrm{0}\mathrm{)}}\mathrm{(}\mathrm{-}\mathrm{\infty }\mathit{,}{v}\mathrm{)}\mathrm{=}\mathrm{0}$. We can write the solution of these equations in the following implicit form: ${\mathit{f}}_{\mathit{N}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{\left(}\mathit{z,}{\mathit{v}}_{\mathrm{\parallel }}\mathit{,}{\mathit{v}}_{\mathrm{\perp }}\mathrm{\right)}\mathrm{=}{\mathit{f}}_{\mathit{N}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{(}\mathrm{-}\mathrm{\infty }\mathit{,}{\mathit{v}}_{\mathrm{\parallel }}\mathit{,}{\mathit{v}}_{\mathrm{\perp }}\mathrm{)}\hspace{0.17em}\exp \left[\mathrm{-}{\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathit{z}}\frac{\mathrm{d}{\mathit{z}}^{\mathrm{\prime }}}{{\mathit{v}}_{\mathrm{\parallel }}}{\mathit{\beta }}_{\mathit{i}}\mathrm{\left(}{\mathit{z}}^{\mathrm{\prime }}\mathit{,}{\mathit{v}}_{\mathrm{\parallel }}\mathit{,}{\mathit{v}}_{\mathrm{\perp }}\mathrm{\right)}\right]\mathit{,}$(5)for k = 0, and ${\mathit{f}}_{\mathit{N}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}\mathrm{=}\mathrm{-}{\mathrm{\int }}_{\mathit{z}}^{\mathrm{+}\mathrm{\infty }}\frac{\mathrm{d}{\mathit{z}}^{\mathrm{\prime }}}{{\mathit{v}}_{\mathrm{\parallel }}}{\mathit{\beta }}_{\mathit{N}}^{\mathrm{(}\mathit{k}\mathrm{-}\mathrm{1}\mathrm{)}}{\mathit{f}}_{\mathit{i}}\exp \left[{\mathrm{\int }}_{\mathit{z}}^{{\mathit{z}}^{\mathrm{\prime }}}\frac{\mathrm{d}{\mathit{z}}^{\mathrm{\prime \prime }}}{{\mathit{v}}_{\mathrm{\parallel }}}{\mathit{\beta }}_{\mathit{i}}\right]\mathit{,}{\mathit{v}}_{\mathrm{\parallel }}\mathit{<}\mathrm{0}\mathit{,}$(6)and ${\mathit{f}}_{\mathit{N}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}\mathrm{=}{\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathit{z}}\frac{\mathrm{d}{\mathit{z}}^{\mathrm{\prime }}}{{\mathit{v}}_{\mathrm{\parallel }}}{\mathit{\beta }}_{\mathit{N}}^{\mathrm{(}\mathit{k}\mathrm{-}\mathrm{1}\mathrm{)}}{\mathit{f}}_{\mathit{i}}\exp \left[{\mathrm{\int }}_{\mathit{z}}^{{\mathit{z}}^{\mathrm{\prime }}}\frac{\mathrm{d}{\mathit{z}}^{\mathrm{\prime \prime }}}{{\mathit{v}}_{\mathrm{\parallel }}}{\mathit{\beta }}_{\mathit{i}}\right]\mathit{,}{\mathit{v}}_{\mathrm{\parallel }}\mathit{>}\mathrm{0}\mathit{,}$(7)for k > 0. Notice that v_|  | = v_z is the velocity component along z, while v_⊥ is orthogonal to z.

A good approximation for the total distribution function of neutrals is obtained by adding a sufficiently large number of ${\mathit{f}}_{\mathit{N}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}$. How many are added is determined by the physical scale of the problem. A rough estimate can be obtained by comparing the length scale for CE with the maximum between the ionization length scale and the precursor length. In general, for typical values of the ambient parameters, the longest scale, hence the one relevant to this comparison, is the size of the CR-induced precursor.

The most expensive part for computation is calculating the interaction rates ${\mathit{\beta }}_{\mathit{N}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}$ from Eq. (2). To reduce the computational time, we used an approximated expression in Paper I, originally proposed by Pauls et al. (1995), which is based on two assumptions: i) the partial distribution functions ${\mathit{f}}_{\mathit{N}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}$ are assumed to be Maxwellian, and ii) the cross section is assumed to be constant (or slowly changing as a function of the relative velocity). Under these assumptions, the interaction rates can be expressed as ${\mathit{\beta }}_{\mathit{N}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}\mathrm{\left(}\mathit{z,}{\mathit{v}}_{\mathrm{\left|}\mathrm{\right|}}\mathit{,}{\mathit{v}}_{\mathrm{\perp }}\mathrm{\right)}\mathrm{=}{\mathit{m}}_{\mathrm{p}}{\mathit{n}}_{\mathit{N}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\hspace{0.17em}\mathit{\sigma }\mathrm{\left(}{\mathit{U}}_{\mathrm{\ast }}\mathrm{\right)}{\mathit{U}}_{\mathrm{\ast }}\mathit{,}$(8)where ${\mathit{U}}_{\mathrm{\ast }}\mathrm{=}\sqrt{\frac{\mathrm{4}}{\mathit{\pi }}{{\mathit{v}}_{\mathrm{th}\mathit{,}\mathrm{N}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}}^{\mathrm{2}}\mathrm{+}\mathrm{(}{\mathit{v}}_{\mathrm{bulk}\mathit{,}\mathrm{N}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}\mathrm{-}{\mathit{v}}_{\mathrm{\parallel }}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}{\mathit{v}}_{\mathrm{\perp }}^{\mathrm{2}}}$, and σ(U^∗) is the relevant cross section calculated at U^∗. Here ${\mathit{n}}_{\mathit{N}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}$ is the numerical density of particles in the distribution function ${\mathit{f}}_{\mathit{N}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}$, while ${\mathit{v}}_{\mathrm{bulk}\mathit{,}\mathrm{N}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}$ and ${\mathit{v}}_{\mathrm{th}\mathit{,}\mathrm{N}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}$ are the moments of the same distribution corresponding to bulk and thermal speed.

In our calculations the ${\mathit{f}}_{\mathit{N}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}$ are neither assumed to be Maxwellian nor found to be Maxwellian, and most of the difference between the results presented here and in previous literature arises exactly from our more realistic treatment of the neutral distribution in phase space. However, we find that the scattering rates can be computed in an accurate way through the following procedure, which allows us to save a large amount of computation time: we compute the moments of the ${\mathit{f}}_{\mathit{N}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}$ and then introduce an equivalent Maxwellian, with the same moments, only with the purpose of calculating the scattering rates. We checked that this procedure does not introduce any critical aspect in the calculations and therefore retained this part of the computation here as well. On the other hand, we relaxed the assumption here that the cross sections are weak functions of the relative velocity. Following Chevalier et al. (1980) the three-dimensional integral in Eq. (2) can be reduced to only one dimension with a change of variable. The final result can be written as ${\mathit{\beta }}_{\mathit{N}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}\mathrm{\left(}\mathit{v,}{\mathit{v}}_{\mathrm{th}}\mathit{,z}\mathrm{\right)}\mathrm{=}{\mathit{n}}_{\mathit{N}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}\mathrm{\left(}\mathit{z}\mathrm{\right)}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\mathrm{d}{\mathit{v}}_{\mathrm{r}}\frac{\mathit{\sigma }\mathrm{\left(}{\mathit{v}}_{\mathrm{r}}\mathrm{\right)}{\mathit{v}}_{\mathrm{r}}^{\mathrm{2}}}{\mathit{v}{\mathit{v}}_{\mathrm{th}}\sqrt{\mathit{\pi }}}\left[{\mathrm{e}}^{\mathrm{-}{\left(\frac{\mathit{v}\mathrm{-}{\mathit{v}}_{\mathrm{r}}}{{\mathit{v}}_{\mathrm{th}}}\right)}^{\mathrm{2}}}\mathrm{-}{\mathrm{e}}^{\mathrm{-}{\left(\frac{\mathit{v}\mathrm{+}{\mathit{v}}_{\mathrm{r}}}{{\mathit{v}}_{\mathrm{th}}}\right)}^{\mathrm{2}}}\right]\mathit{,}$(9)where σ is now an arbitrary function of the relative velocity v_r. The integration in Eq. (9) is performed numerically only once and for each cross section involved in the calculation on a two-dimensional grid in the variables (v,v_th). Using Eq. (9) rather than Eq. (8) the value of ${\mathit{\beta }}_{\mathit{N}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}$ is more accurate by ~10%. At the same time, this improved version of the calculation is also found to be computationally faster by roughly 20%.

3. Results

The most important result of this work is that for shocks moving with speed over ~2500 km s^-1, the neutral distribution function does not reach thermal equilibrium with the ions distribution, not even after many CE reactions. As a consequence, we predict that the widths of the broad Balmer line appreciably smaller than estimated in previous calculations, where the assumption was made that ions and neutrals would thermalize after three CE reactions (vAHMR). The implications for the estimates of the CR acceleration efficiency are rather dramatic and discussed below.

As a benchmark case, we consider a shock with speed V_sh = 4000 km s^-1, total upstream density 1   cm^-3, and ionization fraction x_ion = 50%. With these values for the parameters, the typical interaction length for both CE and ionization is ~ 10¹⁶ cm. For this specific case, we assume, as may be expected for a collisionless shock, that the background electrons and protons are not in thermal equilibrium. We take β_down = 0.01, although this assumption is not crucial for our conclusions. We first consider a case with no particle acceleration.

In Fig. 1 we plot the partial distribution functions ${\mathit{f}}_{\mathit{N}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}$ for k = 1–6 at three different locations downstream, as a function of v_|  | and for v_⊥ = 0, comparing them with the downstream distribution function of ions. Each ${\mathit{f}}_{\mathit{N}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}$ includes three different contributions: a broad component emitted by neutrals that suffered at least one CE downstream, an intermediate component produced by atoms resulting from CE in the neutral return flux precursor upstream (Paper I), and finally a narrow component emitted by truly cold (T ~ 10⁴ K) atoms. As far as the atoms responsible for the broad component are concerned, we clearly see that even for k = 6 (namely after six CE events, basically all in the downstream in this case), their distribution does not approach a Maxwellian, and one can easily judge that its peak is still at a different velocity than for ions. We can make this argument more quantitative by calculating the moments of the distributions.

In Fig. 2 we plot several quantities related to the neutral distribution as a function of the position downstream, for the same benchmark case. Panels (b), (c), and (d) show, for each ${\mathit{f}}_{\mathit{N}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}$, the behavior of density, bulk speed and “equivalent” thermal speed of the broad component (the bulk and thermal velocities are normalized to the same quantities as computed for the ions’ distribution function). The region in which these parameters are most important is where the bulk of Balmer emission is produced. This region is easily identified from panel (a), which shows the numerical density and the Hα emission of both the broad and narrow line emitting components: the bulk of the broad Hα component is produced in the region 10¹⁵ < z < 10¹⁶ cm.

It is easy to see from panel (c) that in the region between 10¹⁵ and 10¹⁶ cm, the bulk speed of all ${\mathit{f}}_{\mathit{N}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}$ is greater than the bulk speed of ions, although v_bulk decreases for higher values of k. We also note that for z ≲ 10¹⁴ cm, v_bulk becomes negative, which is a consequence of the neutral return flux (see Paper I). We see in panel (d) that the thermal speed as calculated for any ${\mathit{f}}_{\mathit{N}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}$ is less than the thermal speed of ions by ~ (25−30)%. These two facts prove that the ionization starts dominating CE before neutrals and ions may reach thermal equilibrium.

Fig. 1 Partial distribution functions of neutrals, ${\mathit{f}}_{\mathit{N}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}$ for k ranging from 1 to 6, as a function of v|  | and for v⊥ = 0. Different panels show the results for different locations downstream: z = 5 × 1014cm (top), 1015 cm (middle), and 5 × 1015 cm (bottom). For comparison, the thick solid lines show the downstream distribution of ions.

Fig. 2 Panel a) total density of broad (dashed) and narrow neutrals (dot-dashed) and the sum of the two (solid line). The thin lines show the Balmer emission of narrow (dashed) and broad (dot-dashed) components (the first divided by 3). Panel b) density of the partial component of broad neutrals normalized to the total density of broad neutrals, for k = 1 to 6. Panel c) bulk speed of the partial component of broad neutrals normalized to the bulk speed of ions. Panel d) thermal speed of the partial component of broad neutrals normalized to the thermal speed of ions. The labels in panel c) also refer to panels b) and d).

Finally, from panel (b) of Fig. 2 we see that for z < 10¹⁵ cm the partial distribution functions with k > 2 contribute about half of the total density of broad neutrals, while for z > 10¹⁵ cm their contribution increases to more than 50%. In other words, the bulk of broad Balmer line emission is dominated by neutrals that have undergone more than two CE events. This is especially important in the comparison of our results with those of van Adelsberg et al. (2008). Like us, they also separate the neutral distribution function into a sum of partial functions, each defined by the number of CE events undergone by the particles. Unlike ours, their partial functions are averaged over the downstream region. There is no description of the upstream region in van Adelsberg et al. (2008). However, what is most important, and makes the biggest difference between our results and those of vAHMR, is their assumption that after two CE events the neutral distribution function acquires the same moments as f_ion, namely, same temperature and same bulk speed.

In the following, we compare our results with those of vAHMR, and show how their assumptions can easily lead to overestimate the CR energy content at the shock. Before proceeding, however, we want caution that a direct comparison is actually difficult, because of the differences in the two methods: vAHMR only account for the dynamics in the downstream, while we also consider the dynamics of ions and neutrals in the upstream, so that our ${\mathit{f}}_{\mathit{N}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}$’s also include particles that undergo one or more CE events in the upstream.

In Fig. 3 we plot the FWHM of the broad Balmer line as a function of the shock speed. Three different sets of curves are shown: the results obtained using our full method and the results of vAHMR. For the sake of comparison of our calculations with those of vAHMR, we also plot the FWHM that would be obtained from our calculations if we made the artificial assumption that all distributions ${\mathit{f}}_{\mathit{N}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}$ with k > 2 Maxwellians with the same moments as the ions’ distribution. For each set of curves, thick and thin lines represent the non equilibration (NE: β_down = 0.01) and full equilibration (FE: β_down = 1) cases between protons and electrons, respectively.

Our results and those of vAHMR coincide for shock speeds ≲2000 km s^-1 for the NE case (or ≲2500 km s^-1 for the FE case), as one may expect based on the fact that ionization becomes important for such low velocities only after a large number of CE path lengths, thereby making the assumption of thermalization between ions and neutrals a better approximation.

Fig. 3 Comparison between the FWHM of the broad Balmer line as calculated in the present work (solid lines) with those calculated by van Adelsberg et al. (2008, dot-dashed lines). The dashed lines represent the FWHM calculated using our technique, but requiring that the partial distribution functions of broad neutrals are Maxwellian for k > 2, with the same temperature and the same bulk velocity as the ions. Dots show the FWHM calculated in the presence of CRs with acceleration efficiency equal to 0.1 and 0.2 and for shock speeds of 4000 and 6000 km s-1.

To provide an easy way to compare our results with measured values of the FWHM of broad Balmer lines, we looked for a fit to the curves shown in Fig. 3. A good fit is provided by the expression, ${\mathit{W}}_{\mathrm{0}}\mathrm{\left(}{\mathit{V}}_{\mathrm{sh}}\mathit{,}{\mathit{\beta }}_{\mathrm{down}}\mathrm{\right)}\mathrm{=}\mathit{A}\sqrt{\frac{\mathrm{3}\ln \mathrm{2}}{\mathrm{2}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{\beta }}_{\mathrm{down}}\mathrm{)}}}\hspace{0.17em}\frac{{\mathit{V}}_{\mathrm{sh}}}{{\mathrm{[}\mathrm{1}\mathrm{+}{\mathrm{(}{\mathit{V}}_{\mathrm{sh}}\mathit{/}{\mathit{V}}_{\mathrm{br}}\mathrm{)}}^{\mathrm{3}}\mathrm{]}}^{\mathrm{1}\mathit{/}\mathrm{5}}}\mathit{,}$(10)which in the limit of low V_sh gives back the linear dependence of line width on shock speed. The term in the square root gives the slope expected if the neutrals have a Maxwellian distribution with exactly the same temperature as ions (see e.g. Heng 2010). The parameter A represents the deviation from this ideal situation, and its value changes only slightly between the NE case (A = 1.126) and the FE case (A = 1.096), while the value of the break velocity V_br changes from V_br = 2104 km s^-1 (for NE) to V_br = 3184 km s^-1 (for FE). Equation (10) provides a good fit, giving an error less than 2%, if the shock speed value is in the range between 500 and 7000 km s^-1.

Aside from the level of equilibration between electrons and protons, the FWHM of the broad Balmer line can also be affected, at least in principle, by the presence of He. The main effect of He is associated to the possibility that fully ionized He by itself can thermalize to a downstream temperature which is m_He/m_p times the proton temperature. If collisionless processes lead to thermal equilibrium of the two species (Coulomb collisions are too slow to have any effect on the relevant time scales), then the presence of He would be reflected in a downstream temperature possibly as high as ~ 1.3T_p. Clearly, in order for this effect to affect Balmer emission, thermalization would have to occur on a spatial scale smaller than the ionization length of hydrogen atoms. On the other hand, if He⁺ and He⁺⁺ were present, they could play an active role in CE and ionization of H atoms, thereby modifying the distribution function relevant for the Balmer emissivity.

While in our calculations we neglect the contribution of He, vAHMR take it into account, under the assumption that it is completely neutral upstream. With this assumption, the effects of He can only appear after it is ionized downstream at least once. This immediately suggests that the differences shown in Fig. 3 between the two calculations (ours and vAHMR’s) cannot be accounted for by the different treatments of He. The results of the two calculations are in fact coincident for slow shocks. But this is exactly the regime for which the effects of He should be most important, because the ratio between its ionization length and that of hydrogen atoms is lower for low shock speeds.

The conclusion we draw is that neglecting the presence of He in our calculations cannot be the main source of the differences between our results and those of vAHMR. At the same time we can conclude, based on vAHMR’s calculations, that He has negligible effects on the system. In fact, their assumption (namely of negligible He ionization upstream) should apply to most (if not all) of the sources of interest, since from Saha’s formula one easily sees that having an appreciable fraction of He⁺/He implies a totally negligible fraction of neutral H, if the overall plasma is in thermal equilibrium. A possible exception to this scenario could be produced by the presence of ionizing radiation which tends to ionize He.

We now turn our attention to the main point we want to address in this paper, namely estimating the CR acceleration efficiency from measurement of the FWHM of the broad Balmer line. It is apparent that our refined treatment of the neutral distribution downstream leads us to expect a FWHM of the broad Balmer line that is smaller than in previous calculations. As a consequence, for a given measured value of the FWHM, we infer a lower CR acceleration efficiency. In fact, in some SNRs, it may turn out to be difficult or even impossible to reach a definite conclusion concerning the presence of CR acceleration.

The theory of Balmer-dominated CR modified shocks has been outlined in our Paper III, where we showed in a quantitative way how the FWHM of the broad Balmer line decreases with increasing acceleration efficiency, as expected on qualitative grounds. This approach has been improved further here with a more accurate calculation of the CE and ionization scattering rates, as discussed in Sect. 2.

Figure 3 represents the FWHM of the broad Balmer line for the cases of NE (β_down = 0.01) and full equilibration (β_down = 1) and for acceleration efficiency of 10% and 20%. Two values of the shock speed have been considered: 4000 and 6000 km s^-1. Here we define the acceleration efficiency as ${\mathit{ϵ}}_{\mathrm{CR}}\mathrm{=}{\mathit{P}}_{\mathrm{CR}}\mathit{/}\mathrm{\left(}{\mathit{\rho }}_{\mathrm{0}\mathit{,}\mathrm{tot}}{\mathit{V}}_{\mathrm{sh}}^{\mathrm{2}}\mathrm{\right)}$, where P_CR is the CR pressure at the shock position, and ρ_0,tot is the total density at upstream infinity (ions plus neutrals). This definition could be somewhat misleading, in that the shock can only accelerate the ionized component. In this sense, the actual efficiency would be ${\mathit{ϵ}}_{\mathrm{CR}}^{\mathrm{\ast }}\mathrm{=}{\mathit{P}}_{\mathrm{CR}}\mathit{/}\mathrm{\left(}{\mathit{\rho }}_{\mathrm{0}\mathit{,}\mathrm{ion}}{\mathit{V}}_{\mathrm{sh}}^{\mathrm{2}}\mathrm{\right)}\mathrm{=}{\mathit{ϵ}}_{\mathrm{CR}}\mathit{/}{\mathit{x}}_{\mathrm{ion}}$. Since in Fig. 3 we assumed an ionization fraction of 50%, we have ${\mathit{ϵ}}_{\mathrm{CR}}^{\mathrm{\ast }}\mathrm{=}\mathrm{2}{\mathit{ϵ}}_{\mathrm{CR}}$.

From Fig. 3, it is clear that, given the strong dependence of the implied efficiency on the assumed level of electron-proton thermalization (parametrized by the value of β_down), it is generally difficult to draw firm conclusions about the CR acceleration efficiency unless additional information on β_down is available. However, if the measured value of the FWHM of the broad line falls below the thin solid line, then one can derive a robust lower limit on the CR acceleration efficiency, provided a reliable measurement of the shock speed is available.

To first approximation, the FWHM of the broad line decreases linearly with ϵ_CR, for ϵ_CR ≲ 20%. A simple linear fit to the FWHM in the presence of CRs reads as $\mathit{W}\mathrm{\left(}{\mathit{V}}_{\mathrm{sh}}\mathit{,}{\mathit{ϵ}}_{\mathrm{CR}}\mathrm{\right)}\mathrm{\approx }{\mathit{W}}_{\mathrm{0}}\mathrm{\left(}{\mathit{V}}_{\mathrm{sh}}\mathrm{\right)}\left[\mathrm{1}\mathrm{-}\mathit{B}\mathit{/}{\mathit{W}}_{\mathrm{0}}\mathrm{\left(}{\mathit{V}}_{\mathrm{sh}}\mathrm{\right)}\hspace{0.17em}{\mathit{ϵ}}_{\mathrm{CR}}\right]\mathit{,}$(11)where W₀(V_sh) is the FWHM in the absence of CR acceleration, and B = 2000 km s^-1. This implies that, in order to obtain a lower limit on the acceleration efficiency, equal, say, to a given value , the accuracy needed in the measurement of W should satisfy the condition . Needless to say, the shock velocity is needed as well with decent accuracy: the lower limit should satisfy . Such constraints can provide hints on the observational requirements needed to determine the presence of CRs in Balmer-dominated shocks.

Finally we note that in the absence of accelerated particles, the FWHM does not depend upon the choice of the density and ionization fraction. What does depend on these two parameters is the absolute value of the length scale where the processes of charge exchange and ionization take place. In the presence of accelerated particles, the results depend upon the gas density very weakly, while the dependence of the inferred acceleration efficiencies on the ionization fraction is more substantial, but it does need to be assesses on a case-by-case basis. This is why we decided to calculate the FWHM for nominal values of the gas density and ionization fraction, since the arguments presented in the paper remain valid in a qualitative way, irrespective of the exact values of parameters.

4. Conclusions

The possibility of using optical telescopes to measure the width of the Balmer lines in SNRs, and then inferring the CR acceleration efficiency, is one of the big hopes for the future of any investigations of the CR origin. The technique is, in a way, calorimetric in that the broad Balmer line contains information about the temperature of the hot proton plasma behind the SNR shock and this temperature is lower when the shock accelerates CRs efficiently. Although this is a simple consequence of energy conservation arguments, a detailed knowledge of the physical phenomena that couple ions and H atoms is needed to achieve a quantitative estimate of the CR acceleration efficiency. A complete theory of the dynamics of a Balmer-dominated, CR-modified shock has recently been presented in Paper III, where we used the mathematical formalism introduced in our Paper I.

In the present paper we have improved on one aspect of that work, by better calculating scattering rates of CE and ionization, the two physical processes responsible for the coupling between ionized and neutral material at the shock. However, the main point we want to stress in this work is that to accurately estimate the CR acceleration efficiency one has to have a space-dependent, kinetic approach to the neutrals, such as the one we put forward in Paper III. To highlight this, we compared our predicted FWHM of the broad Balmer line with the results of vAHMR, which is widely used in the literature to infer information on CR acceleration from optical measurements. In the work by vAHMR, the simplifying assumption was made that neutrals that suffered more than two CE reactions would reach thermal equilibrium with ions. Moreover, their calculation does not provide any space dependence of the distribution functions. In the absence of CRs, at 4000 km s^-1 the difference in the FWHM between the two calculations is ~1000 km s^-1 in the NE case and a few hundred km s^-1 for full equilibration between electrons and protons. This is also the order of magnitude of the reduction in the FWHM expected for CR acceleration, which implies that a proper treatment of the neutral distribution function is absolutely crucial for a reliable quantitative assessment of the CR energy content of SNRs. The assumption of ion-neutral equilibration has been shown to be incorrect for shock speeds above ~2500 km s^-1. The physical reason for this behavior is in the relative importance of CE and ionization: ions and neutrals tend to reach thermal equilibrium through CE, but unless this happens within a length equal to a few times the ionization path length at most, Maxwellian neutrals (with the same temperature as ions) are not what the Balmer line will reflect. For low shock speeds, ionization is inefficient, and thermal equilibrium between ions and neutrals is reached approximately, but for high speed shocks ionization dominates CE and the neutrals do not thermalize.

The application of these calculations has far-reaching implications: the modeling described in Paper III, when including the improvements in terms of accuracy and computation time given in the present paper, provides all the necessary tools to determine the dynamical state of a Balmer-dominated CR-modified shocks and to calculate the shape of the Balmer line (narrow, intermediate, and broad) for any CR acceleration efficiency. The whole computation is fast enough to allow for extensive applications to several SNRs of interest. Future combined observations of the overall shape of the Balmer line, as well as of shock velocity and a multifrequency spectrum from individual remnants, may open the door to a systematic estimate of the CR acceleration efficiency in SNRs. In two forthcoming papers, we will illustrate the results of our investigation of the Balmer emission from the two remnants discussed in the introduction, RCW 86 and SNR 0509–67.5.

Acknowledgments

We are grateful to D. Caprioli and K. Heng for discussions of this topic. This work was partially funded through grant PRIN INAF 2010 and ASTRI.

References

All Figures

Fig. 1 Partial distribution functions of neutrals, ${\mathit{f}}_{\mathit{N}}^{\mathrm{\left(}\mathit{k}\mathrm{\right)}}$ for k ranging from 1 to 6, as a function of v|  | and for v⊥ = 0. Different panels show the results for different locations downstream: z = 5 × 1014cm (top), 1015 cm (middle), and 5 × 1015 cm (bottom). For comparison, the thick solid lines show the downstream distribution of ions.

In the text

Fig. 2 Panel a) total density of broad (dashed) and narrow neutrals (dot-dashed) and the sum of the two (solid line). The thin lines show the Balmer emission of narrow (dashed) and broad (dot-dashed) components (the first divided by 3). Panel b) density of the partial component of broad neutrals normalized to the total density of broad neutrals, for k = 1 to 6. Panel c) bulk speed of the partial component of broad neutrals normalized to the bulk speed of ions. Panel d) thermal speed of the partial component of broad neutrals normalized to the thermal speed of ions. The labels in panel c) also refer to panels b) and d).

In the text

Fig. 3 Comparison between the FWHM of the broad Balmer line as calculated in the present work (solid lines) with those calculated by van Adelsberg et al. (2008, dot-dashed lines). The dashed lines represent the FWHM calculated using our technique, but requiring that the partial distribution functions of broad neutrals are Maxwellian for k > 2, with the same temperature and the same bulk velocity as the ions. Dots show the FWHM calculated in the presence of CRs with acceleration efficiency equal to 0.1 and 0.2 and for shock speeds of 4000 and 6000 km s-1.

In the text