text:
# Hydrodynamic Simulation of Supernova Remnants Including Efficient Particle Acceleration

## 1 Introduction

It is commonly believed that the shocks in supernova remnants (SNRs) produce the majority of galactic cosmic rays (CRs) with energies below $`10^{15}`$ eV via diffusive shock acceleration (e.g., Drury, 1983; Blandford & Eichler, 1987). Convincing support for the production of TeV electrons in SNRs comes from the synchrotron interpretation of nonthermal X-ray emission (e.g., Reynolds, 1998) observed in an increasing number of young SNRs such as SN1006 (Koyama et al., 1995), Cas A (Allen et al., 1997), G347.3$``$0.5 (Slane et al., 1999), and RCW 86 (Borkowski et al., 2001).111As of this writing, there is no unambiguous evidence for the production of TeV ions in SNRs (see Berezhko et al., 2002, for a discussion of SN1006 in this regard). The recent claim that TeV emission from SNR RX J1713.7-3946 (also called G347.3$``$0.5) observed by CANGAROO II is from pion-decay (e.g., Enomoto et al., 2002) is still under debate (see, for example, Reimer & Pohl, 2002; Butt et al., 2002).

In addition to their putative role in accelerating cosmic rays, the shocks in SNRs heat the ambient interstellar medium and ejecta to X-ray emitting temperatures. The interpretation of these X-ray observations leads to inferences for important quantities such as the supernova (SN) explosion energy, ejecta mass and composition, ambient densities, shock speed, and rate of electron and proton equilibration. That CR production and thermal heating in SNRs may be coupled comes from the fact that diffusive shock acceleration is intrinsically efficient in high Mach number shocks if even a small fraction of the shock heated plasma is injected into the acceleration process (e.g., Ellison & Eichler, 1984; Jones & Ellison, 1991; Berezhko & Ellison, 1999; Malkov, 1998; Blasi, 2002).

If strong coupling between particle acceleration and shock heating occurs, the modeling of efficient particle acceleration in SNRs offers the possibility of using the high-quality X-ray and $`\gamma `$-ray data currently being collected by spacecraft (e.g., Chandra, XMM-Newton, INTEGRAL) to address fundamental questions concerning the SNR origin of CRs and the underlying physics of diffusive shock acceleration, particularly the injection of thermal particles into the acceleration process. Furthermore, strong coupling implies that the inferences made from X-ray observations may differ substantially between interpretations which include particle acceleration self-consistently and those that do not.

Despite the expected efficiency of diffusive shock acceleration, X-ray spectra from SNRs have generally been modeled and interpreted assuming that the shocks place an insignificant fraction of their energy in cosmic rays. Exceptions to this include the early works of Chevalier (1983), Heavens (1984), and Boulares & Cox (1988). Chevalier (1983) investigated the effects of cosmic-ray pressure on SNR dynamics using a two-fluid, self-similar solution with an arbitrary fraction of thermal gas (adiabatic index $`\gamma =5/3`$) and relativistic gas ($`\gamma =4/3`$). More recent work has been done by Dorfi & Böhringer (1993) and Dorfi (1994). In our preliminary work (i.e., Decourchelle et al., 2000), we developed a model which coupled the approximate nonlinear (NL) acceleration calculation of Berezhko & Ellison (1999) with an analytic, self-similar description of the SNR hydrodynamics (i.e., Chevalier, 1983) and a non-equilibrium ionization calculation of X-ray emission (e.g., Decourchelle & Ballet, 1994). We illustrated this model by fitting ASCA and RXTE observations of Kepler’s SNR and found adequate fits when acceleration was efficient at the forward shock but inefficient at the reverse shock.

The effects of efficient particle acceleration on SNR hydrodynamics where calculated in a hydro computer simulation of SNRs by Blondin & Ellison (2001). This was done by globally changing the effective ratio of specific heats, $`\gamma _{\mathrm{eff}}`$, from $`5/3`$ to values approaching 1, but did not include coupling between the acceleration and the hydro. As $`\gamma _{\mathrm{eff}}`$ was decreased, the shocked gas became more compressible, the shock compression ratio increased, and the interaction region between the forward and reverse shocks narrowed. Blondin & Ellison (2001) were able to show in two and three-dimensional simulations that if the interaction region was narrow enough, convective instabilities produced Rayleigh-Taylor fingers of dense ejecta material which were able to reach and perturb the forward shock.

Here, we introduce and describe in detail a CR-Hydro model which uses the same NL acceleration calculation of Berezhko & Ellison (1999), but replaces the self-similar description used in Decourchelle et al. (2000) with a 1-D hydro simulation such as that used by Blondin & Ellison (2001). The simulation is more general than the analytic approach used by Decourchelle et al. (2000) since it is not restricted to self-similar evolution and allows for a continuous change in the acceleration efficiency as the SNR evolves. We show, however, that when acceleration efficiency is nearly constant during the self-similar phase, the two models closely correspond, providing an important check on the validity of both models. While we do not calculate X-ray thermal spectra in this paper, we show, with various examples, how the efficient production of CR protons influences SNR evolution and discuss the implications this has on the interpretation of X-ray observations.

## 2 CR-Hydro Model

Our CR-Hydro model couples a spherically symmetric hydrodynamic simulation with a calculation of nonlinear diffusive shock acceleration. The particle acceleration calculation determines how energy is divided between the thermal gas and relativistic particles, and provides the particle distribution over all energies behind the shock, as well as the effective ratio of specific heats, $`\gamma _{\mathrm{eff}}`$, defined below. In future work, we will use the electron and ion distribution functions to calculate the broad-band continuum photon emission from radio to TeV energies (e.g., Ellison et al., 2001), and use the self-consistent thermal properties in a non-equilibrium calculation of X-ray lines (e.g., Decourchelle et al., 2000). For now, we restrict ourselves to calculating CR proton spectra integrated over the SNR lifetime.

### 2.1 Hydrodynamic Simulation

We use a standard hydrodynamic simulation in one dimension to model the effects of a supernova explosion in the interstellar medium (ISM) (see Blondin & Ellison, 2001, and references therein). We are free to choose arbitrary ejecta and ISM mass density profiles, but to facilitate our comparisons discussed below, we adopt the parameters determined for SN1006 by Berezhko et al. (2002), i.e., we assume a constant density, constant temperature ISM, and take the initial ejecta density profile to be $`\rho r^n`$, with $`n=7`$ and a constant density plateau at small radii (i.e., Chevalier, 1982).222At the start of the simulation, the ejecta temperature is assumed to be low enough to be insignificant. While the modeling of particular, young SNRs depends critically on the ejecta and ISM densities (e.g., Decourchelle et al., 2000), the general character of the results we present here are insensitive to these details. Specifically, we begin at some time $`t_{\mathrm{start}}<0.1t_{\mathrm{ch}}`$ with undisturbed ejecta and ISM separated by a contact discontinuity. Here, $`t_{\mathrm{ch}}=R_{\mathrm{ch}}/V_{\mathrm{ch}}`$ is the characteristic age with $`R_{\mathrm{ch}}=[3M_{\mathrm{ej}}/(4\pi \rho _0)]^{1/3}`$, $`V_{\mathrm{ch}}=\sqrt{2E_{\mathrm{sn}}/M_{\mathrm{ej}}}`$, and $`\rho _0=(1+4f_{\mathrm{He}})m_pn_{\mathrm{p0}}`$, where $`E_{\mathrm{sn}}`$ is the explosion energy, $`M_{\mathrm{ej}}`$ is the ejecta mass, $`n_{\mathrm{p0}}`$ is the ISM proton number density, $`f_{\mathrm{He}}`$ is the helium to proton number ratio, and $`m_p`$ is the proton mass.333Throughout this paper the subscript 0 (2) implies values upstream (downstream) from the shock.

The magnetic field, $`B`$, is ignored in our hydrodynamic model (we implicitly assume that $`B^2/8\pi `$ is small compared to the thermal pressure), but it is an important parameter for the particle acceleration discussed next.

### 2.2 Nonlinear Diffusive Shock Acceleration

The full details of the nonlinear acceleration model used here are given in Berezhko & Ellison (1999) and Ellison et al. (2000). This is an approximate, algebraic model of diffusive shock acceleration containing the essential physics of NL acceleration, but which parameterizes important properties of the process such as the injection efficiency and the maximum energy particles achieve. While more complete models of nonlinear shock acceleration exist (e.g., Jones & Ellison, 1991; Berezhko et al., 1996; Malkov & Drury, 2001), our algebraic approximation is easier to include in global models of SNRs. It is computationally fast making it far less time consuming to do parameter searches and to compare model results with observations. In Sect. 2.5 we show by direct comparison that it gives similar results to the more physically complete model of Berezhko et al. (2002).

Briefly, the nonlinear effects in diffusive shock acceleration are: (i) the self-generation of magnetic turbulence by counter-streaming energetic particles. Back streaming particles produce turbulence in the magnetic field which leads to stronger scattering of the particles and hence to more acceleration, quickly leading to saturated turbulence levels near $`\delta B/B1`$ in strong shocks; (ii) the modification (i.e., smoothing) of the shock precursor by the backpressure of energetic particles. The precursor influences the subshock compression, $`r_{\mathrm{sub}}`$, the injection and acceleration efficiencies, and the shape of the accelerated spectrum. Since particle diffusion lengths are generally increasing functions of momentum (e.g., Blandford & Eichler, 1987; Giacalone et al., 1993), high momentum particles sample a broader portion of the flow velocity profile, and hence experience larger effective total compression ratios, $`r_{\mathrm{tot}}`$, than low momentum particles. Consequently, higher momentum particles have a flatter power-law index than those at lower momenta and can dominate the pressure in a NL fashion. The resultant superthermal distribution has a characteristic concave upward curvature until the spectrum turns over at the highest energies from losses (e.g., Ellison & Eichler, 1984; Blasi, 2002; Malkov et al., 2002); and (iii) the increase in $`r_{\mathrm{tot}}`$ from relativistic particle pressure and particle escape. As relativistic particles are produced and contribute significantly to the total pressure, their softer equation of state makes the shocked plasma more compressible ($`\gamma 4/3`$). Even more important, as the highest energy particles escape from strong shocks they drain away energy flux which must be compensated for by ramping up the overall compression ratio to conserve the fluxes. Just as in radiative shocks, this is equivalent to $`\gamma 1`$ and $`r_{\mathrm{tot}}`$ can become arbitrarily large (e.g., Kazanas & Ellison, 1986; Berezhko & Ellison, 1999; Malkov, 1997). As the overall compression increases ($`r_{\mathrm{tot}}>4`$), the subshock compression ratio, $`r_{\mathrm{sub}}`$, which is responsible for heating the gas, must become less than the test-particle (TP) value ($`r_{\mathrm{sub}}<4`$), causing the temperature of the shocked gas to drop below TP values. These changes in shock compression occur simultaneously with changes in the shape of the accelerated particle spectrum, thus linking X-ray heating to cosmic-ray production. For reviews on diffusive shock acceleration see Drury (1983); Blandford & Eichler (1987); Berezhko & Krymsky (1988); Jones & Ellison (1991); and Malkov & Drury (2001).

The most important parameters associated with nonlinear shock acceleration are the Mach numbers (i.e., the shock speed, $`u_0`$, pre-shock hydrogen number density, $`n_{\mathrm{p0}}`$, and preshock magnetic field, $`B_0`$), the injection efficiency, $`\eta _{\mathrm{inj}}`$ (i.e., the fraction of total protons which end up with superthermal energies), and the maximum proton energy produced, $`E_{\mathrm{max}}`$. As described in Berezhko & Ellison (1999), our model includes Alfvén heating in the precursor which reduces the efficiency compared to adiabatic heating and makes the magnetic field strength an important parameter. For given sonic and Alfvén Mach numbers (i.e., given $`M_S=\sqrt{\rho _0u_0^2/(\gamma P_0)}`$ and $`M_A=\sqrt{4\pi \rho _0u_0^2}/B_0`$), and a given shock size and age, $`\eta _{\mathrm{inj}}`$ sets the overall acceleration efficiency and determines the importance of NL effects. With other parameters fixed, Alfvén wave heating causes the acceleration efficiency to decrease with increasing $`B_0`$. In a complete model of diffusive shock acceleration, the injection efficiency would be determined from first principles. However, no current model of diffusive shock acceleration can do this and injection remains dependent on approximations of poorly understood wave-particle interactions.444While in principle plasma simulations, where particles move in response to Newton’s and Maxwell’s equations (particle-in-cell), can fully describe wave-particle interactions and injection, in practice, they have not yet done so because of computational limits. The limits have been insurmountable for two main reasons: (i) The simulations must be performed fully in 3D because 1- or 2-D simulations unphysically prevent cross-field diffusion (Jokipii et al., 1993; Jones et al., 1998). In all cases except strictly parallel shocks (where the upstream magnetic field is parallel to the shock normal), cross-field diffusion will be an essential part of the injection and acceleration process; and (ii) In order for NL effects to become apparent or field amplification to occur on large scales, the simulations must be run long enough in a large enough box with enough particles for a significant population of superthermal particles to be produced. If, in addition, electrons are to be investigated, the simulations must use the short electron time-step yet run for many proton time-scales increasing computation time considerably.  Here, we investigate the effects of efficient acceleration by varying $`\eta _{\mathrm{inj}}`$. A principle aim of future work is to constrain $`\eta _{\mathrm{inj}}`$ from models using X-ray and broad-band observations of particular SNRs.

In order to compare our results directly to those of Berezhko et al. (2002), we assume as they do that the magnetic field is turbulent, adopting the Bohm limit for strong particle scattering, and somewhat arbitrarily take the field downstream from the shock to be the compressed upstream magnetic field i.e., $`B_2=r_{\mathrm{tot}}B_0`$, where $`B_0(B_2)`$ is the upstream (downstream) magnetic field strength. We do not expressly consider shock obliquity, i.e., the angle between the local shock normal and $`B_0`$, even though this may be an important factor for understanding emission around the rims of some SNRs (see Reynolds, 1998, for a discussion of the effects of shock obliquity in a test-particle description of particle acceleration in SNRs). As a crude approximation, we could model the asymmetry seen in many SNRs, including SN 1006, by combining results for different quadrants of the remnant where values of the magnetic field and injection parameter were varied.

In our examples presented here, we take the unshocked ISM field to be $`B_0=20`$ $`\mu `$G to match the value determined by Berezhko et al. (2002) for SN 1006. For simplicity, unless explicitly stated we use the same constant value for the field in the unshocked supernova ejecta even though, in reality, this field is likely to weaken considerably with time because of flux conservation \[in the discussions associated with Figs. 4 (Sect. 2.3) and 8 (Sect. 2.4), we show some effects of a weak ejecta field\].

The maximum energy cosmic rays obtain depends, in part, on the scattering mean free path, $`\lambda `$, which is assumed to be,

$$\lambda =\eta _{\mathrm{mfp}}r_g,$$
(1)
where $`\eta _{\mathrm{mfp}}1`$ is taken to be a constant and $`r_g=p/(qB)`$ is the gyroradius in SI units. Small values of $`\eta _{\mathrm{mfp}}`$ imply strong scattering and allow higher maximum proton energies in a given system. The Bohm limit implies $`\eta _{\mathrm{mfp}}1`$.

The phase-space momentum distributions for protons, $`f(p)`$, are calculated as in Ellison et al. (2000) and consist of a thermal component, a three-component power law at superthermal energies, and a turnover at the highest energies given by

$$\mathrm{exp}\left[\frac{1}{\alpha }\left(\frac{p}{p_{\mathrm{max}}}\right)^\alpha \right].$$
(2)
Here, $`\alpha `$ is a constant and $`p_{\mathrm{max}}=E_{\mathrm{max}}/c`$ is determined by setting the acceleration time equal to the SNR age, $`t_{\mathrm{snr}}`$, or by setting the diffusion length of the highest energy particles equal to some fraction, $`f_{\mathrm{sk}}`$, of the shock radius, $`R_{\mathrm{sk}}`$, whichever gives the lowest $`p_{\mathrm{max}}`$ (see Baring et al., 1999). With these assumptions, $`p_{\mathrm{max}}`$ is proportional to the magnetic field strength. Reynolds (1998) has shown that when synchrotron emission is summed over space in a spherically symmetric SNR model, the resultant photon spectrum tends to be broader than that produced by a single electron distribution falling off as $`\mathrm{exp}(p/p_{\mathrm{max}})`$. This behavior can be approximated in our model by varying $`\alpha `$ and $`\alpha 0.5`$ has been used to fit X-ray and radio spectra in particular SNRs (e.g., Berezhko et al., 1999; Ellison et al., 2001). However, Berezhko et al. (2002) obtain a good fit to SN1006 with a sharp turnover in $`f(p)`$ and to match their results, we take $`\alpha =4`$.

Fig. 1 illustrates the essential differences in particle spectra between TP shock acceleration (dashed line) and NL acceleration with different injection efficiencies (solid and dotted lines). Compared to the TP power law, i.e.,

$$f(p)d^3pp^\sigma d^3p\text{with}\sigma =\frac{3r_{\mathrm{tot}}}{r_{\mathrm{tot}}1},$$
(3)
the high-energy portion ($`p>m_pc`$) of a NL spectrum is flatter, the low-energy portion of the superthermal spectrum (above the thermal peak and below $`m_pc`$) is steeper, and the thermal part is at a lower temperature. In the TP approximation, the normalization of the power-law portion of the spectrum (relative to the thermal peak) is arbitrary as long as it is low enough to contain an insignificant fraction of the total energy. If we arbitrarily set this fraction at 0.01, we can determine, for a given set of shock parameters, the injection efficiency, $`\eta _{\mathrm{TP}}`$, below which the thermal gas is decoupled from the power law, i.e., for any $`\eta _{\mathrm{inj}}\eta _{\mathrm{TP}}`$, the temperature and density of the thermal gas remains constant. For the example shown in Fig. 1, $`\eta _{\mathrm{TP}}2.5\times 10^6`$.

In contrast, the energetic portions of NL spectra are set by the conservation of mass, momentum, and energy fluxes and any change in the relativistic population influences the thermal population. When acceleration is efficient, the relation between the shock velocity, $`u_0`$, and the temperature of the shocked gas, $`T_2`$, is no longer approximated by the commonly used expression for a strong, TP shock. That is,

$$\frac{kT_2}{\mu m_pu_0^2}\frac{3}{16},$$
(4)
where $`\mu m_p`$ is the mean particle mass and $`k`$ is Boltzmann’s constant. If the equality in equation (4) is assumed for a fully ionized plasma with $`n_{\mathrm{He}}/n_{\mathrm{p0}}=0.1`$ ($`n_{\mathrm{He}}`$ is the unshocked number density of helium) and equal downstream temperatures for all species,

$$T_21.3\times 10^7\left(\frac{u_0}{10^3\mathrm{km}\mathrm{s}^1}\right)^2\mathrm{K}.$$
(5)
For young SNRs, this predicts temperatures just behind the blast wave well above those capable of explaining observed X-ray line emission (e.g., $`u_0=3000`$ km s-1 $``$ $`T_21.2\times 10^8`$ K) and suggests that the shocked electron temperature is considerably less than predicted by equation (5) (e.g., Hughes et al., 2000; Hwang et al., 2002; Long et al., 2003).

Normally, it is assumed that electrons obtain a much lower temperature than the ions in the shock layer and take some time to equilibrate downstream through Coulomb collisions (see Decourchelle & Ellison, 2001, for a discussion of electron-ion equilibration in modified shocks). However, the heating process in the shock layer is certainly dominated by collisionless, wave-particle interactions which are poorly understood. The simplest possibility, that all particles upon crossing the shock gain a speed differential $`\mathrm{\Delta }uu_0u_2`$, where $`u_2`$ is the downstream flow speed as seen from the shock, predicts that ions are heated far more than electrons.555That the situation is clearly more complicated than this is evident from the fact that if both species only receive $`\mathrm{\Delta }u`$ upon crossing the shock, energy is not conserved (see Jones & Ellison, 1987). Nevertheless, the possibility exists that electrons equilibrate faster than this simple conjecture implies and the fact that efficient acceleration produces a considerably lower proton temperature (Fig. 1) in strong shocks than the TP prediction of equation (5), suggests that it may not be as obvious that $`T_{\mathrm{e2}}/T_{\mathrm{p2}}1`$ immediately behind the shock.666Estimates of $`T_{\mathrm{e2}}/T_{\mathrm{p2}}`$ based on optical and UV line observations have been obtained just behind a few SNR shocks (e.g., Raymond et al., 1995; Laming et al., 1996; Ghavamian, 1999). These show a wide range, i.e., $`0.2<T_{\mathrm{e2}}/T_{\mathrm{p2}}<0.8`$, with large uncertainties and are generally restricted to regions that are partially neutral where the most efficient particle acceleration may not be occurring (see Drury et al., 2001, for a more complete discussion).

Furthermore, radio synchrotron observations have long confirmed that many SNRs accelerate electrons to relativistic energies, although it is not certain that these electrons are drawn from the shock heated population in all cases. What is certain is that if relatively cold upstream electrons enter the downstream region without being heated and then slowly equilibrate as they advect away downstream, very few if any will be injected into the Fermi acceleration process. In the simplest, kinematic models of injection in Fermi acceleration (e.g., Jones & Ellison, 1991; Baring et al., 1999), the efficient acceleration of thermal electrons is inconsistent with a lack of electron heating in the shock layer.

### 2.3 CR–Hydrodynamic Coupling

When efficient particle acceleration shifts energy from the thermal gas into relativistic particles, the fraction of total energy density in pressure is lowered and the evolution of the SNR is modified. The modified evolution, in turn, modifies the shock parameters which determine the acceleration. We model this coupling by including nonlinear shock acceleration in a standard 1-D hydrodynamic simulation as follows.

Briefly: (i) the hydro simulation is initiated at some time, $`t_{\mathrm{start}}<0.1t_{\mathrm{ch}}`$, after the supernova explosion with a plateau-power-law ejecta distribution and a constant interstellar medium density, $`\rho _{\mathrm{ISM}}`$; (ii) At each time step, the hydro evolves and the radii and Mach numbers of the forward and reverse shocks are determined; (iii) With this shock information, the CR acceleration is calculated using the approximate model of Ellison et al. (2000) with the injection parameter, $`\eta _{\mathrm{inj}}`$, kept constant during the simulation and having the same value at the forward and reverse shocks;777In an actual SNR, of course, the injection efficiency might vary with time, vary over the shock surface, or be different at the forward and reverse shocks (as in our model of Kepler’s SNR; Decourchelle et al., 2000). For the general results presented here, keeping $`\eta _{\mathrm{inj}}`$ constant avoids unnecessary complexity. It is important to note, however, that the acceleration efficiency, i.e., the fraction of ram kinetic energy that is transferred to superthermal particles, depends on the shock parameters (Mach numbers, shock age, etc.) as well as on $`\eta _{\mathrm{inj}}`$ and does vary with age and between the forward and reverse shocks even if $`\eta _{\mathrm{inj}}`$ is constant. (iv) The acceleration calculation provides the overall compression ratio, $`r_{\mathrm{tot}}`$, and the effective ratio of specific heats is determined by

$$\gamma _{\mathrm{eff}}=\frac{M_S^2(r_{\mathrm{tot}}+1)2r_{\mathrm{tot}}}{M_S^2(r_{\mathrm{tot}}1)},$$
(6)
where we note that $`\gamma _{\mathrm{eff}}<(P/\rho _{\mathrm{en}})+1`$ since $`\gamma _{\mathrm{eff}}`$ includes the effects from “escaping” particles.888Performing the acceleration calculation at each time step is computationally expensive. If the full spectrum is not needed and only efficient acceleration is considered, a faster approximation can use $`r_{\mathrm{tot}}=1.3M_A^{3/8}`$ if $`M_A<M_s^2`$ or $`r_{\mathrm{tot}}=1.5M_s^{3/4}`$ otherwise (see Berezhko & Ellison, 1999). Here, $`P`$ ($`\rho _{\mathrm{en}}`$) is the pressure (energy density) of the shocked gas; (v) The coupling between the CR acceleration and the hydrodynamics is accomplished by replacing the ratio of specific heats used in the hydrodynamic equations with $`\gamma _{\mathrm{eff}}`$ in the shells spanning the shock transition. We use a Lagrangian mode in the hydro simulation so the simulation grid moves with the mass. Once a mass element has been shocked and a $`\gamma _{\mathrm{eff}}`$ has been assigned to it using equation (6), that mass element retains that $`\gamma _{\mathrm{eff}}`$ for the rest of the simulation or until it reencounters a shock. Even though energy is leaving the system via the highest energy cosmic rays, we do not explicitly remove energy from the hydro. Lowering $`\gamma _{\mathrm{eff}}`$ causes the gas to be more compressible and mimics the escaping energy in one step; (vi) As the SNR evolves, the proton distribution functions in momentum, $`f(p)`$, associated with each shocked shell undergo adiabatic losses; and finally, (vii) At the end of the simulation, the $`f(p)`$’s from all of the shells are summed to determine the total contribution to the Galactic CR source population.

To start, we show in Fig. 2 a TP case where $`\gamma _{\mathrm{eff}}=5/3`$ always (top panel). For this example and all others presented in this paper, we have taken $`E_{\mathrm{sn}}=3\times 10^{51}`$ erg, $`M_{\mathrm{ej}}=1.4M_{}`$, $`B_0=20`$ $`\mu `$G, and a constant density ISM with $`n_{\mathrm{p0}}=0.3`$ cm-3, and $`f_{\mathrm{He}}=0.1`$. These values give a characteristic age, $`t_{\mathrm{ch}}=R_{\mathrm{ch}}/V_{\mathrm{ch}}210`$ yr. The initial ejecta density distribution is a power law ($`\rho r^7`$) with a flat plateau at small radii. The simulation spans $`6000`$ yr and various quantities are shown, as a function of radius, at different times during the SNR evolution. The temperature shown in panel $`(b)`$ is the average temperature defined as

$$<T_2>=\frac{P}{\rho }\frac{\mu m_p}{k}\times f_{\mathrm{th}},$$
(7)
where $`P`$ is the pressure, $`\rho `$ is the density, $`\mu =(1+4f_{\mathrm{He}})/(2+3f_{\mathrm{He}})`$, $`f_{\mathrm{He}}=0.1`$ is the ratio of helium atoms to protons, and $`f_{\mathrm{th}}`$ is the fraction of pressure in thermal particles. In this plot and all others, we assume the plasma to be fully ionized. When acceleration occurs, superthermal particles contribute to the total hydrodynamic pressure, $`P`$. We define the temperature from that fraction of total pressure in the “thermal” particles, $`f_{\mathrm{th}}`$. This fraction is determined from $`f(p)`$ (see Fig. 1 where the division between thermal and superthermal particles is indicated).

For the particular parameters assumed for the model shown in Fig. 2, the SNR is still in the self-similar phase after 100 yr as the reverse shock moves through the power-law portion of the ejecta. By 300 years, the reverse shock is moving through the plateau region of the ejecta and the forward shock has reached $`R_{\mathrm{FS}}4`$ pc. The density of the unshocked ejecta is now well below the ambient ISM density $`\rho _{\mathrm{ISM}}=1.4m_pn_{\mathrm{p0}}`$ at this stage. By 1000 years, $`R_{\mathrm{FS}}8`$ pc, the RS is now moving back toward the origin, and the energy of the explosion is shifting from being in the kinetic energy of the ejecta to thermal pressure. By $`6000`$ years, the RS has collapsed back through the ejecta to the origin and the SNR is well into the Sedov phase with a FS speed of about $`(4/3)750`$ km s-1.999Note that the speed shown in the bottom panel of Fig. 2 and panel $`(d)`$ of Fig. 4 is the flow speed in the frame of the explosion. The shock speed is given by $`u_0=(r_{\mathrm{tot}}u_{\mathrm{ds}}u_{\mathrm{up}})/(r_{\mathrm{tot}}1)`$, where $`u_{\mathrm{up}}`$ ($`u_{\mathrm{ds}}`$) is the flow speed just upstream (downstream) from the shock and $`r_{\mathrm{tot}}`$ is the shock compression ratio. For the ISM conditions we assume here, $`u_{\mathrm{up}}`$ always equals zero for the forward shock. At $`t_{\mathrm{snr}}=6000`$ yr in Fig. 2 at the FS, $`u_{\mathrm{ds}}750`$ km s-1, $`r_{\mathrm{tot}}4`$, and $`u_01000`$ km s-1.

In Fig. 3 we illustrate how particle acceleration influences the SNR shock evolution by comparing TP results ($`\eta _{\mathrm{inj}}10^6`$) against a simulation with $`\eta _{\mathrm{inj}}=10^3`$. The figure shows the forward and reverse shock radii versus remnant age, $`t_{\mathrm{snr}}`$, for the same parameters as used in Fig. 2.101010We do not follow particle acceleration at the reverse shock after it first collapses to zero radius. The TP results, i.e., those with $`\gamma _{\mathrm{eff}}=5/3`$ everywhere, are shown with solid curves and the NL results are shown with dashed curves. The dotted curve is the TP analytic approximation given by Truelove & McKee (1999) for the FS and this matches the hydro results well over the entire time period. When NL effects are important, the region between the forward and reverse shocks becomes narrower and denser, a result fully consistent with the analytic results of Decourchelle et al. (2000) and the multi-dimensional hydro simulations of Blondin & Ellison (2001).

In Fig. 4 we compare TP profiles to NL ones at $`t_{\mathrm{snr}}=200`$ and 1000 yr. In all panels, the heavy-weight solid curves are the TP results ($`\eta _{\mathrm{inj}}10^6`$), the dotted curves have $`\eta _{\mathrm{inj}}=1.2\times 10^4`$, and the dashed curves have $`\eta _{\mathrm{inj}}=10^3`$. Panels ($`a`$) and ($`e`$) show the reduction in $`\gamma _{\mathrm{eff}}`$ that occurs as shock acceleration shifts energy into relativistic particles. This produces compression ratios greater than 4 and results in denser, narrower regions between the forward and reverse shocks as shown in the density profiles ($`b`$) and ($`f`$). The separation between the forward and reverse shocks varies considerably with $`t_{\mathrm{snr}}`$ and $`\eta _{\mathrm{inj}}`$. At $`t_{\mathrm{snr}}=200`$ yr, this separation is $`0.25`$ pc for $`\eta _{\mathrm{inj}}=10^3`$ and $`0.7`$ pc for the TP case, while at $`t_{\mathrm{snr}}=1000`$ yr it is $`4.2`$ pc for $`\eta _{\mathrm{inj}}=10^3`$ and $`6`$ pc for the TP case. The shift in energy out of the thermal gas also results in a significant decrease in the temperature of the shocked gas in the NL case, as shown in panels ($`c`$) and ($`g`$).

The flow speed is shown in the bottom panel and the forward and reverse shocks are easily identified in all panels. At $`t_{\mathrm{snr}}=1000`$ yr, the flow speeds just behind the forward shocks of the three examples are about equal, i.e., $`u_{\mathrm{ds}}(\mathrm{TP})2600`$ km s-1 and $`u_{\mathrm{ds}}(\eta _{\mathrm{inj}}=10^3)2500`$ km s-1. However, since the shock speed depends on $`r_{\mathrm{tot}}`$, the TP shock is moving faster, i.e., $`u_0(\mathrm{TP})(4/3)26003500`$ km s-1, while $`u_0(\eta _{\mathrm{inj}}=10^3)(7/6)25002900`$ km s-1.

The intermediate NL case with $`\eta _{\mathrm{inj}}=1.2\times 10^4`$ (dotted curves) produces temperatures and compression ratios between the TP and $`\eta _{\mathrm{inj}}=10^3`$ cases, as expected. However, as seen in the variations of $`\gamma _{\mathrm{eff}}`$ (panels $`a`$ and $`e`$ of Fig. 4), the $`\eta _{\mathrm{inj}}=1.2\times 10^4`$ case acts quite differently at the beginning of the simulation. As described in Berezhko & Ellison (1999), unmodified shocks with $`r_{\mathrm{tot}}4`$ can occur for high Mach numbers if $`\eta _{\mathrm{inj}}`$ is small enough. In this case, the pressure in relativistic particles, $`P_{\mathrm{rel}}`$, is small compared to $`\rho _0u_0^2`$ and these particles do not slow the incoming gas enough to produce the nonlinear shock modification. For a given $`\eta _{\mathrm{inj}}`$, as the SNR ages, $`u_0`$ decreases, $`P_{\mathrm{rel}}/(\rho _0u_0^2)`$ increases, and the shock acceleration becomes more efficient and more nonlinear. Thus, for relatively inefficient injection, $`\gamma _{\mathrm{eff}}`$ drops and $`r_{\mathrm{tot}}`$ increases in the early stages of evolution, exactly the opposite behavior as with efficient injection.111111As discussed in Berezhko & Ellison (1999), it is by no means certain that actual SNR shocks will have injection rates low enough for high Mach number, unmodified shocks to occur. Furthermore, there may be large differences in injection rates for parallel and oblique regions of the same SNR blast wave, making it possible that some regions are highly modified and others unmodified.

After the initial transition from a strong, unmodified shock to a strong, modified one, the forward shock will enter the long stage where it slowly weakens, the acceleration becomes less efficient, and $`\gamma _{\mathrm{eff}}`$ increases toward 5/3. However, as shown by the dotted curve in panel $`(e)`$ in Fig. 4, the reverse shock weakens much more quickly and at $`t_{\mathrm{snr}}=1000`$ yr has $`\gamma _{\mathrm{eff}}5/3`$. Since we assume that $`B_0`$ is a constant everywhere, the RS weakens as the ejecta density drops and the magnetic pressure becomes dominant over the gas pressure. In an actual SNR, magnetic flux conservation in the expanding, unshocked ejecta would cause $`B_0`$ to weaken as the ejecta density drops and the RS would remain stronger than with a constant $`B_0`$. For a given unshocked density, a weak $`B_0`$ results in efficient shock acceleration because the Alfvén Mach number increases with decreasing $`B_0`$, and because the amount of energy transferred out of accelerated particles into waves and then into heating of the precursor gets less. Of course, if the magnetic field gets too weak, the gyroradii of accelerated particles will become comparable to the shock radius, or the acceleration time will become comparable to $`t_{\mathrm{snr}}`$, before high particle energies or high efficiencies are obtained. To illustrate the effects of a low ejecta magnetic field, we show in panels (a) and (e) of Fig. 4 $`\gamma _{\mathrm{eff}}`$ produced by our $`\eta _{\mathrm{inj}}=1.2\times 10^4`$ example with the unshocked ejecta field set at 1 $`\mu `$G (light-weight solid curves). As seen clearly in panel (e), $`\gamma _{\mathrm{eff}}`$ remains lower at the end of the simulation than in the $`(B_0)_{\mathrm{RS}}=20`$ $`\mu `$G case, indicating more efficient acceleration.

In Fig. 5 we plot the compression ratio, $`r_{\mathrm{tot}}`$, versus $`t_{\mathrm{snr}}`$ for our TP and two NL examples. There are two sets of curves because there are two ways of determining $`r_{\mathrm{tot}}`$. The solid curves are values obtained from the particle acceleration calculation which are then used in equation (6) to calculate $`\gamma _{\mathrm{eff}}`$. The dashed curves are determined directly from the hydro density ratio taken just downstream and upstream from the shock, $`\rho _{\mathrm{dw}}/\rho _{\mathrm{up}}`$. In principle, these should be identical but in practice differences occur, indicating the inherent errors in our $`\gamma _{\mathrm{eff}}`$ technique. The largest average difference shown here with $`\eta _{\mathrm{inj}}=10^3`$ is $`<10\%`$. The $`\eta _{\mathrm{inj}}=1.2\times 10^4`$ case clearly shows the rapid transition from the strong, unmodified shock with $`r_{\mathrm{tot}}4`$ at early times to the modified shock with $`r_{\mathrm{tot}}>4`$ at intermediate times. The striking difference between the $`\eta _{\mathrm{inj}}=10^3`$ and $`\eta _{\mathrm{inj}}=1.2\times 10^4`$ cases will leave an imprint on the post-shock gas that may be an important diagnostic for determining if efficient Fermi acceleration takes place in young SNRs.

In Fig. 6 we show the postshock proton temperatures for the same three examples. Comparing the heavy-weight solid and dashed curves shows that the RS temperature in the $`\eta _{\mathrm{inj}}=10^3`$ efficient shock is less than $`1/10`$ that of the test-particle shock over most of the time period shown. This is a difference large enough to have profound consequences for the interpretation of X-ray thermal emission in young SNRs. Furthermore, the temperature in the intermediate case ($`\eta _{\mathrm{inj}}=1.2\times 10^4`$) shows a stronger time variation than either the TP or $`\eta _{\mathrm{inj}}=10^3`$ case. This is reflected as a somewhat steeper spatial gradient in temperature as seen in panel $`(g)`$ of Fig. 4, again offering a possible diagnostic to test for efficient particle acceleration.

#### 2.3.1 Comparison with Analytic Model of Decourchelle et al.

In our previous work (Decourchelle et al., 2000) we coupled NL acceleration with an analytic, self-similar description of the SNR hydrodynamics (i.e., Chevalier, 1983) and a non-equilibrium ionization calculation of thermal X-ray emission. That model used the same algebraic acceleration calculation we use here with the same injection parameter, $`\eta _{\mathrm{inj}}`$. The assumptions of the two models differ mainly in the way the hydrodynamics are treated.

While the hydro simulation has a one-fluid approach with an effective ratio of specific heats derived from the compression ratio (equation 6), the analytical, self-similar solutions are based on a two-fluid approach: a nonrelativistic one with $`\gamma =5/3`$ and a relativistic one with $`\gamma =4/3`$. The NL particle acceleration calculation provides the compression ratio and the fraction of total pressure in relativistic gas at both shocks which are used to determine the boundary conditions of the self-similar solution. Imposing the compression ratios, which can be greater than 7, mimics the effects of particle loss. The limitations of the self-similar calculation are that the RS must be propagating in the power-law portion of the ejecta density profile, and that it assumes that the forward and reverse shock compression ratios, as well as the fractions of total pressure in relativistic gas, remain constant throughout the evolution. A limitation of our hydrodynamic simulation is that it assigns at the shocks to an element of gas an effective gamma, which is kept constant during its post-shock evolution.

In the top three panels of Fig. 7 we compare the density structures for our two models at $`t_{\mathrm{snr}}=100`$ yr when the RS is still propagating in the power-law portion of the ejecta density profile. We use the same supernova and ISM parameters as in our previous runs. For the TP (panel $`a`$) and $`\eta _{\mathrm{inj}}=10^3`$ (panel $`b`$) cases, the correspondence between the two models is extremely close with small differences resulting mainly from the finite grid used in the hydro simulation. The reason for this comes from the fact that self-similarity is satisfied in the TP case where $`r_{\mathrm{tot}}4`$ always. It is also approximately satisfied in the $`\eta _{\mathrm{inj}}=10^3`$ case since $`r_{\mathrm{tot}}`$, and therefore $`\gamma _{\mathrm{eff}}`$, remains fairly constant during the 100 yr evolution, as indicated by the dashed curve in panel $`(a)`$ of Fig. 4.

While it was not obvious that the $`\gamma _{\mathrm{eff}}`$ procedure used in the hydro model would be a good approximation, the excellent agreement we see with the analytic model, which only assumes self-similarity, justifies this approach, at least at these early times.

Much larger differences in the two models are seen in panels $`(c)`$ and $`(d)`$ where $`\eta _{\mathrm{inj}}=2\times 10^4`$. In this case, $`r_{\mathrm{tot}}`$ and therefore $`\gamma _{\mathrm{eff}}`$ vary strongly during 100 yr (similar to the dotted curve in panel $`a`$ of Fig. 4), and self-similar conditions no longer apply. Since the analytic model takes $`r_{\mathrm{tot}}`$ at the forward and reverse shocks at the final time and applies these values for the entire evolution, it is effectively using values, at least for the FS in the $`\eta _{\mathrm{inj}}=2\times 10^4`$ case, which are considerably greater than the average values used in the hydro calculation. Using a larger $`r_{\mathrm{tot}}`$ means there is less pressure for a given energy density pushing the shock in the self-similar model than in the hydro model. Therefore, the FS in the hydro calculation travels faster and extends further than in the analytic calculation. The variation in $`\gamma _{\mathrm{eff}}`$ is less in the shocked ejecta (again similar to the dotted curve in panel $`a`$ of Fig. 4) so the difference in the two models is less there as well.

We consider the close correspondence between these two models in the TP and $`\eta _{\mathrm{inj}}=10^3`$ cases to be a clear indication that both adequately model the SNR hydrodynamics during self-similar conditions. The validation of the $`\gamma _{\mathrm{eff}}`$ approach by our two-fluid description of the post-shock flow in the framework of self-similar, cosmic-ray modified Chevalier solutions means the more general hydro simulation can be used with some confidence when self-similar conditions do not apply. This will also allow us to use the analytic model to test the hydro model when X-ray thermal and broad-band photon emission is included.

### 2.4 Energetic Proton Spectra

The proton distribution functions that are calculated at the forward and reverse shocks during the SNR evolution can be summed to determine the total contribution to the cosmic-ray spectrum. The amount of material swept up by each shock is used to weight each spectrum and convert $`f(p)`$ to $`N(p,\mathrm{\Delta }t)`$ where $`4\pi _0^{\mathrm{}}p^2N(p,\mathrm{\Delta }t)𝑑p`$ is the total number of particles overtaken by the shock in the time interval, $`\mathrm{\Delta }t`$. The $`N(p,\mathrm{\Delta }t)`$’s are then added together, with each one adjusted for adiabatic losses during the time from when the spectrum was calculated until the end of the simulation (e.g., Reynolds, 1998), to produce $`N(p)`$.121212We do not include radiation (i.e., synchrotron) losses here since they are always negligible for ions. For electrons, of course, these losses must be considered. Typical proton spectra calculated with $`B_0=20`$ $`\mu `$G at $`t_{\mathrm{snr}}=500`$ yr ($`t_{\mathrm{snr}}/t_{\mathrm{ch}}2.4`$ for $`E_{\mathrm{sn}}=3\times 10^{51}`$ erg, $`M_{\mathrm{ej}}=1.4M_{}`$, $`n_{\mathrm{p0}}=0.3`$ cm-3, and $`f_{\mathrm{He}}=0.1`$) are shown in Fig. 8. The NL examples produce cosmic rays at the expense of heating and show lower temperatures, as indicated by the positions of the thermal peaks. The dashed curve is the result with the most efficient injection ($`\eta _{\mathrm{inj}}=10^3`$) and this spectrum is the flattest at relativistic energies (before the turnover at $`p_{\mathrm{max}}`$). The spectrum with $`\eta _{\mathrm{inj}}=1.2\times 10^4`$ is the most distorted from a Maxwellian at thermal energies due to the summing of spectra which evolve from unmodified ($`r_{\mathrm{tot}}4`$) to strongly modified ($`r_{\mathrm{tot}}>4`$), as discussed above.

To illustrate the effects of a weakening unshocked ejecta magnetic field, we show in the bottom panel of Fig. 8 (thin solid curve) $`N(p)`$ for the case where the unshocked ejecta field is set at $`(B_0)_{\mathrm{RS}}=1`$ $`\mu `$G, a factor of 20 lower than in the other cases. The unshocked ISM field is kept at 20 $`\mu `$G and $`\eta _{\mathrm{inj}}=1.2\times 10^4`$ for both the forward and reverse shocks. The flatness of the superthermal part of the spectrum compared to the dotted curve indicates that the RS is stronger (i.e., $`r_{\mathrm{tot}}`$ is larger) than when $`(B_0)_{\mathrm{RS}}=20`$ $`\mu `$G (see Fig. 4). However, $`p_{\mathrm{max}}`$ is well below that in the $`(B_0)_{\mathrm{RS}}=20`$ $`\mu `$G cases because the weak field produces large gyroradii and high momentum particles cannot be accelerated. The FS spectrum is only slightly changed by changing $`(B_0)_{\mathrm{RS}}`$ and is not shown.

In Fig. 9 we show spectra at 1000, 6000, and $`4\times 10^4`$ yr for the same ISM and SN parameters with $`\eta _{\mathrm{inj}}=10^3`$. As the FS expands and weakens, the acceleration efficiency drops and the relativistic portion of the spectrum steepens. As expected, particles accelerated at the RS decrease in importance in later stages of the SNR. As we show in Fig. 11 below, even moderately efficient injection (i.e., $`\eta _{\mathrm{inj}}10^4`$) results in approximately $`50\%`$ of the explosion energy being put into relativistic particles over the lifetime of the SNR. Another important property shown in the top panel of Fig. 9 is that, even with $`\eta _{\mathrm{inj}}=10^3`$ which produces strongly curved spectra at early times, the total particle distribution function after $`4\times 10^4`$ yr is approximately $`N(p)p^2`$ at relativistic energies. Thus the total cosmic-ray production is fairly consistent with the observed CR flux after energy dependent escape from the galaxy is included.

### 2.5 Comparison with Kinetic Models

One of the most complete models of nonlinear CR production in SNRs is that of Berezhko and co-workers (e.g., Berezhko et al., 1996; Berezhko & Völk, 1997) where the time-dependent CR transport equations are solved together with the gas-dynamic equations in spherical symmetry. This model gives the radial distribution of gas and accelerated CR spectrum at any phase of the SNR evolution and has been used with good success to model a number of young SNRs including SN1006 (Berezhko et al., 2002), Cassiopeia A (Berezhko et al., 2003), and Tycho (Völk et al., 2002). In all of the above cases, the authors conclude that satisfactory fits to the radio and X-ray synchrotron emission, along with $`\gamma `$-rays where observed, are best obtained when the forward shocks are strongly modified by the efficient acceleration of protons in a high magnetic field. While similar results have been presented before (e.g., Reynolds & Ellison, 1992; Ellison et al., 2000, 2001), Berezhko’s model probably contains the most complete physical description of the NL acceleration process at the FS, and the detailed modeling of particular SNRs clearly lends credibility to the suggestion that shocks in young SNRs can be strongly nonlinear with compression ratios greater than 4.

Our CR-hydro model is complementary to that developed by Berezhko and co-workers in that we use a more approximate model of the NL acceleration, but have a more complete model of the hydrodynamics including the reverse shock. In particular, we do not follow the particles through the acceleration process where previously accelerated particles can continue to be re-accelerated at later times. Instead, we assume that all accelerated particles originate as thermal particles in the upstream flow, are accelerated as the shock overtakes them, and then remain in the shell of material where they were accelerated, suffering adiabatic losses as the remnant ages. The re-acceleration of high energy particles with long diffusion lengths can be important as it adds to the maximum energy particles obtain (e.g. Berezhko et al., 1996), but this effect becomes less important as time goes on and new unshocked particles are overtaken and accelerated.

In Fig. 10 we compare our results with those of Berezhko et al. (2002) for the proton distribution at $`t_{\mathrm{snr}}=1000`$ yr they predict for SN1006. We have adjusted the normalization of the FS to obtain a good fit, but otherwise have used the environmental parameters given in Berezhko et al. While both our models require an arbitrary injection parameter defined as the fraction of thermal particles overtaken by the shock that become superthermal, the actual implementation of this parameter may be different in the two models. To fit SN 1006, Berezhko et al. use $`\eta =2\times 10^4`$; we have used $`\eta _{\mathrm{inj}}=1.2\times 10^4`$ to match their results. A distinct advantage of Berezhko’s model is that the maximum CR energy, $`E_{\mathrm{max}}`$, and the shape of $`f(p)`$ near $`E_{\mathrm{max}}`$ are determined self-consistently. In our CR-hydro model, we parameterize these with $`\alpha `$ in equation (2) for the shape and with the fraction of shock radius, $`f_{\mathrm{sk}}`$, set equal to the maximum particle diffusion length to determine $`E_{\mathrm{max}}`$. To match the results given by Berezhko et al. (2002), we use $`\alpha =4`$ and $`f_{\mathrm{sk}}=0.1`$. While fine tuning of our parameters could provide a more exact match than shown in Fig. 10, the particular values are unimportant compared to the fact that both models show a similar level of injection efficiency and shock modification. Considering the fundamental differences in the two models, the fit of superthermal particles accelerated by the FS shown in Fig. 10 is impressive. The concave shape, accentuated by plotting $`p^4N(p)`$, is very similar, although the piecewise nature of our approximation is apparent.

Equally important, the FS compression ratios are very similar, as shown in the top panel of Fig. 11 where the total and subshock compression ratios are compared for this SN1006 model (CR-hydro: solid curves; Berezhko et al.: dashed curves). As we discussed earlier, unmodified strong shocks can exist when $`\eta _{\mathrm{inj}}`$ is relatively low and this is just what is predicted in both our models for $`t_{\mathrm{snr}}200`$ yr in SN 1006. As the FS slows, the acceleration becomes more efficient, $`r_{\mathrm{tot}}`$ increases, and $`r_{\mathrm{sub}}`$ starts to drop below 4 at $`t_{\mathrm{snr}}200`$ yr, indicating that the shocked proton temperature will be less than predicted in the TP case. The dotted curves are CR-hydro results with various $`\eta _{\mathrm{inj}}`$’s and show (curves b and c) that the transition between unmodified and modified shocks can be extremely abrupt. For large enough $`\eta _{\mathrm{inj}}`$ (curve d) unmodified solutions do not occur and $`r_{\mathrm{tot}}`$ is large from the start.

We take the excellent correspondence achieved with Berezhko’s calculation as evidence that our CR-hydro model, with its algebraic NL particle acceleration approximation, is accurate.

As is clear in Fig. 10, Berezhko et al. (2002) do not treat thermal particles explicitly although, in principle, they could convert the thermal pressure and density into a Maxwellian distribution as we do. They also ignore the reverse shock, an approximation that is well justified if only continuum emission from relativistic particles is considered at times $`t_{\mathrm{snr}}t_{\mathrm{ch}}`$. In order to predict thermal X-ray emission in young SNRs, however, thermal electrons must be modeled and particle heating and acceleration at the forward and reverse shocks must be self-consistently determined. We do not show electron spectra here, but our NL acceleration model, with additional parameters, can produce these along with the broad-band continuum emission from radio to TeV $`\gamma `$-rays produced by them (e.g., Ellison et al., 2001). The most important advantage our CR-hydro model has, however, is the ability to model the RS and the potential for including X-ray thermal emission from the shock-heated ejecta calculated using a non-equilibrium ionization calculation (e.g., Decourchelle et al., 2000).

In the bottom panel of Fig. 11 we show the fraction of explosion energy in relativistic particles, $`E_{\mathrm{cr}}/E_{\mathrm{sn}}`$, as a function of $`t_{\mathrm{snr}}`$. As in the top panel, the solid curves are from our CR-hydro model and the dashed curve (FS only) is from Berezhko et al. (2002). The total $`E_{\mathrm{cr}}/E_{\mathrm{sn}}`$ is in good agreement at early times but diverges somewhat at later times, a possible indication that adiabatic losses are treated differently in the two models. Most importantly, both models show that approximately 50% of the explosion energy is placed into cosmic rays during the SNR lifetime.131313Berezhko et al. (2002) make the important point that injection may vary over the surface of the SNR and be significantly less where the magnetic field is highly oblique. They estimate that to supply the galactic cosmic-ray population the overall efficiency need only be $`20`$% of the maximum values shown in Fig. 11.  The dotted curves in the bottom panel show that $`E_{\mathrm{cr}}/E_{\mathrm{sn}}`$ at $`t_{\mathrm{snr}}=6000`$ yr remains in a fairly narrow range ($`40`$ to 65%) for $`\eta _{\mathrm{inj}}`$ ranging from $`8\times 10^5`$ (a), to $`10^3`$ (d).

## 3 CONCLUSIONS

The efficient production of cosmic rays by shocks in SNRs lowers the pressure to energy density ratio in the post-shock gas causing dramatic changes in the thermal properties of the shocked gas and the SNR evolution. We have presented a CR-hydro model that combines a 1-D hydrodynamic simulation of a SNR, including the forward and reverse shocks, with particle acceleration. We explicitly include the effects of particle acceleration on the shock heated ejecta, a critical step in determining how X-ray thermal emission from the hot ejecta is modified by particle acceleration.

In accord with previous results (e.g., Dorfi, 1990; Berezhko et al., 1996, 2002), we find that SNRs can easily transfer $`50`$% of the explosion energy into relativistic particles (bottom panel of Fig. 11). Compared to the situation where acceleration is absent or inefficient, this transfer results in lower FS speeds, larger overall compression ratios, cooler post-shock temperatures, and a smaller and denser interaction region between the forward and reverse shocks. At young ages, compression ratios $`r_{\mathrm{tot}}>6`$ (top panel of Fig. 11) and proton temperatures less than $`1/10`$ the TP value are predicted (Fig. 6) for the particular set of SN and ISM parameters used here to model SN 1006. For some injection efficiencies, an abrupt transition is predicted to occur at early times ($`t_{\mathrm{snr}}<t_{\mathrm{ch}}`$) between high Mach number, unmodified shocks with $`r_{\mathrm{tot}}4`$ and strongly modified shocks with $`r_{\mathrm{tot}}4`$ (top panel of Fig. 11), causing steep spatial gradients in temperature (Fig. 4). We believe the changes in shock speed, density, temperature, and spatial extent and profile are large enough to (i) provide diagnostics for X-ray observations of young SNRs sufficient to place meaningful constraints on the acceleration efficiency, and (ii) to importantly modify the inferred values of $`E_{\mathrm{sn}}`$, $`M_{\mathrm{ej}}`$, and ambient ISM density once a full thermal X-ray emission model is combined with the CR-hydro model.

Despite the strong NL effects expected for young SNRs, the CR spectrum integrated over the age of a remnant (Fig. 9) should have a spectrum not much flatter than $`N(E)E^2`$ at relativistic energies.

Our one-dimensional CR-hydro model uses a computationally fast, approximate calculation of particle acceleration which is coupled to the hydrodynamics by modifying the effective ratio of specific heats. We have verified the accuracy of this approach for parameters applicable to SN1006 by direct comparison with the more physically complete model of acceleration of Berezhko et al. (2002) (Figs. 10 and 11), where the time-dependent CR transport equations are solved self-consistently with the gas-dynamic equations. The main advantage of our model lies in the fact that we can model acceleration at the forward and reverse shocks during all stages of the SNR evolution. In a preliminary work we modeled thermal X-ray emission from Kepler’s SNR with a two-fluid, self-similar description of the SNR hydrodynamics coupled to the same calculation of particle acceleration (Decourchelle et al., 2000). We have demonstrated (Fig. 7) that the two models give essentially identical results when self-similar conditions apply. Our next step will be to include this non-equilibrium calculation of thermal X-ray emission, plus broad-band continuum emission from bremsstrahlung, synchrotron, inverse-Compton, and pion-decay (Baring et al., 1999), in the CR-hydro model.

###### Acknowledgements.

The authors would like to thank J. Blondin for providing his hydrodynamic simulation VH-1 and for help modifying it to include particle acceleration. This work was supported, in part, by a NSF-CNRS grant (NSF INT-0128883) and by a NASA ATP grant (ATP02-0042-0006).


format: text

text:
# Analytical solution of the Gross-Neveu model at finite density

## Abstract

Recent numerical calculations have shown that the ground state of the Gross-Neveu model at finite density is a crystal. Guided by these results, we can now present the analytical solution to this problem in terms of elliptic functions. The scalar potential is the superpotential of the non-relativistic Lamé Hamiltonian. This model can also serve as analytically solvable toy model for a relativistic superconductor in the Larkin-Ovchinnikov-Fulde-Ferrell phase.

In this paper we reconsider the simplest variant of the Gross-Neveu (GN) model, a 1+1 dimensional relativistic field theory with $`N`$ species of fermions interacting via a quartic self-interaction GrossNeveu

$$=\overline{\psi }^{(i)}\mathrm{i}\gamma ^\mu _\mu \psi ^{(i)}+\frac{1}{2}g^2\left(\overline{\psi }^{(i)}\psi ^{(i)}\right)^2.$$
(1)
In a previous work ThiesUrlichs , we have found that the widely accepted phase diagram of this model in the large $`N`$ limit Wolff  needed some revision. The dynamically generated scalar mean field becomes inhomogeneous in a certain region of temperature and chemical potential, a fact which had been overlooked so far. The four-fermion interaction then does not merely lead to mass generation but to formation of a kink-antikink crystal. This in turn reflects the presence of bound baryons in the GN model as can be most clearly seen in the low-density limit. The approach used in ThiesUrlichs  was a numerical implementation of the Dirac-Hartree-Fock method (equivalent to the saddle point method in the functional integral approach). Although one can in principle carry out such calculations to any desired accuracy, it remains somewhat embarrassing to have to rely on numerics when dealing with a supposedly exactly solvable model.

In the meantime we have been able to overcome this deficiency. In the present work, we focus on the $`T=0`$ case and explain how to construct the crystal ground state at any density in closed, analytical form. We expect that the finite temperature calculation can be done along similar lines but have to leave this for future work. Since the results of Ref. ThiesUrlichs  are fully confirmed by our new method we refer to this paper for more details, figures and a discussion of the underlying physics. For a more general introduction into the field of 1+1 dimensional toy models for hot and dense matter, see the review article SchoenThies . Here we concentrate on the technical details of the analytical solution.

We start from the Hartree-Fock Dirac equation,

$$\left(\gamma ^5\frac{1}{\mathrm{i}}\frac{}{x}+\gamma ^0S(x)\right)\psi (x)=\omega \psi (x),$$
(2)
choosing the $`\gamma `$-matrices as follows,

$$\gamma ^0=\sigma _1,\gamma ^1=\mathrm{i}\sigma _3,\gamma ^5=\gamma ^0\gamma ^1=\sigma _2.$$
(3)
In terms of the upper and lower spinor components $`\varphi _\pm `$ the Dirac equation consists of two coupled equations

$`\left({\displaystyle \frac{}{x}}S\right)\varphi _{}`$ $`=`$ $`\omega \varphi _+`$
$`\left({\displaystyle \frac{}{x}}+S\right)\varphi _+`$ $`=`$ $`\omega \varphi _{}.`$ (4)
which can be decoupled by squaring,

$$\left(\frac{^2}{x^2}\frac{S}{x}+S^2\right)\varphi _\pm =\omega ^2\varphi _\pm .$$
(5)
Note that Eqs. (4,5) fall precisely into the pattern of supersymmetric (SUSY) quantum mechanics. Let us now make an ansatz for $`S(x)`$ based on the superpotential of the well-known Lamé potential WhittakerWatson ; Feinberg

$$S(x)=A\kappa ^2\frac{\mathrm{sn}(Ax|\kappa ^2)\mathrm{cn}(Ax|\kappa ^2)}{\mathrm{dn}(Ax|\kappa ^2)}A\stackrel{~}{S}(Ax).$$
(6)
Here, three types of Jacobi elliptic functions with modulus $`\kappa `$ appear Jacobi . Denoting the rescaled space coordinate $`Ax`$ by $`\xi `$, Eq. (5) becomes

$$\left(\frac{^2}{\xi ^2}\frac{\stackrel{~}{S}}{\xi }+\stackrel{~}{S}^2\right)\varphi _\pm =\frac{\omega ^2}{A^2}\varphi _\pm .$$
(7)
The spatial period of $`\stackrel{~}{S}(\xi )`$ is

$$\mathrm{}=2𝐊$$
(8)
where $`𝐊`$ is the complete elliptic integral of the first kind, $`𝐊(\kappa ^2)`$ Jacobi . We shall choose the parameter $`A`$ in such a way that the original potential $`S(x)`$ has the period $`a`$ determined by the mean density ThiesUrlichs ,

$$a=\frac{1}{\rho }=\frac{\pi }{p_f},$$
(9)
hence

$$A=\frac{\mathrm{}}{a}=\frac{2p_f𝐊}{\pi }.$$
(10)
The resulting potential still has one free parameter, $`\kappa `$, which determines both its shape and its size; the period is now fixed by the mean density. Equation (7) for $`\varphi _+`$ can then be converted into

$$\left(\frac{^2}{\xi ^2}+2\kappa ^2\mathrm{sn}^2(\xi |\kappa ^2)\right)\varphi _+=\varphi _+$$
(11)
with

$$=\frac{a^2}{\mathrm{}^2}\omega ^2+\kappa ^2,$$
(12)
which is recognized as the simplest case of the Lamé equation WhittakerWatson . The corresponding equation for $`\varphi _{}`$ differs only by a translation of the potential through half a period and thus yields an identical spectrum Feinberg . Our particular ansatz for the scalar potential was of course designed in such a way as to map the Dirac equation onto a soluble Schrödinger equation with a periodic potential. We can now simply use all the well-known results for the Lamé potential. For our purpose we found Ref. Lietal  particularly useful.

In order to determine the yet unknown parameter $`\kappa `$ we shall minimize the ground state energy density. In the Hartree-Fock approach, this is usually decomposed as

$$E_{\mathrm{g}.\mathrm{s}.}=2_{p_f}^{\mathrm{\Lambda }/2}\frac{\mathrm{d}p}{2\pi }\omega +\frac{1}{2Ng^2a}_0^adxS^2(x)E_1+E_2.$$
(13)
$`E_1`$ is the sum over single particle energies over all filled negative energy states regularized by a cutoff. (For simplicity, we consider antimatter by leaving the valence band in the Dirac sea unoccupied ThiesUrlichs .) $`E_2`$ is the standard correction term for double counting of the interaction energy. Consider $`E_1`$ first, transforming the Bloch momenta $`p`$ and single particle energies $`\omega `$ to the corresponding quantities from the Lamé equation,

$$E_1=2\frac{\mathrm{}^2}{a^2}_{k_{\mathrm{min}}}^{k_{\mathrm{max}}}\frac{\mathrm{d}k}{2\pi }\sqrt{\kappa ^2}$$
(14)
where

$$k_{\mathrm{min}}=\frac{\pi }{\mathrm{}},k_{\mathrm{max}}=\frac{a\mathrm{\Lambda }}{2\mathrm{}}.$$
(15)
It is actually more convenient to integrate over $``$, using Lietal

$$\frac{\mathrm{d}k}{\mathrm{d}}=\frac{𝐄/𝐊+\kappa ^2}{2\sqrt{(1)(\kappa ^2)(1+\kappa ^2)}}.$$
(16)
Here, $`𝐄`$ is the complete elliptic integral of the second kind, $`𝐄(\kappa ^2)`$. We thus have to evaluate

$$E_1=2\frac{\mathrm{}^2}{a^2}_{_{\mathrm{min}}}^{_{\mathrm{max}}}\frac{\mathrm{d}}{2\pi }\left|\frac{\mathrm{d}k}{\mathrm{d}}\right|\sqrt{\kappa ^2}$$
(17)
where now the lower limit is the band edge,

$$_{\mathrm{min}}=1+\kappa ^2,$$
(18)
whereas the upper limit can be inferred from Eq. (16) to be

$$_{\mathrm{max}}=k_{\mathrm{max}}^2+2\left(1\frac{𝐄}{𝐊}\right)$$
(19)
with $`k_{\mathrm{max}}`$ given in Eq. (15). It is necessary to keep the sub-leading term here, since the integral over $`\mathrm{d}`$ is linearly divergent and the divergent part will be subtracted. Performing the integration in Eq. (17) yields

$`E_1`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Lambda }^2}{8\pi }}+{\displaystyle \frac{\mathrm{}^2}{4\pi 𝐊a^2}}\left(4𝐄+(\kappa ^22)𝐊\right)`$ (20)
$`+{\displaystyle \frac{\mathrm{}^2}{2\pi 𝐊a^2}}\left(2𝐄+(\kappa ^22)𝐊\right)\mathrm{ln}\left({\displaystyle \frac{a\mathrm{\Lambda }}{\mathrm{}\kappa }}\right).`$
The term $`\mathrm{\Lambda }^2`$ can be eliminated by subtracting the energy of the trivial vacuum. Now consider the double counting correction $`E_2`$, Eq. (13), in the form

$$E_2=\frac{\mathrm{}}{2Ng^2a^2}_0^{\mathrm{}}d\xi \stackrel{~}{S}^2(\xi ).$$
(21)
Inserting $`\stackrel{~}{S}`$ from Eq. (6) and transforming to the integration variable $`s=\mathrm{sn}\xi `$, the integration can be carried out as follows,

$`{\displaystyle _0^{\mathrm{}}}d\xi {\displaystyle \frac{\mathrm{sn}^2\xi \mathrm{cn}^2\xi }{\mathrm{dn}^2\xi }}`$ $`=`$ $`2{\displaystyle _0^1}ds{\displaystyle \frac{s^2\sqrt{1s^2}}{(1\kappa ^2s^2)^{3/2}}}`$ (22)
$`=`$ $`{\displaystyle \frac{2}{\kappa ^4}}\left(2𝐄+(\kappa ^22)𝐊\right).`$
The coupling constant is related to the cutoff via the (vacuum) gap equation ThiesUrlichs ; Zeitschrift  which reads (in units where the vacuum fermion mass is 1)

$$Ng^2=\frac{\pi }{\mathrm{ln}\mathrm{\Lambda }}.$$
(23)
Combining Eqs. (21)–(23), we find

$$E_2=\frac{\mathrm{}}{\pi a^2}\left(2𝐄+(\kappa ^22)𝐊\right)\mathrm{ln}\mathrm{\Lambda }.$$
(24)
Upon adding $`E_1`$ and $`E_2`$ and recalling Eq. (8), the logarithmically divergent terms cancel and we obtain the finite, renormalized ground state energy density,

$`E_{\mathrm{ren}}`$ $`=`$ $`{\displaystyle \frac{p_f^2𝐊}{\pi ^3}}\left(4𝐄+(\kappa ^22)𝐊\right)`$
$`+{\displaystyle \frac{2p_f^2𝐊}{\pi ^3}}\left(2𝐄+(\kappa ^22)𝐊\right)\mathrm{ln}\left({\displaystyle \frac{\pi }{2p_f\kappa 𝐊}}\right).`$
Let us minimize this expression with respect to the modulus $`\kappa `$, our variational parameter. This yields the simple condition

$$\kappa =\frac{a}{\mathrm{}}=\frac{\pi }{2p_f𝐊},$$
(26)
a transcendental equation for $`\kappa `$. Eliminating $`p_f`$ from $`E_{\mathrm{ren}}`$ with the help of this relation, we finally get the following parametric representation of the ground state energy as a function of density (parameter $`\kappa `$),

$`E_{\mathrm{ren}}`$ $`=`$ $`{\displaystyle \frac{1}{4\pi }}+{\displaystyle \frac{1}{\pi \kappa ^2}}\left({\displaystyle \frac{𝐄}{𝐊}}{\displaystyle \frac{1}{2}}\right),`$ (27)
$`{\displaystyle \frac{p_f}{\pi }}`$ $`=`$ $`{\displaystyle \frac{1}{2\kappa 𝐊}}.`$ (28)
We also give the scalar potential $`S(x)`$ corresponding to the optimal value of the modulus $`\kappa `$,

$$S(x)=\kappa \frac{\mathrm{sn}\left(\frac{x}{\kappa }|\kappa ^2\right)\mathrm{cn}\left(\frac{x}{\kappa }|\kappa ^2\right)}{\mathrm{dn}\left(\frac{x}{\kappa }|\kappa ^2\right)}.$$
(29)
$`S(x)`$ interpolates smoothly between widely spaced kinks and antikinks ($`\pm \mathrm{tanh}x`$) at $`\kappa 1`$ (low-density limit) and the function $`\frac{\kappa }{2}\mathrm{sin}\left(\frac{2x}{\kappa }\right)`$ for $`\kappa 0`$ (high-density limit). Unfortunately, it does not seem possible to express $`E_{\mathrm{ren}}`$ or $`S(x)`$ directly in terms of $`p_f`$ since the relation between $`p_f`$ and $`\kappa `$, Eq. (28), cannot be analytically inverted.

Let us pause here and summarize what has been achieved so far. All we have done may be viewed as a variational calculation of the ground state of baryonic matter in the GN model. Our variational ansatz amounts to generating single particle orbits from a scalar potential $`S(x)`$ and filling all occupied negative energy levels. Guided by severe restrictions from analytical solvability, we choose the one-parameter family of scalar potentials defined in Eqs. (6) and (10). The result of varying the parameter $`\kappa `$ is given in Eqs. (27)–(29). This in itself does not sound very exciting. However, if we now compare the present results with those of Ref. ThiesUrlichs , we discover that the analytical variational solution thus obtained agrees perfectly with the solution of the numerical Hartree-Fock calculation. At all densities considered in ThiesUrlichs , both the values of the ground state energy and the shape and depth of the scalar potentials are indistinguishable if plotted in one graph. In the numerical calculation, no bias about the shape of the potential was put in (except for periodicity), since the Fourier components of $`S(x)`$ were used as independent variational parameters. We therefore conclude that the true ground state happens to lie on the one-parameter trajectory of potentials which can be dealt with analytically. In order to show this without invoking any numerical results, one still has to verify that the ground state expectation value of $`\overline{\psi }\psi `$ is self-consistent, as was done for the single baryon in Zeitschrift . This is indeed possible, but requires more details about the rather involved Lamé wave functions WhittakerWatson  as well as some patience in juggling identities for elliptic functions. In order to keep this paper readable, we have therefore deferred the full analytic proof to the appendix.

Let us now make use of the closed formulae derived above to illustrate certain features of the GN crystal. If we go to the low- or high-density limit, it becomes possible to systematically resolve the transcendental equation relating $`p_f`$ and $`\kappa `$,

$`\kappa `$ $`\begin{array}{c}\\ p_f\mathrm{\hspace{0.17em}0}\end{array}`$ $`18\mathrm{e}^{\pi /p_f}+{\displaystyle \frac{32(\pi +p_f)}{p_f}}\mathrm{e}^{2\pi /p_f}`$ (32)
$`\kappa `$ $`\begin{array}{c}\\ p_f\mathrm{}\end{array}`$ $`{\displaystyle \frac{1}{p_f}}{\displaystyle \frac{1}{4p_f^3}}+{\displaystyle \frac{3}{64p_f^5}}`$ (35)
The non-analytic dependence of $`\kappa `$ on $`p_f`$ at low density is a reflection of tunnelling between the widely separated baryon wells. For the energy as a function of density, one finds

$`E_{\mathrm{ren}}`$ $`\begin{array}{c}\\ p_f\mathrm{\hspace{0.17em}0}\end{array}`$ $`{\displaystyle \frac{1}{4\pi }}+{\displaystyle \frac{2p_f}{\pi ^2}}+{\displaystyle \frac{8p_f}{\pi ^2}}\mathrm{e}^{\pi /p_f}`$ (38)
$`E_{\mathrm{ren}}`$ $`\begin{array}{c}\\ p_f\mathrm{}\end{array}`$ $`{\displaystyle \frac{p_f^2}{2\pi }}{\displaystyle \frac{1}{2^6\pi p_f^2}}+{\displaystyle \frac{3}{2^{14}\pi p_f^6}}`$ (41)
In the low-density limit, the three terms correspond to the vacuum energy density, the contribution from the baryon mass ($`\rho M_B`$ with $`M_B=2/\pi `$) and a term describing the repulsive baryon-baryon interaction. At high densities, we can identify the free massless Fermi gas piece, the leading perturbative correction already given in ThiesUrlichs  and the next term coming from higher order effects, suggesting fast convergence. In view of the comparison with ThiesUrlichs  it is also instructive to determine the Fourier coefficients $`S_n`$ of $`S(x)`$,

$$S(x)=\underset{n}{}S_n\mathrm{e}^{\mathrm{i2}\pi nx/a},(S_n=S_n^{}).$$
(42)
Upon using Eq. (16) on p. 912 of Gradshteyn-Ryzhik Gradshteyn  and correcting a misprint ($`\pi ^2`$ should read $`\pi `$ on the right hand side), we find the following closed expression (only odd $`n`$’s appear),

$$\mathrm{i}S_n=\frac{2p_f}{\mathrm{sinh}(n\pi 𝐊^{}/𝐊)},$$
(43)
where

$$𝐊^{}=𝐊(1\kappa ^2).$$
(44)
Low- and high-density limits of our previous work ThiesUrlichs  can now easily be confirmed, namely

$`\mathrm{i}S_n`$ $`\begin{array}{c}\\ p_f\mathrm{\hspace{0.17em}0}\end{array}`$ $`{\displaystyle \frac{2p_f}{\mathrm{sinh}(n\pi p_f)}}`$ (47)
$`\mathrm{i}S_1`$ $`\begin{array}{c}\\ p_f\mathrm{}\end{array}`$ $`{\displaystyle \frac{1}{4p_f}}`$ (50)
and again one finds excellent agreement with the numerical results at all densities.

Finally, we wish to point out that the GN model at finite density can also serve as a solvable model for a relativistic, inhomogeneous superconductor. Along the lines described in Ref. Thies , one can map the GN Lagrangian onto a “dual” Lagrangian which has quark-quark rather than quark-anti-quark pairing. All one has to do is redefine quarks into anti-quarks for left-handed quarks only. If one works at non-zero chemical potential, a baryonic chemical potential $`\mu `$ in the GN model corresponds to an “axial” chemical potential $`\mu _5`$ in the dual BCS-type model. Left-handed and right-handed fermions have opposite chemical potentials, hence $`\mu _5`$ in 1+1 dimensions acts like a magnetic field in 3+1 dimensions. This favors the appearance of the Larkin-Ovchinnikov-Fulde-Ferrell (LOFF) phase with spatially varying Cooper pair condensate LO ; FF . It is then natural to identify the kink-antikink crystal of the GN model with the LOFF phase of the dual BCS-type model. In this sense, the present study may also be of some use for model studies of relativistic superconductors.

Appendix

Analytical proof of self-consistency of $`S(x)`$
We would like to show that

$$S(x)=Ng^2\underset{\alpha }{\overset{\mathrm{occ}}{}}\overline{\psi }_\alpha (x)\psi _\alpha (x),$$
(51)
where the sum runs over all negative energy levels correponding to the “upper band” of the Lamé spectrum (we are again considering antimatter). First, we have to construct normalized spinor solutions of the Dirac equation (2) out of the known solutions of the Lamé equation (11). We write

$$\psi =𝒩\left(\begin{array}{c}\varphi _+\\ \varphi _{}\end{array}\right)$$
(52)
and choose for $`\varphi _+`$ a solution of the 2nd order differential equation (11). We recall that the relations between Dirac variables and Lamé variables are

$$=\kappa ^2(\omega ^2+1),x=\kappa \xi ,p=k/\kappa ,$$
(53)
where we have used condition (26). Once $`\varphi _+`$ is chosen, $`\varphi _{}`$ follows from Eq. (4),

$$\varphi _{}=\frac{1}{\kappa \omega }\left(\frac{}{\xi }+\stackrel{~}{S}(\xi )\right)\varphi _+.$$
(54)
Let us first compute the normalization factor $`𝒩`$. In a continuum normalization, the (spatially averaged) fermion density is normalized to 1 for each level,

$`1`$ $`=`$ $`{\displaystyle \frac{1}{a}}{\displaystyle _0^a}dx\psi ^{}\psi `$
$`=`$ $`|𝒩|^2{\displaystyle \frac{1}{\mathrm{}}}{\displaystyle _0^{\mathrm{}}}d\xi \left(|\varphi _+|^2+{\displaystyle \frac{1}{\kappa ^2\omega ^2}}|(_\xi +\stackrel{~}{S})\varphi _+|^2\right).`$
The 2nd term can be simplified by partial integration and use of the 2nd order wave equation, Eq. (7). Due to the Bloch theorem the boundary terms vanish. In this way, one finds

$$1=2|𝒩|^2\frac{1}{\mathrm{}}_0^{\mathrm{}}d\xi |\varphi _+|^2.$$
(56)
We now insert the solution $`\varphi _+`$ taken from the literature in terms of Jacobi functions WhittakerWatson ; Feinberg ; Lietal

$`\varphi _+`$ $`=`$ $`{\displaystyle \frac{H(\xi +\alpha )}{\mathrm{\Theta }(\xi )}}\mathrm{e}^{\xi Z(\alpha )}`$ (57)
$`=`$ $`{\displaystyle \frac{\vartheta _1(v+w,q)}{\vartheta _4(v,q)}}\mathrm{e}^{\xi Z(\alpha )}`$
with

$$v=\frac{\pi \xi }{2𝐊},w=\frac{\pi \alpha }{2𝐊},q=\mathrm{nome}(\kappa ).$$
(58)
For the upper band, $`\alpha =\mathrm{i}\eta `$. There is a 2nd solution, $`\varphi _+^{}`$, which will simply be accounted for by a factor of 2 below. For the definitions of the various Jacobi functions, see Jacobi . Using the following addition theorem for $`\vartheta _1`$ Wolfram

$$\vartheta _3^2(0)\vartheta _1(x+y)\vartheta _1(xy)=\vartheta _4^2(x)\vartheta _2^2(y)\vartheta _2^2(x)\vartheta _4^2(y)$$
(59)
together with standard relations between different Jacobi functions Jacobi , we find

$$|\varphi _+|^2=𝒜\left(1\frac{\mathrm{cn}^2(\xi |\kappa ^2)}{\mathrm{cn}^2(\alpha |\kappa ^2)}\right)$$
(60)
with

$$𝒜=\frac{\vartheta _2^2(w,q)}{\vartheta _3^2(0,q)}.$$
(61)
Now the $`\mathrm{d}\xi `$ integration in Eq. (56) can be performed (using the same variable transformation as in Eq. (22) as well as the relation $`\mathrm{}=2𝐊`$) with the result

$$1=\frac{2|𝒩|^2𝒜}{\kappa ^2\mathrm{cn}^2(\alpha |\kappa ^2)}\left(\mathrm{dn}^2(\alpha |\kappa ^2)\frac{𝐄}{𝐊}\right).$$
(62)
This determines the normalization factor $`|𝒩|^2`$. Let us now consider the scalar density. In our representation of the Dirac matrices, Eq. (3),

$`\overline{\psi }\psi `$ $`=`$ $`|𝒩|^2(\varphi _+^{}\varphi _{}+\varphi _{}^{}\varphi _+)`$ (63)
$`=`$ $`{\displaystyle \frac{|𝒩|^2}{\kappa \omega }}(_\xi +2\stackrel{~}{S})|\varphi _+|^2.`$
In the 2nd line, we have used Eq. (54). Inserting the expression Eq. (60) and performing some straightforward calculations yields

$$\overline{\psi }\psi =\frac{2𝒜|𝒩|^2}{\kappa ^3\omega }\frac{\mathrm{dn}^2(\alpha |\kappa ^2)}{\mathrm{cn}^2(\alpha |\kappa ^2)}\stackrel{~}{S}.$$
(64)
With the normalization factor $`|𝒩|^2`$ determined above and expressing $`\stackrel{~}{S}`$ through $`S(x)`$,

$$\overline{\psi }\psi =\frac{1}{\omega }\frac{\mathrm{dn}^2(\alpha |\kappa ^2)}{\mathrm{dn}^2(\alpha |\kappa ^2)𝐄/𝐊}S(x).$$
(65)
This shows that the $`x`$-dependence of $`\overline{\psi }\psi `$ is the same for each orbit, only the prefactors differ. In view of the relations \[see Lietal  and Eq. (12)\]

$`\mathrm{dn}^2(\alpha |\kappa ^2)`$ $`=`$ $`\kappa ^2,`$
$`\kappa \omega `$ $`=`$ $`\pm \sqrt{\kappa ^2},`$ (66)
we get, for negative energy states,

$$\overline{\psi }\psi =\kappa \frac{\sqrt{\kappa ^2}}{\kappa ^2𝐄/𝐊}S(x).$$
(67)
Finally we sum over all filled states. As in the calculation of the ground state energy, we convert the integration over crystal momenta into an integration over $``$, include a factor of 2 for the twofold degeneracy of the orbits and employ the integration limits Eqs. (18-19),

$`{\displaystyle \underset{\alpha }{\overset{\mathrm{occ}}{}}}\overline{\psi }_\alpha \psi _\alpha `$ $`=`$ $`{\displaystyle \frac{2}{\kappa }}{\displaystyle _{_{\mathrm{min}}}^{_{\mathrm{max}}}}{\displaystyle \frac{\mathrm{d}}{2\pi }}\left|{\displaystyle \frac{\mathrm{d}k}{\mathrm{d}}}\right|\overline{\psi }\psi `$ (68)
$`=`$ $`{\displaystyle \frac{1}{2\pi }}S(x){\displaystyle _{_{\mathrm{min}}}^{_{\mathrm{max}}}}d{\displaystyle \frac{1}{\sqrt{(1)(1+\kappa ^2)}}}`$
$`=`$ $`{\displaystyle \frac{1}{\pi }}\mathrm{ln}(\mathrm{\Lambda })S(x).`$
We have dropped terms of order $`1/\mathrm{\Lambda }^2`$ and higher, but of course no finite terms. Inserting the relation between coupling constant and cutoff from the gap equation, Eq. (23), then reproduces Eq. (51) and proves the self-consistency of the scalar potential (29).


text:
# A Black Hole Life Preserver

(August 19, 2003)

## Abstract

Since no one lives forever, all a life preserver can really do is prolong life for longer than would have otherwise been the case. With this rather limited definition in mind we explore in this paper whether in principle you can take a life preserver with you to protect you (for a while at least) against the tidal forces encountered on a trip inside a black hole.

It is well advertised that if you fall into a non-rotating Schwarzschild black hole you will be torn apart by tidal forces. As you free-fall in on a radial trajectory you will experience tidal acceleration given by $`a_z=(2M/r^3)z,a_x=(M/r^3)x,a_y=(M/r^3)y`$, where, using geometrized units, $`M`$ is the mass of the black hole, $`r`$ is the circumferential radius appearing in the Schwarzschild metric, and $`x,y,z,t`$ is a local freely falling coordinate system. If you fall in feet-first there will be a differential tidal acceleration between your head and feet of $`a_{net}=(2M/r^3)h`$, where $`h=1.8`$ meters is your height. It is like being stretched on a rack (by $`a_z`$) and simultaneously crushed in an iron maiden (by $`a_x`$ and $`a_y`$). You can stand this up to $`a_{crit}`$ = 10g = 98 m s<sup>-2</sup>. Fighter pilots withstand 9g turns for example. But beyond 10g’s, the tidal acceleration will cause pain and dismemberment. Thus, the tidal acceleration becomes torturous at $`r_c=(2Mh/a_{crit})^{1/3}`$. (You are still comfortable as you enter the event horizon of the hole only if $`M>13,800`$ solar masses, so $`r_c<2M`$.) What is not generally recognized is that the period of torture is relatively brief. Suppose you fall into the black hole from a large distance (i.e., zero kinetic energy at infinity). The proper time for you to fall from radius $`r_c`$ to $`r=0`$ is $`t_c=(4M/3)(r_c/2M)^{3/2}=(2/3)(h/a_{crit})^{1/2}=0.0904`$ sec, independent of the mass of the black hole. Can you lower $`t_c`$, so that you live longer and are tortured for less time?

You can if you take along a life preserver: a ring of mass $`M_R`$ and radius $`R`$ (see Figure 1). It even looks like a life preserver! The ring produces to first order near the origin $`(x,y,zR)`$ a tidal acceleration $`a_z=(M_R/R^3)z,a_x=(M_R/2R^3)x,a_y=(M_R/2R^3)y`$, where $`M_RM`$, $`Rr`$, and $`M_RR`$ so that the ring’s gravitational field is approximately Newtonian. If $`M_R/R^3=2M/r^3`$, the tidal acceleration by the ring exactly cancels the tidal acceleration by the black hole. If we adjust the forces in the ring so that it is just neutrally stable against collapse due to its own self-gravity, each piece of the ring will fall radially inward keeping $`Rr`$ during the infall so the tidal force of the black hole is countered all the way in. Of course, this eventually fails when $`R`$ becomes smaller than your waist or, more precisely, when $`R`$ is no longer much greater than $`h`$.

A thin ring of radius $`R`$ exerts an acceleration $`a_z`$ along the $`z`$-axis: $`a_z(z)=M_Rz/(z^2+R^2)^{3/2}`$. Consider the ring to be a torus of major radius $`R`$ with cross sectional width $`4\delta `$ and height $`4ϵ`$ where $`\delta ,ϵR`$. It may be approximated as four thin rings, each of mass $`M_R/4`$, located at $`x^2+y^2=(R\pm \delta )^2,z=\pm ϵ`$. This torus produces an acceleration in the $`z`$ direction as a function of $`z`$ along the $`z`$-axis (for $`zR`$) of approximately:

$$a_z(z)\frac{M_Rz}{R^3}\left(1\frac{3}{2}\frac{z^2}{R^2}+6\frac{\delta ^2}{R^2}\frac{9}{2}\frac{ϵ^2}{R^2}\right).$$
(1)
The first term in parentheses exactly cancels the tidal acceleration from the hole leaving a much smaller net acceleration

$$a_{z,net}(z)\frac{M_Rz}{R^3}\left(\frac{3}{2}\frac{z^2}{R^2}+6\frac{\delta ^2}{R^2}\frac{9}{2}\frac{ϵ^2}{R^2}\right).$$
(2)
How might we make the ring neutrally stable against collapse under its own self-gravity? One way would be to use $`Q=M`$ material where the electrostatic repulsion exactly cancels the gravitational attraction (Papapetrou , Majumdar , Bonner and Wickramasuriya ). Surrounding yourself with an electrically conducting, spherical-shell Faraday cage of radius $`h/2`$ would protect you against the large electric fields, while still leaving you subject to the tidal gravitational field of the ring. Since $`m_e,m_pe`$, and $`Q_R=M_Rh/2R`$, the gravitational effects of introducing the Faraday cage and excluding the electric field inside are small. And since $`h/2R`$ the electric quadrupole fields produced by the Faraday cage at the location of the ring are also fractionally small, of order $`(h/2)^5/R^5`$, and so do not affect its dynamics significantly. During free-fall collapse $`Rr,\delta r,ϵr^{1/2}`$ as each part of the ring follows a geodesic trajectory, and we require $`5\text{g}a_{z,net}(z)5\text{g}`$ for $`h/2zh/2`$ up until a time $`t_c`$ from the end. You are neutrally stable at the center to first order because the tidal force of the black hole cancels that of the ring. (For black holes larger than 7 solar masses, gradients in the strength of the hole’s tidal field across your body are small relative to the maximum values of $`a_{z,net}`$ encountered in the situation described above.)

Usefulness of $`Q=M`$ material is limited by electron-positron pair creation, whose rate is given by the relevant limit of Schwinger’s formula (Schwinger ), valid when the electric field $`E`$ varies slowly over an electron Compton wavelength:

$$\frac{dN}{dVdt}=\frac{\alpha E^2}{\pi ^2l_p^2}exp\left(\frac{\pi m_e^2}{|eE|l_p^2}\right),$$
(3)
where, using geometrized units, $`\alpha `$ is the fine structure constant, $`l_p`$ is the Planck length, $`e`$ and $`m_e`$ are the electron charge and mass, and we integrate $`dV`$ over all space. The ring has a total charge $`Q_R=eN_0=M_R`$ and a discharge timescale $`T_c=N_0/(dN/dt)`$ since positrons will be repelled by the positively charged ring and escape while electrons will be trapped and discharge the ring. We estimate $`T_c`$ by noting that, near $`t_c`$, if $`a=2ϵ=2\delta R`$, the discharge timescale for the torus should be approximately equal to the discharge timescale for an infinite cylinder of radius $`a`$ and charge per unit length of $`Q_R/2\pi R`$. Integrating Schwinger’s equation over all space then gives:

$$T_c=\left(\frac{\pi l_p^2R}{4\alpha eM_R}\right)\left[x^4\mathrm{\Gamma }(4,x)+\mathrm{\Gamma }(0,x)\right]^1,$$
(4)
where $`\mathrm{\Gamma }`$ are incomplete gamma functions and $`x=\frac{\pi ^2m_e^2Ra}{el_p^2M_R}`$. As the ring shrinks, the discharge timescale becomes exponentially shorter. Discharge is unimportant as long as $`T_ct_c`$. Using the above equations we establish an optimal solution where at the endpoint of usefulness, $`t_c=3.46\times 10^3`$ sec, $`R=28.47`$ m, $`a=2ϵ=2\delta =1.8`$ m, $`M_R=2.145\times 10^6M_{\text{earth}}`$, and $`T_c=1.3`$ sec $`t_c`$, with 50% of the pairs being created near the surface of the torus (within $`\pm 3\%`$ of the radius $`a=1.8`$ m, the Faraday cage making a negligible correction to the calculation).

Thus, the ring effectively counters black hole tides of 6760 g’s across your body, allowing you to live longer and be tortured for only 1/26$`^{\text{th}}`$ as long as before. If you were killed in 0.0904 seconds, pain signals starting at your waist (where you were being pulled apart) traveling along the fast pathways (6 - 30 m/sec) – those telling you something is happening and where – would not have time to arrive, so it would not “hurt” much. By contrast, being killed in $`3.46\times 10^3`$ sec, you really wouldn’t know what hit you. In this example the number of pairs created prior to $`t_c`$ is large, of order $`1.5\times 10^{25}`$, and the Faraday cage would have to protect you. However, if $`M_R=9.451\times 10^7`$ M$`_{\text{earth}}`$, $`a=2ϵ=2\delta =1.8`$ m, $`R=24.16`$ m, at $`t_c=4.08\times 10^3`$ sec, the pair creation rate would be lowered so that, on average, no pairs would be created prior to $`t_c`$.

Finally, we would note that instead of falling in feet first, a better posture would be to lie in fetal position with the line connecting your shoulders pointing toward the black hole, thus squeezing your body into a cylinder 18 inches high and 36 inches in diameter, taking advantage of the fact that the tidal acceleration is smaller by a factor of 2 in the $`x,y`$ directions than in the $`z`$ direction, giving $`t_c=0.0455`$ sec. A massive ring as described above (at $`t_c`$ with $`R=7.242`$ m, $`a=2ϵ=2\delta =0.4572`$ m, $`M_R=1.393\times 10^7`$ M$`_{\text{earth}}`$ allows you to survive unscathed down to $`t_c=1.75\times 10^3`$ sec.

It is interesting that, using only normal matter, we may in principle counteract tidal forces encountered in extreme situations. This might also find application in trips near neutron stars or small black holes (without falling in) where an adjustable-radius, actively-oriented life preserver might enable you to venture closer than would otherwise have been the case and still return safely home from the adventure. When entering a region of unknown tidal forces one might want to take along a large spherical shell of $`Q=M`$ material. Then, as needed, rings or small spherical masses could be drawn in from the shell to smaller radii to counteract the vacuum tidal fields encountered. A general tidal field produced by an arbitrary distribution of distant matter in the Newtonian limit could be countered by a series of rings, each ring countering each point mass in the distant mass distribution. Strong tidal forces due to a gravitational wave (moving along the $`z`$ axis) could be countered by pulling in six small spherical masses on the $`x`$ and $`y`$ axes out of phase with the gravitational wave. For the case where you fall into an unperturbed Kerr black hole metric along the rotation axis to take advantage of the axisymmetry, you would need to pull in a ring at first as in the Schwarzschild case, but would need to manually adjust its radius as a function of time, finally enlarging it to large radius and replacing it with two small equidistant spherical masses drawn in along the previous ring axis of symmetry upon entering the region where $`r<(3)^{1/2}a`$, where the tide from the black hole changes sign. As you approach the inner Cauchy horizon at $`r_\mathrm{\_}`$, however, other effects would have to be considered, such as infalling photons and gravitons and black hole evaporation. Quantum gravity effects would also have to be considered, as well as limits on the magnitude of the tidal forces you can counter as we have discussed Schwarzschild case. Still, the results presented here prompt one to ask whether masses carried along might allow you to survive longer than would otherwise have been possible.

###### Acknowledgements.

We thank Li-Xin Li and Jeremy Goodman for helpful comments. JRG acknowledges support from his NSF grant AST-9900772.