text:
# Two-dimensional turbulence of dilute polymer solutions

## Abstract

We investigate theoretically and numerically the effect of polymer additives on two-dimensional turbulence by means of a viscoelastic model. We provide compelling evidence that at vanishingly small concentrations, such that the polymers are passively transported, the probability distribution of polymer elongation has a power law tail: its slope is related to the statistics of finite-time Lyapunov exponents of the flow, in quantitative agreement with theoretical predictions. We show that at finite concentrations and sufficiently large elasticity the polymers react on the flow with manifold consequences: velocity fluctuations are drastically depleted, as observed in soap film experiments; the velocity statistics becomes strongly intermittent; the distribution of finite-time Lyapunov exponents shifts to lower values, signalling the reduction of Lagrangian chaos.

Since the discovery of the conspicuous drag reduction obtained by dissolving minute amounts of long chain molecules in a liquid, turbulence of dilute polymer solutions has attracted a lot of attention in view of its industrial applications (see, e.g. Refs. L69 ; NH95 ; SW00 ). The fluid mechanics of polymer solutions is appropriately described by viscoelastic models that are able to reproduce the rheological behavior and many other experimental observations BCAH87 . For example, it has been shown by Sureshkumar et al. that the drag reduction effect can be captured by numerical simulations of the channel flow of viscoelastic fluids SBH97 . Although the parameters used in those simulations do not match the experimental ones, the qualitative agreement is remarkable, and all the hallmarks of the turbulent flow of polymer solutions are recovered in numerical experiments.

Following this premise, it is natural to ask whether a two-dimensional viscoelastic model can reproduce the recent results by Amarouchene and Kellay AK02  showing that the turbulent flow of soap films is spectacularly affected by polymer additives (see also Refs. CBG01 ; PB99 ). Here we show that this is indeed the case, and that the suppression of large-scale velocity fluctuations observed experimentally has a simple theoretical explanation. However, the influence of polymers is not limited to the depletion of mean square velocity, which is a genuinely two-dimensional effect. In the viscoelastic case we observe a strong intermittency, with exponential tails of the velocity probability density. As for the Lagrangian statistics, we show that the values of finite-time Lyapunov exponents lower significantly upon polymer addition, which therefore reduces the chaoticity of the flow. These effects are expected to be independent of the space dimensionality, and thus relevant to three-dimensional turbulence as well.

We also investigate the limit of vanishingly small polymer concentrations, in which the polymer molecules have no influence on the advecting flow. In this case the velocity field evolves according the two-dimensional Navier-Stokes equation with friction, and is therefore smooth at scales smaller than the injection lengthscale KM80 ; NOAG00 . For passive polymers, space dimensionality plays only a minor role, and our system is an instance of a generic random smooth flow to which the theory of passive polymers developed by Chertkov C00  and Balkovsky et al. BFL00 ; BFL01  applies. We check this theory against our numerical results, and find an excellent quantitative agreement.

To describe the dynamics of a dilute polymer solution we adopt the linear viscoelastic model (Oldroyd-B)

$$_t𝒖+(𝒖\mathbf{})𝒖=\mathbf{}p+\nu \mathrm{\Delta }𝒖+\frac{2\eta \nu }{\tau }\mathbf{}𝝈\alpha 𝒖+𝒇$$
(1)
$$_t𝝈+(𝒖\mathbf{})𝝈=(\mathbf{}𝒖)^T𝝈+𝝈(\mathbf{}𝒖)2\frac{(𝝈\mathrm{𝟏})}{\tau }.$$
(2)
The velocity field $`𝒖`$ is incompressible, the symmetric matrix $`𝝈`$ is the conformation tensor of polymer molecules, and its trace $`\text{tr}𝝈`$ is a measure of their elongation note . The parameter $`\tau `$ is the (slowest) polymer relaxation time. The energy source $`𝒇`$ is a large-scale random, zero-mean, statistically homogeneous and isotropic, solenoidal vector field. The pressure term $`\mathbf{}p`$ ensures incompressibility of the velocity field. The matrix of velocity gradients is defined as $`(\mathbf{}𝒖)_{ij}=_iu_j`$ and $`\mathrm{𝟏}`$ is the unit tensor. The solvent viscosity is denoted by $`\nu `$ and $`\eta `$ is the zero-shear contribution of polymers to the total solution viscosity $`\nu _t=\nu (1+\eta )`$. The dissipative term $`\alpha 𝒖`$ models the mechanical friction between the soap film and the surrounding air RW00 , and plays a prominent role in the energy budget of Newtonian two-dimensional turbulence BCMV02 . It should be remarked that a model that describes more accurately the polymer dynamics is the FENE-P model, which accounts for the nonlinear character of polymer elasticity, culminating in a finite molecular extensibility BCAH87 . Although here we limit ourselves to the linear case because it allows a simpler theoretical treatment, our conclusions apply to the nonlinear case as well, provided that the maximal polymer elongation is very large compared to the equilibrium length.

Passive polymers. The effect of polymer concentration $`n`$ is included in Eq. (1) through the parameter $`\eta n`$. In the limit $`\eta =0`$ polymers are passively transported and stretched by Newtonian two-dimensional turbulence. The flow is driven at the largest scales and develops an enstrophy cascade towards the small scales, while the inverse energy flux is immediately halted by friction. The ensuing velocity field is therefore everywhere smooth. We briefly recall that according to Refs.C00 ; BFL00 ; BFL01  the statistics of stretched polymers in random smooth flows is expected to depend critically on the value of the Weissenberg number, here defined as $`Wi=\lambda _N\tau `$, where $`\lambda _N`$ is the maximum Lyapunov exponent of the Newtonian flow. At $`Wi<1`$ the polymer molecules spend most of the time in a coiled state, and stretch occasionally by a considerable amount.

The theory predicts a power law tail for the probability density function of $`\text{tr}𝝈`$, i.e. the square polymer elongation

$$p(\text{tr}𝝈)(\text{tr}𝝈)^{1q}\text{for}\text{tr}𝝈\text{tr}\mathrm{𝟏}.$$
(3)
The exponent $`q`$ is related to the probability of finite-time Lyapunov exponents $`P(\gamma ,t)\mathrm{exp}[tS(\gamma )]`$ via the equation $`L_{2q}=2q/\tau `$, where $`L_{2q}=\mathrm{max}_\gamma [2q\gamma S(\gamma )]`$ is the generalized Lyapunov exponent of order $`2q`$, and $`S(\gamma )`$ is the Cramér rate function (see, e.g., Ref. CPV93 ). The convexity of $`S(\gamma )`$ ensures the positivity of $`q`$ for $`Wi<1`$. Since the distribution of polymer elongations is not accessible experimentally, in order to validate the theory it is necessary to resort to numerical simulations. Eckhardt et al. in Ref. eckhardt  have given the first evidence of a power law tail for the probability distribution function of polymer elongation in three dimensional shear turbulence. As shown in Fig. 1, in our two-dimensional simulations we observe a neat power law as well. In order to check whether the observed exponent coincides with the predicted one, we have also performed direct numerical simulations of particle trajectories, and measured the probability distribution of finite-time Lyapunov exponents, thereby obtaining the expected $`q`$. The numerical result is in close agreement with theory. As the Weissenberg number exceeds unity, the probability distribution of the conformation tensor becomes unstationary and all moments $`(\text{tr}𝝈)^n`$ grow exponentially in time. This “coil-stretch” transition signals the breakdown of linear passive theory. Accounting for the nonlinear elastic modulus of polymer molecules allows to recover a stationary statistics and to develop a consistent theory of passive polymers at all Weissenberg numbers T02 . In the following we do not pursue that approach, but we rather focus on a different mechanism that limits polymer elongation: the feedback of polymers on the advecting flow.

Active polymers. When $`\eta >0`$, polymers can affect significantly the velocity dynamics, provided that they are sufficiently elongated – a condition that is met at $`Wi>1`$. This strong feedback regime is characterized in two dimensions by a suppression of large-scale velocity fluctuations (see Fig. 2), an effect first observed in soap film experiments AK02 . In Fig. 3 we present the time evolution of the total kinetic energy of the system, showing that after polymer injection a drastic depletion of kinetic energy occurs. This should be contrasted with the three-dimensional case where, on the opposite, velocity fluctuations are larger in the viscoelastic case than in the Newtonian one DCBP02 .

The suppression of velocity fluctuations by polymer additives in two-dimensional turbulence can be easily explained in the context of the randomly driven viscoelastic model (1)-(2). Indeed, the average kinetic energy balance in the statistically stationary state reads

$$F=ϵ+\frac{2\eta \nu }{\tau ^2}(\mathrm{tr}𝝈\mathrm{tr}\mathrm{𝟏})+\alpha |𝒖|^2$$
(4)
where $`ϵ=\nu |\mathbf{}𝒖|^2`$ is the viscous dissipation and $`F`$ is the average energy input, which is flow-independent for a Gaussian, $`\delta `$-correlated random forcing $`𝒇`$. To obtain eq. (4) we multiply eq. (1) by $`𝒖`$, add to it the trace of eq. (2) times $`\eta \nu /\tau `$, and average over space and time. Since in two dimensions kinetic energy flows towards large scales, it is mainly drained by friction, and viscous dissipation is vanishingly small in the limit of very large Reynolds numbers BCMV02 . Neglecting $`ϵ`$ and observing that in the Newtonian case ($`\eta =0`$) the balance (4) yields $`F=\alpha |𝒖|^2_N`$, we obtain

$$|𝒖|^2=|𝒖|^2_N\frac{2\eta \nu }{\alpha \tau ^2}(\mathrm{tr}𝝈\mathrm{tr}\mathrm{𝟏}).$$
(5)
Since, as a consequence of incompressibility and chaoticity of the flow, it can be shown from eq. (2) that $`\mathrm{tr}𝝈\mathrm{tr}\mathrm{𝟏}`$, we finally have $`|𝒖|^2|𝒖|^2_N`$, in agreement with numerical results. This simple energy balance argument can be generalized to nonlinear elastic models. As viscosity tends to zero, the average polymer elongation increases so as to compensate for the factor $`\nu `$ in eq. (5), resulting in a finite effect also in the infinite $`Re`$ limit. Since energy is essentially dissipated by linear friction, the depletion of $`|𝒖|^2`$ entails immediately the reduction of energy dissipation. The main difference between two-dimensional “friction reduction” and three-dimensional drag reduction resides in the lengthscales involved in the energy drain – large scales in 2D vs small scales in 3D.

The effect of polymer additives cannot be merely represented by a rescaling of velocity fluctuations by a given factor. In Fig. 4 we show the probability distribution of a velocity component, $`u_x`$. The choice of the $`x`$ direction is immaterial by virtue of statistical isotropy. In the Newtonian case the distribution is remarkably close to the subgaussian density $`N\mathrm{exp}(c|u_x|^3)`$ stemming from the balance between forcing and nonlinear terms in the Navier-Stokes equation, in agreement with the prediction by Falkovich and Lebedev FL97 . On the contrary, the distribution in the viscoelastic case is markedly supergaussian, with approximately exponential tails. This strong intermittency in the velocity dynamics is due to the alternation of quiescent low-velocity phases ruled by polymer feedback and bursting events where inertial nonlinearities take over.

Dilute polymers also alter significantly the distribution of finite-time Lyapunov exponents $`P(\gamma ,t)`$. In Fig. 5 the Cramér rate function $`S(\gamma )t^1\mathrm{ln}P(\gamma ,t)`$ is shown for the Newtonian and for the viscoelastic case. Since in the former situation the Lyapunov exponent $`\lambda _N`$ is greater than $`1/\tau `$, were the polymers passive all moments of elongation would grow exponentially fast. However, the feedback can damp stretching so effectively that after polymer addition $`\lambda `$ lies below $`1/\tau `$. This implies a strong reduction of Lagrangian chaos and a decreased mixing efficiency. Moreover, we find that $`L_n`$ is smaller than $`n/\tau `$ for all $`n`$, a result which guarantees the stationarity of the statistics of $`\text{tr}𝝈`$ in presence of feedback, while imposing a less restrictive condition on $`S(\gamma )`$ than the one proposed in Ref. BFL01 .

Finally, we discuss the influence of polymer concentration on the properties of the flow. As shown in Fig. 5, the Cramér functions at two very different values of $`\eta `$ are practically indistinguishable. The level of velocity fluctuations (not shown) appears to be independent of concentration as well. This property follows from the viscoelastic equations assuming that the term $`\mathrm{𝟏}/\tau `$ in eq. (2) can be neglected if polymers are substantially stretched. In that case the dynamics of the conformation tensor is invariant under rescaling by a constant factor, allowing to absorb $`\eta `$ in the definition of $`𝝈`$, and making the velocity dynamics independent of concentration. This observation poses an interesting question: can there be a concentration-dependent onset of friction or drag reduction within a linear viscoelastic model ? Since polymers are increasingly stretched as $`\eta 0`$ there are two alternatives: either this is a singular limit, and the passive case is not recovered but for $`\eta `$ strictly equal to zero, or the feedback is not uniquely ruled by polymer elongation, and there are other relevant mechanisms for polymer activity, e.g. via the creation of strong gradients of the conformation tensor. The limit of vanishingly small yet finite $`\eta `$ is very demanding at the computational level and its investigation will require further numerical work.

This work has been supported by the EU under the contract HPRN-CT-2002-00300 and by MIUR-Cofin 2001023848. Numerical simulations have been performed at CINECA (INFM parallel computing initiative).


format: text

text:
# Nonlinear self-interaction of plane gravitational waves

## Abstract

Recently Mendonça and Cardoso \[Phys. Rev. D 66, 104009 (2002)\] considered nonlinear gravitational wave packets propagating in flat space-time. They concluded that the evolution equation—to third order in amplitude—takes a similar form to what arises in nonlinear optics. Based on this equation, the authors found that nonlinear gravitational waves exhibit self-phase modulation and high harmonic generation leading to frequency up-shifting and spectral energy dilution of the gravitational wave energy. In this Brief Report we point out the fact—a possibility that seems to have been overlooked by Mendonça and Cardoso—that the nonlinear terms in the evolution equation cancels and, hence, that there is no amplitude evolution of the pulse. Finally we discuss scenarios where these nonlinearities may play a role.

As discussed by Mendonça and Cardoso in 1  the Einstein field equations are highly nonlinear and it is a most interesting question whether gravitational waves share any similarities with other types of more well-known nonlinear waves, for instance those occurring in nonlinear optics. Assuming the presence of a number of wave perturbations

$$h_{ab}(n)=ϵ_{ab}(n)a(n)\mathrm{exp}[iq_c(n)x^c]+\mathrm{c}.\mathrm{c}.,$$
(1)
it was—based on the structure of the vacuum gravitational field equations —argued by Mendonça and Cardoso that the wave amplitude evolution equations for four-wave mixing processes take the form

$$iq_a(1)\frac{}{x_a}a(1)=w(1)a^{}(2)a(3)a(4)\mathrm{exp}[i\mathrm{\Delta }q_bx^b],$$
(2)
where $`(n)`$ labels the wave perturbations, $`ϵ_{ab}(n)`$ is the unit polarization tensor for wave $`(n)`$, $`a(n)`$ is the slowly varying amplitude ($``$ stands for complex conjugate), $`q_c(n)`$ is the wave four-vector, c.c. denotes the complex conjugate of the preceding term, $`\mathrm{\Delta }q_bq_b(3)+q_b(4)q_b(1)q_b(2)`$ and the coupling coefficient $`w(1)`$ is determined by the Ricci tensor cubic in the wave amplitudes (not explicitly computed in 1 ). The background space-time was assumed flat. The authors then considered the implications of Eq. (2), which indeed is similar to what arises in nonlinear optics. The analysis was restricted to the case when all present waves propagate in the same direction. In particular they focused on (i) *self-phase modulation* (when the four waves have identical polarization and wave vector) resulting in that a gravitational wave pulse will be continuously frequency up-shifted as it propagates through space and (ii) *harmonic cascading* through which gravitational wave energy is spectrally diluted.

In this Brief Report we present the results of detailed calculations of the amplitude evolution equation for plane gravitational waves propagating along $`x^3z`$. The perturbed metric $`g_{ab}`$ is

$$g_{ab}=\eta _{ab}+h_{ab}^{(1)}+h_{ab}^{(2)}+h_{ab}^{(3)},$$
(3)
where $`\eta _{ab}=\text{diag}(1,1,1,1)`$, $`h_{ab}^{(2)}`$ and $`h_{ab}^{(3)}`$ are the second and third order nonlinear response to the linear wave perturbation $`h_{ab}^{(1)}`$. As it turns out, there are no third order response terms that couple back to the evolution equation of $`h_{ab}^{(1)}`$, therefore $`h_{ab}^{(3)}`$ is omitted from here on. The vacuum Einstein field equations reads

$$R_{ab}^{(1)}+R_{ab}^{(2)}+R_{ab}^{(3)}=0,$$
(4)
where the Ricci tensor has been divided into $`R_{ab}^{(1)}`$ (linear in $`h_{ab}^{(1)}`$), $`R_{ab}^{(2)}`$ (quadratic in $`h_{ab}^{(1)}`$ and linear in $`h_{ab}^{(2)}`$) and $`R_{ab}^{(3)}`$ (cubic in $`h_{ab}^{(1)}`$ and with terms of order $`h_{ab}^{(1)}h_{ab}^{(2)}`$). Having wave perturbations of the type given by Eq. (1) in mind we note that in $`R_{ab}^{(3)}`$, but not in $`R_{ab}^{(2)}`$, there are terms that are resonant with $`R_{ab}^{(1)}`$. Thus, we identify

$$R_{ab}^{(1)}=R_{ab}^{(3)}$$
(5)
with the amplitude evolution equation, where it should be understood that only the part of $`R_{ab}^{(3)}`$ resonant with $`R_{ab}^{(1)}`$ should be taken into account. The second order nonlinear response terms, $`h_{ab}^{(2)}`$, are related to the $`h_{ab}^{(1)}`$ by

$$R_{ab}^{(2)}=0.$$
(6)
We are free to choose a particular gauge for $`h_{ab}^{(1)}`$ and we choose the TT gauge such that the only nonzero components are $`h_{11}^{(1)}=h_{22}^{(1)}h`$ 2 . With this choice $`R_{ab}^{(2)}`$ reads

$`R_{00}^{(2)}`$ $`={\displaystyle \frac{1}{2}}\left[\lambda _t^2H\right]+{\displaystyle \frac{1}{2}}\left(_th\right)^2+h_t^2h,`$ (7a)
$`R_{03}^{(2)}`$ $`={\displaystyle \frac{1}{2}}_z_tH+{\displaystyle \frac{1}{2}}\left(_th\right)\left(_zh\right)+h_z_th,`$ (7b)
$`R_{33}^{(2)}`$ $`={\displaystyle \frac{1}{2}}\left[\lambda _z^2H\right]+{\displaystyle \frac{1}{2}}\left(_zh\right)^2+h_z^2h,`$ (7c)
and

$`R_{0\alpha }^{(2)}`$ $`={\displaystyle \frac{1}{2}}\left[_z^2h_{0\alpha }^{(2)}+_z_th_{3\alpha }^{(2)}\right],`$ (8a)
$`R_{3\alpha }^{(2)}`$ $`={\displaystyle \frac{1}{2}}\left[_t^2h_{3\alpha }^{(2)}_z_th_{0\alpha }^{(2)}\right],`$ (8b)
$`R_{\alpha \beta }^{(2)}`$ $`={\displaystyle \frac{1}{2}}\left[_t^2_z^2\right]h_{\alpha \beta }^{(2)}{\displaystyle \frac{1}{2}}\left[\left(_th\right)^2\left(_zh\right)^2\right]\delta _{\alpha \beta },`$ (8c)
where $`\alpha ,\beta =1,2`$ and

$`H`$ $`h_{11}^{(2)}+h_{22}^{(2)},`$ (9)
$`\lambda `$ $`_z^2h_{00}^{(2)}2_z_th_{03}^{(2)}+_t^2h_{33}^{(2)}.`$ (10)
When computing $`R_{ab}^{(2)}`$ and $`R_{ab}^{(3)}`$ one should make use of the approximation $`_z_t`$ dz . Applying this to $`R_{ab}^{(2)}`$ it follows that all the components of $`h_{ab}^{(2)}`$ may be set to zero except $`h_{11}^{(2)}`$ and $`h_{22}^{(2)}`$ (the final result does not depend on whether this choice, corresponding to a freedom in gauge, is made or not). Making this choice Eq. (6) becomes

$$_t^2H=\left(_th\right)^2+2_t^2(h)^2.$$
(11)
The nonzero components of $`R_{ab}^{(3)}`$ are

$`R_{11}^{(3)}`$ $`={\displaystyle \frac{1}{4}}(_zh)(_zH){\displaystyle \frac{1}{4}}(_th)(_tH),`$ (12a)
$`R_{22}^{(3)}`$ $`={\displaystyle \frac{1}{4}}(_zh)(_zH)+{\displaystyle \frac{1}{4}}(_th)(_tH),`$ (12b)
and with the approximation $`_z_t`$ we see that $`R_{ab}^{(3)}=0`$. The amplitude evolution equation (5) thus reduces to

$$(_z^2_t^2)h=0,$$
(13)
and we conclude that—up to third order in amplitude—plane gravitational waves propagate with no change in the amplitude modulation, i.e., $`w(1)`$ in Eq. (2) is zero. To some extent this may seem counterintuitive, since a gravitational wave-packet propagating in flat space-time transports energy and momentum as described by the pseudo energy-momentum tensor which is of order $`h^2`$. This produces a local curvature of space-time of the same order, and one could expect a coupling to the original wave perturbation, leading to cubic nonlinearities in the amplitude evolution equation. However, this picture is only partially true. The pseudo energy-momentum tensor does produce curvature quadratic in the amplitude $`h`$ (essentially described by Eq. (11)), and this couples back to the original perturbation leading to cubic nonlinearities in the Riemann tensor. But as follows from Eq. (12a) and Eq. (12b) these terms cancels when computing the Ricci tensor and thus does not enter the amplitude evolution equation.

It is instructive to compare this result to existing exact plane wave solutions of the vacuum Einstein field. One particular exact plane wave solution is given by the metric (see, e.g., Ref. 3 )

$$ds^2=dt^2+a(\xi )^2dx^2+b(\xi )^2dy^2+dz^2,$$
(14)
where $`a`$ and $`b`$ satisfy the equation

$$b_\xi ^2a+a_\xi ^2b=0,$$
(15)
where $`\xi =zt`$. Focusing on weak waves, we take $`a=(1+h+\frac{1}{2}H)^{\frac{1}{2}}`$ and $`b=(1h+\frac{1}{2}H)^{\frac{1}{2}}`$ where $`1h`$ and $`Hh^2`$. Expanding Eq. (15) up to third order in amplitude one find precisely the relation (11) between $`H`$ and $`h`$ and that the third order terms cancels identically. The fact that the metric components of the gravitational wave perturbations in Eq. (1) (with the gauge choices made above) coincide with those of the exact wave solution when expanded in powers of amplitude suggests that Eq. (13) will remain even if one relaxes the approximations made in this Brief Report.

It should be pointed out, however, that gravitational four-wave mixing processes are not ruled out in general. Firstly, we have here only considered parallel propagation, and it is well-known that anti-parallel gravitational waves can interact (see e.g. Ref. 4  and references therein). Secondly, the presence of background curvature and matter modifies the properties of gravitational waves, e.g., giving rise to wave dispersion and new types of nonlinearities. To what extent this increases the possibilities of nonlinear self interaction is an open question. The equations presented here provides a starting point for future investigations of these possibilities.


text:
# The Physics of Supernova Remnant Blast Waves. II. Electron-Ion Equilibration in DEM L71 in the Large Magellanic Cloud

## 1 Introduction

The outer blast wave of DEM L71 provides a valuable laboratory to address the question of electron-ion equilibration at collisionless shocks. These shocks have been traditionally explained by their Rankine-Hugoniot jump conditions (e.g., Spitzer 1978), which quantify the partition of energy into thermal and bulk kinetic forms. However, for collisionless shocks the actual heating is not straightforward because at extremely low densities, the shock front is thinner than the mean free path for collisions. Therefore, the heating of particles at the shock-front must be produced by something else, generally assumed to be plasma waves (for a review see for instance, Draine & McKee 1993; Laming 2000). Since the heating is not collisional there is no guarantee that it will produce a thermal distribution, and in fact both simulations and in-situ observations of solar-wind shocks suggest that it does not. The plasma heating mechanisms, to both thermal and cosmic-ray energies, are as yet poorly understood. These collisionless shocks are ubiquitous in astrophysical situations, but are only spatially resolved in interplanetary shocks (at lower speeds) and SNRs. The partition of energy into thermal and non-thermal populations at collisionless shocks has broad reaching implications for the dynamics of the ISM, shocks from merging galaxies, and hot gas in clusters of galaxies. Our previous work on E0102.2$``$7219 studied the partition of energy into thermal electrons versus cosmic rays (Hughes, Rakowski, & Decourchelle 2000). Here we examine the partition of energy between the two main thermal populations, electrons and protons.

In a small sub-sample of supernova remnants (SNRs), the X-ray and radio emission from the blast wave is accompanied by optical (H$`\alpha `$, Balmer) line emission. The line profiles of these Balmer-dominated spectra provide a direct measurement of the proton temperature at the shock front, thus also placing limits on the shock velocity. The optical spectra of Balmer-dominated SNRs are produced by collisional excitation when partially neutral interstellar hydrogen is overrun by the blast wave (Chevalier & Raymond 1978). Like other non-radiative shocks, post-shock cooling losses are negligible for Balmer-dominated remnants. Each Balmer emission line consists of a narrow and a broad velocity component (Chevalier, Kirshner & Raymond 1980). The narrow component is produced when cold H I, overrun by the shock, is collisionally excited by electrons and protons before being ionized. The cold neutral hydrogen can also charge exchange with fast post-shock protons, producing fast neutrals. These can then be collisionally excited, producing the broad velocity component (Chevalier et al. 1980). The Balmer emission arises in a thin ($`10^{15}`$ cm) ionization zone; therefore, the width of a broad Balmer line yields the proton temperature immediately behind the shock and hence also limits the range of possible shock velocities (Chevalier et al. 1980, Smith et al. 1991). Combining this information with the broad-to-narrow flux ratio allows us to estimate the degree of electron-ion temperature equilibration, $`g_0(T_\mathrm{e}/T_\mathrm{p})_0`$, due to plasma processes at the shock front (Tycho, SN 1006 & RCW 86; Ghavamian et al. 2001).

A more straightforward method of estimating the initial temperature equilibration in a collisionless shock is to combine the proton temperature estimated from the FWHM of the broad Balmer line with the electron temperature measured from X-ray spectra. For this purpose, the Large Magellanic Cloud remnant DEM L71 is particularly appropriate because it is completely encircled by Balmer-dominated filaments (Tuohy et al. 1982). We have utilized Fabry-Perot imaging spectroscopy to measure the width of the H$`\alpha `$ broad component along most of the rim of DEM L71 (Ghavamian et al. 2003; hereafter GRHW03). Here we report results comparing $`T_\mathrm{e}`$ from our Chandra analysis with $`T_\mathrm{p}`$ from our Fabry-Perot analysis to infer the value of $`g_0`$ for multiple apertures around the blast wave of DEM L71 .

Davies, Elliot, & Meaburn (1976) first identified DEM L71 as a SNR candidate in an optical survey of the Large Magellanic Cloud (LMC) because of its shell-like morphology that extended 83$`^{\prime \prime }`$ by 60$`^{\prime \prime }`$ in diameter. The detection of X-ray emission from DEM L71 in the Einstein  survey of the LMC by Long, Helfand and Grabelsky (1981), confirmed its nature as a SNR. DEM L71 is extremely faint in the radio band, where it has only been detected at 843 GHz with a flux of 10 mJy by Mills et al. (1984). In optical follow-up, Tuohy et al. (1982) detected filamentary H$`\alpha `$ emission but only a few faint knots of \[O III\] and \[S II\], and thus categorized it as a Balmer-dominated SNR. Smith et al. (1991) included DEM L71 in their study of 6 Balmer-dominated remnants. Using longslit spectroscopy of the brightest portions of the shock in DEM L71 they estimated shock velocities of 300$``$800 km s-1 from the H$`\alpha `$ broad component widths. In a recent ASCA study of LMC SNRs, Hughes, Hayashi & Koyama (1998) fit the X-ray spectrum of DEM L71 and showed that the data were well described by a non-equilibrium ionization (NEI) Sedov model, particularly for the case of only minimal initial equilibration between electrons and ions.

Further investigations into the nature of DEM L71 utilizing the Chandra observation are given in two companion papers. GRHW03 discusses the details of the Fabry-Perot analysis of the blast wave, and implications for the evolutionary state of the remnant. In Hughes et al. (2003) we investigated general properties of the remnant and the composition of the ejecta. DEM L71 exhibits a prominent double shock morphology with a clear spatial separation between the reverse shock (ejecta) and blast wave (ISM). The mass ($`1.5M_{}`$) and composition of the ejecta (Fe and Si-rich), strengthen the case for a Type Ia SN origin for DEM L71 .

## 2 Observations and Data Reduction

### 2.1 RFP optical observations

Optical spectra of DEM L71 were obtained using the Rutgers Fabry-Perot (RFP) imaging spectrometer as described in GRHW03. For completeness we summarize the relevant analysis here. The 14 RFP scans, centered on the wavelength of $`H\alpha `$ at the systemic velocity of the LMC, provide essentially monochromatic images of the SNR. H$`\alpha `$ line emission spectra can then be extracted from any position on the remnant. We defined 16 apertures around the blast wave, which contained approximately equal numbers of counts and avoided any stars or radiative filaments. Fits of the H$`\alpha `$ line profile determined the broad component widths and broad-to-narrow flux ratios for each aperture. For the purpose of comparison to the X-ray results, larger regions were required to obtain reasonable statistics in the Chandra X-ray spectra. Adjacent regions with consistent broad component widths and broad-to-narrow ratios were combined to provide five final optical apertures (Figure 1, Table 1). The H$`\alpha `$ line profiles from these apertures were then refitted. The proton temperature, $`T_\mathrm{p}`$ , and shock velocity, $`V_s`$, were calculated from the measured broad component width for the full range of $`g_0`$ (GRHW03). The values computed for the extrema of $`g_0=m_\mathrm{e}/m_\mathrm{p}`$ and $`g_0`$=1 are reprinted here in Table 2.

### 2.2 Chandra X-ray observations

We observed DEM L71 using the back-side-illuminated chip (S3) of the Advanced CCD Imaging Spectrometer (ACIS-S) in full-frame timed exposure mode starting on 2000 January 04 for 45.4 ks (OBSID 775). At this time, the ACIS focal plane temperature was $``$110 C. The starting point of the reduction was the events lists from the revision 2 level of the standard processing. The CCD events were corrected for charge transfer inefficiency (CTI) using software developed at Penn State University (Townsley et al. 2000). Once these corrections are applied the entire chip can be modeled with a single response function. Subsequently, the event file was screened for grade (retaining only values 0, 2-4, or 6), bad pixels, high background times and times with incorrect aspect. The final filtered events file corresponds to a livetime-corrected exposure of 34.5 ks.

Using the optical image of DEM L71 as a template, we verified that the Chandra X-ray positions are accurate to $``$0.5′′. This method was used because none of the 16 X-ray point sources detected on the S3 chip were coincident with an optical counterpart from the USNO-A2.0 star catalog.

The spectrum of the blast wave is dominated by emission below 0.8 keV. Therefore, this soft band was used to define the blast wave regions matched to each optical aperture. We adaptively smoothed the soft band image using IMSMOOTH from the Chandra Interactive Analysis of Observations (CIAO) package (with a minimum signal-to-noise ratio of 4), then chose a contour level of 0.8 counts per pixel to define the edge of the blast wave. Three polygonal regions, each four pixels wide ($`2^{\prime \prime }`$), were defined inward of the outer contour towards the local center of curvature (Figure 2). Blast wave spectra were extracted from the filtered events file using these nested regions. The three succeeding regions were constructed in order to determine not only the blast wave electron temperature but also its downstream evolution. The largest temperature gradient is expected close to the current shock front. This and the fact that shortly behind our defined regions other shock fronts and the interior (reverse shock) emission become evident, limited the distance over which we could extract blast wave spectra. In turn, this limited the number of nested regions to three because we required the extraction regions to be wider than the Chandra point-spread-function as well as wide enough to obtain sufficient statistics. Note that the length of each X-ray extraction region was defined to match the full extent of the corresponding optical apertures.

Separate ancillary response files (ARFs) for each spectral extraction region around the blast wave were generated using standard CIAO procedures. Although these were retained in the final fits, we note here that no significant differences between the ARFs were found. For all spectra, background was taken from a large annulus outside the remnant covering radii 1–2. We did not attempt to account for the recently discovered degradation in the low energy quantum efficiency of the CCD.111see the calibration area of the Chandra website for the announcement and details: http://cxc.harvard.edu/cal/Acis/Cal\_prods/qeDeg/index.html  Instead we compensate for this effect by allowing the column density in our fits to be a free parameter.

## 3 X-ray modeling

### 3.1 Theory

To measure the electron temperature at the blast wave, one has to make some assumptions about the evolution of both the temperature and ionization state downstream. These assumptions may bias the measured electron temperature and hence the final derived value $`g_0`$. To address this concern we chose to investigate three very different nonequilibrium ionization (NEI) evolutionary models for the blast wave: (1) a fully equilibrated ($`g_0=1`$), constant velocity, planar shock model, (2) an initially unequilibrated ($`g_0=m_\mathrm{e}/m_\mathrm{p}`$ ) planar shock model with post shock equilibration on Coulomb collisional timescales (Spitzer 1978), and (3) a fully equilibrated ($`g_0=1`$) Sedov model. In all three models there are only two important dynamical parameters, the current average temperature at the shock front, $`T_{ave,s}`$, and the final ionization timescale, $`n_\mathrm{e}t`$ of the innermost shock region. All other temperatures and timescales are constrained to follow the appropriate evolution. The temperatures within a parcel of gas that may evolve as a function of time are the electron temperature $`T_\mathrm{e}`$ , the proton temperature $`T_\mathrm{p}`$ , and their number-density-weighted average temperature, $`T_{ave}`$. Note that throughout this paper, ”average temperature” is used to denote the number-density-weighted average over the particles, and not a spatial or temporal average. All ions were considered to have the same temperature as the protons for the sake of the average temperature, a small effect for the metal-poor LMC interstellar medium.

Planar shocks are a reasonable approximation for (1) straight or small shock segments where the curvature is unimportant for the density profile, and (2) time periods that are short enough that the deceleration of the blast wave due to the accumulation of swept up material is negligible. The planar model has a constant velocity and hence a constant average temperature behind the shock. Initially, we chose to model the two extreme cases for the initial degree of electron-ion equilibration, $`g_0=1`$ and $`g_0=m_\mathrm{e}/m_\mathrm{p}`$ . For $`g_0=1`$, all temperatures ($`T_\mathrm{e}`$ , $`T_{ave,s}`$, and $`T_\mathrm{p}`$ ) are equal and constant throughout the shock and only the ionization timescale varies behind the shock. For $`g_0=m_\mathrm{e}/m_\mathrm{p}`$ , $`T_\mathrm{e}`$ and $`T_\mathrm{p}`$ slowly equilibrate to $`T_{ave,s}`$ via Coulomb collisions (Spitzer 1978). We use the expressions for the time variation of temperature from Itoh (1978) and Cox & Anderson (1982) to describe the variation of $`T_\mathrm{e}`$ as a function of $`n_\mathrm{e}t`$ behind the shock. Note that while this formulation was derived for Sedov evolution, it is actually equally appropriate for planar shocks. For examples of the evolution of planar shocks see section 5.

Considering the advanced age of DEM L71 ($``$4400 years, GRHW03) it was important to include a Sedov solution model. Even under full temperature equilibration, the Sedov solution implies a spatially varying electron temperature behind the shock, which could mimic the temperature variation from Coulomb equilibration in a $`g_0=m_\mathrm{e}/m_\mathrm{p}`$ planar model. The Sedov solution used here assumes a spherical geometry and uniform preshock medium. For simplicity, we only modeled the case of full electron-ion equilibration, i.e., $`T_\mathrm{e}`$ = $`T_{ave}`$ = $`T_\mathrm{p}`$ . There are two mechanisms that determine the evolution of $`T_{ave}`$ and $`n_\mathrm{e}t`$. First, the shock decelerates as it sweeps up interstellar gas, causing $`T_{ave,s}`$ to decrease with time. Second, the volume adiabatically expands as the shocked material moves out at $`(3/4)V_s`$, i.e., the interior density and temperature decrease. Parcels of gas that are shocked at different radii, $`r_i`$, and times $`t_i`$ started off with different $`T_{ave,s}`$ and $`V_s`$. Hence each parcel will follow a unique temperature and density evolution.

The beauty of the Sedov solution lies in its self-similarity. All shock quantities can be related by dimensional analysis, and the interior profiles can be expressed in terms of dimensionless reduced variables. Hence it is possible to construct a generic set of reduced temperature and density histories as a function of current fractional radius (see Figures 3 and 4). We start with the Sedov solutions for the shock quantities and the interior profiles, herein denoted as the shock and interior solutions respectively (see Sedov 1959). The shock solutions can be written to relate the current radius ($`R_s`$), temperature ($`T_{ave,s}`$), velocity ($`V_s`$), and the swept up mass ($`M_{SU}`$) to the explosion energy, the hydrogen number density of the ambient medium ($`n_0`$), and the time since the explosion ($`t`$). Following Cox & Franco (1981), we express the interior solutions in terms of the dimensionless quantity $`\beta `$ which ranges from 0.8 to 1 from the center to the edge of the blast wave and finely samples the outermost radii.222There is a minor typographical error in the Cox & Franco (1981) definition of $`\beta `$: the coefficient three tenths should be inverted to ten thirds, in order to match the solutions given by Sedov (1959). This variable is also convenient because it only depends on nondimensional combinations of the parameters in the problem. These interior solutions describe the profiles of reduced variables, the ratios of physical quantities relative to the immediate post-shock values: radius ($`r`$), density ($`x`$), pressure ($`y`$), initial radius ($`r_i`$), and the mass enclosed ($`\mu `$) within that $`\beta `$. The reduced temperature is just the ratio of the reduced pressure to the reduced density. For clarity, we note that for a given parcel of gas there are three times of importance in the formulation of the temperature and density histories: the time at which the parcel was shocked, $`t_i`$; the time at which the evolution is being evaluated, $`t`$; and the current epoch, $`t_c`$.

To derive the density and temperature histories of a parcel of gas as a function of $`t`$, we need to connect the shock and interior solutions. We make use of the fact that the mass enclosed behind the radius of a given parcel of gas is constant, i.e., that portions of the post-shock material never overtake the gas ahead of them. This allows us to match the equations for the swept up mass at the time the parcel was shocked, $`M_{SU}(t_i)/M_{SU}(t)`$ with those for the fraction of the total swept-up mass that is enclosed within a given interior radius, $`\mu (\beta )=M_{enclosed}/M_{SU}`$ to get $`\beta (t_i/t)`$. This $`\beta `$ history defines where any given parcel was within the remnant as a function of time, as well as its density and temperature relative to the shock values. We can then combine this with the shock solutions which relate the shock quantities at $`t`$ to the current shock quantities. Thus we obtain final solutions for the ratio of the evolving $`T_{\mathrm{ave}}`$ and $`n_\mathrm{e}t`$ of each parcel to the current shock values in terms of the ratio of $`t`$ to $`t_c`$. These generic temperature and ionization timescale histories can be multiplied by any chosen time, density and temperature to give specific solutions.

There are a number of important features of the Sedov solution that can be seen in the generic temperature and ionization timescale histories shown in figures 3 and 4. The three curves shown correspond to the current fractional radii of the three X-ray spectral extraction regions behind the shock in DEM L71 . They represent the evolution of parcels of gas whose final radii are at the innermost edges of these three regions. Note the steep decrease in temperature both in time and current radius. Material that now lies at $`0.85r_s`$ was shocked at about one-fifth the lifetime of the remnant when the shock was almost six times as hot. Even the current temperatures across the three spatial zones vary by more than a factor of two. This leads to a spatial gradient in the electron temperature for the Sedov model that is far larger than the variation in temperature across regions for the minimally equilibrated ($`g_0=m_\mathrm{e}/m_\mathrm{p}`$ ) planar shock model, which is at most 40% (see Table 5).

The ionization timescale of the Sedov model also has a very different evolution than in the planar case. Adiabatic expansion of the post shock gas causes both a decrease in temperature as well as a decrease in the density, reducing the rate of increase of $`n_\mathrm{e}t`$. The global maximum value of $`n_\mathrm{e}t`$ is reached for parcels that are now at $`0.9R_s`$. However, the temperature of the parcel which is currently at 0.85$`r_s`$ was more than three times hotter than the current $`T_{ave,s}`$ for almost half of its evolution. Since the ionization rates depend on the temperature, the final $`n_\mathrm{e}t`$ by itself is not a very accurate measurement of the ionization state profile behind the shock. This complicates a direct comparison of the $`n_\mathrm{e}t`$ values between the planar and Sedov models.

### 3.2 NEI implementation

Other authors have implemented non-equilibrium ionization in Sedov and planar geometries (Shull 1982; Hamilton et al. 1983; Kaastra & Jansen 1993; Hughes et al. 1998; Borkowski, Lyerly & Reynolds 2001) but for the sake of flexibility and the specific needs of the selected regions in DEM L71 , we chose to modify our own NEI code (Hughes & Singh 1994). In all cases we construct one evolving shock model that is then used to fit the three nested regions simultaneously.

For both the $`g_0=m_\mathrm{e}/m_\mathrm{p}`$ and $`g_0=1`$ planar models, we divided the post-shock region into 200 logarithmically spaced regions in $`n_\mathrm{e}t`$ . Since the average shock temperature is constant in the planar model, all parcels of gas will undergo the same electron temperature and ionization state history. At a constant shock velocity, the $`n_\mathrm{e}t`$ steps represent both the time since the material was shocked and the distance behind the shock. This allows us to calculate the emission at each timestep and sum them to find the total emission from the shock. For each step we calculated the ionization state and emission using the electron temperature at that time determined by the solution of Itoh (1978) and advanced the ionization state of the material during that $`n_\mathrm{e}t`$ interval as in Hughes & Singh (1994). The final ionization state from each step was then used as the initial ionization state for the following interval. The emissions were summed, weighting by $`\mathrm{\Delta }n_\mathrm{e}t`$, which for a planar shock is also the width perpendicular to the shock. For the purpose of fitting the nested DEM L71 blast wave spectra, each model was cut into 3 sections linearly spaced in $`n_\mathrm{e}t`$ . We assumed that the current blast wave lies at the outer edge of the outermost aperture. However, the fitted temperatures were not particularly sensitive to the exact choice of blast wave location. Finally, we included a simple spherical projection of the three modeled spatial regions to account for the contribution of the outer emission zones to the inner spectral regions. The projection in this case was done only for the three large data extraction regions and not the individual $`n_\mathrm{e}t`$ intervals, so that the different values of $`R_s`$ for each aperture could be accounted for during the fits.

For the Sedov NEI model, first the generic histories were scaled to a given current shock temperature and product of the ambient density with the remnant age. The ionization state of the material at every radius (parcel) was calculated separately using its $`T_\mathrm{e}`$ and $`n_\mathrm{e}t`$ history, then the emission was computed for the current temperature of that parcel. A grid of 100 points in $`\beta `$ (over its full allowed range of 0.8 to 1.0) was used for both the time evolution histories and the emission as a function of radius. A single set of nominal fractional radii (0.85, 0.89, 0.95, calculated for an 8.8 pc current radius) was used for all regions in DEM L71, despite their slightly different radii (see GRHW03 figures 7 and 8), so that a common grid of models could be applied to all spectra. At these nominal fractional radii spherical projection of every individual segment in $`\beta `$ was used to weigh the contribution of any given parcel to each of the nested regions behind the shock. Note that the Sedov model is inherently spherical, in contrast to the case described in the preceding paragraph where we were simply adding a spherical projection to an otherwise planar model.

## 4 X-ray spectral fits

In the preceding section, we introduced the dynamical variables, $`T_{ave,s}`$ and $`n_\mathrm{e}t`$ that we wish to measure with our X-ray fits. However, the spectra also depend on the composition of the gas and the absorbing column density to the SNR. Ideally, one would like to allow the abundance and column density parameters to be free for each region to allow for variations in the composition of the interstellar medium around the rim of the remnant, and any local absorption. Furthermore, assuming a particular abundance set may bias the temperature measurements, because they are not independent variables. For example, choosing lower abundances will often bias the fits towards higher temperatures. Unfortunately our preliminary fits revealed that the current statistics are insufficient to constrain these additional parameters.

For both planar and Sedov fits the abundances of O, Ne and Fe were at first allowed to be free, with all other species fixed at the LMC abundances of Hughes et al. (1998). The absorbing column density, $`N_\mathrm{H}`$, was also free in these preliminary fits. The abundances included in our column density parameter were kept at solar values. Since the total column densities we found were similar to the minimum absorbing gas column in the Galaxy along the line of sight to the LMC ($`56\times 10^{20}`$ atoms cm-2, Heiles & Cleary 1979), no additional LMC component was deemed necessary. From these preliminary fits we determined the average abundances and column density across all regions and all models. Most regions did not have high enough statistics to constrain the individual abundances. Furthermore, the F-test showed that freeing the abundance and $`N_\mathrm{H}`$ parameters did not significantly improve the fit over using the average values. Hence a common set of abundances for all models and regions was used in order to allow a more straightforward comparison both between different regions and for a single region under different models. In the final fits we fixed the abundances and column density to three sets: the average values and the average values plus or minus their root mean square (RMS) deviations. The final values are listed in Table 3. They are all reasonable for the interstellar medium in the LMC compared to the abundances of Hughes et al. (1998). We fitted the spectra using the RMS deviations in the column density and abundances in order to incorporate the uncertainty in the temperature and timescale due to (1) the uncertainty in our estimate of the average abundances and column and (2) possible variations in abundances and column density around the outer rim of the remnant. This estimate of the errors is quite conservative since the RMS errors are much larger than the errors in any given preliminary fit and include any variations due to the choice of model.

The uncertainty in the abundance and column density parameters prevented us from distinguishing the three shock models on the basis of their $`\chi ^2`$ values. Not only was every model the “best” choice for some region based on $`\chi ^2`$, but the variation in $`\chi ^2`$ across models was always smaller than across the RMS deviations in abundance and column density (not shown).

The nested spectra for each of the five apertures are shown in Figure 5, plotted against the $`g_0=m_\mathrm{e}/m_\mathrm{p}`$ , planar model for comparison. The best-fit parameters of the final fits for all regions and all shock models are given in Table 4. The average shock temperatures predicted by the three shock models are very different because the X-ray emission depends primarily on the current electron temperature. Roughly speaking then, the values for $`T_{ave,s}`$ are just those that give the correct $`T_\mathrm{e}`$ for the “dominant” region in any given aperture under the assumptions of that shock model. The implied average shock temperatures are useful, however, for comparison with the optically determined proton temperatures.

Under the assumption of full equilibration, $`g_0=1`$, the electron, proton and average temperatures are all equal at the shock. Hence if either the planar or the Sedov $`g_0=1`$ models are correct the optically measured $`T_\mathrm{p}`$ should match the X-ray fitted $`T_{ave,s}`$. This is the case only for region X5, the blast wave region with the slowest shock speed. For all other regions $`T_\mathrm{p}`$ is much greater than the $`T_{ave,s}`$ values for full equilibration, implying that the electron temperature at the shock is lower than the proton temperature. The $`g_0=m_\mathrm{e}/m_\mathrm{p}`$ planar fits confirm this result. Under the $`g_0=m_\mathrm{e}/m_\mathrm{p}`$ assumption, $`T_\mathrm{p}(n_\mathrm{e}+n_\mathrm{p})T_{ave}/n_\mathrm{p}`$ or $`T_\mathrm{p}(2.3/1.1)T_{ave}`$. For all regions except X5, the X-ray implied $`T_\mathrm{p}`$ is consistent with the H$`\alpha `$ proton temperature within the 1$`\sigma `$ errors. This agreement occurs on the low side of our allowed $`T_{ave,s}`$ range hinting that partial equilibration may be preferred.

All three NEI models indicate the same general result: by comparing the average shock temperatures to the H$`\alpha `$ proton temperatures, we can already exclude the case of full electron-ion equilibration in four out of five apertures, independent of the global model. In contrast, region X5, which has the slowest shock speed, is consistent with $`g_0=1`$ but not with $`g_0=m_\mathrm{e}/m_\mathrm{p}`$ under the assumption of a constant velocity shock, and is marginally consistent with $`g_0=1`$ under the Sedov model.

In addition to a comparison with the optically measured proton temperature, the parameters of the X-ray fit also imply a certain evolution of the electron temperature and ionization state. Table 5 lists the implied values of $`T_\mathrm{e}`$ and $`n_\mathrm{e}t`$ for each of the nested regions behind the shock in the three NEI models. The different model assumptions lead to contrasting spatial gradients in the electron temperatures. The ionization timescale values also differ strongly from model to model, but this can be understood in terms of their underlying assumptions. For the $`g_0=m_\mathrm{e}/m_\mathrm{p}`$ constant velocity model the fitted $`n_\mathrm{e}t`$ values are larger by factors of 10% to 100% than those for the $`g_0=1`$ constant velocity model. Early on the electron temperatures were lower and hence ionization progressed less rapidly. The fits therefore require a larger $`n_\mathrm{e}t`$ value to reach the same ionization state. For the fully equilibrated Sedov model, the final ionization timescales are much lower than the $`g_0=1`$ planar case, because earlier the material was much hotter than now, so the ionization occurred more rapidly.

Deeper observations, with sufficient statistics to constrain the electron temperature and ionization state in each of the nested regions separately, would clearly be able to discriminate between the three models given here, based on their strong differences in the implied evolutions.

## 5 Estimates of the Initial Electron-Ion Temperature Equilibration

In the previous section we showed that our independent measurements of the electron and proton temperatures led to a consistent picture of the initial temperature equilibration at each blast wave region regardless of which NEI model was used. Here we will assume a constant velocity shock and present an estimate of $`g_0`$ from the independent $`T_\mathrm{e}`$ and $`T_\mathrm{p}`$ measurements, and a direct determination of $`g_0`$ using the proton temperature to constrain our X-ray models.

In figure 6 we plot the electron temperature as a function of ionization timescale for three types of constant velocity shock models, using region X1 for illustration. The diamonds denote the electron temperatures and ionization timescales implied by the fitted values of $`T_{ave,s}`$ and $`n_\mathrm{e}t`$ for the fully equilibrated X-ray model, while the crosses are those for the $`g_0=m_\mathrm{e}/m_\mathrm{p}`$ model. The curves plot the Coulomb evolution of the electron temperature given the best-fit value of the proton temperature, and the full range of possible $`g_0`$ values. All curves eventually equilibrate (flatten out) to the average temperature. The average temperature is not the same for all $`g_0`$ values, because the initial proton temperature is kept constant: $`T_{ave,s}`$ will be lower if the electrons are not fully heated initially. Hence, the top curve represents full initial temperature equilibration, where $`T_\mathrm{e}`$, $`T_\mathrm{p}`$, and $`T_{ave}`$ are equal and remain constant, while the bottom curve represents the minimal initial heating case, $`g_0=m_\mathrm{e}/m_\mathrm{p}`$ .

From Figure 6 it is evident that the electron and proton temperatures are indeed consistent under a planar model with some degree of initial electron heating followed by Coulomb equilibration. The actual estimate of $`g_0`$ will depend on which X-ray model is used to derive the electron temperature to compare to the proton temperature, as well as which region behind the shock is chosen. However, the relative electron temperatures and the ionization timescale differences across the three regions are simply fixed by the model assumptions. For a first estimate of $`g_0`$, we chose to use the $`g_0=m_\mathrm{e}/m_\mathrm{p}`$ planar model to constrain $`T_\mathrm{e}`$ at the outermost region (i.e., shortest ionization timescale, $`n_\mathrm{e}t`$). We compared this to the electron temperature at that timescale implied by the proton temperature as a function of $`g_0`$. These estimates are listed in Table 6. The errors include those due to $`T_\mathrm{p}`$, $`T_\mathrm{e}`$, and abundances summed in quadrature, which were all of approximately equal importance. As seen before from the comparison between $`T_{ave,s}`$ and $`T_\mathrm{p}`$, only X5, the aperture with the slowest shock speed, is consistent with $`g_0=1`$. All others are consistent with $`g_0=m_\mathrm{e}/m_\mathrm{p}`$ . The large errors in both $`T_\mathrm{p}`$ and $`T_\mathrm{e}`$ when added in quadrature lead to large errors in $`g_0`$.

Given that our independent measurements have confirmed that the blast wave temperatures are consistent with a planar shock with some degree of initial electron-ion equilibration, it is now possible to constrain $`g_0`$ by incorporating our knowledge of the proton temperature into the X-ray model fits. We expect this to greatly reduce the errors in our determination of $`g_0`$ for two reasons. Firstly, $`T_\mathrm{e}`$, $`T_\mathrm{p}`$, $`n_\mathrm{e}t`$ and the abundances are all correlated variables and hence summing their errors in quadrature over-predicts the actual error. Secondly, our three nested regions give us some constraint not only on the average electron temperature but also its spatial gradient. Allowing $`g_0`$ to vary gives us access to different temperature evolutions, beyond the three extreme cases tested before.

We adapted the planar model to consider a single $`T_\mathrm{p}`$ and allow for different values of $`g_0`$ from $`m_\mathrm{e}/m_\mathrm{p}`$ to 1. For each aperture we built a grid of models for the best-fit $`T_\mathrm{p}`$ as well as its $`1\sigma `$ error limits. The only free parameters are then $`g_0`$ and $`n_\mathrm{e}t`$ . The average shock temperature and velocity are dependent on $`T_\mathrm{p}`$ and $`g_0`$ (see the above discussion of Figure 6). For any given $`g_0`$, the electron temperature varies with distance behind the shock as the downstream gas undergoes Coulomb collisions. For consistency, the abundances and column density values were fixed to the same values used in the previous fits.

The fits to the variable $`g_0`$ model are as good as the previous three. The resultant $`g_0`$ values (listed in Table 7) are consistent with the previous cruder estimates but with tighter constraints. In many regions $`n_\mathrm{e}t`$ and $`T_\mathrm{p}`$ work together to give very similar $`g_0`$ values within the $`1\sigma `$ range on $`T_\mathrm{p}`$. One aperture, X4, is found to be of intermediate equilibration, excluding both $`g_0=m_\mathrm{e}/m_\mathrm{p}`$ and $`g_0=1`$ (at the $`1\sigma `$ level).

What do these $`g_0`$ values imply for the relationship between velocity and initial equilibration in collisionless shocks? In Figure 7 we compare our DEM L71 results on $`g_0`$ as a function of $`V_s`$ with previous work on other remnants using optical observations of the H$`\alpha `$ line alone (Ghavamian et al. 2001). The results are consistent and indicate a decreasing level of equilibration with increasing shock speed. The method presented in the current paper is considerably different than the previous work modeling the H$`\alpha `$ broad-to-narrow ratios. H$`\alpha `$ emission is only produced extremely close to the shock front, so the $`g_0`$ values modeled by this method truly reflect the initial equilibration, modulo assumptions about the preshock ionization fraction. In contrast the X-ray emitting region is extended. We measure the electron temperature some distance downstream of the shock and then infer the initial equilibration, assuming that only Coulomb collisions have equilibrated the temperatures thereafter. That these two different methods lead to consistent $`g_0`$ values as a function of shock velocity validates both methods. Furthermore it suggests that there are no significant additional heating mechanisms beyond Coulomb collisions that operate downstream of the immediate shock zone. The real test, of course, will come when both methods can be applied to the same shock front. Unfortunately this was not possible here, because the broad-to-narrow ratios were anomalously low, most likely due to precursor contamination of the narrow component of the H$`\alpha `$ emission (see GRHW03 for more details).

## 6 Summary & Concluding Remarks

We have compared the post-shock electron and proton temperatures in DEM L71 to determine the initial electron-ion equilibration, $`g_0(T_\mathrm{e}/T_\mathrm{p})_0`$, over a range of shock speeds. First, we conducted independent measurements of $`T_\mathrm{e}`$ from spatially resolved Chandra ACIS-S X-ray spectra of the outer rim, using three different nonequilibrium ionization models for the post-shock temperature evolution. The first two models assumed a constant velocity shock, one with full initial equilibration, $`g_0=1`$, the other with only minimal initial equilibration, $`g_0=m_\mathrm{e}/m_\mathrm{p}`$ , followed by post-shock Coulomb collisional equilibration. The third model was a fully equilibrated ($`g_0=1`$) NEI Sedov solution. In all cases the ionization state was allowed to progress downstream. We compared the inferred average shock temperatures with the proton temperatures from the H$`\alpha `$ line profiles reported in GRHW03. The slowest shock was consistent with full equilibration in both the constant velocity and Sedov models, but not with minimal heating. In contrast all the faster shocks were consistent with $`g_0=m_\mathrm{e}/m_\mathrm{p}`$ but not $`g_0=1`$.

Our independent measurements of $`T_\mathrm{e}`$ and $`T_\mathrm{p}`$ showed that all regions were consistent with a planar shock model, where the electron and proton temperatures started at some initial degree of equilibration, $`g_0`$, and then slowly equilibrated downstream due to Coulomb collisions. Given this agreement, we constructed a planar X-ray model which incorporated the measured value of $`T_\mathrm{p}`$ and allowed for a variable $`g_0`$. With the additional constraints from the proton temperature we found one region, X4, with intermediate shock velocity, for which both full and minimal equilibration could be ruled out, albeit only at 1$`\sigma `$. Our results support the hypothesis from Ghavamian et al. (2001), that there is an anti-correlation between shock speed and initial equilibration.

The method presented here, of comparing the optically determined $`T_\mathrm{p}`$ to the X-ray measured $`T_\mathrm{e}`$, is a powerful new technique for constraining the initial temperature equilibration of SNR shock fronts. Its strength lies in the fact that neither the proton nor the electron temperature are strongly model dependent, as this work shows. On the other hand, DEM L71 was a particularly good remnant for this study. We do not expect cosmic-ray emission to contaminate our blast wave spectra because DEM L71 has the lowest radio flux of any of the 25 LMC SNRs in Mathewson et al. (1983) and there is no evidence in the X-rays for a high energy continuum from synchrotron emission (Hughes et al. 2003). Additionally, the clear separation between the ejecta and the blast wave allowed us to obtain blast wave spectra that are uncontaminated by SN ejecta. Other Balmer-dominated SNRs, such as Tycho or SN1006, may suffer more from an intermingling of ejecta into the blast wave zone, and complications to the fit from a cosmic ray synchrotron component.

All the shocks we studied in DEM L71 are in the interesting velocity range where intermediate initial equilibrations appear to occur. However, due to the large errors on our current measurements of $`g_0`$, only one region was shown to be of intermediate equilibration, although two other regions were best fit by intermediate values of $`g_0`$. In the current fits these errors were dominated by the fitted error on the electron temperature, by the root-mean-square error on the abundances and column density, and to a lesser extent by the error in the proton temperature. We have been awarded a re-observation of DEM L71 with Chandra to quadruple the total exposure time. With this observation we will constrain the abundances of all elements with prominent lines separately for each aperture and obtain tighter limits on the electron temperature. The increased statistics will also allow us to constrain the temperatures and timescales of each nested region separately. Hence we will be able to directly determine the temperature variation downstream of the shock, and verify whether it follows the evolution given by our shock models. Additional optical spectroscopy of the non-radiative shock regions identified with the Fabry-Perot would be useful to reduce the uncertainties in the proton temperatures, and better sample the broad H$`\alpha `$ component.

We are grateful for the use of Leisa Townsley’s CTI-correction algorithm and associated response matrices. We are pleased to acknowledge useful discussions with A. Sluis and T. B. Williams. Partial support was provided by NASA/Chandra grants GO0-1035X, G01-2052X, and G02-3068X; CER was supported by a NASA Graduate Student Researchers Program Fellowship.


text:
# Information Entropy and Correlations in Prime Numbers

## Abstract

The difference between two consecutive prime numbers is called the distance between the primes. We study the statistical properties of the distances and their increments (the difference between two consecutive distances) for a sequence comprising the first $`5\times 10^7`$ prime numbers. We find that the histogram of the increments follows an exponential distribution with superposed periodic behavior of period three, similar to previously-reported period six oscillations for the distances.

Recent reports indicate that many physical and biological systems exhibit patterns where prime numbers play an important role. Examples range from the periodic orbits of a system in quantum chaos to the life cycles of species BJ2 ; BJ5 ; MP1 ; CMK1 ; SRD1 ; BJ3 ; BJ4 . Recent work reports on a potential for which the quantum energy levels of a particle can be mapped onto the sequence of primes GM1 . Furthermore, it has been shown that a gas of independent bosons with energies equal to the logarithm of consecutive primes possesses a canonical partition function coinciding with the Riemann function BJ1 . The partition function of a system with energies equal to the distances between two consecutive prime numbers behaves like a set of non-interacting harmonic oscillators M1 . Most recently, power-law behavior in the distribution of primes and correlations in prime numbers have been found M2 , along with multifractal features in the distances between consecutive primes M3 . Previous work thus further motivates studies of prime numbers using methods of statistical physics. Here, we focus on the statistical properties of the distances between consecutive prime numbers and the increments in these distances \[Fig. 1\].

Since the distribution of distances is well-studied, we discuss the occurrence frequency of increments between consecutive distances. We find that the distribution of increments \[Fig. 2(a)\] exhibits large peaks for given values of the increments and medium and small peaks for other values, and that these peaks follow period-three oscillation. Specifically, we find that the increments with values of $`6k+2`$ ($`k=0,1,2,3\mathrm{}`$) have the highest occurrence frequency, followed by increments with values of $`6k+4`$. Values of $`6k`$ are relatively rare and correspond to the small peaks in the distribution. This regularity is present for both positive and negative increments and does not depend on the length $`N_p`$ of the sequence. We also find that the occurrence frequency of increments decreases exponentially and that this exponential behavior is well pronounced for both large and small peaks, forming a “double-tent” shape \[Fig. 2(b)\].

We find exponential behavior with superposed periodic behavior with period-three oscillation for the distribution of increments similar to the period-six oscillation for the distribution of distances M1 . Further, we find that the occurrence frequency of a positive increment is almost the same as the occurrence frequency of its negative counterpart for a given sequence length $`N_p`$ \[Fig. 2(c)\].

In summary, we find a new statistical feature in the sequence of increments between consecutive prime distances. We find a period-three oscillation in the distribution of increments and this distribution follows an exponential form.

This empirical observation may be of importance in further understanding the nature of prime numbers as well as those physical and biological processes where prime numbers play a role.

###### Acknowledgements.

We thank M. Wolf, S. Havlin and M. Taqqu for helpful discussions.


text:
# Outbreaks of Hantavirus induced by seasonality

## I Introduction

Hantaviruses are rodent-borne zoonotic agents that may cause diseases in humans such as hemorrhagic fever with renal syndrome and hantavirus pulmonary syndrome mills1 ; hjelle ; hooper . Hantaviruses have been identified at an increasing rate in recent years, and as of now about thirty different ones have been discovered throughout the world. One of these, the Sin Nombre virus, was not isolated until 1993 after an outbreak in the Four Corners Region in the USA mills1 ; hjelle . The host of this particularly dangerous virus is the deer mouse, Peromyscus maniculatus, the most numerous mammal in North America. The virus produces a chronic infection in the mouse population, but it is not lethal to them. It is believed that transmission in the rodent population is horizontal and due to fights, and that the subsequent infection of humans, where the mortality rate can be as high as fifty percent, is produced by their contact with the excreta of infected mice. Moreover, so far there is no vaccine or effective drug to prevent or treat the hantavirus pulmonary syndrome. Therefore, a major effort has been made to understand the population dynamics of deer mice colonies in order to design effective prevention policies mills1 .

It has been noted that environmental conditions are directly connected to outbreaks of Hanta hjelle . For instance, the 1993 and 1998 outbreaks occurring in the Four Corners Region have been associated with the so-called $`El`$ $`Ni\stackrel{~}{n}o`$ Southern Oscillation hjelle . Related to this and other such observations, the effects of seasonality in ecological systems have been a subject of recent interest selgrade ; king . Multi-year oscillations of mammal populations schaffer , prey-predator seasonal dynamics lima , and persistence of parasites in plants between seasons gubbins  are examples that illustrate the importance of seasonality in population dynamics.

Recently, Abramson and Kenkre have proposed a phenomenological model for mice population that successfully reproduces some features of Hantavirus propagation abramson . In particular, that model explores the relation between resources in the medium, carrying capacity, and the spread of the infection in the rodent colony. Herein we study the effects of seasonality in that model. Our motivation is not only to provide more realism to the model, but also to investigate the counterintuitive effects that dynamic alternation may cause in a biological system. Brownian motors astumian  and switching-induced morphogenesis buceta  are examples that show that alternation in time of “uninteresting” dynamics may produce “interesting” outcomes. Along these lines, we will show that alternation of seasons, neither of which by itself fulfills the environmental requirements on the carrying capacity for spreading of the infection, may produce an outbreak of the disease. The mechanism driving this behavior is the interruption of the relaxation process that equilibrates the mouse population from season to season: if the duration of seasons becomes shorter than the relaxation time of the population, the disease spreads.

The paper is organized as follows. In Sec. II we briefly review the model for mouse populations introduced in abramson . In Sec. III we explain how seasonality is introduced in that model and analyze the conditions for Hanta outbreaks to take place due to the alternation of seasons. The exact solution of the model and a particular example that illustrates the phenomenology is given in Sec. IV. Finally, in Sec.V, we summarize the main conclusions and propose some directions for future work.

## II The Model

The model introduced in abramson  for the temporal evolution of a population of mice subjected to the Hantavirus infection reads:

$`{\displaystyle \frac{dM_S}{dt}}`$ $`=`$ $`bMcM_S{\displaystyle \frac{M_SM}{K}}aM_SM_I,`$ (1a)
$`{\displaystyle \frac{dM_I}{dt}}`$ $`=`$ $`cM_I{\displaystyle \frac{M_IM}{K}}+aM_SM_I,`$ (1b)
where $`M_S`$ and $`M_I`$ are the population densities of susceptible and infected mice respectively, $`M=M_S+M_I`$ is the total population of mice, and $`a`$, $`b`$, $`c`$ and $`K`$ are positive constants. The terms on the right hand sides of Eqs. (1a) and (1b) take into account the following processes: births with rate constant $`b`$, depletion by death with rate constant $`c`$, competition for the resources in the medium characterized by the carrying capacity $`K`$, and transmission of the infection with rate coefficient $`a`$. It is worth noting that infected pregnant mice produce Hanta antibodies that keep their fetus free from the infection; that is, all mice are born susceptible mills1 , as indicated by the absence of a birth term in Eq. (1b). Note also the absence of a recovery term in the model since, as mentioned earlier, mice become chronically infected by the virus.

The system of equations (1a) and (1b) has four equilibrium points. Two of them are irrelevant for the analysis: the null state $`M_I=M_S=0`$, which is always unstable if $`b>c`$ (a condition that we will assume throughout this paper), and a meaningless state with $`M_I<0`$ for any value of the parameters. The other two equilibria are

$$M_S=K(bc)M_I=0,$$
(2)
$$M_S=\frac{b}{a}M_I=K(bc)\frac{b}{a},$$
(3)
The stability of the equilibrium points (2) and (3) depends on the value of the carrying capacity. If $`K<K_c=\frac{1}{a}\left(\frac{b}{bc}\right)`$ then (2) is stable and (3) unstable. If $`K>K_c`$ it is the other way around. That is, when the available resources, $`K`$, are below the critical value, $`K_c`$, the infection does not propagate in the colony, the whole population of mice grows healthy, and its size increases proportionally with those resources. As soon as $`K`$ surpasses $`K_c`$ the virus spreads in the colony, the susceptible mouse population saturates, and the fraction of infected mice becomes larger as $`K`$ increases (see Fig. 1).

## III Seasonal Alternation

The Four Corners Region, where an important number of cases of Hantavirus Pulmonary Syndrome have occurred, has a desert climate. The largest climate variations within this region come from the alternation between dry and rainy seasons. We will therefore assume alternation in time of these two seasons. It is important to remark that a two-season assumption is not crucial, and that the analysis with four seasons is also straightforward within the formalism introduced herein. During each of the two seasons under consideration we assume there to be no climate variations, so that each season can be characterized by a set of time-independent parameters $`\{\rho _i\}=\{a_i,b_i,c_i,1/K_i\}`$ where $`i=1,2`$. We implement square-periodic season alternation where the duration of each season is $`T/2`$. Again, this particular alternation pattern is not essential for the mechanism. Other schemes of season alternation, e.g. different duration of the seasons or random switching between seasons, do not qualitatively change the phenomenology.

Any quantity $`\rho (t)`$ alternating in the way described above can be written as:

$$\rho (t)=\rho _++\rho _{}\mu (t),$$
(4)
where $`\rho _\pm =\frac{1}{2}(\rho _1\pm \rho _2)`$ and $`\mu (t)`$ is a periodic square wave

$$\mu (t)=\{\begin{array}{cc}\hfill 1:& 0<t<\frac{T}{2}\hfill \\ \hfill 1:& \frac{T}{2}<t<T\hfill \end{array}.$$

Let us now suppose the following conditions for the sets $`\{\rho _i\}`$ according to seasonality:

$$a_1<a_2,b_1<b_2,c_1<c_2,K_1>K_2,$$
(5)
where 1 stands for the rainy season and 2 for the dry one. The biological motivation for these conditions is the following. The dry season provides less resources for the colony than the rainy season $`\left(K_2<K_1\right)`$, and as a consequence the death rate is higher $`\left(c_2>c_1\right)`$, and the transmission rate is also larger since fights for the available resources are expected to increase $`\left(a_2>a_1\right)`$. However, notice that we consider the birth rate to be larger during the dry season $`\left(b_2>b_1\right)`$. There are two reasons for this. First is the assumption that the colony makes an attempt to maintain its population. Second is an implementation of the maturation process in the model: baby mice do not contribute to the propagation of the disease mills1 , and, therefore, even if births were actually more numerous during the rainy season, the contribution of this fact to the propagation of the infection will only be important in the next (dry) season, when mice have matured and are ready to fight. It will be shown later that these assumptions lead to a situation with high population during the rainy season and low population during the dry one, in agreement with the available data mills1 ; hjelle .

Moreover, we assume that $`K_1<\mathrm{min}(K_{c1},K_{c2})`$, i.e., the resources are all times (during both seasons) below the minimum critical threshold that triggers the propagation of the disease. We will show that nevertheless it is possible for the infection to spread.

The equilibrium populations of the susceptible and infected mice are determined by the set of values $`\{\rho _i\}`$. When switching from season to season, the populations evolve trying to reach a new equilibrium. Therefore, the dynamics is driven by the competition between two characteristic times. On the one hand there is an external time scale determined by the seasonal forcing, $`t_e=T/2`$. On the other hand, the relaxation toward equilibrium after a switching of seasons involves a relaxation time. The latter measures the time required for the mouse colony to relax to the equilibrium state associated with $`\{\rho _j\}`$ after having been driven during the previous season by the conditions $`\{\rho _i\}`$, that is, $`t_r(ij)`$ where $`i,j=1,2`$. The internal time scale is defined as the fastest relaxation process, i.e., $`t_i=\mathrm{min}[t_r(ij)]`$. “Fast” or “slow” seasonal alternation then refers to the comparison between these two time scales. If $`t_et_i`$ the mouse population has enough time to accommodate to the new conditions from season to season and relax to equilibrium. Moreover, since we have imposed the condition that the resources at any time of the year are below the critical thresholds $`K_{ci}`$, there will be no infected mice. In the other limit, $`t_et_i`$, seasonal changes occur too fast, the relaxation process is interrupted, and no equilibrium can be reached from season to season. Then note that an adiabatic elimination can be implemented gardiner , and $`\mu (t)`$ in Eq. (4) can be replaced by its average value, $`\mu (t)=0`$. Therefore, in the limit of fast season alternation the system is driven by the set of averaged values $`\{\rho _+\}=\{a_+,b_+,c_+,(1/K)_+\}`$, and the critical carrying capacity is given by $`K_{c+}=\frac{b_+}{a_+(b_+c_+)}`$. As a consequence, it is possible to find regions of parameters where $`K_{c+}`$ is smaller than the effective value of the carrying capacity associated with the averaged values:

$$\left[\left(\frac{1}{K}\right)_+\right]^1=\left[\frac{1}{2}\left(\frac{1}{K_1}+\frac{1}{K_2}\right)\right]^1=\frac{2K_1K_2}{K_1+K_2},$$

and the infection propagates.

General conditions leading to this behavior can be posed, but the expressions are rather cumbersome. We prefer, for the sake of simplicity, to show a particular typical case. We use the following values for the parameters:

$$a_1=\frac{1}{4},b_1=1,c_1=\frac{1}{3},$$
(6)
$$a_2=4,b_2=\frac{73}{72},c_2=1,$$
(7)
that lead to $`K_{c1}=6`$ and $`K_{c2}=73/4`$ respectively. The dynamics are completely determined once the value of the carrying capacity during each season is specified. According to the previous discussion, these parameters can be chosen such that the following conditions hold:

$$\left[\left(1/K\right)_+\right]^1>K_{c+},K_1>K_2,K_1<\mathrm{min}(K_{c1},K_{c2}).$$

These conditions lead to the points $`(K_1,K_2)`$ that fulfill the seasonal requirements given by Eq. (5), so that slow alternation of seasons leads to infection-free states while fast alternation leads to Hanta outbreaks. This region is plotted in Fig. 2. Notice in particular that the point $`\left(K_1=4,K_2=1\right)`$ lies inside the region and that $`K_i<K_{ci}`$. In the next section we illustrate the seasonality-induced propagation of the disease for this particular point.

## IV The Critical Period

So far we have determined that outbreaks of Hanta induced entirely by seasonal changes can occur if the duration of the seasons are short enough. Now we establish the meaning of “short enough” quantitatively. Since $`K_{c+}`$ is strictly smaller than the effective value of the carrying capacity, there should be a finite value of $`T_c`$ such that for any $`T<T_c`$ the population of infected mice is greater than zero, but for periods above this critical period the infected population goes to zero.

In order to obtain the value of the critical period we need to solve the system of equations (1a) and (1b). In spite of its nonlinearities the system can be solved analytically by means of a reciprocal transformation kenkre  and the following exact solution is obtained:

$`M_I(t,I_0,S_0;\{\rho \})`$ $`=`$ $`{\displaystyle \frac{I_0\mathrm{\Omega }(t)}{S^{aK1}+aI_0_0^t\mathrm{\Omega }(\tau )𝑑\tau }},`$ (8a)
$`M_S(t,I_0,S_0;\{\rho \})`$ $`=`$ $`{\displaystyle \frac{SM_0e^{\frac{S}{K}t}}{(\mathrm{\Omega }(t)e^{ct})^{\frac{1}{aK1}}}}{\displaystyle \frac{I_0\mathrm{\Omega }(t)}{S^{aK1}+aI_0_0^t\mathrm{\Omega }(\tau )𝑑\tau }},`$ (8b)
where $`I_0`$ and $`S_0`$ are the initial conditions for $`M_I`$ and $`M_S`$ respectively, and the following definitions have been introduced,

$$\mathrm{\Omega }(t)=e^{ct}\left(M_0\left(e^{\frac{S}{K}t}1\right)+S\right)^{aK1},$$

$$S=K(bc),M_0=I_0+S_0.$$

Because the external forcing due to the alternation of seasons is periodic, Floquet theory ensures the existence of a periodic solution. The values of $`I_0`$ and $`S_0`$ compatible with the non-equilibrium periodic solution can be obtained by evolving the system during the first half of a period under dynamics $`1`$ and the second half under dynamics $`2`$, and forcing periodity on the solutions after a whole period of evolution, that is,

$`M_I({\displaystyle \frac{T}{2}},M_I({\displaystyle \frac{T}{2}},I_0,S_0;\{\rho _1\}),M_S({\displaystyle \frac{T}{2}},I_0,S_0;\{\rho _1\});\{\rho _2\})`$ $`=`$ $`I_0,`$ (9a)
$`M_S({\displaystyle \frac{T}{2}},M_I({\displaystyle \frac{T}{2}},I_0,S_0;\{\rho _1\}),M_S({\displaystyle \frac{T}{2}},I_0,S_0;\{\rho _1\});\{\rho _2\})`$ $`=`$ $`S_0.`$ (9b)
In order to close the system in the non-equilibrium steady state $`M_I(t,T;\{\rho _{1,2}\})`$ and $`M_S(t,T;\{\rho _{1,2}\})`$, the values of $`I_0`$ and $`S_0`$ that solve that system of equations (9a) and (9b) must be then re-introduced in Eqs. (8a) and (8b). The critical period is then the largest value of $`T`$ satisfying the condition

$$M_I(t,T;\{\rho _{1,2}\})>0.$$

We illustrate the procedure with the example mentioned above where the parameters are given by Eqs. (6) and (7), and with $`K_1=4`$ and $`K_2=1`$. The results are shown in Figs. 3 and 4. The values of $`M_I`$ and $`M_S`$ as a function of the period of the seasons are depicted in Fig. 3, where the populations of the susceptible and infected mice at the end of each season are given. As seen in that figure, the value of the critical period is $`T_c12`$. Note that if the alternation is slow, $`T>T_c`$, all mice grow healthy. On the other hand, if the alternation is faster than the relaxation time required by the colony to accommodate its population from season to season, $`T<T_c`$, the virus spreads and $`M_I>0`$. Notice that in the limit $`T0`$ the dynamics is driven by $`\{\rho _+\}`$ and the populations of susceptible and infected mice are given by Eq. (3) with $`a=a_+`$, $`b=b_+`$, $`c=c_+`$, and $`1/K=(1/K)_+`$. Let us stress again that the carrying capacity is below its critical threshold at any time.

In Fig. 4 we plot, for different period lengths, the exact solutions $`M_I(t,T;\{\rho _{1,2}\})`$ and $`M_S(t,T;\{\rho _{1,2}\})`$ as a function of time through one period of evolution. The first semi-period corresponds to the rainy season and the second to the dry season. When seasons last long (left panel), there are no infected mice and the susceptible population simply oscillates between the two equilibrium states given by Eq. (2). For sufficiently short seasons (right panel), there is propagation of the disease and the values of $`M_I`$ and $`M_S`$ fluctuate around the equilibrium points determined by Eq. (3) and the set of parameters $`\{\rho _+\}`$. Finally, when the period of the seasons is near, but below, the critical period (central panel), the infected population is small and the population $`M_I`$ oscillates in a more pronounced manner.

## V Conclusions

By introducing seasonality in a paradigmatic model for Hantavirus propagation in mice colonies, we have shown that the alternation of seasons may cause outbreaks of the disease. The striking feature of that behavior lays in the fact that neither season satisfies the conditions for the infection to spread in terms of the availability of resources. The mechanism responsible for the phenomenon is the competition between two time scales: an external one, the duration of a season, and an internal one, the relaxation time for the mouse colony to equilibrate its population from season to season. We have shown that if the duration of the seasons is longer than the relaxation time, no propagation of Hantavirus occurs. On the other hand, if the relaxation process is interrupted by a fast seasonal alternation, the disease spreads. We have analyzed the general conditions for which the phenomenon occurrs. Moreover, we have illustrated the mechanism with a particular example that can be solved exactly.

This work may help to clarify the reported relation between climate and propagation of Hanta in mice populations. However, to elucidate whether the proposed phenomenon actually takes place in nature we depend on data that unfortunately are not available in the literature. One can envision further modifications of the model that may improve its features, such as, for example, the inclusion of spatial dependence or of noisy contributions to the dynamics. Finally, we stress that the general idea underlying the mechanism is model-insensitive and can therefore be extended to other systems where seasonality plays a relevant role. Work along these directions is in progress.

## Acknowledgments

The authors gratefully acknowledge fruitful comments from and discussions with V. M. Kenkre during the elaboration of this work. C. Escudero is grateful to the Department of Chemistry and Biochemistry of the University of California, San Diego for its hospitality. This work has been partially supported by the National Science Foundation under grant PHY-9970699, by the Ministerio de Educación y Cultura (Spain) through grants No. AP2001-2598 and EX2001-02880680, and by the Ministerio de Ciencia y Tecnología (Spain), Project No. BFM2001-0291.


text:
# 1 The Dalitz plots of the 𝐷⁰→𝐾⁺𝐾⁻𝜋⁰ decay. Top panel: constructive interference (cos{𝛿_𝐷}=1), 113 events in the square region; bottom panel: destructive interference (cos{𝛿_𝐷}=-1), 4 events in the square region. The total number of events in the bands is 𝑁=1500 in both cases.

## Acknowledgments

We wish to thank D. M. Asner and M. Gronau for helpful correspondence. This work was supported in part by the United States Department of Energy through Grant No. DE FG02 90ER40560.

## Appendix A: Kinematics and decay amplitudes

The first stage of the $`D^0K^{}(K^+\pi ^0)_{K^+}`$ process is the decay of a pseudoscalar meson $`D^0`$ into a pseudoscalar $`K^{}`$ and a (possibly off-shell) vector $`K^+`$. Afterwards, the latter decays into $`K^+`$ and $`\pi ^0`$. From angular momentum conservation the helicity of $`K^+`$ is $`0`$. The corresponding polarization vector is $`ϵ_{K^+}=ϵ^{(\lambda =0)}=(|𝐩_{K^+}|,0,0,E_{K^+})/m_{K^+\pi ^0}`$. ($`m_{K^+\pi ^0}`$ is the invariant mass of $`K^+`$ and the $`z`$ axis is chosen to point in the direction of the $`K^+`$ momentum $`𝐩_{K^+}`$, see Fig. 5).

The amplitude $`A_1(m_{K^+\pi ^0},m_{K^{}\pi ^0})`$ of the $`K^+K^+\pi ^0`$ decay should be Lorentz invariant, i.e., it should contain a product of two 4-vectors. There is only one non-vanishing possibility, $`ϵ_{K^+}(p_{K^+}p_{\pi ^0})`$, since the other, $`ϵ_{K^+}(p_{K^+}+p_{\pi ^0})=ϵ_{K^+}p_{K^+}`$, is identically zero. Then the former can be written in the rest frame of $`K^+`$ and $`\pi ^0`$ as

$$A_1(m_{K^+\pi ^0},m_{K^{}\pi ^0})(0,0,0,1)(E_{K^+}^{}E_{\pi ^0}^{},2𝐩_{K^+}^{})=2|𝐩_{K^+}^{}|\mathrm{cos}\theta ^{},$$
(16)
where $`\theta ^{}`$ is the angle between the negative direction of the $`z`$ axis and the direction of the $`K^+`$ momentum $`𝐩_{K^+}^{}`$ in the rest frame of $`K^+`$ and $`\pi ^0`$. We will keep using the “\*” subscript for quantities determined in this frame.

$`\mathrm{cos}\theta ^{}`$ is given by

$$\mathrm{cos}\theta ^{}=\frac{m_{K^{}\pi ^0}^2m_K^{}^2m_{\pi ^0}^22E_K^{}^{}E_{\pi ^0}^{}}{2|𝐩_K^{}^{}||𝐩_{\pi ^0}^{}|},$$
(17)
so

$$A_1(m_{K^+\pi ^0},m_{K^{}\pi ^0})\frac{m_{K^{}\pi ^0}^2m_K^{}^2m_{\pi ^0}^22E_K^{}^{}E_{\pi ^0}^{}}{|𝐩_K^{}^{}|},$$
(18)
where

$$E_K^{}^{}=(m_{D^0}^2m_K^{}^2m_{K^+\pi ^0}^2)/2m_{K^+\pi ^0},$$
(19)
$$E_{\pi ^0}^{}=(m_{K^+\pi ^0}^2m_{K^+}^2+m_{\pi ^0}^2)/2m_{K^+\pi ^0},$$
(20)
$$|𝐩_K^{}^{}|=\lambda ^{1/2}(m_{D^0}^2,m_K^{}^2,m_{K^+\pi ^0}^2)/2m_{K^+\pi ^0},$$
(21)
$$|𝐩_{\pi ^0}^{}|=|𝐩_{K^+}^{}|,$$
(22)
$$\lambda (x,y,z)x^2+y^2+z^22xy2xz2yz.$$
(23)
Including the finite resonance width $`\mathrm{\Gamma }_{K^+}`$ into the $`K^+`$ propagator, we can write the amplitude of the $`D^0K^{}(K^+\pi ^0)_{K^+}`$ decay as

$$A(D^0K^{}(K^+\pi ^0)_{K^+})\frac{A_1(m_{K^+\pi ^0},m_{K^{}\pi ^0})}{m_{K^+\pi ^0}^2m_{K^+}^2+im_{K^+}\mathrm{\Gamma }}.$$
(24)
As for the $`D^0K^+(K^{}\pi ^0)_K^{}`$ decay, its amplitude can be derived in a similar way and is equal to

$$A(D^0K^+(K^{}\pi ^0)_K^{})r_De^{i\delta _D}\frac{A_2(m_{K^+\pi ^0},m_{K^{}\pi ^0})}{m_{K^{}\pi ^0}^2m_K^{}^2+im_K^{}\mathrm{\Gamma }},$$
(25)
with the kinematic factor $`A_2`$ defined as $`A_2(m_{K^+\pi ^0},m_{K^{}\pi ^0})A_1(m_{K^{}\pi ^0},m_{K^+\pi ^0})`$. The factor $`r_De^{i\delta _D}`$ accounts for possible differences in hadronization as vector particles between quarks arising from the virtual $`W^+`$ and spectator quarks.

### Calculation of the fraction $`ϵ`$ of resonant decays that fall into the square region

For the particular case of an on-shell resonant $`K^+`$ we can neglect the Breit-Wigner denominator of Eq. (24). In this case the amplitude of the $`D^0K^{}(K^+\pi ^0)_{K^+}`$ decay is proportional to $`A_1(m_{K^+},m_{K^{}\pi ^0})`$. The kinematics of the two-body $`D^0K^{}K^+`$ and $`K^+K^+\pi ^0`$ decays determine $`E_K^{}^{}=1.37`$ GeV, $`E_{\pi ^0}^{}=0.32`$ GeV, $`|𝐩_K^{}^{}|=1.27`$ GeV and $`|𝐩_{\pi ^0}^{}|=0.29`$ GeV. As a result, Eq. (17) says:

$$\mathrm{cos}\theta ^{}=1.36(m_{K^{}\pi ^0}^21.135),$$
(26)
where $`m^2(K^{}\pi ^0)`$ is in GeV2. Thus, the amplitude of the $`D^0K^{}(K^+\pi ^0)_{K^+}`$ decays is proportional to $`(m_{K^{}\pi ^0}^21.135)`$. These resonant decays fill the vertical band in a nonuniform way: no decays happen at the middle of the band where $`m_{K^{}\pi ^0}^21.135=0`$. The majority of the events will concentrate near both band ends where $`|m_{K^{}\pi ^0}^21.135|`$ is the largest.

Now we can calculate the fraction of $`D^0K^{}(K^+\pi ^0)_{K^+}`$ decays that fall into the square region,

$$ϵ=\underset{0.75}{\overset{0.84}{}}(x1.135)^2𝑑x/\underset{0.40}{\overset{1.87}{}}(x1.135)^2𝑑x=0.039,$$
(27)
where $`(m_{K^+}\mathrm{\Gamma }/2)^2=0.75`$ GeV2 and $`(m_{K^+}+\mathrm{\Gamma }/2)^2=0.84`$ GeV2 are the boundaries of the square region and $`0.40`$ and $`1.87`$ are the boundaries of the whole band. The latter can be derived from Eq. (26).

This simple calculation implied that the population density of the vertical band is constant along any cross section of the band, i.e., at a fixed $`m_{K^{}\pi ^0}^2`$ it is independent of variations of $`m_{K^+\pi ^0}^2`$ across the band. A more precise discussion involves a simulation of the interference between the Breit-Wigner resonant shapes of Eqs. (24) and (25).

## Appendix B: Influence of scalar resonance $`\kappa `$

The existence of broad scalar resonances below $`1`$ GeV has been a controversial issue for a long time . A few experiments have been able to explore the possibility of their presence as intermediate resonant states in three-body $`D`$ decays. The modes that were studied include $`D^+\pi ^{}\pi ^+\pi ^+`$ , $`D_s^+\pi ^{}\pi ^+\pi ^+`$ , $`D^+K^{}\pi ^+\pi ^+`$  (E791 collaboration), $`D^0K_S^0\pi ^+\pi ^{}`$ , $`D^0K^{}\pi ^+\pi ^0`$  (CLEO), and $`D^0K^0K^{}\pi ^+`$  (BaBar). The first two studies obtained evidence for a light ($`478`$ MeV) $`\sigma `$ resonance and measured the properties of the $`f_0(980)`$. The last four might provide some information on the presence of an intermediate S-wave $`K\pi `$ resonance. Indeed, the E791 analysis of a Dalitz plot found that the best fit to the data is obtained allowing for the presence of an additional scalar resonance $`\kappa (800)^0`$. However, neither CLEO studies found evidence for $`\kappa ^0`$ or its isodoublet partner $`\kappa ^+`$. The preliminary BaBar analysis saw $`\kappa `$ at the level of $`1\sigma `$ which does not allow the confirmation of its presence. Other types of decays could also provide a glimpse of $`\kappa `$. The BES collaboration found $`\kappa ^0`$ as an intermediate state in $`J/\psi \overline{K}^{}(892)^0K^+\pi ^{}`$ decays , while the FOCUS collaboration studied the interference phenomena in $`D^+K^{}\pi ^+\mu ^+\nu `$ decays . Their data can be described by $`\overline{K}^0`$ interference with either a constant amplitude or a broad spin zero resonance.

The $`D^0K^{}K^\pm `$ decays discussed in this note can be affected by the possible presence of $`\kappa ^\pm `$ among the intermediate states. The bands of a broad $`\kappa (800)`$ would cover more than 50% of the Dalitz plot, thereby interfering with the $`K^{}`$ bands and affecting their population. One would expect that in this case the total branching ratio of $`D^0K^+K^{}\pi ^0`$ decays would be considerably larger than the sum of the $`D^0KK^{}`$ modes. Indeed, in $`D^+K^{}\pi ^+\pi ^+`$ an unusually high fraction (over 90%) of decays was found to be non-resonant by previous experiments . That was unusual as the non-resonant (NR) contribution in three-body decays is small in most other cases. That was an indication of a possible broad scalar contribution and motivated the recent searches for it. It was found that the complex structure of the Dalitz plot was best explained when the $`\kappa `$ presence is assumed . Then, intermediate decays through the $`\kappa \pi ^+`$ state account for about 50% of decays while the NR fraction drops to a value of 13% more characteristic of other decays.

The present knowledge of $`D^0K^{}K^\pm `$ decays does not reveal a similar large non-resonant (or broad scalar) contribution. The current data on the resonant  and inclusive  decays comes from CLEO measurements. The inclusive branching ratio is $`(D^0K^+K^{}\pi ^0)=(1.24\pm 0.35)10^3`$. The branching ratios of the $`K\pi `$ resonant decays are $`(D^0K^+(K^{}\pi ^0)_K^{})=\frac{1}{3}(2.0\pm 1.1)10^3=(0.67\pm 0.37)10^3`$ and $`(D^0K^{}(K^+\pi ^0)_{K^+})=\frac{1}{3}(3.8\pm 0.8)10^3=(1.27\pm 0.27)10^3`$. Neglecting the interference between these two channels (it affects just about 4% of these decays; see Appendix A), the two branching ratios add up to $`(1.93\pm 0.45)10^3`$, consistent with the inclusive branching ratio within the current large uncertainties. Basically, there is no room for a broad scalar resonance channel. For example, it cannot negatively interfere with both halves of a $`K^{}`$ band. The phase variation across it would be significant ($`90^{}`$) for a $`K^{}`$ channel and much smaller for a broad $`\kappa `$ one. If this channel is strong enough to cancel half the $`K^{}`$ decays it would contribute many times more than that outside the $`K^{}`$ bands. That would contradict the smallness of the inclusive branching ratio. Thus, we conclude that a broad scalar $`\kappa `$, if present, could only comprise a small fraction of $`D^0K^+K^{}\pi ^0`$ decays and would not significantly affect the estimate of the uncertainty in the strong phase $`\delta _D`$.