text:
# The effects of added salt on the second virial coefficients of the complete proteome of E. coli

## I Introduction

The genomes of a number of organisms are already known and more are being completed at a rate of perhaps one a month. Once a genome nomen  has been sequenced the amino-acid sequences of all its proteins are known. The complete set of proteins of an organism is generally referred to as its proteome. Here we use genome data for E. coli to systematically estimate the change in the interactions of all the proteins of this bacterium when the salt concentration is varied. E. coli can grow in environments with a very wide range of salt concentrations cayley91 , and so its proteins must function in vivo over a wide range of salt concentrations. The potassium ion concentration inside the cell can vary from approximately $`150`$ to $`300`$ mM cayley91 ; potassium is the predominant cation in living cells. Clearly, the proteins must remain soluble over this range, and they should bind to the other proteins which they are required to bind to in order to function, but they should not interact strongly with other proteins. The study of a proteome is often called proteomics. Here, as the physical properties of proteins are studied (as opposed to their chemical properties such as catalytic function) we are doing what may be called physical proteomics.

There has been extensive theoretical work on the salt dependence of the interactions in individual proteins, particularly for the protein lysozyme poon00 ; warren02 ; muschol97 . See Refs. poon00, ; muschol97, ; velev98, ; guo99, ; rosenbaum99, ; piazza00,  for corresponding experimental work. However, as far as the author is aware, this is the first attempt to characterise the interactions of all the proteins of an organism. We will consider the proteins separately, i.e., as single component solutions. Of course inside a bacterium the proteins exist as a mixture of thousands of components. Future work will consider multi-component mixtures of the proteins. We have chosen E. coli as it is a bacterium, and therefore a relatively simple organism, and as it has been extensively studied. However, the distribution of charges on the proteins of almost all organisms is very similar and so our results apply to almost all organisms, including H. sapiens. The only exceptions are some extremophiles runcong01 .

In the next section we use genome data to estimate the charges on the proteins of E. coli. This data is used in the third section where we calculate the variation in their second virial coefficients as the salt concentration is varied. The last section is a conclusion.

## II The charges on proteins of E. coli

E. coli K-12 has a proteome of 4358 proteins. The amino acid sequences of all of them are known from the sequencing of its genome coligen ; ebi . K-12 is the name of a strain of E. coli. Runcong and Mitaku have analysed the charge distributions of a number of organisms using a simple approximate method of estimating the charge on a protein at neutral pH from its amino-acid sequence runcong01 . We will follow their analysis but use a slightly different approximation for the charge on a protein with a given amino-acid sequence.

Of the 20 amino acids, 5 have pK values such that they should be at least partially charged at neutral pH thecell ; stryer . These are two highly acidic amino acids, aspartic acid and glutamic acid, two highly basic amino acids, lysine and arginine and one somewhat basic amino acid, histidine. Aspartic and glutamic acids have pK’s far below 7 and lysine and arginine have pK’s far above 7 and so we assume that all 4 of these amino acids are fully charged at neutral pH. Aspartic and glutamic acids then each contribute $`1`$ to the charge on a protein, and lysine and arginine each contribute $`+1`$. Histidine has a pK of around 6-6.5 (this will depend on the environment of the amino acid). The equation for the fraction $`f`$ of a basic group such as histidine that is charged at a given pH is

$$f=1/\left(1+10^{pHpK}\right),$$
(1)
where pK is the pK value for the basic group. This equation is just the Henderson-Hasselbalch equation stryer  rearranged. Taking pK$`=6.5`$ stryer  and at pH=7 we have that the fraction of histidines charged is $`f=0.24`$. As this is small we assume for simplicity that all the histidine amino acids are uncharged.

Thus, with these assumptions for the charges on these 5 amino acids, our estimate for the net charge $`Q`$ on a protein is simply given by

$$Q=n_K+n_Rn_Dn_E$$
(2)
where $`n_K`$, $`n_R`$, $`n_D`$ and $`n_E`$ are the the protein’s total numbers of lysines, arginines, aspartic acids and glutamic acids, respectively. The subscripts $`K`$, $`R`$ etc. correspond to the standard single letter codes for the amino acids stryer ; thecell . The charge $`Q`$ is in units of $`e`$ where $`e`$ is the elementary charge.

Note that Runcong and Mitaku runcong01  assume that the histidine amino acids contribute $`+1`$ to the charge, that is the only difference between our analysis and that of Runcong and Mitaku. As the histidine amino acid is quite a rare amino acid, approximately 1 in 50 amino acids is a histidine, the difference between the results we obtain and those of Runcong and Mitaku runcong01  is not large but our charges are shifted to more negative values. Using Runcong and Mitaku’s approximation the mean charge on a protein is $`7.11`$ units more positive than the mean charge we find here. As a check on our algorithm, we can compare the prediction of equation (2) for chicken lysozyme to that of a titration experiment to determine the charge. Equation (2) predicts that chicken lysozyme aalys  has a net charge of 8 at neutral pH. Titration experiments on lysozyme give a titratable charge of close to $`8.5`$ at pH=7 lysoq .

Using equation (2) we can obtain estimates for the charges of all 4358 proteins of E. coli coligen ; ebi ; gerstein98 ; mitaku99 . The results are shown in Fig. 1, where we have plotted the number of proteins $`n`$ as a function of net charge $`Q`$. The distribution is centered almost at a net charge $`Q=0`$, and for not-too-large $`|Q|`$ the distribution is roughly symmetric and Gaussian. The mean charge is $`3.15`$. Given the approximate nature of our equation for the charge on a protein, equation (2), the data is probably consistent with a mean charge of 0. The approximation scheme of Mitaku and Runcong runcong01  yields a mean charge of $`+3.96`$. Also, although when $`|Q|`$ is not too large the distribution can be seen to be reasonably symmetric, E. coli has 12 proteins with charges $`<50`$ but none with charges $`>+50`$. Excluding proteins with very large charges, $`|Q|>30`$, the root mean square charge equals $`9.1`$.

A number of other organisms, both other bacteria and eukaryotes such as yeast, have had the charge distribution on their proteomes determined by Runcong and Mitaku runcong01  and by the author unpub . Almost all of them have a roughly Gaussian distribution centered approximately at zero, like the distribution in Fig. 1. The exceptions are some extremophiles. Extremophiles are organisms that live in extreme environments, for example Halobacterium sp. lives in environments with very high levels of salt microbio . The cytosol of Halobacterium sp. contains much higher levels of potassium ions than do other organisms so perhaps it is not a surprise that the distribution of charges on its proteins is different.

We have fitted the Gaussian function

$$n(Q)=\frac{1739}{\sigma }\mathrm{exp}((Q\overline{Q})^2/2\sigma ^2).$$
(3)
to the data for the number of proteins as a function of their charge. It is drawn as the solid curve in Fig. 1. The fit parameters are mean charge $`\overline{Q}=2.16`$ and standard deviation $`\sigma =8.32`$. 1739 is $`4358/(2\pi )^{1/2}`$ and so the distribution is normalised so that its integral gives the total number of proteins. Within a couple of standard deviations of the mean the Gaussian function fits the data well but it underestimates the numbers of proteins with charges such that $`|Q\overline{Q}|`$ is several times the standard deviation.

We also note that there is a correlation between the net charge $`Q`$ on a protein and its size, measured by the number of amino acids $`M`$. Figure 2 is a scatter plot of charge $`Q`$ and number of amino acids $`M`$ for the proteins of E. coli. Although at any particular size $`M`$ there is a wide distribution of charges, on average the more highly charged proteins are larger than average. We expect the volume of a protein to scale with $`M`$.

## III Salt dependence of the second virial coefficients

Consider a dilute solution of a single one of the proteins of E. coli. Apart from water, the only other constituents are a 1:1 salt at a concentration $`c_s`$ and a buffer which controls the pH while making a negligible contribution to the ionic strength. Here we will always assume the pH=7 but other pH’s can be considered if the net charges on the proteins can be calculated. Also, the counterions of the protein are assumed the same as either the anions or cations of the salt, depending on the sign of $`Q`$. The interactions between the protein molecules in the salt solution can be characterised by means of the protein’s second virial coefficient $`B_2`$: a function of temperature, pH and salt concentration.

Proteins are complex molecules and we are unable to calculate from first principles the absolute value of $`B_2`$ for any of the 4358 proteins possessed by E. coli. However, predicting the change in the second virial coefficient when the salt concentration varies is a much easier problem, if we assume that changing the salt concentration changes only the direct electrostatic interaction between the net charges of a protein. This is a strong assumption but studies of the simple protein lysozyme have shown that the variation of its second virial coefficient can be described using a simple model which only includes its net charge poon00 ; warren02 . Here we will follow Warren warren02  and apply his analysis of lysozyme to the complete set of proteins of E. coli. We will discuss which proteins are likely to be less well described by this theory than is lysozyme.

A protein molecule of charge $`Q`$ is surrounded by its counterions and as the concentration of the protein is increased so is the counterion density. This increase in the counterion density decreases the translational entropy of the counterions and this contributes a positive amount to the second virial coefficient. See Warren warren02  and references therein for details. $`B_2`$ has the form warren02

$$B_2=B_2^{(ne)}+\frac{Q^2}{4c_s},$$
(4)
where $`B_2^{(ne)}`$ is an assumed constant term due to excluded volume interactions and other interactions which are insensitive to salt concentration. The second term is from the counterions and the salt. It is quadratic in the charge and so of course is zero for uncharged proteins and is independent of the sign of the net charge on a protein.

As stated above we are unable to calculate the absolute value of $`B_2`$; the values of $`B_2^{(ne)}`$ of the proteins are unknown. However, we can calculate the difference in $`B_2`$ when the salt concentration is changed from $`c_s=c_1`$ to $`c_s=c_2`$. It is

$$\mathrm{\Delta }B_2=\frac{Q^2}{4}\left(\frac{1}{c_2}\frac{1}{c_1}\right).$$
(5)
This is easy to calculate for any protein and in Fig. 3 we have plotted the number of proteins $`n`$ as a function of the change in their second virial coefficient, $`\mathrm{\Delta }B_2`$, when the salt concentration is decreased from 1M to $`0.2`$M. The results are given in units of nm3. For comparison the volume of a typical bacterial protein is about 60nm3 and so if a protein were to interact solely via a hard repulsion it would have a second virial coefficient of about 4 times its volume or about 240nm3.

Results for proteins with $`30Q30`$ are shown. Proteins with larger titratable charges are likely to have an effective charge lower than $`Q`$, see Refs. lobaskin01, ; warren00,  and references therein. From the linear Poisson-Boltzmann equation, the potential (divided by $`e`$) at the surface of a spherical particle with charge $`Q`$ and radius $`a`$ is $`Q\lambda _BkT/((1+\kappa a)a)`$. $`\lambda _B=e^2/(4\pi ϵkT)`$ is the Bjerrum length, and and $`\kappa ^1`$ is the Debye screening length, given by $`\kappa ^2=8\pi \lambda _Bc_s`$. For the dielectric constant of water 80 times that in vacuum and at room temperature, $`\lambda _B=0.7`$nm. Globular proteins are approximately spherical and typically have radii around 2 to 4nm. Taking a protein with a radius of 3nm, in salt at a concentration $`c_s=0.1`$M, we have that for $`Q=30`$, the potential at the surface is about $`2kT`$. Larger charges correspond to larger surface potentials and these large potentials bind oppositely charged ions to the surface reducing the effective charge. On average, this effect will be diminished to a certain extent by the fact that the most highly charged proteins are larger than average. See Fig. 2, where it is clear that the charge and size of a protein are correlated. Recent simulations by Lobaskin et al. lobaskin01  of spheres with radius 2nm and charge $`Q=60`$ in the absence of salt found an effective charge of a little under $`20`$. Thus we restrict ourselves to proteins with charges of magnitude less than or equal to 30. 4300 of the 4358 proteins, or almost 99%, have charges in this range. The mean change in $`B_2`$ of these 4300 proteins when the salt concentration is decreased from 1 to $`0.2`$M is 139nm3 and the standard deviation is 218nm3.

A couple of caveats. The first is that the effect of salt on protein solutions is known to depend not only on whether the salt is a 1:1 salt, a 1:2 salt etc. but also to the nature of ions, whether it is Mg2+ or Ca2+ for example durbin96 . Our generic theory applies only where there are no specific interactions between the salt and the protein. There is good agreement between experiment and theory for lysozyme plus NaCl warren02 ; poon00  and so we may hope that it applies to NaCl and many proteins but it clearly misses potentially important effects for other salts where there are specific protein-salt interactions. The second is that proteins are not simple charged spheres, for example some have large dipole moments. Dipoles exert net attractions which are screened and hence weakened by added salt. Thus proteins with a small charges but large dipole moments are poorly described by the current theory: if the dipole interactions are dominant then the second virial coefficient may even increase when the salt concentration is increased. Velev et al. velev98  discuss this point. Note that although we can estimate the charge on a protein from its amino-acid sequence we cannot estimate its dipole moment without knowing its three-dimensional structure, and so the sequence data from genomics is not adequate to determine dipole moments.

## IV Conclusion

Here we have shown how data from genomics can be used to estimate the charges on the proteins of an organism. We then used these charges to estimate the changes in the second virial coefficients of 4300 (99%) of the proteins of E. coli when the salt concentration is changed. Note that E. coli can survive and multiply in external environments with a very wide range of salt concentrations; Cayley et al. cayley91  studied the growth of E. coli in environments with salt concentrations ranging from very low to $`0.5`$ Molar, corresponding to potassium ion concentrations inside the cell of $`150`$ to $`300`$ mM cayley91 . Thus, studying the change in interactions of proteins with salt concentration is of direct relevance to the in vivo behaviour of proteins. Within molecular biology there is a clear shift of emphasis away from studying the proteins of an organism one or a few at a time, and towards determining the structure and function of large sets of proteins, in particular proteomes. The systematic study of these large sets of proteins is often called proteomics. This work is a first attempt to keep up with this shift by performing a simple theoretical calculation of a solution phase physical property for a complete proteome, rather than for one or a handful of proteins as is usually done. It may be termed physical proteomics.

Future work could consider mixtures of proteins, ultimately aiming to understand the cytosol of a living cell, which is a mixture of of order $`10^3`$ different types of proteins as well as DNA, RNA, ions like ATP and potassium, etc.. This is of course very complex but if in the cytosol the proteins of E. coli are present in amounts which are uncorrelated with their net charge, the mean charge of the proteins will be close to the mean of the distribution of Fig. 1. This is quite small. Neidhardt neidhardt  has taken an inventory of the species inside E. coli, and the results for charged macromolecules are shown in Table 1. The charged macromolecules in a cell are protein, DNA and the various forms of RNA: transfer RNA (tRNA), messenger RNA (mRNA) and the RNA in ribosomes (rRNA). See Refs. stryer, ; thecell,  for an introduction to the proteins, DNA and RNA. Although for every molecule of tRNA molecule there are 10 of protein, for every ribosome there are 100 molecules of protein, and for every mRNA there are 1000 molecules of protein, the contributions of the tRNA, ribosomes, DNA and proteins, to the overall charge density of the macromolecules are very roughly comparable. The ribosomes contribute the largest amount. The macromolecules are negatively charged and this negative charge is balanced by potassium ions thecell . Thus the cytosol resembles a solution of a negatively charged polyelectrolyte, except that there is not one, relatively simple, macromolecular species, but thousands of rather complex and diverse species of macromolecules.

## Acknowledgements

It is a pleasure to acknowledge discussions with J. Cuesta, D. Frenkel and P. Warren.


format: text

text:
# Coexistence of bulk antiferromagnetic order and superconductivity in the QED3 theory of copper oxides

## Abstract

Within the framework of the QED3 theory of cuprate superconductivity it is argued that bulk antiferromagnetic (AF) order can coexist with $`d`$-wave superconductivity in the underdoped region, in agreement with recent experiments. The AF order arises from the phase fluctuating $`d`$-wave superconductor via the mechanism of spontaneous chiral symmetry breaking, provided that fluctuations are sufficiently strong. The phase diagram for this coexistence is mapped out by means of analytical and numerical solutions of the underlying Dyson-Schwinger equation in the large $`N`$ limit.

Interplay between the AF order and superconductivity remains one of the central themes in the physics of the high-$`T_c`$ cuprates anderson1 ; zhang1 ; lee1 ; sachdev1 . Recently, an appealing new connection between AF and $`d`$-wave superconducting ($`d`$SC) orders has been pioneered, based on the ideas originally articulated by Emery and Kivelson emery1 . When long range $`d`$SC order is destroyed by thermal or quantum vortex-antivortex fluctuations corson1 ; ong1 ; fischer1 ; balents1 ; ft1 , the resulting state can be either a symmetric algebraic Fermi liquid ftv1  or, if the fluctuations are sufficiently strong, an incommensurate antiferromagnet herbut1 ; ft2 ; herbut2  with ordering vector $`𝐐`$ illustrated in Fig. 1. In the latter case AF order arises through an inherent dynamical instability of the underlying effective low energy theory of a phase fluctuating $`d`$-wave superconductor, a (2+1) dimensional quantum electrodynamics, QED3 ft1 ; ftv1 . This instability is known as the spontaneous chiral symmetry breaking (CSB) pisarski1 ; appelquist1 ; nash1  and is a well studied phenomenon in the particle physics.

In this communication we investigate the possibility that bulk AF order can set in within the superconducting phase, resulting in the region of coexistence of AF and $`d`$SC orders in the phase diagram shown in Fig. 1, within the QED3 framework. Other approaches to this problem have been considered previously; see e.g. Ref. chen1 . Recently, experiments have found tantalizing hints of such coexistence in zero applied magnetic field in Y and La based cuprates brewer1 ; nieder1 ; sonier1 ; mook1 ; sidis1 . It has been shown previously fst1  that within the QED3 framework such coexistence indeed can occur locally in the vicinity of fluctuating field-induced vortices, as found in neutron scattering lake1 ; lake2 , $`\mu `$SR miller1  and STM davis1  experiments.

Within the QED3 theory ft1  the dynamical agent responsible for the emergence of AF order is a noncompact U(1) Berry gauge field $`a_\mu `$ which encodes the topological frustration encountered by the nodal fermions as they propagate on the background of the fluctuating vortex-antivortex plasma. In the non-superconducting phase Berry gauge field is massless ft1 ; ftv1 ; herbut1 . Its quanta, “berryons”, have the same effect as photons in ordinary QED3: they mediate long range interactions between fermions and precipitate the chiral instability if the number of fermion species $`N`$ (pairs of Dirac nodes per unit cell in a $`d`$SC) is less than a critical value $`N_c3.2`$ pisarski1 ; appelquist1 . In the $`d`$SC phase, as vortices bind to finite pairs or loops, the Berry gauge field becomes massive ft1 ; ftv1 ; herbut1 . Nominaly, one would expect that coupling to such massive gauge field should become irrelevant for the low energy physics. Here we show that this is not necessarily so. Based on the approximate analytical and full numerical solutions of the underlying large-$`N`$ Dyson-Schwinger equation for the fermion self energy we find that solutions with broken chiral symmetry and finite fermion mass $`m_D`$ can be found even when the gauge field acquires small mass $`m_a`$. This opens up a possibility for the coexistence of bulk AF and $`d`$SC orders in the QED3 theory of underdoped cuprates.

As shown in Refs. ft1 ; ftv1  the low-energy dynamics of fermionic quasiparticles in a $`d`$-wave superconductor coupled to fluctuating vortices is governed by the QED3 Euclidean action $`S=d^3x_D`$ with

$$_D\underset{l=1}{\overset{N}{}}\overline{\mathrm{\Psi }}_l(x)\gamma _\mu (i_\mu a_\mu )\mathrm{\Psi }_l(x)+_B[a(x)].$$
(1)
Here, $`\mathrm{\Psi }_l(x)`$ is a four component Dirac spinor representing the “topological”fermion excitations associated with a pair of antipodal nodes, $`x=(\tau ,𝐫)`$ denotes the space-time coordinate, and $`\gamma _\mu `$ ($`\mu =0,1,2`$) are the gamma matrices satisfying $`\{\gamma _\mu ,\gamma _\nu \}=2\delta _{\mu \nu }`$. The number $`N`$ of

fermion species is equal to $`2n_{\mathrm{CuO}}`$, with $`n_{\mathrm{CuO}}`$ denoting the number of CuO planes per unit cell ftv1 . Lagrangian $`_B`$ encodes the dynamics of the gauge field $`a_\mu `$ and is given by $`_B[a]=\frac{1}{2}\mathrm{\Pi }_{\mu \nu }(q)a_\mu (q)a_\nu (q)`$ with

$$\mathrm{\Pi }_{\mu \nu }(q)=\frac{1}{e^2}\left(m_a^2+\alpha |q|+q^2\right)\left(\delta _{\mu \nu }\frac{q_\mu q_\nu }{q^2}\right).$$
(2)
The mass and Maxwell terms in Eq. (2) originate from the bare vortex action; $`m_a=0`$ in the non-superconducting phase, reflecting the unbound nature of vortex-antivortex fluctuations, while $`m_a>0`$ in the superconductor ft1 ; ftv1 ; herbut1 . The $`q`$-linear term, with $`\alpha =Ne^2/8`$, reflects the partial screening of the gauge field by the medium of topological fermions, expressed to one loop order.

We analyze the CSB phenomenon in the Lagrangian (1,2) by solving the Dyson-Schwinger (DS) equation for the fermion self energy $`\mathrm{\Sigma }(p)`$ to leading order in $`1/N`$ expansion. This technique is known to be robust against the higher order $`1/N`$ corrections even for small $`N2`$ nash1 ; gusynin1 , and the results appear to be insensitive to vertex corrections maris1 .

With the above definitions we can write the DS equation for the fermion self energy as appelquist1

$$\mathrm{\Sigma }(p)=\frac{8\alpha }{N}\frac{d^3k}{(2\pi )^3}\frac{\gamma _\mu D_{\mu \nu }(pk)\mathrm{\Sigma }(k)\gamma _\nu }{k^2+\mathrm{\Sigma }^2(k)},$$
(3)
where $`D_{\mu \nu }(q)=\mathrm{\Pi }_{\mu \nu }^1(q)`$ is the gauge field propagator in the Landau gauge. Following Pisarski pisarski1 , the simplest way to establish the spontaneous mass generation for fermions is to assume that, crudely, $`\mathrm{\Sigma }(p)`$ is constant over most of the integration range, $`\mathrm{\Sigma }(p)=m_D`$, and demand self consistency for the constant value $`m_D`$:

$$m_D=\frac{8\alpha }{N}\frac{d^3k}{(2\pi )^3}\frac{\gamma _\mu D_{\mu \nu }(k)m_D\gamma _\nu }{k^2+m_D^2}.$$
(4)
The above DS equation always has a trivial solution $`m_D=0`$, corresponding to the chirally symmetric state. A non-trivial, symmetry broken solution $`m_D>0`$ satisfies

$$1=\frac{8\alpha }{\pi ^2N}_0^{\mathrm{}}\frac{k^2dk}{(k^2+m_D^2)(m_a^2+\alpha k+k^2)},$$
(5)
where we have taken the trace of both sides and performed the angular integrals. For small gauge field mass, $`m_a^2\alpha m_D`$, we observe that the fermion mass effectively cuts off the integral in the infrared while $`\alpha `$ plays the role of the ultraviolet cutoff. We may thus approximate Eq. (5) by

$$1=\frac{8\alpha }{\pi ^2N}_{m_D}^\alpha \frac{k^2dk}{k^2(\alpha k+m_a^2)}=\frac{8}{\pi ^2N}\mathrm{ln}\left(\frac{\alpha ^2}{\alpha m_D+m_a^2}\right).$$

In this approximation the fermion mass is given by

$$m_D=\alpha e^{\pi ^2N/8}\frac{m_a^2}{\alpha }.$$
(6)
This is just the classic result of Pisarski pisarski1  with a small correction due to the gauge field being massive. Eq. (6) thus suggests that CSB in QED3 could survive in the presence of small gauge field mass. Translated into the language of condensed matter physics Eq. (6) suggests that even bound vortex-antivortex fluctuations could give rise to AF instability within the state with the true $`d`$SC long range order.

One can also perform the integral appearing in Eq. (5) exactly and confirm that the solution (6) indeed emerges in the limit $`m_a^2\alpha m_D`$ (up to an unimportant prefactor). The approximation employed to obtain Eq. (6) is instructive in that it reveals the fundamental reason for the insensitivity of CSB to small gauge field mass: the largest contribution to the RHS of Eq. (5) comes from the region $`m_D<k<\alpha `$. Therefore, modifying the gauge boson propagator at momenta $`km_D`$ by introducing the mass term has very little effect on the system.

The above solution is highly suggestive but because of the crude approximations involved it is deficient in at least two ways. First, for $`m_a=0`$, it fails to reproduce the result of a more refined treatment of DS Eq. (3) appelquist1 ; nash1  that no CSB occurs for $`N>N_c`$, with $`N_c=32/\pi ^2`$ to leading order in $`1/N`$ appelquist1 . Second, Eq. (6) is only valid for $`m_a^2\alpha m_D`$. In what follows we present a more careful treatment of Eq. (3) which is free of the above deficiencies. We use these solutions to map out the phase diagram of CSB as a function of $`N`$ and $`m_a`$.

To this end we must relax our assumption of constant $`\mathrm{\Sigma }(p)`$ in Eq. (3) and treat its full functional dependence on the three-momentum $`p`$. It is still possible to perform the angular integral to obtain

$`\mathrm{\Sigma }(p)`$ $`=`$ $`{\displaystyle \frac{4\alpha }{\pi ^2Np}}{\displaystyle \frac{1}{\sqrt{\alpha ^24m_a^2}}}{\displaystyle _0^{\mathrm{}}}𝑑k{\displaystyle \frac{k\mathrm{\Sigma }(k)}{k^2+\mathrm{\Sigma }^2(k)}}`$
$`\times `$ $`\left[\theta _1\mathrm{ln}\left({\displaystyle \frac{k+p+\theta _1}{|kp|+\theta _1}}\right)\theta _2\mathrm{ln}\left({\displaystyle \frac{k+p+\theta _2}{|kp|+\theta _2}}\right)\right]`$
where $`\theta _{1,2}=\frac{1}{2}(\alpha \pm \sqrt{\alpha ^24m_a^2})`$. In the limit $`\alpha ^2m_a^2`$ we recover the equation studied in Ref. appelquist1 ,

$$\mathrm{\Sigma }(p)=\frac{8}{\pi ^2Np}_{\frac{m_a^2}{\alpha }}^{\mathrm{}}𝑑k\frac{k\mathrm{\Sigma }(k)}{k^2+\mathrm{\Sigma }^2(k)}\mathrm{min}(k,p),$$
(8)
but with modified lower bound on the integral. In the massless theory ($`m_a=0`$), Eq. (8) is solved by linearizing the integrand in $`\mathrm{\Sigma }`$ and seeking the solution in the form of a power law, $`\mathrm{\Sigma }(p)p^b`$. Only solutions with complex $`b`$ turn out to be physically admissible placing a condition on the number of fermion species $`N<N_c=32/\pi ^2`$. The chiral mass is then given by appelquist1

$$m_D\mathrm{\Sigma }(0)=c\alpha e^{2\pi /\sqrt{N_c/N1}},$$
(9)
with $`c`$ a numerical constant.

Returning to the theory with massive gauge field, we observe that a simple power law ansatz no longer solves the linearized integral equation (8), essentially because of the appearance of the nonzero lower bound of the integral. Consequently, it is a non-trivial issue to extend such an analysis to our case of interest with $`m_a>0`$ gusynin2 . Fortunately, it is relatively straightforward to analyze the full nonlinear, angle integrated, DS Eq. (Coexistence of bulk antiferromagnetic order and superconductivity in the QED3 theory of copper oxides) numerically. We start from a constant $`\mathrm{\Sigma }(k)`$ and numerically evaluate the RHS at $`n`$ discrete points $`p_i`$. We then iterate this procedure until the solution for $`\mathrm{\Sigma }(p)`$ no longer changes appreciably between iterations. In this we adopt an upper cutoff, $`\mathrm{\Lambda }`$, for the integral. We find that when $`\mathrm{\Lambda }\alpha `$, our solutions no longer depend on $`\mathrm{\Lambda }`$.

Figures 2 and 3 summarize our findings.

Fig. 2 displays the dependence of chiral mass $`m_D`$ on the gauge field mass $`m_a`$ for selected values of $`N<N_c`$. A notable feature here is the surprising resilience of the chiral mass: $`m_D`$ persists until very large values of $`m_a`$. A hint of why this may be so is offered by a simple power counting in Eq. (5). This reveals that the relevant scale actually is not $`m_a`$ itself but $`m_a^2/\alpha `$, with $`m_a/\alpha 1`$. The latter quantity also enters our estimate for $`m_D`$, Eq. (6), which Fig. 2 confirms to be valid for small $`m_a`$.

When scaled by $`m_D^{}`$ and $`m_a^{}`$, defined as endpoints of the curves in Fig. 2, the points for different $`N`$ fall on the same universal curve (see inset). One way to understand this universality is from Eq. (5): its exact solution is given by the implicit equation

$$f(\frac{m_D}{m_D^{}},\frac{m_a^2}{m_a^2})=0,$$
(10)
with $`f(x,y)=x^2\mathrm{ln}x+y^2\mathrm{ln}y+\frac{\pi }{2}xy`$, independent of $`N`$. This solution is displayed as a dashed line in the inset of Fig. (2). We find that the scaling holds very well up to $`N2.9`$; closer to $`N_c^0N_c(m_a=0)=32/\pi ^2`$ the agreement gets progressively worse. This comparison indicates that the approximation of constant self energy, leading to Eq. (5), is very accurate for determining $`m_D`$, as long as $`N`$ is not too close to $`N_c^0`$.

Another way of viewing the above data is to consider $`N_c`$ as a function of $`m_a`$. Fig. 3 shows this dependence as determined by our numerical solution of Eq. (Coexistence of bulk antiferromagnetic order and superconductivity in the QED3 theory of copper oxides). As expected $`N_c`$ is a decreasing function of $`m_a`$. A striking feature is that $`N_c`$ remains relatively large for $`m_a`$ up to significant fraction of $`\alpha `$. The following simple argument leads to a highly accurate heuristic formula for $`N_c(m_a)`$ At low energies there are only two scales in the problem: $`m_D`$ and $`m_a`$. Increasing $`m_a`$ leads to reduction of $`m_D`$. It is then natural to expect that $`m_D`$ vanishes when the two scales cross, i.e. $`m_am_D^{}`$. Here $`m_D^{}`$ is the fermion mass for given $`N`$ but with massless gauge field, given by Eq. (9). Given our previous insight that the relevant quantity is $`m_a^2/\alpha `$ we take our criterion to read

$$(m_a/\alpha )^2=C(m_D^{}/\alpha ),$$
(11)
with $`C`$ is a numerical constant. Inserting $`m_D^{}`$ from Eq. (9) and solving for $`N`$ we obtain

$$N_c(m_a)=\frac{32}{\pi ^2}\left[1+\left(\pi /\mathrm{ln}\frac{𝒞\alpha }{m_a}\right)^2\right]^1,$$
(12)
where $`𝒞=\sqrt{cC}`$. As seen from the fit in Fig. 3 for $`𝒞2.7`$ this leads to an excellent agreement with our numerical results. We note that the same expression, with $`𝒞=e^27.4`$, was obtained by Gusynin et al. gusynin2  from the linearized version of the DS equation with an infrared cutoff.

Our results presented above firmly establish the existence of nontrivial, chiral symmetry broken solution of the QED3 Dyson-Schwinger equation in the presence of a small gauge field mass. The question arises whether this represents a true ground state of the system, as is known to be the case for the massless-photon QED3 appelquist1 . We believe this remains true when $`m_a>0`$. Imagine turning on the gauge field mass from zero deep inside the symmetry broken phase, where there is large energy difference between the ground state and the “false vacuum” represented by chirally symmetric phase. It is then clear that introduction of arbitrarily small gauge field mass cannot raise the energy of the ground state above that of the false vacuum: the chirally broken phase must remain a ground state of the system for some range of values $`m_a<m_a^{}`$, as indicated by our solutions (6) and Fig. 2.

Detailed nature of the criticality in the presence of gauge field mass appears to be a nontrivial issue. The transition in the ordinary massless-photon QED3 is thought to be of a special “conformal” variety, i.e. a continuous transition of infinite order appelquist2 ; miransky1 . This is due to the long range nature of interactions mediated by the massless gauge field. One would therefore expect that in the massive case the transition becomes a conventional second order gusynin2 ; gusynin3 .

Recently, powerful field theoretic duality arguments have been advanced cohen1  leading to a proposal for an upper bound $`N_c^0\frac{3}{2}`$. This suggestion finds some support in numerical simulations of noncompact lattice QED3 kogut1  which found no decisive signal for chiral mass generation at $`N=2`$. If these results are correct this would suggest that the conventional state of the art analysis based on the Dyson-Schwinger equation (3) overestimates $`N_c^0`$ by more than a factor of 2. Our analysis, based on this same technique, would then also be quantitatively inaccurate. However, we expect its qualitative features to remain valid.

If we stipulate that chiral symmetry breaking occurs in the massive-photon QED3, does this necessarily imply coexistence of the AF and $`d`$SC orders in cuprates? Answering this question involves additional subtlety that is related to the detailed nature of the criticality at the transition to the $`d`$SC state. According to the standard phenomenology ft1 ; ftv1 ; herbut1 ; ft2 , approaching the transition from the above (i.e. from the symmetric pseudogap phase) berryon is massless but the charge $`e`$ tends to zero, transforming eventually the Maxwell term into the mass term in the $`d`$SC phase. Since according to Eq. (9) $`m_D`$ is proportional to $`e^2`$ (through $`\alpha `$) this would suggest that chiral mass should vanish at the transition. However, there is presumably higher order $`(q^4)`$ term in the bare berryon action whose prefactor also diverges at the transition giving rise to a Maxwell term in the $`d`$SC phase. The chiral mass would then be proportional to this new “charge” in the $`d`$SC phase. Our calculations above strongly suggest that such a chirally broken phase can be a globally stable ground state of a system with massive berryons. Again, while this work might have resolved the nature of the phases on both sides of the transition, the detailed understanding of the criticality itself appears to be a complex problem awaiting future resolution.

QED3 theory as formulated in Refs. ft1 ; ftv1  therefore appears to support a locus of AF/$`d`$SC coexistence in the low doping region of the phase diagram Fig. 1. In this region the system remains superconducting while fermionic excitations become fully gapped. This small gap should be observable in thermodynamic and transport measurements. In particular the superfluid density should become exponentially activated in very clean samples at low temperatures.

The authors are indebted to I. F. Herbut, D. E. Sheehy and Z. Tešanović for many penetrating discussions. This work was supported by NSERC; in addition M.F. acknowledges the support of the A. P. Sloan Foundation.


text:
# Onset of Delocalization in Quasi-1D Waveguides with Correlated Surface Disorder

## Abstract

We present first analytical results on transport properties of many-mode waveguides with rough surfaces having long-range correlations. We show that propagation of waves through such waveguides reveals a quite unexpected phenomena of a complete transparency for a subset of propagating modes. These modes do not interact with each other and effectively can be described by the theory of 1D transport with correlated disorder. We also found that with a proper choice of model parameters one can arrange a perfect transparency of waveguides inside a given window of energy of incoming waves. The results may be important in view of experimental realizations of a selective transport in application to both waveguides and electron/optic nanodevices.

During last few years there is a burst of interest to solid state models with correlated disorder (see, e.g.,  and references therein). This fact is mainly due to the possibility to observe an anomalous transport in 1D models with random correlated potentials. In particular, it was shown  that specific long-range correlations in potentials give rise to a complete transparency of electron waves for given energy windows. Experimental realization of such potentials for single-mode waveguides with delta-like scatters  has confirmed theoretical predictions.

The subject of wave propagation through surface-disordered waveguides (see, e.g., and references therein) is important both from the theoretical viewpoint and for experimental applications such as optic fibers, remote sensing, radio wave propagation, shallow water waves, etc. On the other hand, it describes electronic transport in mesoscopic conducting channels as well. So far, main results in this field are obtained for random surfaces with fast-decaying correlations along a scattering surface. It is now of great importance to understand main properties of multiple scattering from rough surfaces with long-range correlations.

In this contribution we theoretically analyze transport properties of quasi-1D waveguides with rough surfaces. Our main interest is in exploring the possibility of constructing such surface profiles that result in frequency windows of a perfect propagation of waves.

In what follows we consider a plane waveguide (or conducting wire) of average width $`d`$ and length $`L`$, stretched along the $`x`$-axis. One (lower) surface of the waveguide is a rough profile $`z=\xi (x)`$, and the other (upper) surface, for simplicity, is taken to be flat, $`z=d`$. The function $`\xi (x)`$ describes the boundary roughness and is assumed to be statistically homogeneous and isotropic, with standard statistical properties,

$$\xi (x)=0,\xi (x)\xi (x^{})=\sigma ^2𝒲(|xx^{}|).$$
(1)
Here the angular brackets stand for a statistical average over different realizations of a profile $`\xi (x)`$. The binary correlator $`𝒲(|xx^{}|)`$ is assumed to decrease on a characteristic scale $`k_0^1`$, and is normalized to its maximal value, $`𝒲(0)=1`$.

Both the amplitude and gradient roughness of the corrugated surface is taken to be small, $`\sigma d`$ and $`k_0\sigma 1`$. These limitations are common in surface scattering theories that are based on an appropriate perturbative approach, see for example, . We should also add that zero (Dirichlet) boundary conditions are applied at both surfaces of the waveguide.

Let us first consider a single-mode waveguide in which only one normal mode (channel) has real lengthwise wave number $`k_1=\sqrt{k^2(\pi /d)^2}`$ (for details see ). It is known  that in this case the problem of surface scattering is reduced to a purely 1D model with the potential determined by the random profile $`\xi (x)`$. The corresponding wave equation takes the form,

$$\left(\frac{d^2}{dx^2}+k_1^2\right)\mathrm{\Psi }(x)\frac{2}{\pi }\left(\frac{\pi }{d}\right)^3\xi (x)\mathrm{\Psi }(x)=0.$$
(2)
Since one channel case corresponds to inequality $`1<kd/\pi <2`$, the value of $`k_1`$ is restricted by the relation $`0<k_1d/\pi <\sqrt{3}`$. Note also that the weak scattering assumption $`\sigma d`$ leads to the condition, $`k_1\sigma 1`$.

The above equation (2) demonstrates an exact correspondence between the surface scattering in one-mode waveguides, and bulk scattering for electrons in 1D models of the Anderson (or Kronig-Penney) type. Specifically, the wave number $`k`$ is equal to $`\omega /c`$ for a classical wave of frequency $`\omega `$, and to the Fermi wave number for electrons.

In order to solve the wave equation (2) and to find the averaged transmittance, or, the same, the dimensionless conductance (in units of $`e^2/\pi \mathrm{}`$), one can use methods developed in the transport theory of one-dimensional disordered systems, such as the perturbative diagrammatic method of Berezinski , invariant imbedding method , or the two-scale approach . All the methods allow to take adequately into account the effects of coherent multiple scattering of waves (electrons) from the corrugated surface.

Main theoretical result is that the average transmittance $`T`$ is described by a universal expression which depends on one parameter $`\mathrm{\Lambda }=L_{loc}/L`$ only (one-parameter scaling). Here the quantity $`L_{loc}/4`$ is the backscattering length in a 1D disordered system. Its inverse value $`L_{loc}^1`$ is equal to the Lyapunov exponent that can be found in the transfer matrix approach, so that the quantity $`L_{loc}`$ is known as the localization length .

The expression for $`T`$ is fully consistent with the theory of 1D Anderson localization. Specifically, the transmittance $`T`$ exhibits the ballistic behavior for large localization length, $`L_{loc}/L1`$,

$$T(L/L_{loc})14L/L_{loc},$$
(3)
therefore, in this case the waveguide is practically transparent. On the contrary, the transmittance is exponentially small when the localization length is much less than the length of the waveguide, $`L_{loc}/L1`$,

$$T(L/L_{loc})\frac{\pi ^{5/2}}{16}\left(L/L_{loc}\right)^{3/2}\mathrm{exp}\left(L/L_{loc}\right).$$
(4)
This implies strong wave localization. Since in this case the transmittance is not a self-averaged quantity, for completeness of presentation one should refer to the expression for its logarithm, $`\mathrm{ln}T(L/L_{loc})=4L/L_{loc}`$.

Thus, in order to describe transport properties of a disordered waveguide, one needs to know the value of the localization length $`L_{loc}`$. The latter is entirely determined by the type of disorder. For the surface disorder, from the solution of the equation (2) one can obtain ,

$$\frac{1}{L_{loc}}=\frac{\sigma ^2}{\pi ^2}\left(\frac{\pi }{d}\right)^6\frac{W(2k_1)}{(2k_1)^2}.$$
(5)
This formula contains the Fourier transform $`W(k_x)`$ of the binary correlator $`𝒲(|xx^{}|)`$, which is called the roughness power spectrum,

$$𝒲(|xx^{}|)=_{\mathrm{}}^{\mathrm{}}\frac{dk_x}{2\pi }\mathrm{exp}\{ik_x(xx^{})\}W(k_x).$$
(6)
One can see that the localization length and, correspondingly, transport properties of a single-mode waveguide are completely defined by the dependence $`W(k_x)`$. In particular, if $`W(2k_1)`$ abruptly vanishes within some interval of wave number $`k_1`$, then the localization length $`L_{loc}`$ diverges and the waveguide of any length will be fully transparent. Therefore, an interesting problem arises: how to construct such (random) surface profiles that result in a complete transparency of waveguides within a predefined part of the allowed region for the wave vector $`k_1`$? The answer to this question was done in Ref. where the results obtained for discrete Anderson-like disordered models  have been extended to the surface scattering.

The recipe is as follows . Having a desirable dependence for the Fourier transform $`W(k_x)`$, one obtains the function $`\beta (x)`$,

$$\beta (x)=_{\mathrm{}}^{\mathrm{}}\frac{dk_x}{2\pi }\mathrm{exp}\left(ik_xx\right)W^{1/2}(k_x).$$
(7)
Then, the surface profile can be constructed by the convolution of the delta-correlated random process (white noise $`Z(x)`$) with the function $`\beta (x)`$,

$$\xi (x)=\sigma _{\mathrm{}}^{\mathrm{}}𝑑x^{}Z(xx^{})\beta (x^{}).$$
(8)
The proposed method allows us to practically solve the inverse scattering problem of a reconstruction of rough surfaces from their power spectrum. We would like to emphasize that by a proper choice of the binary correlator, the transition between localized and delocalized wave/electron transport can occur abruptly at a given point inside the allowed interval for the wave number $`k_1`$ (see details and examples in ). This means that the roughness power spectrum $`W(k_x)`$ is discontinuous at this point. Consequently, the correlator $`𝒲(|xx^{}|)`$ of the corrugated surface $`\xi (x)`$ should be of specific form with long-range correlations along the waveguide.

It is important to stress that surfaces with the above discussed properties (when $`W(k_x)`$ vanishes in some regions of $`k_x`$) are not exotic. They have been recently fabricated in experimental studies of an enhanced backscattering .

Let us now consider many-mode waveguides with correlated surface disorder. In the spirit of the Landauer’s concept , the total transmittance $`T`$ of any quasi-1D structure is expressed as a sum of partial transmittances $`T_n`$ that describe the transport for every $`n`$-th propagating normal mode, $`T=_{n=1}^{N_d}T_n`$. Here $`N_d=[kd/\pi ]`$ is the total number of propagating modes (open channels) determined by the integer part $`[\mathrm{}]`$ of the ratio $`kd/\pi `$.

From general theory of scattering systems  it follows that the wave propagation through any $`n`$-th open channel is determined by two attenuation lengths, the length $`L_n^{(f)}`$ of forward scattering and the backscattering length $`L_n^{(b)}`$. For quasi-1D waveguides with surface disorder the inverse scattering lengths are given by,

$$\frac{1}{L_n^{(f)}}=\sigma ^2\frac{(\pi n/d)^2}{k_nd}\underset{n^{}=1}{\overset{N_d}{}}\frac{(\pi n^{}/d)^2}{k_n^{}d}W(k_nk_n^{}),$$
(9)
$$\frac{1}{L_n^{(b)}}=\sigma ^2\frac{(\pi n/d)^2}{k_nd}\underset{n^{}=1}{\overset{N_d}{}}\frac{(\pi n^{}/d)^2}{k_n^{}d}W(k_n+k_n^{}),$$
(10)
where $`k_n=\sqrt{k^2(\pi n/d)^2}`$.

The results (9) and (10) can be obtained by the diagrammatic Green’s function approach  as well as with the use of the technique developed in Ref. . Also, these expressions directly follow from the invariant imbedding method extended to quasi-1D structures (see Ref. ). Note that in a single-mode waveguide with $`N_d=1`$ the sum over $`n^{}`$ contains only one term with $`n^{}=n=1`$. Therefore, the backscattering length $`L_1^{(b)}`$ is four times less then the 1D localization length $`L_{loc}`$, see Eq. (5).

The expressions (9) and (10) manifest that, in general, both attenuation lengths are formed by scattering of a given $`n`$-th propagating mode into all other modes. This is the case when, for example, a surface profile is either delta-correlated random function (“white noise”) with constant power spectrum $`W(k_x)`$, or the function with fast decreasing binary correlator (or, the same, with a slow decrease of its Fourier components).

One can see that the larger the mode number $`n`$ is, the stronger is the scattering of this mode into others. As was shown in Ref. , in this case a very interesting phenomenon of the coexistence of ballistic, diffusive, and localized transport arises, which seems to be generic for propagation through waveguides with random surfaces. Specifically, while lowest modes can be in the ballistic regime, the intermediate and highest modes exhibit the diffusive and localized behavior, respectively.

Now we are ready to show that the situation is fundamentally different when the surface disorder has specific long-range correlations. To demonstrate this, let us take the surface profile $`\xi (x)`$ with the binary correlator in the following form, $`𝒲(|xx^{}|)=\mathrm{sin}\{k_0(xx^{})\}/k_0(xx^{})`$. In this case the Fourier transform $`W(k_x)`$ is a ”window function”,

$$W(k_x)=(\pi /k_0)\mathrm{\Theta }(k_0|k_x|),k_0>0,$$
(11)
where $`\mathrm{\Theta }(z)`$ stands for the unit-step function.

In accordance with Eqs. (7), (8) rough surfaces that have such correlation properties can be constructed as follows,

$$\xi (x)=\frac{\sigma }{\sqrt{\pi k_0}}_{\mathrm{}}^{\mathrm{}}𝑑x^{}Z(xx^{})\frac{\mathrm{sin}(k_0x^{})}{x^{}}.$$
(12)
Now one can see that in the case under consideration the number of modes into which a given $`n`$-th mode is scattered, is entirely determined by the width $`k_0`$ of the rectangular spectrum (11). It is clear that if the distance $`|k_nk_{n\pm 1}|`$ between neighboring values of $`k_n`$ is larger than the width $`k_0`$,

$$|k_nk_{n\pm 1}|>k_0,$$
(13)
then the transitions between all modes are forbidden. In this case the sum over $`n^{}`$ in Eq. (9) contains only one term with $`n^{}=n`$ which describes scattering inside the channels. Moreover, each term in the sum of Eq. (10) is equal to zero, and the scattering lengths are,

$$\frac{1}{L_n^{(f)}}=\frac{\pi \sigma ^2}{k_0}\frac{(\pi n/d)^4}{(k_nd)^2},\frac{1}{L_n^{(b)}}=0.$$
(14)
The above analysis shows that a remarkable phenomenon can arise. Specifically, all the propagating modes with index $`n`$, that satisfy the condition (13), are independent of other waveguide modes, in spite of the interaction with a rough surface. In other words, they represent a coset of non-interacting channels. Since the transport of any 1D disordered system is determined by the backscattering length only, and the latter is infinite according to Eq.(14), the transmittance for such channels is $`T_n=1`$. This means that according to the Landauer’s formula for $`T`$, the transmittance of the coset of such independent transparent waveguide modes is simply equal to their total number.

As for other propagating modes with index $`n`$, that do not satisfy the condition (13), they remain to be mixed by surface scattering because the roughness power spectrum (11) is non-zero, $`W(k_nk_n^{})=\pi /k_0`$. Since these mixed modes have finite backscattering lengths, for large enough size $`L`$ they do not contribute in the total transmittance and the latter is equal to the number of independent transparent modes.

For large number $`N_d1`$ of open channels one can perform further analytical analysis. In this case the inequality (13) is analogous to the requirement $`|k_n/n|>k_0`$ which can be written in the following form,

$$n>N_{k_0}\left[\frac{(kd/\pi )}{\sqrt{1+(k_0d/\pi )^2}}\right].$$
(15)
We remind again that square brackets stand for the integer part of the inner expression.

The expression (15) determines the total number $`N_{k_0}`$ of mixed non-transparent modes, total number $`N_t=N_dN_{k_0}`$ of independent transparent modes, and the critical value of $`n`$ that divides these two groups. All propagating modes with $`n>N_{k_0}`$ are independent and fully transparent, otherwise, they are mixed and characterized by finite backscattering lengths $`L_n^{(b)}`$.

The numbers $`N_{k_0}`$ and $`N_t`$ of mixed non-transparent and independent transparent modes are determined by two parameters, $`kd/\pi `$ and $`k_0d/\pi `$. In the case of “weak” correlations, $`k_0d/\pi 1`$, the number of mixed modes $`N_{k_0}`$ is of the order of $`N_d`$,

$$N_{k_0}\left[\left(\frac{kd}{\pi }\right)\frac{1}{2}\left(\frac{kd}{\pi }\right)\left(\frac{k_0d}{\pi }\right)^2\right].$$
(16)
Consequently, in this case the number of independent modes $`N_dN_{k_0}`$ is small, or there are no such modes at all. If the parameter $`k_0d/\pi `$ tends to infinity, $`k_0d/\pi \mathrm{}`$, the rough surface profile becomes white-noise-like and, naturally, $`N_{k_0}N_d`$.

The most interesting case is when surface profiles are strongly correlated, which is specified by the condition $`k_0d/\pi 1`$. Then the characteristic integer $`N_{k_0}`$ is much less than the total number of modes $`N_d`$,

$$N_{k_0}\left[\left(\frac{kd}{\pi }\right)\left(\frac{k_0d}{\pi }\right)\right]N_d,$$
(17)
and the number of transparent modes $`N_t`$ is large. When $`k_0d/\pi `$ decreases and becomes small, $`k_0d/\pi <(kd/\pi )^11`$, the number $`N_{k_0}`$ vanishes and all modes become independent and fully transparent. Evidently, if $`k_0d/\pi `$ tends to zero, the roughness power spectrum $`W(k_x)`$, see (11), becomes delta-function-like and, as a consequence, $`N_{k_0}=0`$. In this case the correlated disorder results in a perfect transmission of waves, in spite of their scattering from a rough surface.

We also would like to note that the scattering through many-mode waveguides with correlated surfaces of the above kind reveals a step-wise dependence of the transmittance on the value of the parameter $`\gamma _k=kd/\pi `$. This is clear from the explicit expression for the transmittance,

$$T=[\gamma _k][\alpha \gamma _k],\alpha =\{1+(k_0d/\pi )^2\}^{1/2}.$$
(18)
The effect is similar to that known to occur for the conductance of quasi-1D ballistic structures (see, e.g., ). Interesting enough, in our model the step-wise dependence of the transmittance arises both for integer and non-integer values of $`\gamma _k`$, due to the presence of the second term in Eq.(18).

In conclusion, we have studied transport properties of quasi-1D waveguides with random correlated surfaces. It was analytically shown that for a specific choice of long-range correlations a subset of propagating modes emerges that are independent and fully transparent. These modes do not interact with each other and the theory of 1D transport can be rigorously applied for them. The rest of propagating modes creates the second set of mixed and non-transparent modes. With a proper choice of model parameters, one can construct such waveguides for which all modes turn out to be completely delocalized. This means a perfect transparency of waveguides in spite of random character of surface scattering. These results may be used for experimental fabrication of waveguides and supperlattices with selective windows of a complete transparancy in dependence of the wave number of incoming waves.

This research was supported by Consejo Nacional de Ciencia y Tecnología (CONACYT, México) grant 34668-E, and by the Universidad Autónoma de Puebla (BUAP, México) under the grant II-63G01.


text:
# Neutron-fed afterglows of gamma-ray bursts

## 1 Introduction

The presence of neutrons in gamma-ray bursts (GRBs) and their possible role for the observed emission was realized in the past few years (Derishev, Kocharovsky, & Kocharovsky 1999a,b; Bahcall & Mészáros 2000; Mészáros & Rees 2000; Fuller, Pruet, & Abazajian 2000; Pruet & Dalal 2002). In an accompanying paper, we study in detail the nuclear composition of GRB fireballs and show that the presence of neutrons among the ejected baryons is practically inevitable (Beloborodov 2003, Paper 1). One implication is an observable multi-GeV neutrino emission from inelastic neutron-ion collisions in the fireball. Here, we focus on a different aspect. We show that the neutrons have a dramatic impact on the explosion dynamics at radii as large as $`10^{17}`$ cm and propose a novel mechanism for GRB afterglow emission.

Let us remind what happens in a relativistic explosion without neutrons (see Mészáros 2002 for a review). The ejected fireball with mass $`M_{\mathrm{ej}}`$ and Lorentz factor $`\mathrm{\Gamma }_{\mathrm{ej}}`$ sweeps up an ambient medium and gradually dissipates its kinetic energy. The dissipation rate peaks at a characteristic “deceleration” radius $`R_{\mathrm{dec}}`$ where half of the initial energy is dissipated. This radius corresponds to swept-up mass $`m_{\mathrm{dec}}=M_{\mathrm{ej}}/\mathrm{\Gamma }_{\mathrm{ej}}`$. Further dynamics is described by the self-similar blast wave model of Blandford & McKee (1976). How does this picture change in the presence of neutrons?

## 2 Neutron front and its trail

Neutrons develop a Lorentz factor $`\mathrm{\Gamma }_n=10^210^3`$ at the very beginning of the explosion when the fireball is accelerated by radiation pressure (Derishev et al. 1999b). They are collisionally coupled to the ions in the early dense fireball, and decouple close to the end of the acceleration stage. Then the neutrons coast and gradually decay with a mean lifetime $`\tau _\beta 900`$ s and a mean decay radius $`R_\beta =c\tau _\beta \mathrm{\Gamma }_n`$,

$$R_\beta =0.8\times 10^{16}\left(\frac{\mathrm{\Gamma }_n}{300}\right)\mathrm{cm}.$$
(1)
At radii under consideration, $`R>10^{15}`$ cm, the ejected fireball is a shell of thickness $`\mathrm{\Delta }R`$. In contrast to neutrons, the ion component of the fireball is aware of the external medium and its Lorentz factor $`\mathrm{\Gamma }`$ decreases. As $`\mathrm{\Gamma }`$ decreases below $`\mathrm{\Gamma }_n`$, the ions fall behind and separate from the neutrons. Thus the fireball splits into two relativistic shells which we name N (neutrons) and I (ions); the N-shell will also be called “neutron front”. The I-shell lags behind by a distance $`l`$,

$$\frac{l}{R}\beta _n\beta _i\frac{1}{2\mathrm{\Gamma }^2}\frac{1}{2\mathrm{\Gamma }_n^2},$$
(2)
where $`\beta _i`$ and $`\beta _n`$ are the shell velocities in units of $`c`$. For simplicity, we hereafter focus on the stage of complete separation $`l>\mathrm{\Delta }`$ (it sets in right after the beginning of the I-shell deceleration if $`\mathrm{\Delta }10^{11}`$ cm, and it covers the whole explosion if $`\mathrm{\Delta }0`$).

The mass of the N-shell gradually decreases because of the $`\beta `$-decay,

$$M_n(R)=M_n^0\mathrm{exp}\left(\frac{R}{R_\beta }\right).$$
(3)
However, even at $`R>R_\beta `$, the N-shell energy $`E_n=\mathrm{\Gamma }_nM_nc^2`$ is huge compared to the ambient rest mass $`mc^2`$. For example, at $`R=R_{\mathrm{dec}}`$, $`E_n/m_{\mathrm{dec}}c^2=X_n\mathrm{\Gamma }_n\mathrm{\Gamma }_{\mathrm{ej}}\mathrm{exp}(R_{\mathrm{dec}}/R_\beta )`$ where $`X_n=M_n^0/M_{\mathrm{ej}}`$.

The neutron decay products $`p`$ and $`e^{}`$ share immediately their huge momentum with ambient particles due to two-stream instability (the instability timescale is set by the ion plasma frequency $`\omega _p`$ and it is the shortest timescale in the problem). Thus, the N-shell leaves behind a mixed trail with a relativistic bulk velocity $`\beta <\beta _n`$ which we calculate now.

Let $`\mathrm{d}m=(\mathrm{d}m/\mathrm{d}R)\mathrm{d}R`$ be ambient mass overtaken by the N-shell as it passes $`\mathrm{d}R`$ and $`\mathrm{d}M_n=(M_n/R_\beta )\mathrm{d}R`$ be mass of decayed neutrons. The $`\mathrm{d}m`$ and $`\mathrm{d}M_n`$ share momentum and form a trail element with proper mass $`\mathrm{d}m_{}=\mathrm{d}M_n+\mathrm{d}m+\mathrm{d}m_{\mathrm{heat}}`$ which includes heat dissipated in the inelastic $`\mathrm{d}m`$-$`\mathrm{d}M_n`$ collision. The laws of energy and momentum conservation read

$`\mathrm{\Gamma }_n\mathrm{d}M_n+\mathrm{d}m=\gamma \mathrm{d}m_{},\beta _n\mathrm{\Gamma }_n\mathrm{d}M_n=\beta \gamma \mathrm{d}m_{},`$ (4)
where $`\gamma =(1\beta ^2)^{1/2}`$ is the trail Lorentz factor. Let us denote

$$\zeta (R)=\frac{\mathrm{d}M_n}{\mathrm{d}m}=\frac{M_n}{R_\beta }\left(\frac{\mathrm{d}m}{\mathrm{d}R}\right)^1.$$
(5)
From equations (4) we find

$$\beta (R)=\frac{\beta _n}{1+(\mathrm{\Gamma }_n\zeta )^1},\gamma (R)=\frac{\mathrm{\Gamma }_n\zeta +1}{(\zeta ^2+2\mathrm{\Gamma }_n\zeta +1)^{1/2}}.$$
(6)
It gives $`\gamma 1`$ until essentially all neutrons have decayed. For illustration, let us specialize to a power-law density profile of the ambient medium. Then

$$m(R)=m_\beta \left(\frac{R}{R_\beta }\right)^k,\zeta (R)=\frac{M_n}{km_\beta }\left(\frac{R}{R_\beta }\right)^{1k},$$
(7)
and $`\zeta `$ evolves from $`\zeta =M_n^0e^1/km_\beta 1`$ at $`R=R_\beta `$ to $`\zeta 1`$ at $`RR_\beta `$ as $`M_n`$ decays exponentially. There exists a characteristic radius $`R_{\mathrm{trail}}`$ where the trail becomes nonrelativistic ($`\beta =0.5`$). It is defined by condition $`\zeta =\mathrm{\Gamma }_n^1`$ (see eq. 6), which requires about 10 e-folds of the decay (for a typical $`m_\beta m_{\mathrm{dec}}10^5M_n^0`$). Thus,

$$R_{\mathrm{trail}}10R_\beta =0.8\times 10^{17}\left(\frac{\mathrm{\Gamma }_n}{300}\right)\mathrm{cm}.$$
(8)
$`R_{\mathrm{trail}}`$ depends very weakly (logarithmically) on the ambient density and $`X_n`$.

Now that we know $`\beta `$, the trail density $`n`$ can be easily calculated. The neutron front has a small thickness $`\mathrm{\Delta }`$ and the ambient particles cross it almost instantaneously (on timescale $`\gamma ^2\mathrm{\Delta }/cR/c`$). Measured in the N-shell frame, the flux of ambient particles is $`\beta _n\mathrm{\Gamma }_nn_0=(\beta _n\beta )\mathrm{\Gamma }_n\gamma n_{\mathrm{amb}}`$ where $`n_0`$ and $`n_{\mathrm{amb}}`$ are proper densities of the ambient particles ahead and behind the N-shell, respectively. This gives the medium compression as it is accelerated in the neutron front: $`n_{\mathrm{amb}}/n_0=[\beta _n/\gamma (\beta _n\beta )]`$. The total density of the trail includes the neutron decay products, so $`n=(1+\zeta )n_{\mathrm{amb}}`$. Using equation (6), we get

$$\frac{n}{n_0}=\frac{\beta _n(1+\zeta )}{\gamma (\beta _n\beta )}=(1+\zeta )\left(\zeta ^2+2\mathrm{\Gamma }_n\zeta +1\right)^{1/2}.$$
(9)
The energy dissipated in the neutron front is given by

$$\mathrm{d}E_{\mathrm{n}.\mathrm{f}.}=\gamma \left(\mathrm{d}m_{}\mathrm{d}m\mathrm{d}M_n\right)c^2=\left(\mathrm{\Gamma }_n\gamma \right)\mathrm{d}M_nc^2(\gamma 1)\mathrm{d}mc^2.$$
(10)
After simple algebra (using eqs. ), we find

$$\frac{\mathrm{d}E_{\mathrm{n}.\mathrm{f}.}}{\mathrm{d}R}=\left[\frac{\mathrm{\Gamma }_n\beta _n\gamma \beta }{\mathrm{\Gamma }_n\gamma (\beta _n\beta )}1\right]\gamma \frac{\mathrm{d}m}{\mathrm{d}R}c^2.$$
(11)
It simplifies to $`\mathrm{d}E_{\mathrm{n}.\mathrm{f}.}=2\gamma ^2\mathrm{d}mc^2`$ for the interesting regime $`1\gamma \mathrm{\Gamma }_n`$ ($`\mathrm{\Gamma }_n^1<\zeta <\mathrm{\Gamma }_n`$).

The trail is formed very hot. This can be seen when comparing inertial mass $`\mathrm{d}m_{}`$ (that includes heat) with rest mass $`\mathrm{d}M_b=\mathrm{d}m+\mathrm{d}M_n=(1+\zeta )\mathrm{d}m`$, which gives dimensionless relativistic enthalpy

$$\mu =\frac{\mathrm{d}m_{}}{\mathrm{d}M_b}=\frac{(\zeta ^2+2\mathrm{\Gamma }_n\zeta +1)^{1/2}}{1+\zeta }.$$
(12)
For $`\mathrm{\Gamma }_n^1<\zeta <\mathrm{\Gamma }_n`$ we find $`\mu 1`$, i.e. internal energy of the trail far exceeds its rest-mass energy. The trail parameters are summarized in Table 1.

It is instructive to view the dissipation process in the rest frame of the trail. Here, the elements $`\mathrm{d}m`$ and $`\mathrm{d}M_n`$ have initial Lorentz factors $`\gamma `$ and $`\stackrel{~}{\gamma }=\mathrm{\Gamma }_n\gamma (1\beta \beta _n)`$, share their opposite momenta, and come at rest. This is achieved via a plasma instability that isotropizes the particle momentum distribution. It may result in two isotropic ion populations in the trail: ambient ions with mean Lorentz factor $`\gamma `$ and $`\beta `$-decay protons with mean Lorentz factor $`\stackrel{~}{\gamma }`$. Both populations are relativistically hot, $`\gamma 1`$ and $`\stackrel{~}{\gamma }1`$, as long as $`\mathrm{\Gamma }_n^1<\zeta <\mathrm{\Gamma }_n`$. They have equal energies $`\stackrel{~}{\gamma }\mathrm{d}M_n=\gamma \mathrm{d}m`$ (which correspond to equal relativistic bulk momenta of $`\mathrm{d}M_n`$ and $`\mathrm{d}m`$ before they come at rest in the trail frame). The ratio of their temperatures is the reciprocal of their density ratio.

A detailed model of dissipation in the neutron front is an interesting plasma physics problem, which we defer to a future work. We emphasize here that it is different from dissipation in collisionless shocks. The thickness of a shock, $`\delta c/\omega _p`$, is set by the timescale of two-stream instability, and it is a discontinuity in the hydrodynamical sense. By contrast, the N-shell has thickness $`\mathrm{\Delta }10^{11}10^{12}`$ cm which is orders of magnitude larger than $`\delta `$. The neutrons decay and cause dissipation everywhere in the N-shell. The ion medium velocity grows smoothly from 0 at the leading edge of the neutron front to $`\beta `$ at its trailing edge. This can be called “volume dissipation”.

## 3 Shock wave

The ion fireball follows the neutron front and collects the trail. As a result, (1) the ion Lorentz factor $`\mathrm{\Gamma }`$ decreases and (2) a shock wave propagates in the trail material. The shock has a Lorentz factor $`\mathrm{\Gamma }<\mathrm{\Gamma }_{\mathrm{sh}}<\mathrm{\Gamma }_n`$ and cannot catch up with the neutron front (unless $`n_0[R]`$ falls off steeper than $`R^3`$). The shock lags behind the neutron front by distance $`l_{\mathrm{sh}}R(\beta _n\beta _{\mathrm{sh}})`$ and its time lag is

$$t_{\mathrm{sh}}=\frac{l_{\mathrm{sh}}}{c(\beta _{\mathrm{sh}}\beta )}\frac{\gamma ^2}{(\mathrm{\Gamma }_{\mathrm{sh}}^2\gamma ^2)}\frac{R}{c}.$$
(13)
It is much shorter than $`R/c`$ as long as $`\mathrm{\Gamma }_{\mathrm{sh}}\gamma `$. Hereafter we use the approximation $`t_{\mathrm{sh}}<R/c`$ and assume that the trail is quickly picked up by the shock, before it could expand and change its density or velocity. We allow, however, for rapid radiative losses of the trail: it may cool on a timescale $`t_{\mathrm{cool}}R/c`$ and possibly $`t_{\mathrm{cool}}<t_{\mathrm{sh}}`$. Below we consider two extreme regimes: $`t_{\mathrm{cool}}t_{\mathrm{sh}}`$ (adiabatic) and $`t_{\mathrm{cool}}t_{\mathrm{sh}}`$ (radiative with $`\mu =1`$).

In the radiative regime, the energy dissipated by the shock is

$$\frac{\mathrm{d}E_{\mathrm{sh}}}{\mathrm{d}R}=\mathrm{\Gamma }(\mathrm{\Gamma }_{\mathrm{rel}}1)(1+\zeta )\frac{\mathrm{d}m}{\mathrm{d}R}c^2,$$
(14)
where $`\mathrm{\Gamma }_{\mathrm{rel}}=\mathrm{\Gamma }\gamma (1\beta _i\beta )`$ is the trail Lorentz factor with respect to the I-shell. For $`\mathrm{\Gamma }_n^1\zeta \mathrm{\Gamma }_n`$ one can use approximate expressions $`\zeta =2\gamma ^2/\mathrm{\Gamma }_n`$ and $`\mathrm{\Gamma }_{\mathrm{rel}}=(1/2)(\gamma /\mathrm{\Gamma }+\mathrm{\Gamma }/\gamma )`$ to get $`\mathrm{d}E_{\mathrm{sh}}=(1/2\gamma +\gamma /\mathrm{\Gamma }_n)(\mathrm{\Gamma }\gamma )^2\mathrm{d}mc^2`$. Note that the dissipation is smaller than it would be in the absence of the neutron front, $`\mathrm{d}E_{\mathrm{sh}}=\mathrm{\Gamma }(\mathrm{\Gamma }1)\mathrm{d}mc^2`$. Hence, in the radiative regime, the neutrons delay the deceleration of the ion fireball.

In the adiabatic regime, the preshock material has a high enthalpy $`\mu `$ given by equation (12). The postshock heat, measured in the lab frame, is $`\mathrm{d}E_{\mathrm{heat}}=\mathrm{\Gamma }(\mathrm{\Gamma }_{\mathrm{rel}}\mu 1)(1+\zeta )\mathrm{d}m`$. It includes heat deposited into the preshock medium by the neutron front (eq. 10), which should be subtracted to get the energy dissipated in the shock itself. Hence,

$$\frac{\mathrm{d}E_{\mathrm{sh}}}{\mathrm{d}R}=\left[\mathrm{\Gamma }(\mathrm{\Gamma }_{\mathrm{rel}}\mu 1)\gamma (\mu 1)\right](1+\zeta )\frac{\mathrm{d}m}{\mathrm{d}R}c^2.$$
(15)
For $`\mathrm{\Gamma }_n^1\zeta \mathrm{\Gamma }_n`$ we get $`\mathrm{d}E_{\mathrm{n}.\mathrm{f}.}=2\gamma ^2\mathrm{d}mc^2`$, $`\mathrm{d}E_{\mathrm{sh}}=(\mathrm{\Gamma }^2\gamma ^2)\mathrm{d}mc^2`$, and the total dissipated energy $`\mathrm{d}E_{\mathrm{heat}}=(\mathrm{\Gamma }^2+\gamma ^2)\mathrm{d}mc^2`$. As long as $`\mathrm{\Gamma }\gamma `$, the bulk of energy is dissipated in the shock rather than in the neutron front. Moreover, the dissipated energy is the same as in the absence of a neutron front, and hence the fireball deceleration is the same.

In both radiative and adiabatic regimes, the shock dissipation is suppressed when the blast wave decelerates to $`\mathrm{\Gamma }\gamma `$. At $`\mathrm{\Gamma }=\gamma `$ the shock would disappear. Thus, $`\mathrm{\Gamma }`$ is bound from below by $`\gamma `$.

## 4 Numerical examples

Let us consider the simple radiative regime. The I-shell and the shocked part of the trail can be treated as a single shell with a growing mass $`M(R)=M_{\mathrm{ej}}+mM_n`$ (a known function of radius) and a decreasing $`\mathrm{\Gamma }`$. Mass gain $`\mathrm{d}M`$ causes deceleration $`\mathrm{d}\mathrm{\Gamma }`$ that is found from the energy and momentum conservation: $`\mathrm{d}(\mathrm{\Gamma }M)=\gamma \mathrm{d}M\mathrm{d}E_{\mathrm{rad}}/c^2`$ and $`\mathrm{d}(\beta _i\mathrm{\Gamma }M)=\beta \gamma \mathrm{d}M(\mathrm{d}E_{\mathrm{rad}}/c^2)\beta _i`$ where $`\mathrm{d}E_{\mathrm{rad}}=\mathrm{d}E_{\mathrm{sh}}`$ is the radiated energy. Excluding $`\mathrm{d}E_{\mathrm{rad}}`$, we get the dynamic equation,

$$M\frac{\mathrm{d}\mathrm{\Gamma }}{\mathrm{d}R}=\mathrm{\Gamma }^2\gamma \beta _i(\beta _i\beta )(1+\zeta )\frac{\mathrm{d}m}{\mathrm{d}R}.$$
(16)
Using equations (6) and (7) we can solve numerically equation (16) for $`\mathrm{\Gamma }(R)`$ with an initial condition $`\mathrm{\Gamma }_{\mathrm{ej}}`$. The results depend on the type of the ambient medium. Figure 1 shows two examples: a standard ISM and a wind from a Wolf-Rayet progenitor.

The radiative explosion has two separate emission regions: behind the neutron front and behind the shock front. Their luminosities equal the corresponding dissipation rates $`(c\mathrm{d}E_{\mathrm{n}.\mathrm{f}.}/\mathrm{d}R)`$ and $`(c\mathrm{d}E_{\mathrm{sh}}/\mathrm{d}R)`$ (Fig. 1). The neutron front dissipation peaks at $`R_\beta `$. The shock dissipation can have two peaks (if $`m_\beta M_{\mathrm{ej}}/\mathrm{\Gamma }_{\mathrm{ej}}`$ as in the wind example in Fig. 1). The first peak marks the beginning of the I-shell deceleration, and it is followed by a minimum when $`\mathrm{\Gamma }`$ approaches $`\gamma `$. At $`R>2R_\beta `$, $`\gamma `$ falls down steeply and the shock becomes strong again (high $`\mathrm{\Gamma }_{\mathrm{rel}}`$), leading to the second peak at $`RR_{\mathrm{trail}}`$.

## 5 Discussion

Dynamics of neutron-fed explosions is a clean physical problem: existence of neutrons in GRBs and their mean lifetime of 15 min is all we needed to construct the model. We focused here on the simple case where all neutrons have equal Lorentz factor $`\mathrm{\Gamma }_n`$, and data may require a multi-shell picture with variable $`\mathrm{\Gamma }_n`$. This extension could affect our model in case of high variations $`\mathrm{\Delta }\mathrm{\Gamma }_n>\mathrm{\Gamma }_n`$ (though at large $`R`$, where neutrons die out exponentially, it is sufficient to consider a highest $`\mathrm{\Gamma }_n^{\mathrm{max}}`$ as slower neutrons have shorter lifetimes and their population is smaller by $`\mathrm{exp}[\mathrm{\Gamma }_n^{\mathrm{max}}/\mathrm{\Gamma }_n]`$). Another simple extension is beamed ejecta. Then at late stages the ion fireball spreads laterally while the neutron beaming remains constant.

In contrast to dynamics, emission is really complicated. In general, radiation from collisionless blast waves is not derived from first principles, and neutrons do not make the problem simpler. The afterglow emission is believed to be synchrotron, and it depends on poorly understood generation of magnetic field and electron acceleration. The standard model without neutrons relies on the field generation in the shock by the two-stream instability (Sagdeev 1966, Medvedev & Loeb 1999, Gruzinov 2001). The shock front is, however, extremely thin ($`\delta c/\omega _pR/\mathrm{\Gamma }^2`$), and the postshock field decays quickly. The model needs a significant remnant field in an extended layer behind the shock, which is uncertain. This problem is alleviated in neutron-fed explosions. Here, the leading neutron front is an additional dissipation region maintained in a turbulent state with generated magnetic fields. The N-shell has thickness $`\mathrm{\Delta }`$ comparable to $`R/\mathrm{\Gamma }_n^2`$. It can produce a significant synchrotron radiation by itself111The N-shell emission decays exponentially on observed timescale $`R_\beta /c\mathrm{\Gamma }_n^2=\tau _\beta /\mathrm{\Gamma }_n110`$ s and it can be related to smooth (fast-rise-exponential-decay) pulses in prompt GRBs. and leave behind remnant fields and hot plasma to be used by the ensuing shock wave for afterglow emission. We emphasize differences from a customary external shock: the neutron-trail shock propagates in a relativistically moving, dense, hot, and possibly magnetized medium. The fate of magnetic fields in such shocks and the resulting emission is an interesting issue for a future study.

The neutron impact ceases at $`R_{\mathrm{trail}}10^{17}`$ cm, which can leave an imprint on the observed afterglow. For example, the shock dissipation can have a bump (Fig. 1), and a spectral transition is also possible. The arrival time of radiation emitted at $`R_{\mathrm{trail}}`$ is $`R_{\mathrm{trail}}/2\mathrm{\Gamma }^2c`$ (counted from the arrival of first $`\gamma `$-rays). It may be as long as 30 days or as short as a few seconds, depending on the fireball Lorentz factor $`\mathrm{\Gamma }(R_{\mathrm{trail}})`$. Remarkably, $`R_{\mathrm{trail}}`$ is almost independent of the ambient medium, and its observational signature would give information on $`\mathrm{\Gamma }(R_{\mathrm{trail}})`$. Recent early observation of a GRB afterglow (GRB 021004) discovered an interesting re-brightening at $`10^3`$ s which would correspond to $`\mathrm{\Gamma }(R_{\mathrm{trail}})30(\mathrm{\Gamma }_n/300)^{1/2}`$. Also, we do not exclude a possible relevance of neutrons to the 20 day bumps observed in a few GRBs, as the time coincides with $`R_{\mathrm{trail}}/c`$.

Neutron signatures should be absent if the fireball is strongly dominated by a Poynting flux and has extremely low baryon loading. Then the neutron component decouples early, with a modest Lorentz factor $`\mathrm{\Gamma }_n`$, and decays at small radii. The upper bound on $`\mathrm{\Gamma }_n`$ due to decoupling is $`\mathrm{\Gamma }_n300(\dot{M}_\mathrm{\Omega }/10^{26})^{1/3}`$ where $`\dot{M}_\mathrm{\Omega }`$ \[g/s\] is the mass outflow rate per unit solid angle of the fireball (see § 4.1 in Paper 1).

We focused here on the neutron front and did not account for the $`\gamma `$-ray precursor that impacts the blast wave dynamics at $`R<R_{\mathrm{acc}}=0.7\times 10^{16}(E_\gamma /10^{53})^{1/2}`$ cm, where $`E_\gamma `$ \[erg\] is the isotropic energy of the GRB (see Beloborodov 2002 and refs. therein). The analysis in this Letter is strictly valid for afterglows emitted at $`R>R_{\mathrm{acc}}`$. Then the radiation-front effects, including the gap opening at $`R<0.3R_{\mathrm{acc}}`$, occur at smaller radii, and apply to an earlier afterglow. For a dense medium, where $`R_{\mathrm{dec}}<R_{\mathrm{acc}}`$, effects of the neutron and $`\gamma `$-ray fronts should be studied together.

I am grateful to P. Goldreich, R. Blandford, and TAPIR group for their hospitality during my stay at Caltech, where part of this work was done. This research was supported by NSERC and RFBR grant 00-02-16135.


text:
# Decay of scalar turbulence revisited

## Abstract

We demonstrate that at long times the rate of passive scalar decay in a turbulent, or simply chaotic, flow is dominated by regions (in real space or in inverse space) where mixing is less efficient. We examine two situations. The first is of a spatially homogeneous stationary turbulent flow with both viscous and inertial scales present. It is shown that at large times scalar fluctuations decay algebraically in time at all spatial scales (particularly in the viscous range, where the velocity is smooth). The second example explains chaotic stationary flow in a disk/pipe. The boundary region of the flow controls the long-time decay, which is algebraic at some transient times, but becomes exponential, with the decay rate dependent on the scalar diffusion coefficient, at longer times.

Study of advection of a steadily supplied scalar in smooth flow, pioneered by Batchelor 59Bat  for the case of a slowly changing strain and extended by Kraichnan to the opposite limit of short-correlated flow 67Kra , has developed into a universal theory, applicable to any kind of statistics and temporal correlations of the smooth chaotic flow 94SS ; 95CFKLa . (See also the reviews 00SS ; 01FGV .) Extension of this theory to the decay problem 99Son ; 99BF  suggested a fast (exponential) decay of the passive scalar in the case of spatially smooth (that is approximated by linear velocity profiles at relevant scales) and unbounded random flow. (We will call this case the “pure Batchelor” one.) In parallel, laboratory experiments 96MD ; 97WMG ; 00JCT ; 01GSb  and numerical simulations 93Pie ; 94HS ; 95AFO ; 97BDY  were conducted to verify the Batchelor-Kraichnan theory. However, a comparison between theory, experiment and numerical simulations was not conclusive. Different experiments and simulations showed different results, inconsistent with each other and with the theory. In this letter we explain a possible source of the discrepancy. Our major point is that the Batchelor-Kraichnan theory does not apply to the late-time stage of the scalar decay measured in the decay laboratory and numerical experiments 93Pie ; 96MD ; 97WMG ; 01GSb . This is because the long-time decay is dominated by the regions of the flow, where the mixing is not as efficient as in the pure Batchelor case.

Before constructing a quantitative theory, let us draw a qualitative physical picture. We start from the first example of a multi-scale turbulent flow with the following hierarchy of scales assumed: $`L_v\eta r_d`$, where $`L_v`$ is the energy containing scale of the flow, $`\eta `$ is the Kolmogorov (viscous) scale and $`r_d`$ is the diffusive scale of the scalar. It is known that the passive scalar mixing in the inertial range of turbulent flow is slower at larger scales, since the turnover time grows with scale. This explains the algebraic-in-time decay of the scalar for the inertial range of scales, $`L_vr\eta `$ 00EX ; 01CEFV . If we now take into account that turbulent velocity is smooth at scales smaller than $`\eta `$, it becomes important to understand what happens with the scalar fluctuations at these small scales. We show in the letter that scalar decay is accompanied by both upscale and downscale transport, however, the situation is asymmetric. The decay of scalar fluctuations in time is due to inertial interval physics, while the spatial correlations of the scalar at the smallest scales are mainly controlled by the viscous part of the flow. (Here, the spatial correlations of the scalar are logarithmic, like in the stationary case 95CFKLa .) An analogous situation is realized if the scalar is advected by a chaotic (spatially smooth) flow in a box or a pipe. Advection of the scalar slows down in the boundary region, while mixing continues to be efficient in the bulk, i.e. far from the boundaries. The result is that the boundary region supplies the scalar into the bulk. Plumes (blobs) of the scalar, injected from the boundary domain into the box center, cascade in a fast way down to $`r_d`$ where diffusion smears it out. Therefore, the temporal decay of the scalar is mainly controlled by the rate of scalar injection from the boundary layer. This causes an essential slow down of the decay in comparison with the pure Batchelor case, since in this case there are no stagnation (boundary) regions. The decay of the scalar is algebraic during the transient stage, when the scalar is supplied to the bulk from the boundary region which is wider than the diffusive boundary layer. The boundary region width decreases with time. Whenever the boundary layer width becomes of the order of the diffusive layer width the decrease stops and the algebraic temporal decay turns into an exponential one. The exponential decay rate appears to be parametrically smaller than in the pure Batchelor case.

In all the cases discussed, we consider decay of a passive scalar, $`\theta `$, described by the equation

$`_t\theta +(𝒗)\theta =\kappa ^2\theta ,`$ (1)
where $`𝒗`$ is the flow velocity and $`\kappa `$ is the diffusion coefficient. Our goal is, starting from the equation (1) and assuming some statistical properties of the velocity field $`𝒗`$, to describe statistics of the passive scalar decay. In order to understand the qualitative features of the passive scalar decay, we choose (in the spirit of Kraichnan 94Kra ) the simplest possible case of a short-correlated in time velocity, possessing Gaussian statistics. This model enables us to obtain analytically a time evolution of the passive scalar correlation functions. Two situations are discussed. The first is of a statistically homogeneous flow (infinite volume) with a multi-scale velocity which is smooth at small scales, and possesses a scaling behavior at larger scales. The Reynolds number $`Re=(L_v/\eta )^{4/3}`$ and the Prandtl number $`Pr=(\eta /r_d)^2`$ are both assumed to be large. The second situation is $`2d`$ chaotic flow confined within a disk. We also discuss an extension to the case of a more realistic pipe flow. The flow is assumed to be smooth and obeying the standard viscous behavior in the neighborhood of the boundary: the longitudinal and transverse components of the velocities tend to zero quadratically and linearly (respectively) with the distance to the boundary.

We begin with an analysis of the homogeneous unbounded case. The statistics of a short-correlated in time Gaussian velocity is completely characterized by its pair correlation function

$`v_\alpha (t_1,𝒓)v_\beta (t_2,0)=\delta (t_1t_2)\left[V_0\delta _{\alpha \beta }𝒦_{\alpha \beta }(𝒓)\right],`$ (2)
$`𝒦_{\alpha \beta }{\displaystyle \frac{r}{2}}W^{}(r)\left(\delta _{\alpha \beta }{\displaystyle \frac{r_\alpha r_\beta }{r^2}}\right)+{\displaystyle \frac{d1}{2}}W(r)\delta _{\alpha \beta },`$ (3)
where $`d`$ is dimensionality of space, $`V_0`$ is a constant characterizing fluctuations of the velocity at the integral scale, $`L_v`$, and $`W(r)`$ is a function, which determines the scale-dependence of velocity fluctuations. To describe both viscous and inertial ranges, we introduce

$`W=Dr^2/[1+(r/\eta )^\gamma ],`$ (4)
where $`D`$ stands for the amplitude of the velocity fluctuations and $`\gamma `$ characterizes the velocity scaling in the inertial interval. The Kolmogorov scale $`\eta `$ separates the inertial and viscous intervals. The diffusive length is now expressed as $`r_d=\sqrt{\kappa /D}`$. Measuring time, $`t`$, in the units of $`1/D`$, one replaces $`D`$ by unity in all the forthcoming formulas. The object we examine here is the simultaneous pair correlation function of the passive scalar, $`F(t,𝒓)=\theta (t,𝒓)\theta (t,0)`$, which is the simplest non-zero correlation function in the homogeneous case. Under the assumption of isotropy one derives from Eqs. (1-3)

$`r^{1d}_r\left\{r^{d1}\left[W(r)+\kappa \right]_rF_\lambda \right\}=\lambda F_\lambda ,`$ (5)
where $`F_\lambda (r)`$ is the Laplace transform of $`F(t,r)`$ with respect to $`t`$. The operator on the left-hand side of Eq. (5) is self-adjoint and non-negative (with respect to the measure $`\text{d}rr^{d1}`$ and under condition that $`F(r)`$ is smooth at the origin and $`_rF(r=0)=0`$).

We assume that $`\text{d}𝒓\theta =0`$ (the Corrsin invariant is zero), and also, that initially the scalar is correlated at the scale $`r_0`$ from the viscous-convective range, $`r_dr_0\eta `$, i.e. $`F(0,r)`$ is $`F(0,0)`$ at $`r<r_0`$ and decays fast enough at $`r>r_0`$. Then some part of the initial evolution of $`F(t,r)`$ is not distinguishable from the decay in the pure Batchelor case, described by the limit $`\eta \mathrm{}`$ in Eq.(4). Let us briefly describe this initial stage of the dacay. The solutions of Eq. (5) are power-like, $`F_\lambda r^{\pm \sqrt{d^2/4\lambda }d/2}`$, outside the diffusive range, $`rr_d`$. The functions are normalizable only if $`\lambda >d^2/4`$. Therefore there is a gap in the spectrum ($`0<\lambda <d^2/4`$ are forbidden). The existence of this gap guarantees an exponential decay of $`F(t,r)`$ with time. This feature is, of course, in agreement with the previous analysis of the pure Batchelor case 99Son ; 99BF . A complete set of functions at $`rr_d`$ is $`F_\lambda =r^{d/2}\mathrm{cos}[\alpha _k+k\mathrm{ln}(r/r_d)]`$, where $`k`$ is a positive number related to $`\lambda `$ via $`\lambda =d^2/4+k^2`$, and $`\alpha _k`$ are phases, determined by matching at $`rr_d`$. Taking into account the orthogonality condition for $`F_\lambda `$, $`\text{d}rr^{d1}F_\lambda (r)F_\lambda ^{}(r)=\pi \delta (kk^{})/2`$, one can express $`F(t,r)`$ via $`F(0,r)`$. Analysis of this expression leads to the following picture. The correlation function, $`F(t,r)`$, does not change at $`r<r_{}=r_0\mathrm{exp}(td)`$. At the scales larger than $`r_{}`$, $`F`$ decays exponentially, $`\mathrm{exp}[td^2/4]`$. Note, that $`r_{}`$, also, marks the position of the maximum of the spatial derivative of $`F(t,r)`$, $`_rF(t,r)`$. Therefore, $`r_{}(t)`$ describes the front running downscale (from $`r_0`$). The front reaches $`r_d`$ at $`t_d=\mathrm{ln}(r_0/r_d)/d`$, so that at $`t>t_d`$, $`F(t,r)`$ decays exponentially at all the small scales. Complimentary to the front running downscale, there exists another front running upscale. Indeed, the integral quantity, $`_0^s\text{d}rr^{d1}F(t,r)`$, which determines the overall amount of the scalar at all the scales smaller than $`s`$, achieves its maximum around $`r_+=r_0\mathrm{exp}(td)`$.

At $`t>t_\eta \mathrm{ln}[\eta /r_0]/D`$ the pure Batchelor description ceases to be valid and one does need to account for complete multiscale form of eddy-diffusivity function $`W(r)`$ given by Eq. (4). Our further analysis is devoted to this general case. The spatial decay of the functions $`F_\lambda `$ in the multi-scale case is not as steep as in its pure Batchelor one. It is convenient to change to a new field, $`\mathrm{\Phi }_\lambda `$, $`F_\lambda =r^{1d}_r[r^{d/2}\mathrm{\Phi }_\lambda (z)]`$, where $`z=2\sqrt{\lambda }r^{\gamma /2}/\gamma `$. The general solution of the resulting equation for $`\mathrm{\Phi }_\lambda (z)`$ is a linear combination of the Bessel functions, $`\mathrm{\Phi }_\lambda =J_{\pm \nu }(z)`$, where $`\nu =2\sqrt{d^2/4\lambda }/\gamma `$. If $`\lambda <d^2/4`$, the positive root for $`\mathrm{\Phi }_\lambda `$ has to be chosen to satisfy the matching conditions at $`r_d`$. The eigen function is normalizable for any positive $`\lambda `$, so that even the smallest positive $`\lambda `$ are not forbidden. The lack of the gap in the spectrum of the multi-scale model means an algebraic in time decay of $`F(t,r)`$, i.e. much slower decay than the one found for the Batchelor-Kraichnan model. We present here the expression for the long-time, $`tt_\eta `$, asymptotic behavior of $`F(t,r)`$:

$`F(t,r)\{\begin{array}{cc}\left[td+\mathrm{ln}(r/\eta )+(r/\eta )^\gamma /\gamma \right]^{1d/\gamma },& rr_+;\\ t^{5/4d/2\gamma }r^{d/2+\gamma /4},& rr_+\end{array}`$ (8)
Here $`r_+=(\gamma t)^{1/\gamma }`$ stands for position of the front running upscale, $`r_+`$ lies in the scaling region $`r_+\eta `$, where the multi-scale model turns into the so-called Kraichnan model 94Kra , $`Wr^{2\gamma }`$. (Therefore, it is not surprising, that our result (8) is consistent at $`rr_+`$ with the expression for the pair correlation function derived before for the Kraichnan model 00EX ; 01CEFV  in the regime of zero Corrsin invariant.) Note, that when $`tt_\eta `$ is yet moderate in value, the small-scale part of the asymptotic expression (8) is correct only at $`rR_{}(t)=\eta e^{td}`$, where $`R_{}`$ is the position of the front initiated at $`tt_\eta `$ and running downscale from $`\eta `$. At $`tt_\eta +\mathrm{ln}(\eta /r_d)/d`$, the downscale front reaches the dissipative scale, and afterwards the expression (8) is correct for any scales $`rr_d`$. Eq. (8) shows that just like in the pumped case 59Bat ; 67Kra ; 95CFKLa  the scalar structure function is logarithmic, $`[\theta (t,𝒓)\theta (t,0)]^2t^{2d/\gamma }\mathrm{ln}[r/r_d]`$, at $`r_dr\eta `$ and large times, where, therefore, the time-dependent factor at the logarithm can be interpreted as a scalar flux from the inertial range, $`r\eta `$, downscale to the viscous range, $`r\eta `$.

To conclude, the existence of the inertial region in velocity affects drastically the spatio-temporal distribution of the scalar in the viscous range (where the velocity can be approximated by linear profiles). The inertial interval, where the scalar is mainly concentrated, serves as a kind of reservoir for the smaller scales. This explains why the decay of scalar correlation function in the viscous is much slower (algebraic) than in the pure Batchelor case (when it would be exponential).

Now we proceed to the situation of spatially bounded flows. Aiming to demonstrate the main qualitative features of the confined geometry, that is strong sensitivity of the scalar decay to the peripheral (close to the boundary) part of the flow, we examine a model case of $`2d`$ chaotic (i.e. consisting of only few spatial harmonics) flow inside a disk. (It is, actually, clear that the picture of temporal decay and spatial distribution of the scalar, described below, is universal, i.e. it applies to a typical chaotic flow confined in a close box, particularly in $`3d`$.) Incompressible flow in $`2d`$ can be characterized by a stream function $`\psi `$, then the radial and azimutal components of the velocity are $`v_r=_\phi \psi /r`$ and $`v_\phi =_r\psi `$. Our model of $`\psi `$ is

$`\psi ={\displaystyle \frac{\xi _1}{2}}r^2U(r)\mathrm{sin}(2\phi )+{\displaystyle \frac{\xi _2}{2}}r^2U(r)\mathrm{cos}(2\phi ),`$ (9)
$`\xi _i(t_1)\xi _j(t_2)=2D\delta _{ij}\delta (t_1t_2),`$ (10)
where $`U(r)`$ is a function of $`r`$, finite at the origin, $`U(0)=1`$, and becoming zero, together with its first derivative, at the disk boundary, $`r=1`$, and $`\xi _{1,2}(t)`$ are zero mean short correlated random Gaussian functions. The value of the passive scalar, $`\theta `$, averaged over the statistics of $`\xi `$, $`\theta (t,𝒓)`$, is the object of our interest here. One examines how the average concentration of the scalar evolves with time at different locations $`𝐫`$ within the disk. The short-correlated feature of the velocity field allows a closed description for $`\theta `$ in terms of a partial differential equation. Considering the spherically symmetric part of $`\theta `$ only (asymptotically, at large times only this $`varphi`$-independent part remains essential), and passing from $`\theta `$ and $`r`$ to $`\mathrm{{\rm Y}}`$ and $`q`$, respectively, where $`\theta =r^1_r[r^2\mathrm{{\rm Y}}(t,q)]`$ and $`q=\mathrm{ln}r`$, one finds that the Laplace transform of $`\mathrm{{\rm Y}}`$ with respect to $`t`$ (measured, again, in $`D^1`$ units) satisfies

$`(U^2+\kappa e^{2q})(_q^2\mathrm{{\rm Y}}_\lambda 2_q\mathrm{{\rm Y}}_\lambda )+\lambda \mathrm{{\rm Y}}_\lambda =0.`$ (11)
The diffusion term is important only at the center of the disk and near the boundaries. Analysis of Eq. (11) is straightforward but bulky. Below we will present only selected details of the analysis, aiming to describe the general picture of the phenomenon. (The complete account for the derivation details will be published elsewhere 02CLT .)

In the bounded flow the decay of passive scalar splits into three distinct stages. The major effect dominating the first stage (just as in the pure Batchelor case, explained by Eq. (11) with $`U1`$) is formation of elongated structures (stripes) of the scalar in the bulk region of the flow. The stripes are getting thinner with time, i.e. inhomogeneities of smaller and smaller scales are produced. Once the width of the stripes decreases down to the dissipative scale, $`r_d\sqrt{\kappa }`$, the stripes are smeared out by diffusion. The stretching-contraction process is exponential in time, so that the initial stage when the stripes are formed lasts for $`\tau _1\mathrm{ln}[1/r_d]`$. By the end of this first stage the scalar is exhausted in the central region of the flow. The stretching rate, however, is smaller in the peripheral domain than in the bulk. Thus the scalar literally remains longer in the peripheral domain. This defines the second, transient, stage. At $`\kappa ^{1/4}q1r\lambda ^{1/4}1`$, the solution of Eq. (11) is $`\mathrm{{\rm Y}}_\lambda q\mathrm{sin}\left(\alpha _\lambda +\sqrt{\lambda }q/U\right)`$. In the other asympotitc region at $`q\lambda ^{1/4}`$ one can drop the second derivative (with respect to $`q`$) and diffusive terms in Eq. (11). One finds that, $`_q\mathrm{ln}\mathrm{{\rm Y}}_\lambda 1`$ at $`q\lambda ^{1/4}`$. Matching those two asymptotic regions at $`q\lambda ^{1/4}`$ we come to the $`\alpha _\lambda =0`$ condition. Once the asymptotic form of the eigen-function $`\mathrm{{\rm Y}}_\lambda `$ is known, it is straightforward to restore the behavior of $`\mathrm{{\rm Y}}(t,q)`$. One finds that $`\mathrm{{\rm Y}}(t,q)`$ (and, therefore, $`\theta `$) is concentrated in the $`\delta (t)1/\sqrt{t}`$-small (and shrinking with time) vicinity of the boundary. In the domain outside of the shrinking layer the decay is algebraic, $`\theta t^{3/2}q^3`$. The second stage lasts for $`\tau _2\kappa ^{1/2}`$, i.e. until $`\delta (t)`$ shrinks down to $`r_{bl}=\kappa ^{1/4}`$, which is the width of the diffusive boundary layer. Account for diffusivity in the $`q1`$ analysis gives yet another matching condition for $`\mathrm{{\rm Y}}_\lambda `$, at $`rr_{bl}`$, resulting in formation of a discrete spectrum for $`\lambda `$. The smallest $`\lambda `$ (and therefore the value of the level spacing in the discrete spectrum) is estimated by $`\sqrt{\kappa }`$. (Using the periodic character of the $`\mathrm{sin}`$-function, one finds that the $`n`$-th level eigenvalue $`\lambda _n`$ is estimated by $`n^2\sqrt{\kappa }`$ at large $`n`$.) Therefore, at $`t\tau _2`$, $`\theta `$ decays exponentially, $`\mathrm{exp}(t/t_d)`$, where $`t_d\kappa ^{1/2}`$. The smallest $`\lambda `$ eigen-function is localized at $`qr_{bl}`$. Thus the third (final) stage of evolution is characterized by the majority of the scalar remaining in the diffusive boundary layer.

To conclude, the temporal behavior of the passive scalar correlations, predicted by the bounded flow theory, is complicated. It involves two different exponential regimes separated by an algebraic one, so that a special accuracy is required in order to quantify the theory against various experiments in the field.

The bounded flow theory can be applied to the experiment of Groisman and Steinberg 01GSb , where the passive scalar is advected through a pipe by dilute-polymer-solution flow. The polymer-related elastic instability makes the flow chaotic. The mean flow along the pipe is also essential, so that according to the standard Taylor hypothesis, measurements of scalar concentration at different positions along the pipe are interpreted as correspondent to different times in an artificial decay problem. Exponential decay of the scalar was reported in 01GSb  at long times (long pipes). The theoretical picture explaining the experiment is similar to the one proposed above for the disk. The width of the diffusive boundary layer near the pipe boundaries is estimated by $`RPe^{1/4}`$, where $`R`$ is the pipe radius and $`Pe`$ is the Peclet number, $`Pe=R^2\sigma /\kappa `$ and $`\sigma `$ is a typical value of the velocity gradient ($`Pe`$ is $`10^4`$ in the conditions of 01GSb ). The characteristic decay time (fixed by the diffusive boundary layer width) is estimated by $`t_d\sigma ^1Pe^{1/2}`$. However, one should be careful when converting this time to the decrement of the passive scalar decay along the direction of the flow (pipe). The average velocity of the flow tends to zero near the boundary of the pipe, making the advection near the boundary less efficient than in the bulk. Assuming a linear profile of the average velocity near the boundary, one gets the following estimate for the law of scalar fluctuation decay along the pipe, $`\mathrm{exp}(\gamma z/u_0)`$, where $`z`$ is the coordinate in the direction of the pipe, $`u_0`$ is the average velocity at the center of the pipe, and $`\gamma \sigma Pe^{1/4}`$. The factor $`Pe^{1/4}`$ here is a manifestation of a slow down of the passive scalar decay in comparison with the pure Batchelor case. The same factor characterizes decay of higher order correlation functions of the scalar as well. For example, one can consider the passive scalar pair correlation function $`F(r)`$ near the center of the disk. One finds that $`F(r)r^\alpha `$, where $`\alpha `$ is small and can be estimated as $`\alpha \text{Pe}^{1/4}`$. (Note, that this behavior is almost indistinguishable from the logarithmic one, found for the stationary case 59Bat ; 67Kra ; 95CFKLa .) This estimations agree with the experimental data of 01GSb .

Let us now briefly discuss two other available numerical and experimental results. Pierrehumbert 93Pie  reported exponential decay of the scalar correlations in his numerical experiment with chaotic map velocity in a periodic box. This observation is consistent with our results, since the periodic boundary conditions are less restrictive than the zero condition for velocity at the box boundary, and, therefore, should leave a finite gap in the $`\lambda `$-spectrum even in the limit $`\kappa 0`$. In the experiment of Jullien, Castiglione and Tabeling on $`2d`$ stationary flow steered by magnets in a finite volume beneath the $`2d`$ layer 00JCT , the decay rate of the scalar seems slower than exponential (Fig. 2 of 00JCT ). This observation is in agreement with the absence of the gap in the $`\lambda `$-spectrum (at $`\kappa 0`$) we found for the finite box model.

The brevity of this letter does not allow us to discuss effects of intermittency (which manifests in high-order moments of the scalar), of higher order angular harmonics (in the disk or pipe geometry), and details of the spatial distribution of the scalar in a variety of other inhomogeneous cases (e.g. periodic flow). This detailed discussion is postponed for a longer paper to be published elsewhere 02CLT  (where we will also present results of direct numerical simulations of the passive scalar decay for these cases).

We thank B. Daniel, R. Ecke, G. Eyink, G. Falkovich, A. Fouxon, A. Groisman, I. Kolokolov, M. Riviera, V. Steinberg, P. Tabeling and Z. Toroczkai for helpful discussions.