# Collisionless shocks in a partially ionized medium: III.

Efficient cosmic ray acceleration

###### Abstract

In this paper we present the first formulation of the theory of non-linear particle acceleration in collisionless shocks in the presence of neutral hydrogen in the acceleration region. The dynamical reaction of the accelerated particles, the magnetic field amplification and the magnetic dynamical effects on the shock are also included. The main new aspect consists however in accounting for charge exchange and ionization of neutral hydrogen, which profoundly change the structure of the shock, as discussed in our previous work. This important dynamical effect of neutrals is mainly associated to the so-called neutral return flux, namely the return of hot neutrals from the downstream region to the upstream, where they deposit energy and momentum through charge exchange and ionization. We also present the self-consistent calculation of Balmer line emission from the shock region and discuss how to use measurements of the anomalous width of the different components of the Balmer line to infer the cosmic ray acceleration efficiency in supernova remnants showing Balmer emission: the broad Balmer line, which is due to charge exchange of hydrogen atoms with hot ions downstream of the shock, is shown to become narrower as a result of the energy drainage into cosmic rays, while the narrow Balmer line, due to charge exchange in the cosmic-ray–induced precursor, is shown to become broader. In addition to these two well-known components, the neutral return flux leads to the formation of a third component with intermediate width: this too contains information on ongoing processes at the shock.

###### Subject headings:

acceleration of particles – Cosmic Rays – Balmer emission – ISM: supernova remnants## 1. Introduction

In the context of the so-called Supernova Remnant (SNR) paradigm for the origin of Galactic Cosmic Rays (CRs), the bulk of CRs are accelerated at SNR shocks through diffusive shock acceleration (DSA). Under this assumption, in order to explain the observed fluxes at Earth and, at the same time, the B/C ratio (sensitive to the grammage traversed by CRs during propagation) one is forced to require an acceleration efficiency of per SNR, depending on details of the propagation in the Galaxy (e.g. Blasi & Amato, 2012a, b). In these conditions, the non-linear effects induced by the dynamical reaction of accelerated particles on the SNR shock cannot be neglected and reflect on the shape of the spectrum (e.g. Jones & Ellison, 1991; Malkov & Drury, 2001, for a review) and on the temperature of the ionized plasma downstream of the shock Decourchelle et al. (2000); Drury et al. (2009); Patnaude, Ellison & Slane (2009); Castro et al. (2010). Moreover, the streaming instability induced by accelerated particles ahead of the shock is thought to be responsible for the magnetic field amplification that is a crucial ingredient to reach maximum energies close to the knee Blasi, Amato & Caprioli (2007) and to explain the thin non-thermal filaments observed in X-rays (e.g. Vink, 2012).

The non-linear theory of diffusive shock acceleration (NLDSA) has been developed by several authors through several different approaches: two-fluid models (e.g. Drury & Voelk, 1981), kinetic models (e.g. Malkov, 1987; Malkov, Diamond & Voelk, 2000; Blasi, 2002) and numerical approaches, both Monte Carlo and other simulation procedures (e.g. Ellison, Moebius & Paschmann, 1990; Ellison, Baring & Jones, 1995, 1996; Kang, Jones & Gieseler, 2002; Kang & Jones, 2005; Vladimirov, Ellison & Bykov, 2006).

The description of modified shocks we adopt in the following, however, is based on a semi-analytical model presented in a series of papers by our group Amato & Blasi (2005, 2006); Caprioli et al. (2009). For fully ionized plasmas, this model gives a complete description of the full non-linear chain of processes related to efficient particle acceleration, including magnetic field amplification due to streaming instabilities, and the backreaction of the amplified field on the system. In terms of amplifying the field, both resonant (e.g. Bell, 1978) and non-resonant modes Bell (2004) of the streaming instability can contribute. However the non-resonant mode is typically important only in the case of very fast shocks Amato & Blasi (2009), corresponding to few hundred years old remnants, for a typical evolution scenario. Moreover, while the magnetic fields inferred in several SNRs seem to be in agreement with the saturation levels predicted by Bell (2004) (e.g. Voelk, Berezhko & Ksenofontov, 2005), it is not clear at all that the non resonant mode can enhance particle scattering, which is what is required for more efficient acceleration: this is because this mode has very short wavelengths compared with the typical particles’ gyroradius, and very efficient inverse cascading is needed before it can resonantly scatter the particles.

Both test-particle DSA and NLDSA lead to predict spectra which are rather hard at high energies, with a power-law index in momentum . This implies a steep dependence on energy of the Galactic diffusion coefficient Berezhko & Völk (2007), if the CR spectrum at Earth is to be reproduced. Such energy dependence appears to be at odds with the observed anisotropy of Galactic CRs Ptuskin et al. (2006); Blasi & Amato (2012b). Accounting for the velocity of scattering centers in the amplified magnetic field might lead to somewhat steeper injection spectra thereby mitigating the problem Caprioli (2012).

In general, anything that affects the structure of the shock also leads to changes in the spectrum of accelerated particles. One important source of shock modification is represented by the presence of neutral atoms in the acceleration region: as discussed by Blasi et al. (2012, hereafter Paper I), the reactions of charge exchange (CE) and ionization produce important effects in the shock region. The most important of these effects is induced by the so-called neutral return flux (NRF): a hydrogen atom crossing a collisionless shock front can undergo a charge-exchange with a hot ion downstream of the shock, thereby leading to the formation of a cold ion and a hot neutral. It may happen that the latter is produced with velocity in the direction of the shock; being neutral, it can recross the shock and suffer CE or ionization upstream, hence releasing energy and momentum there.

The existence of a NRF was highlighted for the first time by Lim & Raga (1996). These authors, however, adopted a very low fraction of neutrals, so that the effect of the NRF on the system dynamics was not very important (in addition the neutrals were treated as maxwellian, which is not a good approximation and may lead to an incorrect estimate of the NRF). The dynamical role of the returning neutrals was worked out under general conditions in Paper I, where we showed how these lead to heating of the upstream ionized plasma and to a decrease of the Mach number. Even in the context of test-particle DSA it is easy to show that this phenomenon leads to a steepening of the spectrum of particles accelerated in a SNR where appreciable neutral hydrogen is present (Paper I). The NRF affects the upstream plasma on a scale of a few interaction lengths of CE.

The presence of neutral hydrogen also represents a unique diagnostic tool of the acceleration process, through the Balmer emission associated with CE to excited states: in a Balmer dominated shock Chevalier & Raymond (1978), the Balmer line has a narrow and a broad component (see Heng, 2009, for a review), the first associated with hydrogen atoms that did not suffer CE with higher temperature ions, and the second that reflects the temperature of shocked ions, due to CE reactions downstream van Adelsberg et al. (2008); Chevalier, Kirshner & Raymond (1980). In the presence of CR acceleration the ions’ temperature downstream is lower, and this should reflect in decreased Full-Width at Half-Maximum (FWHM) of the broad Balmer line. If the CR pressure is large enough, a CR-induced precursor is formed in front of the shock Malkov & Drury (2001) so that CE can also take place upstream and eventually lead to a broadening of the narrow Balmer line. In addition, it was shown in Morlino et al. (2012, hereafter Paper II), that the shape of the Balmer line is also affected by the NRF discussed above through the formation of an intermediate component of the line.

Remarkably, anomalous widths of both broad and narrow lines have been observed in Balmer-dominated shocks. For instance, Helder et al. (2009) combined proper-motion measurements of the shock and broad H line width for the remnant RCW 86 to demonstrate that the inferred low temperature behind the shock may be interpreted as a signal that a sizable fraction of the energy is being channeled into CRs. A similar conclusion has been reached by Helder et al. (2010) by analizing the SNR 0509-67.5 in the LMC, even if in this case the shock velocity has been estimated using an evolution model with information derived from the line width of elements in the shocked ejecta (wich can be used to infer the speed of reverse shock).

Anomalous width of narrow Balmer lines have been also observed in several SNRs (see, e.g. Sollerman et al., 2003). The width of such lines is in the 30–50 km s range, implying a pre-shock temperature around 25,000–50,000 K. If this were the ISM equilibrium temperature there would be no atomic hydrogen, implying that the pre-shock hydrogen is heated by some form of shock precursor in a region that is sufficiently thin so as to make collisional ionization equilibrium before the shock unfeasible. The CR precursor is the most plausible candidate to explain such a broadening of the narrow line width.

A previous attempt to include neutral particles in the shock acceleration theory has been done by Morlino et al. (2010) using a fluid approach. Nevertheless, as discussed in Paper I, neutrals cannot be treated as a fluid because the length scales involved are much smaller than the equilibration scale. More recently Raymond et al. (2011) tried to describe the interaction of neutrals with the CR precursor using a phenomenological approach where both the CR spectrum and spatial distribution are assumed rather than calculated.

An accurate description of the CR acceleration process in the presence of neutrals can only be achieved by formulating a NLDSA theory where the dynamical action of neutrals is taken into account. This is a needed preliminary step, if one wants to make use of measurements of Balmer line emission from SNRs to quantitatively assess their efficiency as particle accelerators. In this paper we present the first treatment of NLDSA in partially ionized media and discuss how to use the theory to calculate the shape of the Balmer line and infer information on particle acceleration at the shock.

The paper is structured as follows: in Section 2 we describe the calculations used to describe neutrals, ions and the acceleration process and we explain in detail how to assemble the parts together to describe NLDSA in the presence of neutrals. The formalism previously introduced in Paper I — as far as neutrals are concerned — and in Amato & Blasi (2005, 2006); Caprioli et al. (2009) — as far as NLDSA is concerned — will be heavily used. In Section 3 we illustrate our results in terms of shock modification, spectrum of accelerated particles and shape of the Balmer line. We summarize and conclude in Section 4.

## 2. Calculations

We consider a plane-parallel shock wave propagating in a partially ionized proton-electron plasma with velocity along the direction. The fraction of neutral hydrogen is fixed at upstream infinity 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. Both CRs and neutrals profoundly change the shock structure, especially upstream where both create a precursor: the CR-induced precursor reflects the diffusion properties of accelerated particles and has a typical spatial scale of the order of the diffusion length of the highest energy particles. The neutral-induced precursor develops on a spatial scale comparable with a few interaction lengths of the dominant process between CE and ionization. The downstream region is also affected by the presence of both CRs and neutrals and the velocity gradients that arise from ionization have a direct influence on the spectrum of accelerated particles.

A self consistent description of the overall system requires to consider four mutually interacting components: thermal particles (protons and electrons), neutrals (hydrogen), accelerated protons (CRs) and turbulent magnetic field. We neglect the presence of helium and heavier chemical elements. In the following subsections we write down the basic equations used to describe each component, including the interaction terms. We are looking for stationary solutions, hence all equations are written as time-independent. The interaction terms make the system highly non linear and the solution is found using a iterative scheme similar to the one we introduced in some of our previous work (Blasi, 2002; Amato & Blasi, 2005, Paper I). The details on how to generalize such scheme to the case in which neutrals are present are discussed in § 2.5.

### 2.1. Neutral distribution

The behavior of neutrals has been extensively discussed in Paper I. Neutral hydrogen interacts with protons through CE and ionization and with electrons through ionization only. The hydrogen distribution function, , can be described using the stationary Boltzmann equation

(1) |

where is the distance from the shock (which is located at the origin), is the velocity component along the axis and the electron and proton distribution functions, and , are assumed to be Maxwellian at each position. The collision terms, , describe the interaction (due to CE and/or ionization) between the species and . The interaction rate is formally written as

(2) |

where and is the cross section for the relevant interaction process. More precisely, is the rate of CE of an ion that becomes a neutral, is the rate of CE plus ionization of a neutral due to collisions with protons, while 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 Paper II.

A novel method to solve Eq. (1) has been developed in Paper I, where we showed how to calculate starting from the distribution of protons and electrons. The method is based on decomposing the distribution function in the sum of the neutrals that have suffered 0, 1, 2, … , processes of charge exchange. Each distribution function is named , and clearly . Using this formalism, Eq. (1) can be rewritten as a set of equations, one for each component:

(3) |

(4) |

Each one of these equations is linear, and the solution for can be obtained from . The boundary conditions are imposed at upstream infinity where ions and neutrals are assumed to start with the same bulk velocity and temperature, therefore charge exchange occurs at equilibrium and the distributions do not change. As a consequence, the boundary conditions are and . A good approximation to the total distribution function of neutrals is obtained when a sufficient number of is taken into account. This number is determined by the physical scale of the problem. A rough estimate can be obtained comparing the lengthscale for CE whith the maximum between the ionization lengthscale and the precursor length. In general, for reasonable values of the ambient parameters, the longest scale, and the one relevant for comparison, is that of the CR-induced precursor.

### 2.2. CR distribution

The isotropic distribution function of CRs satisfies the following transport equation in the reference frame of the shock:

(5) |

The -axis is oriented from upstream infinity to downstream infinity with the shock located at . We assume that the injection occurs only at the shock position and is monoenergetic at , hence the injection term can be written as where . Here is the number density of ions immediately upstream of the subshock, and is the fraction of particles that is going to take part in the acceleration process. Following Blasi et al. (2005), can be related to the subshock compression factor as:

(6) |

Here is defined by the relation , where is the momentum of the thermal particles downstream. parametrizes the poorly known microphysics of the injection process and is taken as a free parameter with a typical value between 2 and 4.

The diffusion properties of particles are described by the diffusion coefficient . We assume Bohm diffusion in the local amplified magnetic field:

(7) |

where is the Larmor radius in the amplified magnetic field. The calculation of is described in § 2.3.

The solution of Eq. (5) can be easily obtained by using standard techniques. Here we generalize the approach of Blasi (2002) and Amato & Blasi (2005) in order to relax the assumption of null gradient of the distribution function downstream (homogeneity), which does not apply to the general case in which CRs and neutrals are present.

The solution of Eq. (5) in the upstream is the same found by Amato & Blasi (2005). Caprioli et al. (2010) showed that an excellent approximation to this solution is:

(8) |

where is the distribution function at the shock position. The fact that the expression in Eq. 8 provides a very good approximation to the solution of Eq. 5 is easy to understand: Eq. 8 is the exact solution of Eq. 5 when the third term can be neglected. This condition is fulfilled almost everywhere, since the gradient in the precursor is very localized (at a distance determined by the diffusion length of the particles that carry most of the energy).

Through the same approach used by Amato & Blasi (2005) we can express the solution in the downstream in the following implicit form:

(9) |

We found that an excellent approximation to the expression in Eq. (9) is given by:

(10) |

where is the slope of the distribution function at the shock position. Both Equations (9) and (10) give when the downstream speed is constant.

The distribution function at the shock is obtained, as usual, by integrating Eq. (5) across the subshock:

(11) |

Here and in the rest of the paper, labels “1” and “2” refer to quantities immediately upstream and immediately downstream of the subshock, respectively. The value of can be obtained by integrating Eq. (5) between upstream infinity () and position 1 Blasi (2002). In the standard stationary solution the downstream distribution is homogeneous, hence . This condition does not hold in our case because of ionization of neutrals downstream. The general expression can be derived from Eq. (9) and reads

(12) |

Substituting the last expression in Eq. (11) we obtain for :

(13) |

where we defined and can be interpreted as an estimate of the mean plasma velocity that a particle with momentum experiences in the precursor region (see Equation (8) in Amato & Blasi, 2005). Eq. (13) differs from the standard solution only for the term including , where

(14) |

The mean plasma speed experienced in the downstream by particles that are able to return to the subshock is . When the downstream plasma has a constant velocity, and Eq. (13) reduces to the standard solution (see Equation (7) in Amato & Blasi, 2005).

### 2.3. Transport equation for waves

Following Amato & Blasi (2006) we describe the scattering of accelerated particles in the acceleration region as due to the waves that the particles generate through resonant streaming instability. These waves are also damped due to several processes. In particular, when the plasma is not fully ionized, the presence of neutrals can damp Alfvèn waves via ion-neutral damping.

The equation for transport of waves can be written as:

(15) |

where and are, respectively, the energy flux and the pressure per unit logarithmic bandwidth of waves with wavenumber . is the growth rate of magnetic turbulence, while is the damping rate. For resonant wave amplification the growth rate of Alfvén waves is Bell (1978):

(16) |

where is the resonant momentum. The damping of the waves is mainly due to non-linear Landau damping and ion-neutral damping. For the sake of simplicity here we adopt a phenomenological approach in which the damping results in a generic Alfvén heating at a rate . This expression assumes that a fraction of the power in amplified waves is locally damped and results in heating of the background plasma Berezhko & Ellison (1999).

Following Caprioli et al. (2009), Eq. (15) can be solved by integrating in -space and using the relation between energy density and pressure of Alfvèn waves, namely , which holds upstream when . This leads to the following equation:

(17) |

where we used the normalized quantities , , . Eq. (17) can be used to derive in the upstream region.

The jump conditions for the magnetic field are calculated following Caprioli et al. (2009):

(18) |

(19) |

In the downstream section we solve the same transport Eq. (15) but we assume the relation between energy flux and pressure that derives from Eqs. (18) and (19), namely . Moreover, in the downstream region we assume for simplicity that damping occurs on scales much larger than the typical scale length for CE and ionization (this assumption is not strictly needed, it simply makes the results of the calculations easier to interpret).

### 2.4. Dynamics of the ionized gas

The dynamics of the background plasma is affected by the presence of accelerated particles and by CE and ionization of neutrals. Protons and electrons in the plasma are assumed to share the same local density, , but not necessarily the same temperature, i.e., may be different from . The equations describing the conservation of mass, momentum and energy taking into account the interactions of the plasma fluid with CRs are:

(20) |

(21) |

(22) |

Here , and are respectively the fluxes of mass, momentum and energy of neutrals along the direction (labelled as ). They can be easily computed once the neutral distribution function is known. and are the pressure and energy flux of waves as defined in section § 2.3, while is the CR pressure computed from the CR distribution function:

(23) |

The dynamical role of electrons in the conservation equations is usually neglected due to their small mass. However, collective plasma processes could contribute to equilibrate electron and proton temperatures, at least partially. If the equilibration occurs in a very efficient manner, the electron pressure cannot be neglected and the total gas pressure needs to include both the proton and electron contributions, namely , where is the electron to proton temperature ratio and is taken here as a free parameter. While it is well established that electron-ion equilibration in the downstream might be only partial Ghavamian et al. (2001, 2007), in the presence of a precursor (either induced by the CRs or by the NRF), also upstream of the shock the level of equilibration becomes an unknown.

We note that the electron-ion equilibration level can be different between upstream and downstream because the plasma conditions in these two regions are different. Hence we distinguish between upstream and downstream using two separate parameters, and , respectively. On the right side of Eq. (22) we include the heating produced by the compression of CR acting on the plasma plus the heating due to the damping of magnetic waves.

The computation strategy we adopt is as follows: we first solve Eqs. (20)-(22) for the upstream dynamics starting from the boundary conditions at upstream infinity, i.e. and . In order to determine the conditions immediately downstream of the subshock we use the jump conditions derived from the conservation equations:

(24) |

(25) |

(26) |

where the brackets indicate the difference between quantities upstream and downstream of the shock (). The terms describing neutrals and CRs are absent because of the continuity of and across the subshock. Finally we solve the conservation equations, Eqs. (20)-(22) in the downstream region, using the values of the quantities at position 2 as boundary conditions.

### 2.5. Iterative method for the solution of the non-linear problem

In order to solve the set of non-linear equations involving neutrals, ions, CRs and magnetic field, we adopt an iterative method that we illustrate in Fig. (1) in the form of a flow diagram. The input quantities are the values of the shock velocity and all environmental quantities at upstream infinity, where the distribution function of neutrals is assumed to be Maxwellian at the same temperature as that of ions. The method proceeds through two nested cycles of iteration: in the first cycle, we start with a distribution function of ions that corresponds to an unperturbed distribution that one would obtain in the absence of both neutrals and CRs and we calculate the neutral distribution function from Eq. (1). With the given we iterate to obtain the ion distribution, the CR distribution and the wave energy density. The outer loop of iteration allows us to update the neutral distribution function using the newly available ion and CR distribution functions. The iteration ends when the predefined accuracy is achieved. At this point the distribution functions are used to calculate the structure of the Balmer line, following the method illustrated in Paper II.

The typical computation time required to complete the calculations for a given set of environment parameters is of the order of one hour on a standard laptop. The most time-consuming part is the solution of Eq. (1) which increases linearly with the number of partial functions, , required for a correct description of the neutral distribution function. In particular the number of increases linearly with the ratio between the precursor length and the CE interaction length.

## 3. Results

The formalism presented above represents the first theory of particle acceleration at collisionless shocks in the presence of neutral atoms and allows us to calculate the shock structure, the spectrum of accelerated particles and the Balmer emission from the shock region. In the subsections below we describe the main results concerning these three topics.

### 3.1. Shock structure and spectra of accelerated particles

We adopt a test case consisting of a shock moving at in a medium with temperature K, magnetic field G, and densities (the neutral fraction is 50%). CRs are injected with an efficiency of (in fact we fix the injection parameter ) and the efficiency of Turbulent Heating (TH hereafter) is fixed as (see Eq. 17). The shock structure is illustrated in Fig. 2: the top panel shows the pressure of gas (continuous solid line), magnetic field (dotted line) and CRs (dashed line) as functions of , plotted here in logarithmic scale for both the upstream and downstream regions. All pressure components are normalized to the ram pressure at upstream infinity. Clearly the CR pressure is continuous across the shock front. The self-generated magnetic field acquires a pressure that is of the same order of magnitude of the thermal pressure, which implies that its dynamical effect is non-negligible, as already discussed in Caprioli et al. (2008, 2009).

The background ions are substantially heated in the precursor due to TH, as it is shown by the increase in the gas temperature on a spatial scale comparable with , i.e., the diffusion length of particles at (mid panel of Fig. 2). The upstream profile of the fluid velocity, showed in the bottom panel of Fig. 2 as normalized to , is a distinctive signature of the CR-induced precursor.

Finally, one last thing to notice in this Figure is the increase of both the gas pressure and velocity at some distance from the shock in the downstream. This is a result of progressive ionization of fast cold neutrals, and indeed occurs on scales which are of the order of the ionization length. We will see that the situation is different, with the plasma velocity decreasing, rather than increasing, towards downstream infinity, when the shock velocity is lower and the NRF starts to play an important role (see Fig. 4 below). Spatial variations of the plasma velocity are important for determining the accelerated particle spectrum at the shock, whose slope at each given energy always reflects the mean compression ratio experienced by particles of that energy while moving away from the shock.

The spectrum of accelerated particles at different locations upstream and downstream of the shock is shown in Fig. 3 for a neutral fraction . The solid lines refer to the downstream section of the plasma, the thick dot-dashed line is the spectrum at the subshock location, the dashed lines are the spectra at different locations upstream of the shock.

One can see that, moving away from the shock in the upstream direction, the spectrum is more and more depleted of the lowest energy particles: at any position , only particles energetic enough to have a diffusion length are present. The transport of low-energy particles in the downstream, is rather dominated by advection, and the spectrum is never truncated: at downstream infinity, the low-energy spectrum has the same slope as at the shock, but with a higher normalization. This is due to the compression that occurs around , where ionization takes place. At the highest energies, instead, diffusion dominates over advection and the CR spectrum far downstream is expected (and found) to be exactly the same as at the shock location. As a result, at large distances downstream of the shock the spectrum turns out to be slightly steeper than at the shock.

As already illustrated in Paper I in the test-particle limit, the effect of the NRF is to create an upstream precursor that steepens the spectra of accelerated particles. For this reason the spectra of accelerated particles are less concave than they would be in the presence of the CR reaction alone. This result is qualitatively general, though the strength of the effect is quantitatively dependent upon the specific environment that the shock propagates into, as discussed below.

In Fig. 4 we show the normalized plasma velocity as a function of the distance from the shock for different values of the neutral fraction (left panel), and the slope of the spectrum of accelerated particles as a function of momentum (right panel). All the calculations refer to a shock speed , a total gas density at upstream infinity , an unperturbed magnetic field strength G and an injection parameter . The solid lines refer to the case with fully ionized plasma (). In this case the shock is appreciably modified by the CR pressure, the acceleration efficiency being , and the spectrum is quite concave. The spectrum is steeper than only for energies below GeV.

Due to the relatively-low shock velocity, the effects of CE and ionization are quite relevant; therefore, we expect the results to change appreciably for non vanishing values of , due to the NRF. This is illustrated by the long-dashed, dashed and dotted lines in Fig. 4. For (long-dashed line) and (dashed line) the formation of the precursor upstream is dominated by the dynamical reaction of CRs on the plasma, although the heating of the upstream due to the NRF is very evident. In these two cases the acceleration efficiency is (for ) and (for ). The spectra in the same cases are correspondingly less modified by the CR pressure (less concave) and, on average, somewhat steeper: for the spectrum is steeper than for energies GeV, while for it remains steeper than even at GeV.

The case with large neutral fraction, is rather extreme: the precursor becomes very short, as it becomes dominated by the NRF, and the acceleration efficiency is still . The spectrum is quite steeper than at all energies, so that the bulk of the pressure is carried by relatively low energy particles. This implies that the CR pressure can slow down the incoming plasma only very close to the subshock (see left panel in Fig. 4).

As discussed in Paper I, the effects of the NRF are most important for shock speeds ; for larger shock speeds, the typical relative velocity between neutrals and ions is such that the cross section for CE is so small that atoms get ionized before suffering a CE reaction. As a consequence, the NRF is suppressed. It is useful to explore the dependence of the spectrum on the shock velocity in a case in which the neutral fraction is . In the left panel Fig. 5 we plot the normalized plasma velocity as a function of distance from the shock for shock velocities (solid line), (dashed line) and (dotted line). The other parameters are the same as in Fig. 4. The right panel shows the spectral slope for the same values of the shock velocity. For comparison, the thick lines refer to , while the thin lines of the same type refer to completely ionized gas ().

The left panel qualitatively confirms our expectations based on the test-particle results: when increasing the shock speed from to , the effect of neutrals becomes smaller, namely the shock becomes increasingly more modified by the reaction of CRs. This is even more evident in the right panel of Fig. 5: for a neutral fraction of 50% is sufficient to make the spectrum of accelerated particles steeper than , at all energies (thick solid line). With the same shock speed, but for a completely ionized plasma, the spectrum is instead harder than above GeV (, thin solid line).

### 3.2. Balmer emission

The Balmer line is affected by the presence of CRs in at least two ways. First of all, the width of the broad Balmer line is expected to be narrower if efficient CR acceleration takes place. This is a straightforward consequence of energy conservation: energy is channelled into accelerated particles, therefore the plasma temperature downstream is lower. Moreover, the formation of a CR-induced precursor leads to a non-negligible heating upstream of the subshock, so that the narrow Balmer line may become broader. In this section we investigate these possibilities in a more quantitative way, by using the formalism introduced in Paper II. The different components of the line are identified as discussed in Paper II. The global shape of the Balmer line is fitted with three Gaussian components: a narrow line, an intermediate line (reflecting directly the effect of the NRF) and a broad line.

In Fig. 6 we show the temperature throughout the shock region for different values of the efficiency of TH, , for a shock propagating at in a medium with density . The case without CR acceleration is shown as a thick solid (black) line. In this case the precursor is solely due to the NRF and is limited to the region very close to the subshock.

When the acceleration process is turned on, several phenomena become visible. In the absence of TH (, thin solid line), the upstream fluid is mildly heated by the adiabatic compression due to the CR-induced precursor only. Nevertheless, the temperature of the downstream plasma drops by almost a factor 2 as a consequence of the energy channelled into the accelerated particles. Increasing the value of , the heating in the precursor becomes more and more marked, and it extends over larger and larger distances. However, the downstream temperature is not appreciably affected by a change in the efficiency of TH.

The shape of the global Balmer line emission is plotted in the left panel of Fig. 7, while the right panel shows a zoom in the region of the narrow Balmer line. The curves refer to the same cases as in Fig. 6, labelled as indicated. The width of the broad component of the Balmer line is appreciably reduced when CR acceleration is turned on, but the results are not sensitive to the amount of TH. On the other hand, as one can notice in the right panel, the narrow Balmer line becomes broader when CR acceleration is efficient; the effect is more pronounced when non-negligible TH is taken into account. Hence, the distribution function of neutrals becomes broader mainly because of the scattering with a warmer ion distribution in the far precursor (i.e., where TH is more effective), rather than because of the NRF, which operates only within a few CE interaction lengths from the subshock.

In Fig. 8 we show the FWHM of the broad Balmer line as a function of the CR acceleration efficiency (left panel) and as a function of the shock velocity (right panel). In the left panel, the two curves illustrate the change in the FWHM obtained by assuming full or partial temperature equilibration between electrons and protons in the downstream ( and , with the ratio between electron and ion temperature), for and .

The FWHM decreases when the energy density in CRs increases, confirming the very important notion that the width of the broad Balmer line can be used to measure the CR acceleration efficiency. There is however an important caveat to keep in mind: for a given value of , the uncertainty in the electron-proton equilibration translates into a spread in the FWHM of . Hence, unequivocal evidence of efficient CR acceleration is likely to be achievable only if the measured FWHM is below the line corresponding to full equilibration. Even in this case the value of is likely to be quite uncertain because of the unknown level of equilibration. For instance, if a FWHM of were measured for a shock with , the result could be interpreted as a measurement of in the case of full equilibration. In the more general case of partial equilibration, the above value of can be taken as a lower limit. On the other hand, for the same shock velocity, if one measured a FWHM of , it would be hardly possible to say anything at all about the efficiency of CR acceleration, without knowing the level of electron-proton thermalization. In order to have an independent estimate of , it would be very helpful to pair optical observations with the emission of SNR shocks in the (thermal) X-rays as well.

In the right panel of Fig. 8 we fix the injection parameter and we calculate the FWHM of the broad Balmer line as a function of the shock velocity. The two thin (red) lines refer to the cases with no CR acceleration with (upper curve) and (upper curve). The dashed () and dash-dotted lines () refer to the cases with CR acceleration (efficiencies in parentheses), again for (upper curves) and (lower curves). From the plot it is clear that the TH does not affect appreciably the width of the broad line. As pointed out above, if the electron-proton thermalization is inefficient, it may be difficult to draw any strong conclusion about the CR acceleration efficiency unless additional information is available to constrain the value of . Nevertheless, measuring a FWHM of the broad line below the solid thick line allows one to put a lower limit on the CR acceleration efficiency.

Additional discrimination power might derive from a combined analysis of the broad and narrow Balmer line components. The FWHM of the narrow Balmer line as a function of the maximum momentum is plotted in Fig. 9 for (top left panel), (top right panel) and (bottom panel). The three curves in each panel are obtained for . The shock velocity and total density are set to and , which approximately correspond to 0.4, 0.2 and 0.1 (weakly dependent upon ) for the three values of . The maximum momentum, determines the spatial extent of the CR-induced precursor. Larger values of imply that there is more time (space) for depositing heat in the upstream, and the width of the narrow Balmer line broadens correspondingly. The effect becomes more pronounced for larger values of the parameter . It is worth noticing that the behavior of the curves in Fig. 9 is the same in all cases and reflects the slight broadening of the narrow Balmer line as due to the NRF. In fact, one can estimate the momentum at which the diffusion length equals the size of the precursor induced by the NRF, which is of order cm for the parameters used here. By imposing one gets TeV/c (from our calculations the amplified magnetic field close to the subshock is G). This interpretation of the modest increase in the FWHM for low values of is further confirmed by the fact that the width is very weakly dependent on both and for TeV/c.

For completeness, Figure 10 shows the width of the intermediate component of the Balmer line as a function of the CR acceleration efficiency , for , and neutral fraction . The two curves refer again to (upper curve) and (lower curve). In the upstream, energy is assumed to be deposited only in the protons, which therefore become warmer than the electrons. We assume that they do not equilibrate very effectively since the Coulomb collision timescale is typically longer than the advection one.

The reason why the width of the intermediate Balmer line depends on the CR acceleration efficiency is that when increases, the amount of TH close to the subshock increases as well, leading to a broader distribution of neutrals in the region affected by the NRF. The dependence of the width of the intermediate line on lies in the fact that the NRF carries upstream the information about the ion temperature behind the shock. We further notice that the width of the intermediate component may also depend on the electron-ion equilibration upstream, . In Figure 10 we assume no equilibration upstream, namely while ions are heated by NRF and TH, electrons always have K. Hence one has to be careful in using the intermediate line to infer shock properties, because this line is a function of several parameters.

## 4. Summary and Conclusions

In this paper we have presented the first formulation of the non-linear theory of diffusive shock acceleration at SNR shocks propagating in a partially ionized medium. At the same time, we have put forward the first fully self-consistent calculation of how the shape of Balmer line is affected by the presence of accelerated particles. This development is crucial if one wants to use the width of the various components of the Balmer line to infer information on the efficiency of CR acceleration in SNRs.

The three most important physical points discussed here in a quantitative way are the following: 1) neutrals produced in CE events with hot ions downstream can return upstream and deposit energy and momentum there. This phenomenon, that we refer to as the neutral return flux (NRF), has an important dynamical effect on the structure of collisionless shocks as well as on the spectrum of accelerated particles (see also Paper I). 2) The structure of a collisionless shock is affected in a non-trivial way by the combined action of neutrals and CRs. 3) Also the shape of the Balmer line reflects the effects of the two phenomena mentioned above.

The NRF was not predicted in most of the previous calculations of the structure of SNR shocks, mainly because of their intrinsic limitation to fluid approaches (that cannot account for this kinetic effect) or restriction to the downstream plasma alone. Both CRs and the NRF produce a precursor upstream of the shock: the CR-induced precursor develops on a scale of the diffusion length of particles at the maximum energy , while the neutral-induced precursor develops on a scale of a few CE interaction lengths. Both precursors cause heating of the ion plasma upstream. The heating associated with the NRF is mediated by CE and ionization upstream. The CR-induced heating is most effective only in the presence of turbulent (namely non adiabatic) heating.

For typical values of the parameters relevant for SNRs, the spatial extent of the neutral-induced precursor is such as to overlap with the diffusion range of particles with TeV/c. In this energy region the spectrum of accelerated particles is steepened by the dynamical effect of the NRF, as already found in Paper I, in the context of test-particle calculations. The calculations presented here are a generalization of previous work on NLDSA Amato & Blasi (2005, 2006); Caprioli et al. (2009) and on the dynamical effects of neutrals on the structure of collisionless shocks Blasi et al. (2012). The technique adopted for the solution of the non-linear equations that describe the whole problem is basically an iterative method in two stages: in a first stage, we calculated iteratively the distribution of neutrals and ions (see Sections 2.1, 2.4 and Paper I). The distribution function at a generic location is not Maxwellian therefore the calculation of the neutral distribution is kinetic, while ions (protons and electrons) are treated as fluids. In a second stage, the updated profile of neutrals and ions is used to calculate the spatial distribution and spectrum of accelerated particles and waves (Sections 2.2 and 2.3). The procedure is repeated until convergence is reached.

On average, as expected, the presence of neutrals leads to less concave and steeper spectra of accelerated particles than in a completely ionized medium. These effects are quite substantial for shocks with . For faster shocks, the NRF becomes dynamically less important because neutrals are ionized downstream before they suffer a CE reaction. In these cases the shock modification is mainly induced by CRs.

Besides the general interest in a theory of NLDSA in the presence of neutrals, we think that the present work is especially important in terms of providing an essential tool to interpret the shape of the line and translate it into quantitative information on the presence of CRs in SNR shocks. The widths of the narrow and broad component of the Balmer line in SNR shocks have long been known to bear information on the temperature of the ions upstream and downstream of the shock respectively (see Chevalier & Raymond, 1978; van Adelsberg et al., 2008; Smith et al., 1991; Ghavamian et al., 2001; Heng, 2007, and the recent high quality results of Nikolić et al. (2013)). Both quantities are sensitive, in turn, to the CR acceleration efficiency: if CRs are efficiently accelerated, the plasma upstream is heated by the CR-induced precursor and the downstream plasma is heated less than in the absence of CRs because part of the ram pressure is channelled into CR acceleration. In terms of Balmer line, one expects the narrow component to be somewhat broader and the broad component to be somewhat narrower whenever CR acceleration is efficient. A phenomenological approach to the problem was recently put forward by Raymond et al. (2011), where however the CR spectrum and spatial distribution are both assumed rather than calculated. As described in the present paper, the CR properties are profoundly affected by the presence of neutrals. Moreover Raymond et al. (2011) do not treat the downstream plasma, therefore there is no NRF in their calculation.

In Fig. 6 we can see the temperature of the plasma calculated self-consistently. Both the heating in the precursor (for different values of the efficiency of TH) and the reduced temperature in the downstream (as a result of CR acceleration) are clearly visible.

These effects reflect in the shape of the Balmer line, which we calculated by using the formalism introduced in Paper II. The width of the narrow Balmer line is found to depend appreciably upon the level of TH and the value of maximum momentum . The dependence on can be understood, because increasing leads to larger CR precursors and therefore neutrals incoming from upstream infinity have a longer time to experience CE and adapt to the ion temperature. Larger TH reflects in an increased width of the narrow Balmer line.

The level of TH does not affect appreciably the plasma temperature downstream, therefore the width of the broad Balmer line mainly reflects the CR acceleration efficiency. However, the broad line is sensitive to the level of equilibration between protons and electrons downstream: if electrons and protons reach thermal equilibrium in a time short compared with CE (clearly this cannot happen through collisions) then the ion temperature is lower, thereby mimicking the presence o CRs. The estimate of the CR acceleration efficiency might be less dependent upon this ambiguity if the width of the broad Balmer line is quite small. These cases however, suggesting a very high acceleration efficiency, are very difficult to obtain in the presence of TH and magnetic reaction on the shock Caprioli et al. (2008, 2009). Moreover these cases would correspond to very concave spectra of accelerated particles and of the associated gamma rays, that so far have never been observed. The observation of the structure of the Balmer line (intensity and width of its components) might also help breaking the degeneracy between CR acceleration efficiency and electron-proton equilibration. By the same token, observation of an X-ray spectrum of thermal origin may independently provide a measurement of the electron temperature.

Finally, we stress that the existence of the NRF leads to the formation of an intermediate component of the Balmer line, with a width that depends on the temperature of protons downstream, on the electron-proton equilibration and on the level of TH. The practical possibility to infer from data on individual SNR will need to be assessed on the case-by-case basis.

## Aknowledgments

We are grateful to the anonymous referee for several comments that helped us improve the paper. This work was partially funded through grant ASI-INAF I/088/06/0 and PRIN INAF 2010. The research work of D.C. was partially supported by NSF grant AST-0807381.

## References

- van Adelsberg et al. (2008) van Adelsberg, M., Heng, K., McCray, R., and Raymond, J. C. 2008, ApJ, 689, 1089
- Amato & Blasi (2005) Amato, E. & Blasi, P., 2005, MNRAS, 364, L76
- Amato & Blasi (2006) Amato, E. & Blasi, P., 2006, MNRAS, 371, 1251
- Amato & Blasi (2009) Amato, E. & Blasi, P., 2009, MNRAS, 392, 1591
- Bell (1978) Bell, A.R., 1978, MNRAS, 182, 147
- Bell (2004) Bell, A. R., 2004, MNRAS, 353, 550
- Berezhko & Ellison (1999) Berezhko, E.G., and Ellison, D.C., 1999, ApJ, 526, 385
- Berezhko & Völk (2007) Berezhko, E.G. and Völk, H.J., 2007, ApJ, 661, L175
- Blasi (2002) Blasi, P., 2002, Astropart. Phys., 16, 429
- Blasi & Amato (2012a) Blasi, P., and Amato, E., 2012, JCAP, 01, 10
- Blasi & Amato (2012b) Blasi, P., and Amato, E., 2012, JCAP, 01, 11
- Blasi, Amato & Caprioli (2007) Blasi, P., Amato, A., and Caprioli, D., 2007, MNRAS, 375, 1471
- Blasi et al. (2005) Blasi, P., Gabici, S., Vannoni, G., 2005, MNRAS, 361, 907
- Blasi et al. (2012) Blasi, P., Morlino G., Bandiera R., Amato, E. & Caprioli, D. 2012, ApJ, 755, 121 [Paper I]
- Caprioli (2012) Caprioli, D., 2012, JCAP, 7, 38
- Caprioli et al. (2008) Caprioli, D., Blasi, P., Amato, E. & Vietri, M. 2008, ApJ, 679, L139
- Caprioli et al. (2009) Caprioli, D., Blasi, P., Amato, E. & Vietri, M. 2009, MNRAS, 395, 895
- Caprioli et al. (2010) Caprioli, D., Amato, E., Blasi, P., 2010, Astropart. Phys., 33, 307
- Castro et al. (2010) Castro, D., Slane, P. O., Gaensler, B. M., Hughes, J. P., Patnaude, D. J. 2011, ApJ, 734, 85
- Chevalier, Kirshner & Raymond (1980) Chevalier, R. A., Kirshner R. P., and Raymond, J. C. 1980, ApJ, 235, 186
- Chevalier & Raymond (1978) Chevalier, R. A., and Raymond, J. C. 1978, ApJ, 225, L27
- Decourchelle et al. (2000) Decourchelle, A., Ellison, D. C. & Ballet, J. 2000, ApJ, 543, L57
- Drury & Voelk (1981) Drury, L. O.’C., Voelk, H. J. 1981, ApJ, 248, 344
- Drury et al. (2009) Drury, L. O.’C., Aharonian, F. A., Malyshev, D., Gabici, S. 2009, A&A, 496, 1
- Ellison, Moebius & Paschmann (1990) Ellison D. C., M¬obius E., Paschmann G., 1990, ApJ, 352, 376
- Ellison, Baring & Jones (1995) Ellison D. C., Baring M. G., Jones F. C., 1995, ApJ, 453, 873
- Ellison, Baring & Jones (1996) Ellison D. C., Baring M. G., Jones F. C., 1996, ApJ, 473, 1029
- Ghavamian et al. (2001) Ghavamian, P., Raymond, J., Smith, R. C., Hartigan, P., 2001, ApJ, 47, 995
- Ghavamian et al. (2007) Ghavamian, P., Laming, J. M., Rakowski, C. E., 2007, ApJ, 654, 69
- Helder et al. (2009) Helder, E. et al. 2009, Science, 325, 719
- Helder et al. (2010) Helder, E. A., Kosenko, D., Vink, J. 2010, ApJ, 719, L140
- Heng (2007) Heng, K., van Adelsberg, M., McCray, R., Raymond, J. C., 2007, ApJ, 668, 275
- Heng (2009) Heng, K., 2009, PASA, 27, 23
- Jones & Ellison (1991) Jones, F. C. & Ellison, D. C., 1991, Space Science Reviews, 58, 259
- Kang, Jones & Gieseler (2002) Kang H., Jones T. W., Gieseler U. D. J., 2002, ApJ, 579, 337
- Kang & Jones (2005) Kang H., Jones T. W., 2005, ApJ, 620, 44
- Lim & Raga (1996) Lim, A. J., Raga, A. C., 1996, MNRAS, 280, 103
- Malkov (1987) Malkov, M. A., 1997, ApJ, 485, 638
- Malkov, Diamond & Voelk (2000) Malkov, M. A., Diamond, P. H., Voelk, H. J., 2000, ApJ, 533,171
- Malkov & Drury (2001) Malkov, M.A. and Drury, L. O’C., 2001, Rep. Prog. Phys., 64, 429
- Morlino et al. (2012) Morlino G., Bandiera R., Blasi, P. & Amato, E. 2012, 760, 137 [Paper II]
- Morlino et al. (2010) Morlino, G., Amato, E., Blasi, P., Caprioli, D. 2010, proceedings of the 12th ICATPP Conference
- Nikolić et al. (2013) Nikolić, S., van de Ven, G., Heng, K., Kupko, D., Husemann, B., Raymond, J. C., Hughes, J. P., Falcón-Barroso, J. 2013, arXiv:1302.4328
- Patnaude, Ellison & Slane (2009) Patnaude, D. J., Ellison, D. C., Slane P. 2009, ApJ, 696, 1956
- Ptuskin et al. (2006) Ptuskin, V.S., Jones, F.C., Seo, E.S., and Sina, R., 2006, Adv. in Sp. Res., 37, 1909
- Raymond et al. (2011) Raymond, J.C., Vink, J., Helder, E.A, and de Laat, A, 2011, ApJ, 731, L14
- 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
- Vink (2012) Vink, J., 2012, Astr. Astrophys. Rev, 20, 49
- Vladimirov, Ellison & Bykov (2006) Vladimirov, A., Ellison, D. C., Bykov, A., 2006, ApJ, 652, 1246
- Voelk, Berezhko & Ksenofontov (2005) Voelk H. J., Berezhko E. G., Ksenofontov L. T., 2005, A&A, 433, 229