text:
# Cosmic Ray History and its Implications for Galactic Magnetic Fields

## 1 Introduction

Superthermal particles, or cosmic rays, are a major constituent of the interstellar medium in galaxies. They are dynamically coupled to the thermal gas, and collectively provide pressure support, can drive outflows, and are thought to modify the structure of shocks. Cosmic rays are the dominant source of ionization, and a major source of heating, in regions of high visual extinction. Spallation reactions between cosmic rays and ambient interstellar nuclei are the primary mechanism for synthesizing light elements.

The energy spectrum, directional anisotropy, and chemical composition of cosmic rays are measured directly through in situ observations on the Earth and in the heliosphere. The global properties of cosmic rays in the Galaxy are probed through their radio frequency synchrotron emission and $`\gamma `$-ray line and continuum emission. A recent summary of the observations is given by Schlickeiser (2002).

Although synchrotron and $`\gamma `$-ray radiation from normal galaxies have up to now been directly detected only within the local universe, there is evidence that cosmic rays existed in galaxies even at early times. The light elements present in metal poor halo stars are thought to originate by spallation reactions between energetic particles and the material from which these stars formed \[Duncan, Lambert, & Lemke (1992), Gilmore et al. (1992), Duncan et al. (1998), Boesgaard et al. (1999) and references therein\]. Further evidence is provided by the spectrum of the diffuse $`\gamma `$-ray background, which is best fit by a superposition of active and normal galaxy $`\gamma `$-ray spectra at moderate redshift \[Vasiliki & Fields (2002)\].

The energy densities in the Galactic magnetic field and cosmic rays are now roughly equal. If this relationship is universal and necessary, then the presence of cosmic rays in young galaxies implies the existence of equipartition magnetic fields. This would be an important constraint on theories of the history of magnetic fields in galaxies, a subject on which evidence is sparse. At present, the only direct observations of magnetic fields associated with young galaxies are the Faraday rotation measurements in damped Ly$`\alpha `$ systems at redshift $`z1`$ reported by Oren & Wolfe (1995). These observations suggest magnetic fields of $`1\mu `$G strength which are spatially coherent over several kpc, similar to the fields in the Milky Way and other spiral galaxies. Inferences based on cosmic rays could probe galactic magnetic field evolution at earlier epochs, and, we will see, could be directly relevant to conditions in star forming regions. This is significant because of the important role magnetic fields can play in star formation and their possible effect on the Initial Mass Function.

The leading theories of cosmic ray acceleration and confinement depend on the strength of the magnetic field. In this paper we use this dependence to estimate the minimum magnetic field at which a population of relativistic cosmic rays with energy density comparable to the present Galactic population can be accelerated and confined. We find that a field several orders of magnitude weaker than the present Galactic field is sufficent for both acceleration and confinement, provided that the coherence length of the field exceeds the cosmic ray gyroradius.

The plan of this paper is as follows. In §2 we state our assumptions about those conditions in early galaxies which affect the acceleration and propagation of cosmic rays. In §3 we briefly review the theory of energetic particle transport and its application to cosmic rays. In §4 we derive constraints which arise by assuming that the cosmic rays diffuse at the minimum possible rate, which at a given fieldstrength gives the maximum acceleration and confinement. In §5 we discuss hydromagnetic turbulence in a weak magnetic field and find that it is drastically affected by even a small population of cosmic rays. Fluctuations grow rapidly in the presence of cosmic ray anisotropy, and, in contrast to fluctuations driven by cosmic ray anisotropy in the present Galactic environment, are not Alfvénic in character. At small amplitudes, the damping rate is much lower than the excitation rate. Therefore, we predict that the waves grow to nonlinear amplitudes. This may have implications for the nature of turbulence in other environments in which the magnetic fields is relatively weak and a population of energetic particles is present. Section 6 is a summary and conclusion. Most of the technical material is in Appendices A - C.

We use Gaussian cgs units throughout, except in expressing the energies of cosmic rays, where we follow standard practice and use eV, and also in certain formulae, in which we use “natural” units such as parsecs or years.

## 2 Conditions in Early Galaxies

We assume that the density and temperature of interstellar gas in young galaxies are in the range found in contemporary galaxies. At lower metallicity the cooling rate is reduced, and the period in which supernova remnants are nonradiative is extended. This moderately enhances the efficiency of cosmic ray acceleration. We assume that Type II supernovae are clustered, as they are now, and may occur at a higher rate than at present. We consider the overall galactic size scale to be similar to present galaxy sizes, or perhaps somewhat smaller. The values we choose are appropriate to spheroidal systems with radii comparable to the present thickness of galactic disks.

The most important assumptions concern the magnetic field. A recent review of the origin of cosmological magnetic fields is given by Widrow (2002), and discussions of magnetic fields in galaxies are given in Beck et al. (1996) and Zweibel & Heiles (1997). Briefly, theories of the origin of galactic magnetic fields can be classified as top down or bottom up. In top down theories, the fields arise through large scale coherent processes and have large scale structure, although they are generally very weak. Examples of top down theories include magnetogenesis during inflation (Turner & Widrow, 1988), by the operation of the Biermann battery in cosmological shock fronts or other vortex structures \[Pudritz & Silk (1989), Kulsrud et al. (1997), Davies & Widrow (2000)\], in cosmological ionization fronts \[Subramanian, Narasimha, & Chitre (1994), Gnedin, Ferrara, & Zweibel (2000)\], and through the manyfold expansion of plasma lobes associated with active galactic nuclei \[Furlanetto & Loeb (2001), Kronberg et al. (2001)\]. Although only the first of these theories produces fields which are truly coherent over cosmological scales, all three produce fields which have some coherence over galactic scales.

In bottom up theories, there are many independent, small scale sources of field. The magnetized plasma from these sources expands and diffuses to fill the interstellar medium, eventually merging to create pervasive fields with largely random structure. Examples of the sources include plerion supernova remnants and jets or winds from accretion disks surrounding stellar mass compact objects \[Zeldovich, Ruzmaikin, & Sokoloff (1983), Rees (1987)\]. The magnetic field is inhomogeneously distributed until it has had time to diffuse. For example, if magnetic fields were seeded by massive stars, there could have been fields in OB associations before there was a large scale galactic field.

Both top down and bottom up theories predict magnetic fields which are much weaker than present galactic fields, and must be amplified by several orders of magnitude or more to resemble the fields of several $`\mu `$G which are observed in galaxies. It is thought that amplification is due to the action of a hydromagnetic dynamo, which stretches the field by fluid motions while dissipating it at the resistive scale. Since resistive effects typically operate at an AU or less, while the coherence lengths of galactic magnetic fields are observed to be several kpc or more, a longstanding issue is how the fieldlines can be lengthened, within a fixed volume, without the extensive folding and tangling that characterizes a predominantly small scale field \[see e.g. Kulsrud & Anderson (1992) and Schekochihin et al. (2002)\]

For our purposes, the chief distinctions between top down and bottom up theories are that the former predict a large scale, spatially homogeneous field while the latter predict a small scale, inhomogeneous field. Partly in recognition of bottom up theories, we will discuss cosmic ray acceleration by supernova shocks and propagation in the general galactic field more or less independently of one another.

## 3 Transport Theory with Application to Galactic Cosmic Rays

There is a well developed theory for interactions between energetic particles and low frequency electromagnetic fluctuations $`𝑩_1`$, $`𝑬_1`$ superimposed on a mean magnetic field $`𝑩_0`$. Reviews are given by Jokipii (1971), Wentzel (1974), Toptygin (1985), and Schlickeiser (2002). This theory has been reasonably successful in explaining the $`10^7`$ yr confinement time of cosmic rays, their near isotropy, and their energy density, the latter provided that approximately 10% of supernova energy is converted to cosmic rays. In the interests of keeping the paper self contained, we provide a brief treatment here.

A particle with velocity $`𝒗`$ and cyclotron frequency $`\omega _{cr}ZeB_0/mc\gamma `$ interacts primarily with the Fourier components of $`𝑩_1`$ and $`𝑬_1`$ which satisfy the cyclotron resonance condition

$$\omega k_{}v_{}\pm \omega _{cr}=0,$$
(1)
where $`\omega `$ and $`k`$ are the frequency and wavenumber of the Fourier mode and the subscript “$``$” denotes the projection along $`𝑩_0`$.

If the spectrum of waves is incoherent, if a particle interacts with a wave for one gyroperiod, and if the effect of any single interaction is small, then the the wave-particle interaction produces diffusion in momentum space. Although the particles diffuse in both the magnitude and direction of their momentum $`𝒑`$, the diffusion in angle is dominant as long as the phase velocity of the resonant waves is much less than the speed of light. The rate of angular scattering as a function of the variable $`\mu p_{}/p`$ is

$$\nu (\mu )\frac{\pi }{4}\omega _{cr}\left(k_R(k_R)\right)^1,$$
(2)
where $``$ is the power spectrum of the magnetic field fluctuations, normalized to the large scale magnetic field, and $`k_R`$ is the value of $`k_{}`$ which satisfies eqn. (1).

The resonance condition eqn. (1) involves only the parallel component of $`𝒌`$. However, if $`k_{}k_{}`$, the oscillations of the wave electromagnetic fields over the gyrorbit diminish the response of the particle, and the scattering frequency is reduced (Chandran, 2000a).

Locally, the fluctuations define a frame with velocity $`𝒖_f`$, which is the mean fluctuation velocity measured in the rest frame of the observer and weighted by the intensity spectrum. If the fluctations propagate isotropically, or if the velocity $`𝒖`$ of the fluid is much larger than the relative velocities of the fluctuations and the fluid, then $`𝒖_f`$ can be approximated by $`𝒖`$ itself. If the scattering time $`\nu ^1`$ and associated mean free path $`\lambda _{}\mu v/\nu `$ are short compared to the global lengthscales and timescales of interest, the particles are well coupled to the fluctuation frame, with a small amount of diffusion caused by scattering. Under these conditions, the phase space distribution function $`f(𝒙,𝒑,t)`$ of the cosmic rays is nearly isotropic with respect to $`𝒑`$. The isotropic part of $`f`$, averaged over the ensemble of fluctuations, is governed by a transport equation which includes advection and compression by the fluctuations, spatial diffusion with diffusion tensor $`𝑫`$, and a source $`S_0(𝒙,𝒑,t)`$

$$\frac{f_0}{t}+𝒖_f\frac{f_0}{𝒙}\frac{1}{3}\left(\mathbf{}𝒖_f\right)p\frac{f_0}{p}\frac{}{𝒙}𝑫\frac{f_0}{𝒙}=S_0$$
(3)
(Skilling, 1975). Particle drifts caused by large scale curvature or gradients in $`𝑩_0`$ can be included explicitly in eqn. (3), in which case they appear in the advection term (Forman, Jokipii, & Owens, 1974).

For our purposes, it suffices to consider a diagonal diffusion tensor with components $`D_{}`$ and $`D_{}`$. The parallel diffusivity is related to the scattering frequency $`\nu `$ by

$$D_{}\frac{v^2}{4}_1^1𝑑\mu \frac{1\mu ^2}{\nu (\mu )}$$
(4)
According to eqns. (2) and (4), $`D_{}`$ can be approximated by

$$D_{}\frac{1}{3}\frac{v^2}{\omega _c}\left(k_R(k_R)\right)^1.$$
(5)
The perpendicular diffusivity $`D_{}`$ is given approximately in terms of $`D_{}`$ by (Jokipii, 1987)

$$D_{}\frac{D_{}}{1+\left(k_R(k_R)\right)^2}.$$
(6)
Equation (6) shows that cross field diffusion is weak if the fluctuation amplitude is small, but becomes nearly isotropic if the fluctuations are nonlinear ($`k_R(k_R)1`$). According to eqn. (2), if the fluctuation amplitude is nonlinear, the scattering frequency is approximately the cyclotron frequency. This is also the timescale on which a particle interacts with a single wave.

Equation (3) is the foundation of the leading theories of cosmic ray acceleration and confinement. It gives the confinement time $`\tau _c`$ to a region of size $`L`$ as

$$\tau _c\frac{L^2}{D},$$
(7)
where $`D`$ is the largest component of the diffusion tensor, generally $`D_{}`$. Equation (7) is consistent with the observed confinement time of Galactic cosmic rays if $`k_R(k_R)10^610^7`$ (Cesarsky, 1980). Equation (7) is an upper limit in the sense that it does not include transport due to migration of the fieldlines themselves, which may be considerable if the lines of force have a large stochastic component (Berezinskii et al, 1990).

Equation (3) predicts that in a steady state in which diffusion of cosmic rays is balanced by sources, with little advection or compression, the relationship between the cosmic ray energy density $`U_{cr}`$ and the power density of the source $`𝒫_{cr}`$ is given approximately by

$$U_{cr}\tau _c𝒫_{cr}𝒫_{cr}\frac{L^2}{D}.$$
(8)
In the Milky Way, $`U_{cr}10^9\mathrm{GeV}\mathrm{cm}^3`$, and $`\tau _c`$ is known from the isotopic ratios of light elements to be $`10^7`$ yr. This implies that $`𝒫_{cr}5\times 10^{27}`$ erg cm$`{}_{}{}^{3}s_{}^{1}`$, which is about 5% of the power density in supernovae.

If the scattering is strong and the distribution function $`f`$ is nearly isotropic, the cosmic rays can be treated as a fluid, and the dynamic and energetic coupling of the cosmic ray and thermal fluids, which is mediated by scattering, can be expressed in hydrodynamic form. The rate of momentum transfer $`𝑭_{cr}`$ takes the form of a pressure gradient force

$$𝑭_{cr}=\mathbf{}P_{cr}$$
(9)
acting on the thermal fluid. The rate of energy transfer to the thermal fluid $`\dot{E}_{cr}`$ can be written as

$$\dot{E}_{cr}=\gamma _{cr}(𝒖_f𝒖)\mathbf{}P_{cr}$$
(10)
(Völk, Drury, & McKenzie, 1984), where $`\gamma _{cr}`$, the polytropic exponent of the cosmic rays, is usually taken to be 4/3. The rate of energy transfer is proportional to the frictional force between the cosmic rays and thermal fluid. The dynamical effects represented by eqns. (9) and (10) are responsible for hydrostatic support, acceleration of outflows, and modification of shock structure in the thermal interstellar gas.

According to the leading theory of cosmic ray origin \[original papers by Axford, Leer, & Skadron (1977), Bell (1978a), Bell (1978b), Blandford & Ostriker (1978); reviews by Blandford & Eichler (1987), Jones & Ellison (1991)\], the particles are accelerated by the first order Fermi process operating in strong interstellar shocks driven by supernovae. Acceleration requires hydromagnetic turbulence both upstream and downstream of the shock. In a frame in which the shock is steady, the upstream turbulence approaches the shock at nearly the shock speed $`u_S`$ (the wave speed being negligible compared to the shock speed), and the downstream turbulence recedes from the shock at speed $`u_S/r`$, where $`r`$ is the compression ratio. Particles which propagate upstream away from the shock, are reflected back across the shock, and are then scattered back upstream by the postshock turbulence gain an angle averaged net momentum of $`4m\gamma u_S(r1)/3r`$, much as if they were trapped between converging walls.

The steady state spectrum of particles accelerated by a parallel shock is a power law in momentum space, $`f(p)p^q`$, where the spectral index $`q`$ is a function only of $`r`$: $`q=3r/(r1)`$. The power spectrum is realized only within a range $`p_{min}<p<p_{max}`$. The lower limit $`p_{min}`$ is the minimum momentum at which the conditions of shock acceleration theory are fulfilled: the particle speed $`v`$ must exceed the shock speed $`u_S`$ and Coulomb losses must be insignificant. The upper limit $`p_{max}`$ is the maximum momentum a particle can have reached since the onset of acceleration. In §4, we will use the expression for $`p_{max}`$ derived by Lagage & Cesarsky (1983) to derive a lower limit on the magnetic fieldstrength such that shocks can accelerate particles to relativistic energies.

The original theory of shock acceleration predicts the shape of the energetic particle spectrum, but not its normalization. The latter depends on the mechanisms by which particles are injected into the shock. This so-called injection problem is closely related to the issue of how the energetic particles modify the structure of the shock itself. We will parameterize the injection rate by an efficiency factor representing the fraction of supernova energy which is converted to cosmic rays. Estimates based on theory and/or numerical experiment suggest that the efficiency is likely to be at least 10% \[Ellison & Eichler (1984), Blandford & Eichler (1987) Jones & Ellison (1991), Berezhko & Ellison (1999), Malkov, Diamond, & Völk (2000), Ellison (2002)\]

## 4 Constraints at the Minimum Diffusion Rate

In this section, we assume that cosmic rays are scattered by nonlinear turbulence: $`k_R(k_R)1`$. According to eqns. (4) and (6), the spatial diffusion coefficients $`D_{}`$ and $`D_{}`$ are approximately

$$D_{}D_{}\frac{v^2}{3\omega _{cr}}3.1\times 10^{16}\frac{\gamma \beta ^2}{B}\frac{A}{Z},$$
(11)
for particles with charge $`Z`$, atomic number $`A`$, speed $`\beta c`$, and Lorentz factor $`\gamma `$. Equation (11) corresponds to scattering by one gyroradius $`r_g`$ per gyroorbit, which is generally thought to be the minimum rate at which cosmic rays can diffuse. We defer the question of whether such a high level of turbulence can be maintained to §5. Here, we merely find the minimum value of $`B`$ which is consistent with the theories of cosmic ray acceleration and propagation at maximum efficiency.

In this section and the following one, we focus on mildly relativistic cosmic rays with energies of about 1 GeV/nucleon. This is slightly more than the energy required for spallation reactions (a few hundred MeV/nucleon), and the energy required to generate $`\pi `$ mesons from nuclear collisions ($``$280 MeV/nucleon). That is, we focus on the energy range required to create light elements and the $`\gamma `$-ray continuum resulting from $`\pi ^0`$ decay. We do not discuss the acceleration of electrons. They are underabundant in cosmic rays, probably because they are injected at lower efficiency, but the injection problem is beyond the scope of the present paper.

The acceleration and propagation of the MeV cosmic rays responsible for ionization and heating the ISM is not, in our view, a plasma physics problem. The energy of these ions is comparable to the energy per nucleon of supernova ejecta, and their range against inelastic collisions in neutral hydrogen is 1020 \- 1021 cm-2. Therefore, these particles probably originate as supernova debris, and decay within their first traversal of the galaxy rather than being confined for times much longer than their transit time. Interstellar gas has probably been heated and ionized by low energy cosmic rays since the first generation of supernova explosions.

### 4.1 Acceleration

Lagage & Cesarsky (1983) have derived the rate of particle acceleration by a strong parallel shock when the fluctuation amplitude is unity. Their result can be written as

$$\frac{dE}{dt}=\frac{3}{20}\frac{ZeB}{c}u_S^2=1.5\times 10^{18}ZBu_S^2\mathrm{GeV}\mathrm{s}^1,$$
(12)
where $`B`$ and $`u_S`$ are expressed in G and $`\mathrm{cm}\mathrm{s}^1`$, respectively.

Equation (12) can be integrated in time for any model of the shock velocity $`u_S(t)`$. LC found that most of the energy gain occurs during the blast phase, when the remnant is essentially undecelerated and $`u_S`$ is constant, and the energy conserving Sedov phase, when $`u_St^{3/5}`$. The duration $`t_b`$ of the blast phase is

$$t_b=6.2\times 10^9\frac{M_e^{5/6}}{n^{1/3}E_{51}^{1/2}}\mathrm{s},$$
(13)
where $`M_e`$ is the mass of the ejecta in $`M_{}`$, $`n`$ is the ambient interstellar number density in cm-3, and $`E_{51}`$ is the energy released in the supernova explosion, in units of 1051 erg. From eqn. (12),

$$E(t_b)=9.1\times 10^9\frac{E_{51}^{1/2}}{M_e^{1/6}n^{1/3}}ZB\mathrm{GeV}.$$
(14)
The energy gained between the beginning of the Sedov phase and its termination by radiative cooling at time $`t_r`$ can be written in terms of the shock radii $`R_b`$ and $`R_r`$ at times $`t_b`$ and $`t_r`$

$$\mathrm{\Delta }E_S=1.3\times 10^{10}\frac{E_{51}^{1/2}ZB}{M_e^{1/6}n^{1/3}}\left(1\left(\frac{R_b}{R_r}\right)^{1/2}\right)\mathrm{GeV}.$$
(15)
If the metallicity of the ambient medium is low, $`R_r`$ will be somewhat larger than its contemporary value, somewhat increasing the energy gain. In any case, Equations (14) and (15) show that in order to accelerate protons to relativistic energies, $`E\mathrm{GeV}`$, $`B`$ must be at least $`10^{10}`$G.

Since Type II supernovae are observed to be clustered, we also consider a shock driven by expansion of a superbubble driven by a constant (radiative plus mechanical) luminosity $`L_0`$ \[Shull & Saken (1995)\]. In this case, $`u_St^{2/5}`$ and the bubble radius $`R_St^{3/5}`$. The particle energy $`E_{sb}`$ increases steadily with time, and can be written in terms of the bubble radius $`R_S`$ as

$$E_{sb}=8.2\times 10^8\left(\frac{L_{38}R_{pc}}{n}\right)^{1/3}ZB\mathrm{GeV},$$
(16)
where $`L_{38}`$ and $`R_{pc}`$ are $`L_0`$ in units of 1038erg s-1 and the bubble radius in pc, respectively. Equation (16) shows that particle acceleration by superbubbles is generally weaker than acceleration driven by single supernovae. A similar conclusion was reached by Bykov & Fleishman (1992) on the basis of a somewhat different acceleration model.

If the upstream magnetic field lies nearly in the plane of the shock, instead of nearly perpendicular to it, particles drift along the shock face while scattering back and forth across it. The maximum energy $`E_d`$ which can be acquired by drifting a distance $`R`$ is of order (Jokipii, 1987)

$$E_d\frac{ZeB}{c}u_SR=30ZBu_SR_{pc}\mathrm{GeV}.$$
(17)
If $`u_S10^9\mathrm{cm}\mathrm{s}^1`$, and $`R`$ 1 pc, $`B`$ must again be about 10-10G in order to accelerate protons to 1 GeV.

The requirement that the particle gyroradius must be less than the radius of curvature of the shock also imposes a lower limit on the magnetic fieldstrength. However, the gyroradius of a GeV proton in a 10-10G magnetic field is $``$0.01 pc, so the constraint arising from the acceleration rate is more severe.

### 4.2 Confinement

The constraint that the cosmic ray gyroradius must be less than the typical magnetic lengthscale $`L^B`$ in the galaxy is thought to set an upper bound on the energies of cosmic ray nuclei in contemporary galaxies, and has led to the inference that the highest energy cosmic rays are extragalactic in origin. This limit arises because if $`r_g>L^B`$, the cross field drift velocity becomes comparable to the particle velocity itself. The constraint derived from setting $`r_g=L^B`$ is

$$B>10^{15}\frac{A}{Z}\frac{\sqrt{\gamma ^21}}{L_{kpc}^B},$$
(18)
where $`L^B`$ is expressed in kpc.

The limit on $`B`$ derived from eqn. (18) is well below the value derived from the acceleration time \[eqns. (14), (15), (16), or (17)\], provided that $`L^B`$ is comparable to the scale of the galaxy itself. However, it is a nontrivial constraint on those dynamo theories which predict that much of the power in the magnetic field is on small scales, especially when the field is still very weak \[Kulsrud & Anderson (1992), Schekochihin et al. (2002)\]. If the power spectrum in the magnetic field peaked at the Ohmic scale, or even the viscous scale, a much stronger field would be required to confine cosmic rays.

If we require that cosmic rays are confined within superbubbles long enough to produce light elements by spallation reactions, the resulting condition on $`B`$ is more stringent 111Unraveling the history of light elements from measurements of their relative abundances in the oldest stars is a complex matter which may require more than one type of spallation source; see Parizot & Drury (1999), Parizot & Drury (2000), Suzuki & Inoue (2002). We provide only a rough estimate here.. Using eqns. (11) and (7), the resulting constraint on $`B`$ for particles with Lorentz factor $`\gamma `$ and speed $`v=\beta c`$ is

$$B>1.0\times 10^7\frac{A}{Z}\frac{\gamma \beta ^2\tau _6}{R_{pc}^2}G,$$
(19)
where $`\tau _6`$ is $`\tau _c`$ expressed in units of 106 yr. Equation (19) shows that a field of 10-11 G is sufficient to confine GeV cosmic rays for 106 yr in a 100 pc superbubble, provided that the turbulence is strong. This is comparable to, but slightly weaker than, the field required for acceleration of such cosmic rays.

Equations (8) and (11) lead to an expression for $`B`$ in terms of the energy density in cosmic rays $`U_{cr}`$ and the power density of the cosmic ray source $`𝒫_{cr}`$, which we write as a fraction $`ϵ`$ of the supernova power density $`𝒫_{sn}`$. Taking $`\gamma \beta ^2=1`$ as average, setting $`A=Z=1`$ because most cosmic rays are protons, normalizing $`U_{cr}`$ and $`P_{sn}`$ by their contemporary values: $`U_{cr}10^9u_{cr}\mathrm{GeV}\mathrm{cm}^3`$, $`𝒫_{cr}10^{22}p_{cr}ϵ\mathrm{GeV}\mathrm{cm}^3\mathrm{s}`$, and expressing the confinement lengthscale $`L`$ in kpc gives

$$B=3.3\times 10^{14}\frac{u_{cr}}{ϵL_{kpc}^2p_{cr}}G.$$
(20)
Nonlinear models of cosmic ray acceleration in shocks suggest that $`ϵ`$ 0.1. If all the other parameters take their fiducial values, a field of about 10-12 \- 10-13G would sustain a cosmic ray energy density comparable to that found in contemporary galaxies. An even lower fieldstrength would have sufficed if the supernova rate was much higher at these early times than it is now.

## 5 Wave Propagation, Excitation, and Damping in Weak Fields

In this section we consider the propagation, excitation, and damping of hydromagnetic fluctuations at very low magnetic fieldstrength, and attempt to estimate whether early galaxies could have sustained a nonlinear level of magnetic fluctuations.

At the magnetic fieldstrengths considered here, the frequencies of Alfvén waves at the resonant wavenumbers are extremely low. The Alfvén speed $`v_{\text{Ai}}`$ defined with respect to the ionized mass density is $`v_{\text{Ai}}B/(4\pi m_in_i)^{1/2}=2.2\times 10^{11}B/n_i^{1/2}`$cm s-1. Motivated by the resonance condition eqn. (1), we define a fundamental wavenumber $`k_0`$ in terms of the nonrelativistic proton gyrofrequency $`\omega _{cp}`$

$$k_0\frac{\omega _{cp}}{c}=3.5\times 10^7B\mathrm{cm}^1.$$
(21)
The Alfvén wave frequency is

$$\omega k_0v_{\text{Ai}}=7.3\times 10^4\frac{B^2}{n_{i}^{}{}_{}{}^{1/2}}\mathrm{s}^1.$$
(22)
For example, if $`B10^{10}G`$ and $`n_i1`$, the wave period is more than 108 yr.

The implication of the very long Alfvén wave period is that the waves are strongly affected by cosmic rays. Under contemporary interstellar conditions, cosmic rays influence hydromagnetic waves only through the resonant interaction, described by eqn. (1). This means that only the imaginary part of the wave frequency is affected by the cosmic rays. But when the magnetic field is weak, the electromagnetic response of the plasma is so feeble that even a low density of cosmic rays has a significant effect on both the real and imaginary parts of $`\omega `$. A similar situation prevails if a very high flux of cosmic rays is present \[Kulsrud & Zweibel (1975), Zweibel (1979)\].

At low magnetic fieldstrength, Alfvén waves are also affected by the thermal plasma. Foote & Kulsrud (1979) have shown that the the large thermal ion gyroradius causes the ions to deviate from the $`𝑬\times 𝑩`$ drift which nearly exactly characterizes their motion in a cold plasma, leading to modifications to the Alvén wave dispersion relation. Thermal plasma also causes efficient damping outside a small cone of propagation angles centered on $`𝑩_0`$. It turns out that with the lower bounds on $`B`$ imposed in the previous section, the thermal modifications to the real part of $`\omega `$ are not as important as the cosmic ray modifications. We justify this statement below. Collisionless damping is potentially important, however, and we consider it in §5.2. Because the collisionless damping increases with the propagation angle $`\theta `$ as $`\theta ^2`$, while the excitation rate due to cosmic ray streaming is relatively insensitive to $`\theta `$ for small $`\theta `$, we consider only excitation of waves with $`\theta =0`$.

Much of the interstellar gas in early galaxies may have been only weakly ionized, leading to damping of the waves by ion-neutral friction. We consider weakly ionized gases in §5.3.

In §2 we mentioned the possibility that young galaxies do not have large scale magnetic fields. In discussing wave propagation, we tacitly assume that there is a well ordered field on a scale at least several times larger than the wavelength of the resonant wave; $`\lambda 6\times 10^{12}B^1`$ pc. Thus, a field of 10-12G would have to be coherent over tens of pc in order to support waves. This is much less than the current coherence length of the Galactic magnetic field, which is at least a few kpc. However, the possibility remains that even this degree of coherence is lacking. In this case, it seems unlikely that the particles can be confined at all; see eqn. (18).

### 5.1 Dispersion relation with cosmic rays

We assume the cosmic rays are protons. As the presumptive dominant species, they control collective effects. The dispersion relation for hydromagnetic waves in a medium with cold protons and cosmic ray protons drifting relative to each other with speed $`v_D`$, and with electrons cancelling the charge and current of each population of protons, is derived in Appendix A \[eqn. (A10)\]. The dispersion relation is modified by cosmic rays because of their drift and by the same mechanism that modifies the dispersion relation in a hot plasma: the electron motion in the wave is nearly the $`𝑬\times 𝑩`$ drift, but the cosmic rays deviate from the $`𝑬\times 𝑩`$ drift because of their large gyroradii. Therefore, the currents associated with perturbations of the two species do not come close to cancelling, and even at low density the cosmic rays have a strong effect on the waves.

In the rest frame of the thermal plasma, the dispersion relation is

$$\omega ^2+\omega _{cp}\frac{n_{cr}}{n_i}\zeta _{lr}(k)\left(\omega kv_D\right)k^2v_{\text{Ai}}^2=0.$$
(23)
The complex functions $`\zeta _{lr}(k)`$ (the subscript denotes the sign of circular polarization) are defined in Appendix A and plotted for a sample cosmic ray distribution function in Figure 1. The real and imaginary parts of $`\zeta `$ represent the contributions of the nonresonant and resonant cosmic rays, respectively, and are of similar magnitude for cosmic rays near the mean momentum of the distribution. Because the $`\zeta _{lr}`$ cannot in general be easily evaluated, we treat $`\zeta _{lr}`$ as a parameter of order unity in what follows, except in computing the group velocity \[see eqn. (39)\]. If $`v_D=v_{\text{Ai}}`$, $`\omega `$ is real, and the waves are marginally stable.

It is convenient to scale the wavenumber by

$$k=yk_0,$$
(24)
where $`k_0`$ is defined in eqn. (21), and to absorb the ratio of cosmic ray to thermal plasma density into the $`\zeta _{lr}`$ by introducing

$$ϵ_{lr}\zeta _{lr}\frac{n_{cr}}{n_i}.$$
(25)
With these definitions, the general solution of eqn. (23) is

$$\frac{\omega }{\omega _{cp}}=\frac{1}{2}\left[ϵ_{lr}\pm \left(ϵ_{lr}^2+4y^2\frac{v_{\text{Ai}}^2}{c^2}+4ϵ_{lr}y\frac{v_D}{c}\right)^{1/2}\right].$$
(26)
The properties of eqn. (26) depends on the parameter

$$X\frac{ϵ_{lr}}{y}\frac{c}{v_{\text{Ai}}}\frac{v_D}{v_{\text{Ai}}}1.5\times 10^{10}\left(\frac{10^9n_{cr}\zeta _{lr}}{yn_i^{1/2}B}\right)\left(\frac{v_D}{v_{\text{Ai}}}\right)$$
(27)
which is the ratio of the third term under the radical to the second term. This parameter has a more physical interpretation: if we assume that most of the cosmic ray pressure is provided by mildly relativistic particles ($`\gamma 1`$) and use eqns. (25) and (27), then the ratio of the pressure in the mean magnetic field to the cosmic ray pressure can be written approximately as

$$\frac{B^2}{8\pi P_{cr}}\frac{1}{2X}\frac{v_D}{c}.$$
(28)
Equation (28) shows that if $`X`$ is large, the magnetic field must be far below equipartition with the cosmic rays.

Under contemporary interstellar conditions, $`X1`$ unless the cosmic ray flux is much larger than average (Zweibel, 1979), and eqn. (26) reduces to

$$\frac{\omega }{\omega _{cp}}=\pm y\frac{v_{\text{Ai}}}{c}+\frac{ϵ_{lr}}{2}\left(1\pm \frac{v_D}{v_{\text{Ai}}}\right).$$
(29)
Equation (29) agrees with the expression for the growth rate given in Kulsrud & Cesarsky (1971).

The drift velocity $`v_D`$, diffusion coefficient $`D`$, and cosmic ray gradient lengthscale $`L^c`$ are related by $`v_DD/L^c`$. Setting $`D`$ equal to is minimum value \[eqn. (11)\] gives the minimum possible value of $`v_D/v_{\text{Ai}}`$

$$\frac{v_D}{v_{\text{Ai}}}\frac{\lambda }{L^c}\frac{c}{v_{\text{Ai}}}4.7\times 10^{17}\frac{n_i^{1/2}}{L_{kpc}^cB^2},$$
(30)
where in the last step $`L^c`$ is expressed in kpc. At the fieldstrengths of interest here, $`B10^{10}`$G, the anisotropy produced by cosmic ray sources drives the instability far above threshold. Under these conditions, the solution of eqn. (26) is approximately

$$\frac{\omega }{\omega _{cp}}=\pm \left(ϵ_{lr}y\frac{v_D}{c}\right)^{1/2}.$$
(31)
The real and imaginary parts of eqn. (31) are comparable in size, and the wave periods are much shorter than those implied by eqn. (22). Evaluating eqn. (31) numerically gives

$$\omega =3.2\times 10^1\frac{B}{n_i^{1/2}}\left(\frac{v_D}{c}\right)^{1/2}\left(10^9n_{cr}y\zeta _{lr}\right)^{1/2}.$$
(32)
Although the turbulence required for diffusive shock acceleration might always be present in the ambient medium, it is of interest to investigate the circumstances under which the streaming instability grows quickly enough for the cosmic rays to confine themselves. A rough estimate of $`n_{cr}/n_i`$ follows from the efficiency assumed for shock acceleration and the ratio of the mean cosmic ray energy to the postshock thermal energy: if 10% of the shock energy is converted to cosmic rays and the energy ratio is 105, then $`n_{cr}/n_i10^6`$. If $`v_Du_Sc/30`$, eqn. (32) gives $`\omega 1.9(y\zeta _{lr})^{1/2}B`$. If $`B10^{10}`$ G, which is approximately the limit implied by eqns. (14), (15), and (17), the growth time of the waves is a few hundred years. Thus, it is not implausible that cosmic rays in the vicinity of a strong shock can trap themselves.

In order to investigate confinement on the galactic scale, we express $`v_D`$ in terms of the confinement time $`\tau _7`$, the confinement time $`\tau _c`$ in units of 107 yr, and $`L_{kpc}`$, the confinement lengthscale in units of kpc. If $`ϵ_{lr}10^9`$, $`\omega \tau _c10^{12}B`$, which exceeds unity as long as $`B>10^{12}`$G. This is comparable to, although slightly larger than, the field required to maintain cosmic rays at their present energy density \[eqn. (20)\].

Finally, we return to the effect of thermal plasma on the wave frequency. According to the results of Foote & Kulsrud (1979), the ion gyroradius correction is represented by adding the term

$$\pm \omega \omega _{cp}\frac{y^2}{2}\frac{v_i^2}{c^2}$$
(33)
to eqn. (23), where the $`\pm `$ signs again represent left and right circular polarization, respectively, and $`v_i^22k_BT/m_i=1.6\times 10^8T`$. The relative effect of this thermal term on the dispersion relation is measured by the quantity

$$\frac{y^3}{16ϵ_{lr}}\left(\frac{v_i}{c}\right)^4\frac{c}{v_D}3.4\times 10^{26}T^2\left(\frac{n_i}{n_{cr}}\right)\left(\frac{c}{v_D}\right).$$
(34)
The right hand side of eqn. (34) is typically much less than unity, so we neglect thermal corrections to $`Re(\omega )`$.

### 5.2 Landau damping

Unless the angle $`\theta `$ between the wave vector and $`𝑩_0`$ is exactly zero, the wave electric field has a component $`E_{}`$ parallel to $`B_0`$. Thermal ions which satisfy the Landau resonance condition $`\omega =k_{}v_{}`$ are steadily accelerated by $`E_{}`$ and absorb energy from the wave. This process is known as Landau damping, and is the primary linear damping mechanism in fully ionized regions.

As we argue below, the condition $`\theta 0`$ is unlikely to be fulfilled because of inhomogeneity in the background density and magnetic field. Therefore, $`\theta 0`$ is the generic case, and we must consider damping.

Foote & Kulsrud (1979) evaluated the rate of Landau damping in a hot thermal plasma, with $`v_i^2/v_{\text{Ai}}^21`$. We show in Appendix B that their result is unchanged by the presence of cosmic rays. The damping rate $`\mathrm{\Gamma }_L`$ is given to within a factor of order unity by

$$\mathrm{\Gamma }_L=\frac{\sqrt{\pi }}{2}k_{}v_i\theta ^2,$$
(35)
where we have approximated $`\mathrm{tan}\theta `$ by $`\theta `$ \[see eqn. (B15)\] within the range of interest.

We calculate the wave damping time in two ways. First, we assume that the waves are excited at $`\theta =0`$, but $`\theta `$ grows as the waves propagate through the inhomogeneous background. We define the Landau damping time $`t_L`$ by the condition

$$_0^{t_L}\mathrm{\Gamma }_L𝑑t=1,$$
(36)
where $`t`$ is the time along a ray path.

Let $`L^b`$ denote the gradient lengthscale of the background, and $`v_g`$ the group velocity of the wave; $`v_gd\omega /dk`$. The evolution equation for $`\theta `$ is

$$\frac{d\theta }{dt}=\frac{v_g}{L^b},$$
(37)
with the solution $`\theta =v_gt/L`$ starting from zero initial condition. The solution of eqn. (36) is

$$t_L=\left(\frac{6}{\sqrt{\pi }k_{}v_i}\frac{(L^b)^2}{v_g^2}\right)^{1/3}.$$
(38)
According to eqns. (31) and (25),

$$v_g=\left(\chi _{lr}\frac{n_{cr}}{n_i}\frac{\omega _{cp}v_D}{4k_{}}\right)^{1/2},$$
(39)
where $`\chi _{lr}^{1/2}Re[\zeta _{lr}^{1/2}d(k\zeta _{lr})/dk]`$. Using eqn. (39), we rewrite eqn. (38) as

$$t_L=\left(\frac{24(L^b)^2}{\sqrt{\pi }v_Dv_i\chi _{lr}\omega _{cp}}\right)^{1/3}3.0\times 10^9\left(\frac{(L_{pc}^b)^2n_i}{10^9n_{cr}\chi _{lr}T^{1/2}B}\frac{c}{v_D}\right)^{1/3}s,$$
(40)
where $`L_{pc}^b`$ is $`L^b`$ expressed in pc.

Damping and excitation balance when $`\omega t_L1`$. Using eqns. (32) and (40) yields the condition

$$1.1\times 10^9\left(BL_{pc}^b\right)^{2/3}\frac{\left(\zeta _{lr}y\right)^{1/2}}{\chi _{lr}^{1/3}}\left(\frac{10^9n_{cr}}{n_iT}\frac{v_D}{c}\right)^{1/6}1.$$
(41)
Evaluating the right hand side of eqn. (41) with $`v_D`$ set equal to its minimum value given in eqn. (30) yields

$$\omega t_L|_{min}=2.9\times 10^6B^{1/2}\frac{(L_{pc}^b)^{2/3}}{(L_{kpc}^c)^{1/6}}\left(\frac{10^9n_{cr}\zeta _{lr}}{n_iT}\right)^{1/6}.$$
(42)
If $`\omega t_L|_{min}`$ exceeds unity, Landau damping is always weaker than excitation by cosmic ray anisotropy. This is likely to be the case for $`B10^{11}10^{10}`$G. Presumably, the wave amplitude is limited to $`B_1/B_01`$ by a nonlinear mechanism.

As an alternative method of calculating the damping, we suppose that the initial condition $`\theta =0`$ can never be realized in an inhomogeneous system. If instead of solving eqn. (37) we simply set $`\theta =(k_{}L)^1`$ and define the damping time $`t_L^{}`$ as $`\mathrm{\Gamma }_L^1`$, with $`\mathrm{\Gamma }_L`$ taken from eqn. (35), we recover the condition on $`v_D`$ given by eqn. (41) with the factor of 1.1. reduced by about a factor of 2. The insensitivity of the damping rate to the initial conditions justifies our approximate treatment of the problem, and reinforces the conclusion that Landau damping by the thermal plasma is weak compared with excitation by cosmic rays.

### 5.3 Ion-neutral damping

Large volumes of the interstellar gas in young galaxies may be only weakly ionized. In such regions, hydromagnetic waves are damped by ion-neutral friction. Under present conditions, frictional damping is so strong that cosmic ray stream through H I regions virtually without scattering (Kulsrud & Cesarsky, 1971).

The rate $`\mathrm{\Gamma }_{in}`$ at which Alfvén waves are damped by friction is given in Kulsrud & Pearce (1969). The modification of the waves by cosmic rays significantly reduces the damping rate. As shown in Appendix C, most of the energy in the waves is in the cosmic rays themselves, with a smaller fraction in electromagnetic fields, and the energy carried by the thermal gas a distant third. Since the ions are coupled to the wave by electromagnetic forces, while collisions with the neutrals directly tap only the small amount of energy carried by the ions, the overall effect of ion-neutral friction is doubly weakened.

The damping rate is calculated by an energy method in Appendix C. The result can be written in terms of the ion-neutral collision frequency $`\nu _{in}\rho _n\sigma v/(m_i+m_n)`$ and the phase velocity $`v_\varphi \omega /k`$ as

$$\mathrm{\Gamma }_{in}=\nu _{in}\frac{v_{\text{Ai}}^4}{v_\varphi ^2\left(v_\varphi ^2+v_{\text{Ai}}^2\right)}.$$
(43)
We assumed in deriving eqn. (43) that $`\omega >\nu _{in}`$. In the parameter regime of interest, this appears to be well justified.

If the waves are Alfvèn waves, $`v_\varphi =v_{\text{Ai}}`$, and eqn. (43) reduces to $`\nu _{in}/2`$. This is the usual expression for the damping rate of high frequency Alfvén waves in a weakly ionized gas (Kulsrud & Pearce, 1969). According to eqns. (27) and (31),

$$\frac{v_\varphi }{v_{\text{Ai}}}X^{1/2}1$$
(44)
in the regime of interest. Not only is the frictional damping rate a factor of $`X^2`$ lower than it is for Alfvèn waves, it actually decreases as $`v_D/c`$ increases. Therefore, the final streaming rate cannot be determined by balancing ion-neutral friction against excitation by cosmic ray anisotropy, and the ionization state of the medium has little relevance to cosmic ray acceleration or confinement.

## 6 Summary and Conclusions

The light elements present in the oldest Galactic halo stars, and the spectrum of the diffuse $`\gamma `$-ray background, are indirect evidence that cosmic rays were present in young galaxies. Cosmic rays cannot be accelerated or confined without some level of magnetization. In this paper, we investigated just how strong a magnetic field is required, and what the presence of cosmic rays implies about the evolution of galactic magnetic fields.

Cosmic rays can be confined in appreciable numbers only if the magnetic coherence length $`L^B`$ is larger than their gyroradii. This is the standard argument for the extragalactic origin of the cosmic rays with energies above $`10^9`$ GeV. Generalizing the argument to early times leads to a lower limit on the product $`BL^B`$ \[eqn. (18)\]. If $`L^B`$ were as large as a kpc, typical of the gradient lengthscale in the galaxy itself, then $`B`$ must have exceeded 10-15G if GeV protons were confined. If the magnetic energy density was concentrated at a much smaller scale such as the the Ohmic dissipation length, cosmic rays probably could not have been confined at all.

Under contemporary conditions, cosmic rays are scattered by gyroresonant interactions with small amplitude ($`B_1/B_010^310^4`$) hydromagnetic turbulence \[eqn. (1)\], and propagate diffusively. Although the turbulence can in principle be excited by the streaming anisotropy of the cosmic rays themselves, it is also strongly damped, and it appears that supplemental sources of turbulence are required222If the turbulence is highly anisotropic, so that the efficiency of scattering is reduced, then processes other than gyroresonant scattering, such as magnetic mirroring, may be necessary to explain confinement (Chandran, 2000b). The excitation of waves by cosmic rays is, however, important in the vicinity of strong shock waves, where it is believed that the bulk of galactic cosmic rays are accelerated.

If cosmic rays in early galaxies propagated diffusively, the limits on the fieldstrength are more stringent. We derived lower bounds on the fieldstrengths necessary for diffusive shock acceleration and for confinement by assuming that the scattering frequency is the gyrofrequency. This corresponds to a nonlinear level of turbulence \[$`B_1/B_01`$; eqn. (2)\]. We used the upper bound on the acceleration rate derived by Lagage & Cesarsky (1983) to estimate the minimum magnetic fieldstrength required to accelerate protons to relativistic energies \[eqns. (14), (15), (16), and (17)\] in supernova and superbubble driven shocks, and found it to be about 10-10 G. This turns out to be similar to the minimum fieldstrength required to confine cosmic rays to superbubbles long enough for them to undergo nuclear collisions with ambient material and synthesize the light elements observed in the oldest stars \[eqn. (19)\]. At this fieldstrength, the growth rate of the streaming instability is fast enough that the cosmic rays may be able to trap themselves in the vicinity of the shock front \[see discussion following eqn. (32)\].

The magnetic fieldstrengths in regions of active star formation may not have been representative of the global galactic fieldstrength (see the discussion in §2). Therefore we derived a relationship between the average interstellar energy density in cosmic rays, the power in cosmic ray sources, and the global confinement time, assuming the maximum scattering rate. If the intensity of cosmic ray sources and the energy density in cosmic rays were comparable to their current values, the magnetic field could have been as low as 10-12 \- 10-13 G. The fieldstrength scales inversely with the luminosity of the sources, so an enhanced supernova rate would have permitted an even smaller magnetic field. A field of this strength is consistent with growth of the waves on a timescale less than the nominal confinement time of 107 yr, provided that the drift velocity $`v_D`$ exceeds $``$100 km s-1.

Estimates of the fieldstrength based on cosmic ray diffusion theory scale with the turbulent amplitude as $`(B_1/B_0)^2`$ \[see eqn. (2)\]. There are two reasons why the turbulence could have been much stronger at early times than it is today. The sources extraneous to cosmic rays range from stellar radiative and mechanical luminosity to large scale dynamical effects associated with differential rotation and self gravity \[Sellwood & Balbus (1999), Kim & Ostriker (2000)\]. Under present conditions, the gyroradius of a GeV cosmic ray is about 1 AU \[see eqn. (21)\], much smaller than the scale at which turbulent energy is injected. But at fieldstrengths of 10-12 G, the gyroradii were at parsec scales, closer to the energy injection scale, at which the turbulent amplitude is significantly larger.

The strength of these turbulent sources, however, can only be speculated upon. Turbulence is driven by the anisotropy of the cosmic rays themselves. The anisotropy is an inevitable consequence of the discrete and inhomogeneous distribution of cosmic ray sources, and the finite size of the galaxy. In §5 we considered the propagation, excitation, and damping of low frequency electromagnetic waves in a medium consisting of cosmic rays, thermal plasma, and a weak magnetic field. We found that the circularly polarized, parallel propagating Alfvén waves which efficiently scatter cosmic rays in contemporary galaxies, and which are excited by cosmic ray streaming anisotropy, are dramatically modified because the electrons $`𝑬\times 𝑩`$ drift in the wave electromagnetic fields but the cosmic rays do not \[Appendix A and eqn. (27)\]. These modifications are only significant when the pressure in the mean field is much less than the cosmic ray pressure \[see eqn. (27)\]. We found that if the drift velocity $`v_D`$ is far above the instability threshold $`v_{\text{Ai}}`$, the real and imaginary parts of the wave frequency are comparable, and much larger than the Alfvén frequency \[see eqn. (31)\]. In contrast to the current situation, neither Landau damping in fully ionized regions (§5.2 and Appendix B) nor frictional damping in weakly ionized regions (§5.3 and Appendix C) appears capable of competing with the excitation rate. Therefore, we predict that the waves grow to nonlinear amplitudes, and the cosmic rays are scattered at nearly the maximum possible rate. The limits on $`B`$ obtained under the assumption of scattering at the maximum rate are probably close to correct.

Cosmic rays will interact dynamically with the thermal gas as long as they are strongly scattered \[see eqns. (9) and (10)\]. The rate of energy transfer is usually assumed to be proportional to the Alfvén speed under standard conditions, but is much larger in the present case because the turbulence which scatters the cosmic rays is highly super-Alfvénic \[see eqn. (44)\]. Cosmic rays could potentially play an important role in driving outflows from young galaxies.

The fact remains that magnetic fields and cosmic rays are in approximate equipartition in our own Galaxy at the present time. We suggest that this is so because they share a common energy source: supernovae. Given enough time, magnetic fields, cosmic rays, and for that matter turbulent bulk motions will probably achieve a state near equipartition, as currently observed in the Galaxy. However, as we have shown in this paper, this is not required by anything fundamental to the physics of cosmic ray acceleration or confinement. Why the magnetic field and cosmic ray energy density have saturated at values close to the turbulent energy density, instead of increasing beyond it, is not completely clear. The reason may be that if the magnetic field or cosmic ray energy density were appreciably larger, the vertical stratification of the interstellar medium would be unstable to buoyancy driven instabilities (Parker, 1966), which would probably lead to separation of the thermal and nonthermal components of the interstellar medium, and possibly facilitiate magnetic field and cosmic ray escape. We have briefly considered the role of these buoyancy instabilities at early times, and find it to be ambiguous. If the ratio of cosmic ray pressure to thermal gas pressure is denoted by $`\alpha `$, the effective polytropic exponents of the thermal gas and cosmic rays are $`\gamma _g`$ and $`\gamma _{cr}`$, and magnetic pressure is negligible, then the system is unstable if $`\gamma _g1+\alpha (\gamma _{cr}1)<0`$ (Zweibel & Kulsrud, 1975). Thus, if the cosmic rays are well coupled by scattering and $`\gamma _{cr}=4/3`$, they can stabilize the system, but if they diffuse rapidly, they are destabilizing. While a detailed investigation of this problem is beyond the scope of this paper, it is safe to say that if the instability occurs, its nonlinear development differs from that of the instability in a strong magnetic field.

If cosmic rays can be accelerated and confined in the presence of magnetic fields which are several orders of magnitude weaker than the fields in contemporary galaxies then evidence for cosmic rays in young galaxies is not evidence for magnetic fields at the equipartition level. The fieldstrengths required are several orders of magnitude larger than the seed fields produced by mechanisms operating on large scales, but are representative of the fields expected from a superposition of many small sources, such as plerion supernova remnants. The dominance of small scale fields predicted by some dynamo theories is inconsistent with the presence of cosmic rays, unless the field is strong enough that the cosmic ray gyroradius is less than the size of the system. Thus, the lower bounds on $`B`$ implied by the presence of cosmic rays impose meaningful constraints on theories of the origin and evolution of galactic magnetic fields (§2).

Finally, the drastic changes in the properties of hydromagnetic turbulence brought about by cosmic rays may be of some importance beyond the present application, although they fall outside its immediate domain. It has been pointed out elsewhere that cosmic rays can be a significant source of hydromagnetic fluctuations at the scale of their gyroradius, and that this may enhance the rate at which they are accelerated in shocks \[Eichler (1985), Lucek & Bell (2000), Bell & Lucek (2001)\]. Because the gyroradius scale falls between the global galactic scale and the resistive scale, fluctuations driven by cosmic rays may affect the operation of the galactic dynamo. The novel dispersion relation of these fluctuations will affect their nonlinear behavior, and will modify the properties of a turbulent cascade. These effects could operate in any environment in which the magnetic field is weak and energetic particles are present.

This paper was begun at the conference “Making Light of Gravity” held in honor of Sir Martin Rees. I am grateful to the organizers for the opportunity to speak there. Fabian Heitsch, J. Michael Shull, and especially an anonymous referee provided useful comments on the manuscript. Material support was provided by NSF Grant 0098701 to the University of Colorado.

## Appendix A. Dispersion Relation with Cosmic Rays

The dispersion relation for plane electromagnetic waves in a plasma follows from the equation for the first order electric field $`𝑬_1`$

$$𝒌\times \left(𝒌\times 𝑬_1\right)+\frac{\omega ^2}{c^2}𝑲𝑬_1=0,$$
(A1)
where the response of the plasma is contained in the dielectric tensor $`𝑲`$. Equation (A1) and general expressions for the components of $`𝑲`$ are given in texts such as Krall & Trivelpiece (1973).

We consider circularly polarized waves propagating parallel to the ambient magnetic field $`\widehat{z}B_0`$. In this case, $`K_{yy}=K_{xx}`$ and $`K_{yx}=K_{xy}`$, and eqn. (A1) reduces to

$$\left[K_{xx}\pm iK_{xy}\frac{c^2k^2}{\omega ^2}\right]\left(E_{1x}iE_{1y}\right)=0.$$
(A2)
The $`\pm `$ and $``$ signs denote the sign of circular polarization. The dispersion relation follows by setting the term in square brackets equal to zero.

The dielectric terms are

$$K_{xx}\pm iK_{xy}=1+\underset{\alpha }{}\frac{2\pi q_\alpha ^2}{\omega }_0^{\mathrm{}}_1^1\frac{p^2v\left(1\mu ^2\right)}{\omega kv\mu \pm \omega _{c\alpha }}\left[\frac{f_{0\alpha }}{p}+\left(\frac{kv}{\omega }\mu \right)\frac{1}{p}\frac{f_{0\alpha }}{\mu }\right]𝑑p𝑑\mu ,$$
(A3)
where the summation index $`\alpha `$ represents particle species and the $`f_{0\alpha }`$ are the zero order phase space distribution functions averaged over gyrophase. The $`\omega _{c\alpha }`$ are the relativistic gyrofrequencies.

We evaluate eqn. (A3) for a system with four populations of particles: cold (zero temperature) protons with number density $`n_i`$ drifting with speed $`\widehat{z}v_D`$, cold electrons with the same number density and drift speed as the ions, proton cosmic rays which are isotropic in velocity space and have number density $`n_{cr}`$, and an isotropic population of cold electrons, also with number density $`n_{cr}`$. This choice guarantees that to zero order the system is electrically neutral and there is no net current flow, consistent with the assumption of a uniform magnetic field. We have chosen to work in the frame in which the cosmic rays are isotropic, and transform the resulting dispersion relation to the rest frame of the plasma, because that is the simplest way to do the calculation.

Foote & Kulsrud (1979) considered the effect of finite temperature (but did not include cosmic rays). As we argue in §5.2, Landau damping is the most important thermal effect. It requires wave propagation at an angle to $`𝑩_0`$, and is considered in Appendix B.

The distribution functions of the cold, drifting species are

$$f_{0\alpha }=\frac{n_i}{p^2}\delta \left(pp_{D\alpha }\right)\delta \left(\mu +1\right),$$
(A4)
where $`p_{D\alpha }m_\alpha v_D`$, and $`\alpha `$ represents protons or electrons. The distribution function of the cold, stationary electrons is

$$f_{0\alpha }=\frac{1}{2}\frac{n_{cr}}{p^2}\delta (p),$$
(A5)
and we write the distribution function of the cosmic rays as

$$f_{0cr}=\frac{1}{2}n_{cr}\varphi (p),$$
(A6)
where $`\varphi (p)`$ is normalized to unity. Substituting eqns. (A4), (A5), and (A6) into eqn. (A3) and carrying out the integrations gives the dispersion relation

$`{\displaystyle \frac{c^2k^2}{\omega ^2}}`$ $`=`$ $`1{\displaystyle \frac{4\pi n_ie^2}{\omega ^2}}\left(\omega +kv_D\right)\left[{\displaystyle \frac{1}{m_e}}{\displaystyle \frac{1}{(\omega _{ce}\pm \omega \pm kv_D)}}+{\displaystyle \frac{1}{m_i}}{\displaystyle \frac{1}{(\omega _{cp}\pm \omega \pm kv_D)}}\right]`$ (A7)
$`{\displaystyle \frac{4\pi n_{cr}e^2}{\omega }}[{\displaystyle \frac{1}{m_e}}{\displaystyle \frac{1}{(\omega _{ce}\pm \omega )}}\pm {\displaystyle \frac{i\pi }{4}}{\displaystyle \frac{p_1}{m_p\omega _{cp}}}{\displaystyle _{p_1}^{\mathrm{}}}(p^2p_1^2){\displaystyle \frac{d\varphi }{dp}}dp`$
$`\pm {\displaystyle \frac{1}{m_p\omega _{cp}}}{\displaystyle \frac{𝒫}{4}}{\displaystyle _0^{\mathrm{}}}{\displaystyle _1^1}{\displaystyle \frac{p_1p^2\left(1\mu ^2\right)}{\mu p_1/p}}{\displaystyle \frac{d\varphi }{dp}}dpd\mu ],`$
where $`𝒫`$ denotes the principal part of the integral,

$$p_1\frac{m_p\omega _{cp}}{k}$$
(A8)
is the minimum momentum at which the resonance condition eqn. (1) can be satisfied, and we have dropped $`\omega `$ compared with the other terms in the denominator of the last integrals. The right hand side of eqn. (A7) is $`K_{xx}\pm iK_{xy}`$.

The “1” on the right hand side of eqn. (A7) represents the displacement current. The first square bracket represents the thermal plasma. Since the gyrofrequencies $`\omega _{cp}`$, $`\omega _{ce}`$ much exceed $`\omega `$ and $`kv_D`$, we approximate the electron term as $`1/m_e\omega _{ce}`$ and the proton term as $`1/m_p\omega _{cp}(\omega +kv_D)/m_p\omega _{cp}^2`$. The dominant terms cancel: $`\omega _{ce}m_e=\omega _{cp}m_p`$. The cancellation occurs because in the low frequency waves considered here, the electron and ion motion is nearly the $`𝑬\times 𝑩`$ drift, which gives no net current.

The second square bracket represents the cosmic rays and the cold electrons that cancel their current. The first cosmic ray integral represents the resonant cosmic rays. It is this term which leads to instability if the relative drift of the ions and thermal plasma exceeds $`v_{\text{Ai}}`$. The second term represents the nonresonant cosmic rays. If we revert to the cold plasma approximation by taking the limit $`p_1/p1`$, this term reduces to $`1/m_p\omega _{cp}`$, which cancels the electron term. The condition $`p_1/p1`$ is equivalent to the condition that the wavelength of the wave is much larger than the particle gyroradius. If this condition is not fulfilled, the particles do not $`𝑬\times 𝑩`$ drift, and their current does not cancel the cold electron current. This can lead to a relatively large effect on the dispersion relation even when $`n_{cr}/n_i1`$.

We simplify the cosmic ray term somewhat by integrating the resonant term by parts and performing the integration over $`\mu `$ in the nonresonant term. The cosmic ray integrals in the second square bracket of eqn. (A7) then reduce to

$$\frac{i\pi }{2}\frac{p_1}{m_p\omega _{cp}}_{p_1}^{\mathrm{}}p\varphi 𝑑p\pm \frac{p_1}{m_p\omega _{cp}}\frac{𝒫}{4}_0^{\mathrm{}}\left[\left(p^2p_1^2\right)\mathrm{ln}\left|\frac{1p/p_1}{1\pm p/p_1}\right|2pp_1\right]\frac{d\varphi }{dp}𝑑p.$$
(A9)
We then incorporate the electrons directly into the nonresonant cosmic ray term by replacing $`1/m_e(\omega _{ce}\pm \omega )`$ by $`1/m_p\omega _{cp}`$.

Finally, we drop the displacement current and multiply eqn. (A7) by $`\omega ^2v_{\text{Ai}}^2/c^2`$. The resulting dispersion relation is

$`k^2v_{\text{Ai}}^2`$ $`=`$ $`(\omega +kv_D)^2+\omega \omega _{cp}{\displaystyle \frac{n_{cr}}{n_i}}\times `$ (A10)
$`\left\{{\displaystyle \frac{i\pi }{2}}{\displaystyle _{p_1}^{\mathrm{}}}p_1p\varphi 𝑑p{\displaystyle \frac{p_1𝒫}{4}}{\displaystyle _0^{\mathrm{}}}\left[\left(p^2p_1^2\right)\mathrm{ln}\left|{\displaystyle \frac{1p/p_1}{1\pm p/p_1}}\right|2pp_1\pm {\displaystyle \frac{4}{3}}{\displaystyle \frac{p^3}{p_1}}\right]{\displaystyle \frac{d\varphi }{dp}}𝑑p\right\}.`$
If $`v_D0`$ and $`p/p_1`$ is small but finite for the bulk of the particles, eqn. (A10) reverts to the dispersion relation for hydromagnetic waves in a hot plasma (Foote & Kulsrud, 1979), with a factor of $`n_{cr}/n_i`$ multiplying the thermal correction term.

The quantity in $`\{\}`$ is the function $`\zeta _{lr}`$ introduced in eqn. (25). It is clear from eqn. (A10) that $`\zeta _l=\zeta _r^{}`$. The behavior of $`\zeta `$ is illustrated in Figure 1 for the normalized distribution function

$$\varphi (p)=\frac{4}{\pi p_0^3}\left(1+\frac{p^2}{p_0^2}\right)^2,$$
(A11)
which is close to the Galactic cosmic ray spectrum for $`p/p_01`$ (the observed power law index is approximately 4.6 rather than 4).

Equation (23) follows from eqn. (A10) by shifting back to the rest frame of the thermal plasma.

## Appendix B. Landau Damping Rate

We calculate the Landau damping rate using the dispersion relation derived by Foote & Kulsrud (1979), but replacing the thermal terms entering into the real part of $`\omega `$ with cosmic ray terms. The dispersion relation for small propagation angle $`\theta `$ can be written as

$$\left(K_{xx}+iK_{xy}\frac{c^2k^2}{\omega ^2}\right)\left(K_{xx}iK_{xy}\frac{c^2k^2}{\omega ^2}\right)+i\pi ^{1/2}\frac{kv_i}{\omega }\left(K_{xx}\frac{c^2k^2}{\omega ^2}\right)\frac{c^2}{v_{\text{Ai}}^2}\mathrm{tan}^2\theta =0,$$
(B12)
where $`v_i(2k_BT/m_p)^{1/2}`$, and the components of the dielectric tensor are given in Appendix A. For $`\theta =0`$, eqn. (B12) reduces to the dispersion relation eqn. (A2).

We write $`\omega =\omega _0+ϵ\omega _1`$, where $`\omega _0`$ is a solution of eqn. (B12) for $`\theta =0`$ and assume $`ϵ𝒪(\theta ^2)1`$. Expanding eqn. (B12) to first order in $`ϵ`$ and using the dispersion relation at $`\theta =0`$ gives

$$ϵ\omega _1=i\frac{\pi ^{1/2}}{2}\frac{kv_i}{\omega }\frac{c^2}{v_{\text{Ai}}^2}𝒟(\omega _0)^1\mathrm{tan}^2\theta $$
(B13)
where

$$D(\omega _0)\frac{}{\omega }\left(K_{xx}\pm iK_{xy}\right)|_{\omega _0}.$$
(B14)
If the waves are Alfvén waves, $`\omega _0^2=k^2v_{\text{Ai}}^2`$, $`D=2c^2k^2/\omega ^3`$, and

$$ϵ\omega _1=i\frac{\pi ^{1/2}}{2}\frac{kv_i}{\omega }\mathrm{tan}^2\theta ,$$
(B15)
which to order unity is the result of Foote & Kulsrud (1979). If the parameter $`X`$ defined in eqn. (27) is large, $`\omega (ϵ_{lr}\omega _{cp}kv_D)^{1/2}`$ \[see eqn. (31)\], and $`D2c^2/\omega v_{\text{Ai}}^2`$. Substituting this result into eqn. (B13) shows that $`ϵ\omega _1`$ is still given by eqn. (B15).

Foote & Kulsrud (1979) found that the waves are critically damped - the real and imaginary parts of $`\omega `$ are comparable - when $`\theta (v_{\text{Ai}}/v_i)^{1/2}`$. When the wave frequency is determined by cosmic rays, the angle for critical damping is increased to $`(cv_Dϵ_{lr}/v_i^2)^{1/4}10(v_D/c)^{1/4}(n_iT)^{1/4}`$. The damping is weakened because the thermal plasma energy is a relatively small part of the wave energy, as shown in Appendix C.

## Appendix C. Damping Rate from Ion-Neutral Friction

We calculate the damping rate $`\mathrm{\Gamma }_{in}`$ using an energy method, rather than directly solving the dispersion relation Let the energy density in waves be $`U_w`$ and the rate per unit volume at which energy is dissipated be $`\dot{U}_w`$. These quantities are related to $`\mathrm{\Gamma }_{in}`$ by

$$\dot{U}_w=2\mathrm{\Gamma }_{in}U_w.$$
(C16)
The dissipation rate can be written in terms of the first order thermal ion and neutral velocities $`𝒗_{1i}`$ and $`𝒗_{1n}`$ as

$$\dot{U}_w=\rho _i\nu _{in}𝒗_{1i}𝒗_{1n}^2.$$
(C17)
The wave energy density includes contributions from the electric and magnetic fields, thermal plasma, and cosmic rays (the contribution from the neutrals is much smaller than the plasma contribution because the wave frequency is so high that the neutrals remain nearly at rest), and can be written as (Bekefi, 1966)

$$U_w=\frac{1}{8\pi }𝑩_1^{}𝑩_1+\frac{1}{8\pi }𝑬_1^{}\frac{}{\omega }\left(\omega 𝑲\right)𝑬_1,$$
(C18)
where $`𝑲`$ is the dielectric tensor (see Appendix A). The first term in eqn. (C18) represents magnetic energy density, and the second accounts for the energy density in the electric field, thermal plasma, and cosmic rays.

We assume the waves are circularly polarized and propagate parallel to $`𝑩_0`$, which we take to define the $`\widehat{z}`$ direction, and that first order quantities vary with $`z`$ and $`t`$ as $`\mathrm{exp}i(\omega tkz)`$. The neutral velocity is determined solely by friction with the ions, and can be written in terms of $`𝒗_{1i}`$ as

$$i\omega 𝒗_{1i}=\nu _{ni}\left(𝒗_{1n}𝒗_{1i}\right),$$
(C19)
where the neutral-ion collision frequency $`\nu _{ni}=\rho _i\sigma v/(m_i+m_n)`$. The ion velocity is determined by both frictional and magnetic forces (it suffices to use a one fluid treatment of the plasma, $`𝒗_{1i}𝒗_{1e}`$, because the thermal electron and ion velocities are essentially the $`𝑩\times 𝑩`$ drift. The ion equation of motion is

$$i\omega 𝒗_{1i}=\frac{i}{4\pi \rho _i}kB_0𝑩_1+\nu _{in}\left(𝒗_{1n}𝒗_{1i}\right).$$
(C20)
Solving eqn. (C19) for $`𝒗_{1n}`$ in terms of $`𝒗_{1i}`$ and expressing $`𝒗_{1i}`$ in terms of $`𝑩_1`$, we rewrite eqn. (C17) as

$$\dot{U}_w=\nu _{in}\frac{v_{\text{Ai}}^2}{v_\varphi ^2}\frac{1}{4\pi }𝑩_1^{}𝑩_1,$$
(C21)
where we have assumed $`\omega /\nu _{in}1`$.

The second term on the right hand side of eqn. (C18) can be written as

$$\frac{1}{4\pi }E_{1x}E_{1x}^{}\frac{}{\omega }\omega \left(K_{xx}\pm iK_{xy}\right).$$
(C22)
From eqns. (A2) and (31),

$$K_{xx}\pm iK_{xy}ϵ_{lr}\frac{c^2}{v_{\text{Ai}}^2}\frac{\omega _{cp}kv_D}{\omega ^2}.$$
(C23)
It is straightforward to evaluate eqn. (C22) using eqn. (C23). We replace $`E_{1x}E_{1x}^{}`$ by $`𝑩_1𝑩_1^{}v_\varphi ^2/2c^2`$ using Faraday’s law. Equation (C18) becomes

$$U_w=\frac{1}{8\pi }𝑩_1𝑩_1^{}\left(1+\frac{v_\varphi ^2}{v_{\text{Ai}}^2}\right).$$
(C24)
The first term in eqn. (C24) represents the magnetic energy density, while the second is primarily cosmic ray energy, with a small contribution from the first order electric field.

Equation (43) follows from eqn. (C16) if we substitute for $`U_w`$ and $`\dot{U}_w`$ from eqns. (C21) and (C24).


format: text

text:
# Electric Conduction in Short DNA Wires

## Abstract

A strict method is used to calculate the current-voltage characteristics of a double-stranded DNA. A more reliable model considering the electrostatic potential drop along an individual DNA molecular wire between the contacts is considered and the corresponding Green’s Function is obtained analytically using Generating Function method, which avoids difficult numerical evaluations. The obtained results indicate that the electrostatic drop along the wire always increases the conductor beyond the threshold than without considering it, which is in agreement with recent experiments. The present method can also be used to calculate the current-voltage characteristics for other molecular wires of arbitrary length.

Recently, the progress of molecular wires A  attracts much attention on their transport behaviors. A number of experimental groups have reported measurements of the current-voltage (I-V) characteristics of molecules, especially Deoxyribonucleic Acid (DNA)AR ; KH . DNA has a special double-helix structure with complementary nucleotide base-pairs stacking onto each other, which could possibly be a candidate for one-dimensional electronic transport10 . Intense experimental investigations have already been made on the transport properties of DNAC . Using scanning tunneling microscope technique, Dunlap et al. D  found that DNA is an insulator. Fink and schoenenberger F  measured I-V characteristics of $`\lambda `$-DNA ropes consisting of a few double-stranded DNA (dsDNA) molecules by low energy electron bombardment technique, and found a linear current-voltage relation. Kasumov et al. K  measured small number of DNA molecules, and observed proximity-induced superconductivity. In particular, Porath et al. P  directly measured electrical transport through individual DNA molecules, and suggested DNA to be a semiconductor with a voltage gap. To explain experimental disputes, several theoretical models ranging from the Hükel model 6  to the density functional theory den ; Mattias  have been developed.

In the present paper, we investigate theoretically Porath’s experiment P  and calculate the current between two electrodes which are connected by a DNA molecular wire. To use a simple model illustrating the basic physics, we follow the algorithm proposed by Mujica et al. 6 . The model Hamiltonian can be written as the sum of two terms,

$$H=H^0+H^{^{}},$$
(1)
where $`H^0`$ is the non-interacting Hamiltonian for the electrodes (electron reservoirs) and DNA molecule, and $`H^{^{}}`$ represents the coupling of contacts and DNA molecule. Using extended Hükel model, $`H^0`$ and $`H^{^{}}`$ can be expressed as 6

$`H^0`$ $`=`$ $`{\displaystyle \underset{i}{}}E_i^0ii+{\displaystyle \underset{\alpha =1}{\overset{N}{}}}E_\alpha ^0\alpha \alpha +{\displaystyle \underset{f}{}}E_f^0ff+`$
$`{\displaystyle \underset{\alpha =1}{\overset{N1}{}}}\beta \alpha \alpha +1+h.c.`$
$`H^{^{}}`$ $`=`$ $`{\displaystyle \underset{i}{}}V_{i1}i1+{\displaystyle \underset{f}{}}V_{Nf}Nf+h.c.,`$ (3)
where the sum on $`i`$ ($`f`$) runs over the state in the continuum of left (right) reservoirs. A DNA molecular wire consists of N sites, with one state per site, which is denoted by $`\alpha `$. $`V_{i1}`$ ( $`V_{Nf}`$) is the tunneling parameter from the left (right) reservoir to the single electron state $`1(N)`$ in the molecule. T-matrix formalism of scattering theory T  gives differential conductance 6 :

$$g=\frac{2e^2}{\pi \mathrm{}}G_{1N}^2\mathrm{\Delta }_A(E_f)\mathrm{\Delta }_B(E_f),$$
(4)
where $`G_{1N}`$ is the $`(1N)`$ element of the matrix formula Green’s Function G, it is related to the molecule’s site $`1`$ and $`N`$. $`\mathrm{\Delta }_1`$ and $`\mathrm{\Delta }_N`$ are semielliptical reservoir spectral densities defined through Newns’ chemisorption theory N :

$$\mathrm{\Delta }_K(E)=\{\begin{array}{cc}\frac{V_K^2}{\gamma }\sqrt{1(E/2\gamma )^2},& E/2\gamma <1,\\ 0,& E/2\gamma >1,\end{array}$$
(5)
where E is measured from the center of reservoir energy band caused by chemisorption in the surface of the electrode. $`V_K`$ $`(K=1,N)`$ describe the strength of the chemisorption coupling between the DNA molecule and the reservoirs, and $`4\gamma `$ is the bandwidth of the reservoir. Through Löwdin’s matrix partition technique 15 , the electrode Hamiltonian can be replaced by a self-energy:

$`H_{eff}`$ $`=`$ $`\mathrm{\Sigma }_111+{\displaystyle \underset{\alpha =1}{\overset{N}{}}}E_\alpha ^0\alpha \alpha +{\displaystyle \underset{\alpha =1}{\overset{N1}{}}}\beta \alpha \alpha +1`$
$`+\mathrm{\Sigma }_NNN+h.c.,`$
where $`\mathrm{\Sigma }_K`$ $`(K=1,N)`$ are respectively the self-energy resulting from the coupling of the molecule to the left (right) electrode. The Green’s Function is expressed as

$$G=\frac{1}{(zIH_{eff})},$$
(7)
where z is a complex number, whose real part $`E`$ is the energy of the transfer electrons.

Mujica et al. 6  ignored the electrostatic drop along the molecule. Here we consider the linear voltage drop along the molecule, i.e., the electric potential between the electrodes varying linearly with distance. This assumption is a good approximation to the computed potential profile through the molecule between the two electrodes RE . Thus the energies of the sites are function of bias V:

$`E_\alpha ^0`$ $`=`$ $`E_bqV_0(\alpha 1),\alpha =(1,\mathrm{},N),`$ (8)
$`V_0`$ $`=`$ $`V/(N1),`$ (9)
The effective Hamiltonian can be expressed in the following matrix form:

$$H_{eff}=\left[\begin{array}{ccccc}E_b+\mathrm{\Sigma }_1\hfill & \beta \hfill & 0& \mathrm{}\hfill & \mathrm{}\hfill \\ \beta \hfill & E_bqV_0\hfill & \beta & 0\hfill & \mathrm{}\hfill \\ 0\hfill & \mathrm{}\hfill & \mathrm{}& \mathrm{}\hfill & \beta \hfill \\ 0\hfill & \mathrm{}\hfill & 0& \beta \hfill & E_b(N1)qV_0+\mathrm{\Sigma }_N\hfill \end{array}\right],$$
(10)
where $`q`$ is the average effective charge on each site. $`G_{1N}`$ can be obtained:

$$G_{1N}=\frac{(1)^{N1}\beta ^{N1}}{|zIH_{eff}|},$$
(11)
where $`|zIH_{eff}|`$ denotes determinant. Then our task is to deduce the expression of $`|zIH_{eff}|`$. For convenience, we define:

$$A_n=\left[\begin{array}{ccccc}a\hfill & \beta \hfill & 0& \mathrm{}\hfill & \mathrm{}\hfill \\ \beta \hfill & aqV_0\hfill & \beta & 0\hfill & \mathrm{}\hfill \\ 0\hfill & \mathrm{}\hfill & \mathrm{}& \mathrm{}\hfill & \beta \hfill \\ 0\hfill & \mathrm{}\hfill & 0& \beta \hfill & a(n1)qV_0\hfill \end{array}\right],$$
(12)
and from linear algebra knowledge,

$$G_{1N}=\frac{(1)^{N1}\beta ^{N1}}{D_{1,N}D_{1,N1}\mathrm{\Sigma }_ND_{2,N}\mathrm{\Sigma }_1+D_{2,N1}\mathrm{\Sigma }_1\mathrm{\Sigma }_N}$$
(13)
where $`a=E_bE_F`$, i.e. the site’s energy measured from the Fermi energy of the reservoirs without the electric field. $`D_{PQ}(P=1,2;Q=N1,N)`$ is the determinant of the matrix obtained from A by taking rows and columns only in the range from P to Q 6 .

We use Generating Function method to obtain $`D_{1N}`$. Setting $`A_n=D_{1n}`$, and $`A_N=D_{1N}`$, we find that $`A_n`$ satisfies the following recursion relation:

$$A_n=[a(n1)qV_0]A_{n1}\beta ^2A_{n2}.$$
(14)
If we define

$$F(x)=\underset{n=1}{\overset{\mathrm{}}{}}A_nx^n,$$
(15)
then we obtain that

$$F(x)[1ax+\beta ^2x^2]+(qV_0)x^2F^{}(x)=ax\beta ^2x^2.$$
(16)
Since $`qV_0`$ is small, we can use a perturbation method, and regard $`qV_0`$ as a perturbation parameter. Up to the first order approximation, we have:

$`F(x)`$ $`=`$ $`F_0(x)+(qV_0)F_1(x),`$ (17)
$`F_0(x)`$ $`=`$ $`1+{\displaystyle \frac{1}{1ax+\beta ^2x^2}},`$ (18)
$`F_1(x)`$ $`=`$ $`{\displaystyle \frac{x^2F_0^{}(x)}{(1ax+\beta ^2x^2)}}.`$ (19)
Then we decompose $`F(x)`$ into a sum of several fractions, and expand them to a power series, and obtain the $`A_n`$:

$`A_n`$ $`=`$ $`A_n^{(0)}+(qV_0)A_n^{(1)},`$
$`A_n^{(0)}`$ $`=`$ $`{\displaystyle \frac{1}{2^{n+1}}}{\displaystyle \frac{(a+\sqrt{a^24\beta ^2})^{n+1}(a\sqrt{a^24\beta ^2})^{n+1}}{\sqrt{a^24\beta ^2}}},`$
$`A_n^{(1)}`$ $`=`$ $`{\displaystyle \frac{\beta ^{2n1}}{8\mu ^3}}(n1)(\lambda _1^n\lambda _2^n)[1+n(1\alpha ^2+\alpha \mu {\displaystyle \frac{\lambda _1^n+\lambda _2^n}{\lambda _1^n\lambda _2^n}})],`$
$`\alpha `$ $`=`$ $`{\displaystyle \frac{a}{2\beta }},`$ (20)
$`\mu `$ $`=`$ $`\sqrt{\alpha ^21},`$
$`\lambda _1`$ $`=`$ $`{\displaystyle \frac{\alpha \sqrt{\alpha ^21}}{\beta }},`$
$`\lambda _2`$ $`=`$ $`{\displaystyle \frac{\alpha +\sqrt{\alpha ^21}}{\beta }}.`$
In terms of our convention, $`A_N`$ is the determinant $`D_{1,N}`$. Using similar steps, we can obtain $`D_{1,N1}`$, $`D_{2,N}`$, $`D_{2,N1}`$, and then obtain $`G_{1N}`$. It is straightforward to perform similar calculations for high order perturbation.

Then, we can calculate the current through the single molecular dsDNA using wide energy band approximation. We assume $`\mathrm{\Delta }_K(K=1,N`$) to be energy independent. For convenience, we consider that the two electrodes have identical Fermi energies, which are set to zero. For DNA molecules in equilibrium, the bases of DNA are neutral, while in transport process, there are charges introduced by the contacts, and charges in DNA will be redistributed. Since the average charge $`q`$ on each site is small, the derivation of the actual potential drop along DNA molecule from the assumed linear drop can be ignored.

In previous efforts, many papers ignored the situation that the sites’ energies will vary along the molecular wire because of the voltage drop in high-intensity electric field. Some papers did consider this effect (e.g., Ref. 7 ), but they did not obtain an analytic expression of the differential conductance. The present paper achieves this task. In general, the modified Hamiltonian is more complex. Since the Green’s Function must be obtained by calculating the inverse matrix, it is not a trivial task, especially when the matrix is large.

Using the linear voltage drop approximation, we first obtain the explicit expression of $`G_{1N}`$ up to the first order perturbation, and get a more reasonable current-voltage relation. We show the I-V curve in Fig.1 for a 4-nucleotide DNA, and Fig. 2 for a 30-nucleotide DNA. We find when considering the voltage drop in the molecule, the current-voltage curves always become much steeper near the threshold. This is in accordance with the experiment of Porath et al. P : some I-V curves in their paper are much steeper near the threshold. Our result is intuitional, because voltage drop effect makes the electron transport more easily in strong electric field. Moreover, Fig. 3a of Porath’s paper contains two different I-V characteristics of the same 30-base pair DNA sample. The difference between these two curves may be the result of a sudden change in the DNA wire, possibly a conformational change. In transport process, the charge distribution along DNA may depend on DNA’s chemical nature, for example, the molecule’s conformation. Although there was a similar voltage drop in the sample of both cases, the observed electrostatic drop effect are different. Finally, our work can be used to calculate other molecular wires. Since the computational complexity of our analytical deduction is not sensitive to the number of sites, so we can calculate long chains.

An issue of interest is to investigate theoretically the transport property of DNA molecules made of heterogeneous base sequences. We are now working on this possibility.

In summary, we consider the electrostatic potential across the DNA, which may change the site energy along the DNA with the applied bias. This made the computation more difficult. We have proposed an analytical computational method to deal with this difficulty. It can be used to calculate the I-V characteristic of homogeneous DNA wires of arbitrary length. We can use the potential drop effect to explain recent single molecular DNA experiments.

The authors are grateful to Guangping Gao for helpful discussions and to Haijun Zhou for reading earlier versions of the manuscript. We finally think Prof. J. Kong for modifying the manuscript.