text:
# Incoherent Noise and Quantum Information Processing

## I Introduction

One of the biggest challenges in Quantum Information Processing (QIP) is the precise control of quantum systems. Errors in the control are conveniently classified as coherent, decoherent, and incoherent PraviaRFI . Coherent errors are systematic and differ from the desired operation by a unitary operation. Decoherent errors can be expressed by completely positive (CP) superoperators Alicki  and can be counteracted by techniques such as Quantum Error Correction (QEC) Shor ; NielsenChuang . An incoherent process can also be described by a completely positive superoperator PraviaRFI ; BoulantEntSwap . The apparent non-unitary behavior of the incoherent process arises due to a distribution over external experimental parameters. The incoherent process is described by a superoperator $`S`$ acting on Liouville space which can be written, when acting on columnized density matrices $`|\rho `$ obtained by stacking their columns on top of each other from left to right Havel:03 , as

$`S`$ $`=`$ $`{\displaystyle p(\alpha )\overline{U}(\alpha )U(\alpha )𝑑\alpha },`$ (1)
where $`p(\alpha )`$ is a probability density, i.e. the fraction of quantum systems within an ensemble that sees a given $`U(\alpha )`$ within an interval $`d\alpha `$, $`p(\alpha )𝑑\alpha =1`$, and $`\overline{U}`$ denotes the complex conjugate of $`U`$. This decomposition of a CP map into unitary Kraus operators is sometimes called a random unitary decomposition (RUD) LeungThesis . A RUD exists for an incoherent process, but such a decomposition is sometimes possible even for a very general decoherent process Streater  when there is no correspondence between the unitary operators in the decomposition and some actual distribution of associated experimental control parameters $`\alpha `$. The distinction between the two therefore is practical, and depends primarily on the correlation time of the variation of experimental parameters. If the latter quantity is longer than the inverse of the typical modulation frequency, the process falls into the class of incoherent noise Carr ; Hahn . The point of making this distinction is that, whilst correcting for decoherent errors requires the full power of QEC, in practice incoherent noise effects are often reduced directly through the design of the time-dependence of the control fields. This is possible since the operators underlying the incoherence $`U(\alpha )`$ are assumed to be time-independent over the length of the expectation value measurement. Common approaches for instance in Nuclear Magnetic Resonance (NMR) include composite and adiabatic pulses Tycko ; Shaka ; Levitt ; Jones ; Pines ; Silver . Furthermore, the work done by Tycko Tycko  and Jones Jones2 ; Jones3  on composite pulses finds a great application in QIP since the schemes proposed are universal and therefore work regarless of the input state.

In this paper, we demonstrate how incoherent errors introduce particular limitations to Quantum Process Tomography (QPT) NielsenChuang ; Childs ; BoulantQPT ; Cirac  due to the correlations they introduce with an environment in the QPT input states. Prior work has been devoted to the study of the implications of such correlations in the system’s reduced dynamics Buzek ; BuzekErra ; Pechukas . However, to our knowledge they have not explicitly been studied within the context of QPT to show that the method may output non-completely positive (NCP) maps. If the existence or origin of such noise is unknown, it is shown in PraviaRFI  by means of explicit examples how one can eventually infer qualitative information, say its symmetry, about the probability distribution underlying the incoherent process from superoperator eigenvalue spectra. In this paper we tackle the inverse problem of extracting an effective probability distribution $`p(\alpha )`$ representing the incoherent superoperator, given a model for the source of the incoherent noise in the system. Such information is crucial for counteracting the incoherent errors PraviaRFI  which, due to their slow variation would otherwise persist during the experiment and quickly accumulate.

## II Incoherent Noise and Quantum Process Tomography

QPT measures the experimental map associated with the implementation of a desired quantum operation by passing a complete set of input states through the gate and measuring the corresponding output states (see Fig. 1).

This procedure is important for the experimental study of noise processes and for the design of quantum error correcting codes Shor . If incoherent noise is present in the preliminary step of QPT, then the prepared (input) states will be classically correlated with the control parameter $`\alpha `$ characterizing the incoherence (Eq 1). Furthermore, if the correlation time of the noise in the control parameter is long compared to the coarse-grained time at which the evolution of the system is monitored, then the subsequent dynamics is non-Markovian Cohen-Tannoudji . More specifically, the gate applied to the input state (i.e., the gate being characterized by QPT) will be correlated with the same (slowly-varying) control parameter, and therefore also correlated with the input state to which it is applied. In such cases, the measured dynamics are not guaranteed to be completely positive, and need not correspond to a linear map NielsenChuang ; Buzek ; BuzekErra ; Pechukas .

More generally, any correlations arising from non-Markovian dynamics, whether quantum or classical, can lead to incorrect interpretations of the data obtained from QPT. The “environment” which we assume from the outset is correlated with the system is in general defined by the degrees of freedom that are not part of the system. For example, the different spatial locations of individual qubits in a NMR ensemble, or the bosonic bath producing the fluctuations of the gate charge in a superconducting qubit Weiss .

The basic issues can be seen by exploring the QPT of a spin-$`1/2`$ system $`A`$, coupled to a second spin-$`1/2`$ system $`B`$ as its environment. Borrowing the example used in Buzek , the initial density matrix of the total system may be written as

$`\rho _{AB}={\displaystyle \frac{1}{4}}(I_{AB}+\alpha _i\sigma _iI_B+\beta _jI_A\sigma _j+\gamma _{ij}\sigma _i\sigma _j),`$ (2)
where $`I`$ and $`\sigma _i`$ denote the identity and Pauli matrices respectively. The density matrix of the system $`A`$ is then obtained by tracing over the environment $`B`$:

$`\rho _A={\displaystyle \frac{1}{2}}(I_A+\alpha _i\sigma _i).`$ (3)
The dynamics of the whole system $`\rho _{AB}`$ is a unitary evolution $`U_{AB}`$, so that the resulting density matrix of $`A`$ is

$`\rho _A^{}`$ $`=`$ $`{\displaystyle \underset{\mu }{}}\mu |U_{AB}(\rho _A\rho _B)U_{AB}^{}|\mu `$
$`+{\displaystyle \underset{\mu }{}}\mu |U_{AB}\gamma _{ij}^{^{}}\sigma _i\sigma _jU_{AB}^{}|\mu ,`$
where $`\gamma _{ij}^{^{}}=(\gamma _{ij}\alpha _i\beta _j)/4`$ and $`\rho _B=\mathrm{Tr}_A(\rho _{AB})`$. Let $`\rho _B=_\nu p_\nu |\nu \nu |`$ and $`M_{\mu \nu }=\mu |U_{AB}|\nu `$. Then the above expression can be reexpressed as

$`\rho _A^{}`$ $`=`$ $`{\displaystyle \underset{\mu \nu }{}}M_{\mu \nu }\rho _AM_{\mu \nu }^{}`$
$`+{\displaystyle \underset{\mu }{}}\mu |U_{AB}\gamma _{ij}^{^{}}\sigma _i\sigma _jU_{AB}^{}|\mu `$
The first line therefore corresponds to the Kraus operator sum form Kraus2  of the evolution when initial correlations are not present, while the second line represents the contribution from these correlations. We can easily see that there exist $`U_{AB}`$, e.g., the swap gate, for which these initial correlations are not observable. In general, however, initial correlations cause the map to be non-linear or NCP.

Within the context of QPT, we now investigate an explicit example of how NCP superoperators can arise. We take a set of $`4`$ initial density matrices $`\rho _{AB}`$ such that $`\rho _B`$ is the same in each case, and where the input states $`\rho _A`$ span the Hilbert space of the system $`A`$ (as required by the QPT procedure)

$`\rho _{AB}^1`$ $`=`$ $`(I_{AB}+\beta I\sigma _z)/4\rho _{A,in}^1=I/2`$
$`\rho _{AB}^2`$ $`=`$ $`(I_{AB}+\alpha \sigma _xI+\beta I\sigma _z+\gamma \sigma _x\sigma _z)/4\rho _{A,in}^2=(I+\alpha \sigma _x)/2`$
$`\rho _{AB}^3`$ $`=`$ $`(I_{AB}+\alpha \sigma _yI+\beta I\sigma _z+\gamma \sigma _y\sigma _z)/4\rho _{A,in}^3=(I+\alpha \sigma _y)/2`$
$`\rho _{AB}^4`$ $`=`$ $`(I_{AB}+\alpha \sigma _zI+\beta I\sigma _z+\gamma \sigma _z\sigma _z)/4\rho _{A,in}^4=(I+\alpha \sigma _z)/2`$
where in each case $`\rho _B=(I+\beta \sigma _z)/2`$. With the example $`U_{AB}=e^{i\frac{\pi }{4}\sigma _z\sigma _z}`$, the corresponding $`4`$ outputs are $`\stackrel{~}{\rho }_A^1=I/2`$, $`\stackrel{~}{\rho }_A^2=(I+\gamma \sigma _y)/2`$, $`\stackrel{~}{\rho }_A^3=(I\gamma \sigma _x)/2`$ and $`\stackrel{~}{\rho }_A^4=(I+\alpha \sigma _z)/2`$. We write the density matrices $`\rho `$ as vectors in Liouville space in the Zeeman basis Havel:03 , which are obtained by first writing the density matrix in the Zeeman basis and then stacking their columns on top of each other from left to right. We refer to the resulting vector simply as the “columnized density matrix”, and will denote it as a ket $`|\rho `$. If we set $`\alpha =\beta =0.5`$ and $`\gamma =0.6`$, the map $`S`$ is

$$\begin{array}{c}\hfill \begin{array}{ccccc}& \text{Output 1}& \text{Output 2}& \text{Output 3}& \text{Output 4}\end{array}\\ \hfill \left[\begin{array}{ccccc}& 0.5& 0.5& 0.5& 0.75\\ & 0& 0.3i& 0.3& 0\\ & 0& 0.3i& 0.3& 0\\ & 0.5& 0.5& 0.5& 0.25\end{array}\right]\end{array}\begin{array}{cc}& \begin{array}{cccc}\text{Input 1}& \text{Input 2}& \text{Input 3}& \text{Input 4}\end{array}\hfill \\ & =S\left[\begin{array}{ccccc}& 0.5& 0.5& 0.5& 0.75\\ & 0& 0.25& 0.25i& 0\\ & 0& 0.25& 0.25i& 0\\ & 0.5& 0.5& 0.5& 0.25\end{array}\right]\hfill \end{array}\begin{array}{cc}& \\ \hfill S=\left[\begin{array}{ccccc}& 1& 0& 0& 0\\ & 0& 1.2i& 0& 0\\ & 0& 0& 1.2i& 0\\ & 0& 0& 0& 1\end{array}\right]\end{array}$$
(4)
which is in general non-linear, since one can no longer predict the output for an arbitrary input state given the action of the gate on these four specific input states alone. However, the map can be considered to act linearly on the system $`A`$ Hilbert space, i.e. on the linear combinations of input states which contain the right correlations with the environment. For instance, the action of the gate on the input state $`(\rho _A^2+\rho _A^3)/2=(I+\frac{\alpha }{2}(\sigma _x+\sigma _y))/2`$ can be computed by using the above matrix expression if the total input state is $`\rho _{AB}=(I_{AB}+\frac{\alpha }{2}(\sigma _x+\sigma _y)I+\beta I\sigma _z+\frac{\gamma }{2}(\sigma _x+\sigma _y)\sigma _z)/4`$. This result conveniently allows one to treat the map as linear. If treated as linear, the Choi matrix Choi  $`𝒞=_{i,j=0}^{N1}(E_{ij}I)S(IE_{ij})`$, where $`N`$ is the dimension of the system’s Hilbert space and $`E_{ij}`$ is the $`N\times N`$ elementary matrix (with a ”1” in the $`ijth`$ position and zeros elsewhere), corresponding to the superoperator $`S`$ is not positive semidefinite and consequently $`S`$ can not be CP Havel:03 .

It is suggested in Havel:03  how the NCP part of the superoperator can be removed, namely, by removing the negative eigenvalues of the Choi matrix and then renormalizing so that the trace is equal to the dimension of the Hilbert space, $`N`$. We shall call this method CP-filtering. The Choi matrix corresponding to $`S`$ in (4) has two non-zero eigenvalues $`(2.2,0.2)`$. In this example, the CP-filtering procedure replaces the negative eigenvalue by $`0`$ and then renormalizes the new Choi matrix to trace $`N=2`$. The CP-filtering method outputs one unitary Kraus operator. On the other hand, if no initial correlations were present, the superoperator would be equal to $`S=\mathrm{diag}(1,0.5i,0.5i,1)`$, and the corresponding Choi matrix would have two positive eigenvalues $`(0.5,1.5)`$, yielding two Kraus operators. This superoperator is therefore completely positive. Unsurprisingly, we thus observe that the superoperator obtained from the CP-filtering procedure is not equal to the superoperator obtained by removing the initial correlations. Therefore, unless the initial correlations are very small, CP-filtering is a fairly uncontrolled procedure, giving CP superoperators that may significantly misrepresent the true quantum dynamics.

We now take the case $`\alpha =\beta =\gamma =0.5`$ with the same $`U_{AB}`$ as previously. Although initial correlations are still present, the superoperator obtained by the above QPT method, without CP-filtering, is CP with one Kraus operator. In contrast, if the correlations in the initial states are removed whilst keeping all other things equal, the superoperator obtained is CP with two Kraus operators, not one. So even when initial correlations are present, the superoperator obtained via QPT may be completely positive. Therefore, one cannot rule out the presence of initial correlations merely by the existence of a valid Kraus operator sum form via QPT. Initial correlations can masquerade as CP maps, and in reality the process may not be linear with respect to arbitrary input states.

To summarize, the results of this section are: (i) incoherent errors introduce correlations between the system and the environment in the QPT input states which can persist during the implemented transformation, (ii) these correlations can yield non-completely positive superoperators or non-linear maps, (iii) the CP-filtering method suggested in Havel:03  is not equivalent to removing these initial correlations, and (iiii) initial correlations can masquerade as CP maps which misrepresent what is in reality a non-linear process with respect to the input states. See Table 1 for a summary. This simple analysis explains the apparent NCP behavior measured in experiments reported in BoulantQPT ; WeinsteinQFT . This motivates the need to characterize the incoherent noise, and to provide ways to correctly interpret QPT data. If the noise can be successfully characterized, we may use this information to better counteract the noise in the first place.

## III Extracting Probability Distributions from Superoperator Eigenvalue Spectra

By applying a first order perturbation theory analysis of the eigenvalues of superoperators, we now present a method to extract the probability distribution profile $`p(\alpha )`$ of unitary matrices present in incoherent processes. For an incoherent noise to be refocused Hahn , knowledge about the noise is *a priori* required. If qualitative information about the inhomogeneity in the Hamiltonian is known, spectroscopic techniques can be used to obtain the missing quantitative information. In the following analysis, we assume that the physical origin of the incoherent noise is unknown or hidden due to the complexity of the system-apparatus interactions, but that a mathematical model is presumed.

As presented in the introduction, an incoherent process implies a random unitary distribution. Incoherent processes are thus unital, which means that the maximally mixed density matrix is left unchanged. A linear, completely positive, trace preserving and unital map is called a doubly stochastic map Alberti . Although a single necessary and sufficient condition for a doubly stochastic map to possess a RUD has not been found to our knowledge, examples of doubly stochastic maps which do not possess a RUD are reported in LeungThesis ; Streater . However, since many decoherent unital processes can be modelled by a stochastic Hamiltonian, i.e. semiclassically, we believe that many instances of decoherent processes can have a RUD. This belief is supported by the following two facts. Any two density operators $`(\rho ,\rho ^{})`$ connected by a doubly stochastic map, $`\rho ^{}=\mathrm{\Lambda }(\rho )`$, can always be related by a transformation of the form $`\rho ^{}=_ip_iU_i\rho U_i^{}`$, where $`_ip_i=1`$ and $`\{U_i\}`$ is a set of unitary operators Alberti . Furthermore, all unital maps for a two-level quantum system always have a RUD Streater .

### III.1 Perturbation Theory Analysis of the Eigenvalue Spectra

In what follows, we take examples from NMR Ernst  physics where the main source of incoherence comes from Radio Frequency (RF) power inhomogeneity. Due to the spatial extent of the sample, individual spins during the course of a RF field see different powers PraviaRFI  and evolve according to different unitary evolutions with different characteristic frequencies. Note that identical spins can have different resonance frequencies due to inhomogeneity in the static magnetic field within the ensemble, which is another source of incoherence. However, as shown in BoulantEntSwap , the non-unitary features arising from this static external field inhomogeneity are usually much smaller than those arising from RF inhomogeneity and will be therefore ignored in this example. Finally while the distribution of RF fields can be easily measured via a nutation experiment on a single spin, the method presented here is quite general and can for instance account for the correlation between multiple sources of incoherence (several RF fields, DC field etc…).

Let $`n`$ be the number of spin-1/2 particles in the ensemble, and $`U_k`$ denote the unitary operator generated by the RF field in the $`k`$th frequency interval of the RF amplitude profile. The eigenvalues of the superoperator $`S=\overline{U}_kU_k`$ are products of the eigenvalues of $`U_k`$ with those of $`\overline{U}_k`$. This yields $`2^n`$ eigenvalues that are equal to unity and $`(2^{2n1}2^{n1})`$ pairs of eigenvalues $`(\lambda ,\overline{\lambda })`$. In general, the eigenvalues of CP superoperators come in conjugate pairs, but only in the case of unitary superoperators do all the eigenvalues lie on the complex unit circle.

The incoherent process resulting from an inhomogeneous distribution of $`U_k`$ processes is given by the superoperator $`S=_kp_k\overline{U}_kU_k`$, where $`p_k`$ is the fraction of spins that sees the unitary evolution $`U_k`$. The more broadly the $`\{p_k\}`$ are distributed, the larger the degree of inhomogeneity in the evolution, and the more incoherent noise enters into the evolution. Estimates of the actual eigenvalues of $`S=_kp_k\overline{U}_kU_k`$ will now be obtained using non-degenerate first-order perturbation theory. Because the RF pulses are not perfect, even in the absence of RF field inhomogeneity, we may assume that the unperturbed eigenvalues are generically non-degenerate. The unitary operator $`U_k`$ generated by the RF field acting at position $`k`$ may be written in exponential form as $`U_k=e^{iH_kt}`$ where $`H_k`$ represents the effective Hamiltonian of the evolution over the period $`t`$ for which the pulse is applied ($`\mathrm{}`$ has been set equal to $`1`$). Defining $`H_0`$ to be the unperturbed (and desired) Hamiltonian, the eigenvalues $`\varphi _j`$ and eigenstates $`|\varphi _j`$ of $`H_0`$ satisfy the eigenvalue equation

$$U_0|\varphi _j=e^{i\varphi _jt}|\varphi _j,$$
(5)
where $`U_0=\mathrm{exp}(iH_0t)`$. The Hamiltonian of a particular $`U_k`$ is assumed to be a perturbation of the desired, homogeneous Hamiltonian

$$H_k=H_0+K_k,$$
(6)
where $`K_k`$ is the perturbation. To first order, the new eigenvalues of $`H_k`$ are

$$\stackrel{~}{\varphi }_{j,k}=\varphi _j+\varphi _j|K_k|\varphi _j,$$
(7)
and the corresponding eigenvalues of $`U_k`$ are

$$U_k|\stackrel{~}{\varphi }_{j,k}=e^{i\stackrel{~}{\varphi }_{j,k}t}|\stackrel{~}{\varphi }_{j,k}.$$
(8)
To first order, the spectral decomposition of $`S`$ is

$$S=\underset{k,m,j}{}p_k\left(e^{i\stackrel{~}{\varphi }_{m,k}t}\overline{|\varphi _m}\overline{\varphi _m|}e^{i\stackrel{~}{\varphi }_{j,k}t}|\varphi _j\varphi _j|\right)$$
(9)
and the eigenvalues of $`S`$ are given approximately by

$`\lambda _{jm}`$ $`=`$ $`{\displaystyle \underset{k}{}}p_ke^{i(\stackrel{~}{\varphi }_{j,k}\stackrel{~}{\varphi }_{m,k})t}`$
$`=`$ $`e^{i(\varphi _j\varphi _m)t}{\displaystyle \underset{k}{}}p_ke^{i\left(\varphi _j|K_k|\varphi _j\varphi _m|K_k|\varphi _m\right)t}.`$
We now imagine the scenario where $`K_k`$ is given by $`K_kt=\frac{\omega _k\omega _0}{\omega _0}K`$. This result would in fact be exact for one spin on resonance. In this case, $`\frac{\omega _k\omega _0}{\omega _0}`$ is the parameter $`\alpha `$ defined in the introduction (which parameterizes the inhomogeneity) and represents the normalized RF power deviation from the desired power $`\omega _0`$. Defining $`\mathrm{\Delta }\omega =\frac{\omega _k\omega _0}{\omega _0}`$ and $`K_{jm}=\varphi _j|K|\varphi _j\varphi _m|K|\varphi _m`$, the previous equation in the continuous limit becomes

$$\lambda _{jm}=e^{i(\varphi _j\varphi _m)t}p(\mathrm{\Delta }\omega )e^{iK_{jm}\mathrm{\Delta }\omega }𝑑\mathrm{\Delta }\omega .$$
(10)
We see in this case that to first order the eigenvalue $`\lambda _{jm}`$ is just the unperturbed eigenvalue $`e^{i(\varphi _j\varphi _m)t}`$ times the Fourier transform of the RF distribution profile evaluated at $`K_{jm}`$. This result demonstrates that the probability distribution profile $`p(\mathrm{\Delta }\omega )`$ of an incoherent process can be determined, within some degree of approximation, from the eigenvalue structure $`\{\lambda _{jm}\}`$ of an experimental superoperator, given some model for the incoherence $`K`$. Knowing $`K`$ would indeed allow one to build the correspondence between $`\lambda _{jm}`$ and $`K_{jm}`$, and then to determine $`p(\mathrm{\Delta }\omega )`$ by performing an inverse Fourier transform. Of course, this result holds for general $`K`$ only when the perturbation is in the first order regime, and when the unperturbed eigenvalues are non-degenerate. But if $`K`$ approximately commutes with $`H_0`$, then the first-order perturbation in the eigenvalues is close to an exact correction, and the above analysis gives a very accurate description of the incoherent process.

### III.2 Recovery of the Profile

We now demonstrate via a numerical example how one can recover the profile $`p(\mathrm{\Delta }\omega )`$ from the eigenvalue spectrum of a measured superoperator. In the theory derived in the previous subsection, $`\lambda _{jm}`$ is the data, i.e., the eigenvalues from the measured superoperator. A model $`K`$ is needed for the perturbation, while $`\varphi _j`$ and $`|\varphi _j`$ are known through the knowledge of $`H_0`$. Formally solving for $`p(\mathrm{\Delta }\omega )`$ from (10),

$$p(\mathrm{\Delta }\omega )=\frac{1}{2\pi }\lambda _{jm}e^{i(\varphi _j\varphi _m)t}e^{iK_{jm}\mathrm{\Delta }\omega }𝑑K_{jm}.$$
(11)
The different eigenvalues $`\lambda _{jm}`$ multiplied by $`e^{i(\varphi _j\varphi _m)t}`$ therefore allow us to obtain the complex function of $`K_{jm}`$ which we shall call $`f`$, corresponding to the Fourier transform of the distribution profile.

For the numerical demonstration of this technique we take a $`3`$-qubit system. We choose $`H_0`$ and $`K`$ such that $`|\varphi _l|K|\varphi _n/(\varphi _l\varphi _n)|0.1`$ for $`nl`$ and $`[H_0,K]0`$ so that first order perturbation theory can be used Sakurai . We then use a measured RF inhomogeneity profile shown in Fig. 2 to construct the following superoperator acting on Liouville space:

$`S={\displaystyle \underset{\mathrm{\Delta }\omega }{}}p(\mathrm{\Delta }\omega )\overline{U}(\mathrm{\Delta }\omega )U(\mathrm{\Delta }\omega ).`$ (12)
Provided with this superoperator, we compute its eigenvalues and plot them on the Argand diagram (see Fig. 3). The perturbation $`K`$ is such that the first order limit condition is fulfilled, but that the size of its diagonal elements in the $`H_0`$ eigenvectors basis is large enough to generate significant dephasing and attenuation in the eigenvalues. The correspondence between $`\lambda _{jm}`$ and $`K_{jm}`$ is needed to recover the distribution profile (11).

The correspondence can be established by first computing $`\overline{\varphi _j},\varphi _m|S|\overline{\varphi _j},\varphi _m`$, then searching for the eigenvalue of $`S`$ closest to it. This allows us to make the correspondence between one unperturbed eigenvalue with eigenvector $`|\overline{\varphi _j},\varphi _m=|\overline{\varphi _j}|\varphi _m`$ (obtained from the knowledge of $`H_0`$) and one eigenvalue of $`S`$. The function $`f=p(\mathrm{\Delta }\omega )e^{i\mathrm{\Delta }\omega K_{jm}}𝑑\mathrm{\Delta }\omega `$ with respect to $`K_{jm}`$ can then be constructed. The real and imaginary parts of that function are plotted in Fig. 4.

Note that we ignore the degenerate points at $`f=1`$ because they do not provide any information about $`p(\mathrm{\Delta }\omega )`$ other than normalization. The $`64`$ eigenvalues of the superoperator minus the $`8`$ degenerate ones equal to $`1`$ (at $`K_{jm}=0`$), plus $`1`$ eigenvalue added at $`K_{jm}=0`$ to avoid a DC offset in the reciprocal domain, yield a complex function $`f`$ made of $`57`$ unequally spaced points. The function $`f`$ is conjugate symmetric with respect to $`0`$, which is consistent with the fact that we are supposed to recover a probability distribution, i.e. a real function, after computing the inverse Fourier transform. To perform the inverse Fourier transform of a function sampled at unequally spaced points, we used an algorithm prescribed in MathPaper . The result is shown in Fig .5. The width of the probability distribution and its skewness are recovered to a good extent, the discrepancy being due to the lack of information about the function $`f`$. It is worth mentioning that with $`57`$ sample points, the window of $`K_{jm}`$ values should be large enough to allow low frequency components of the profile to be reliably extracted. If the incoherent perturbations were very small, there will not be as many large values of $`K_{jm}`$, and therefore less low frequency information would be available.

However, the perturbations can be made larger without changing the mathematical model, by simply repeating the control sequence several times, provided other noise mechanisms do not play a significant role. In addition, more points could be used to get a better sampling resolution, with, for instance, a $`4`$-qubit superoperator yielding $`241`$ points. One cannot have arbitrarily many points, however, because the correspondence between the eigenvalues of $`S`$ and those of $`S_0`$ could quickly become impossible to establish unless a very good knowledge of the perturbation $`K`$ is available. The density of points in the Argand diagram becomes so large that eigenvalues can easily become confused. Here, the $`3`$-qubit superoperator was enough to recover the essential features of the probability distribution.

If $`H_0`$ is not exactly known, and in fact a constant offset Hamiltonian which is proportional to $`K`$ is present, then a different function $`f`$ is obtained :

$$f=p(\mathrm{\Delta }\omega )e^{iK_{jm}\mathrm{\Delta }\omega }e^{i\beta K_{jm}}𝑑\mathrm{\Delta }\omega ,$$
(13)
where $`\beta `$ is a constant real number. Taking the inverse Fourier transform of $`f`$ would reveal a distribution $`p(\mathrm{\Delta }\omega )`$ centered around $`\beta `$ rather than $`0`$, indicating that $`H_0+\beta K`$ is in fact the unperturbed Hamiltonian. Perfect knowledge about $`H_0`$ as a result is not required provided the offset is approximately proportional to $`K`$.

It is important for this method to work that the model $`K`$ chosen *a priori* is a reasonably faithful one, and that it approximately commutes with $`H_0`$. This ensures that the applied first-order perturbation theory is valid. In the extreme case, in which $`K`$ anticommutes with $`H_0`$, $`\varphi _j|K|\varphi _j=0j`$ and no “data” would be available for analysis. It is also worth mentioning that this method is, needless to say, not scalable. However, in many instances, as in NMR, the inhomogeneity features are apparatus dependent, so that reasonably small physical systems can be used to probe them. The scalability of the method is not necessarily a requirement. The distribution of some control parameters, once obtained, can be valuable in designing robust control sequences PraviaRFI  for larger and more complex systems.

## IV Conclusion

Here we reviewed that when incoherence is present during the preparation of the input states for QPT, the resulting correlations between the system and the ”environment” can play an important role on the subsequent system’s dynamics. The map obtained by right multiplying the matrix of output states by the inversion of the matrix of input states still has a meaning, but a correct interpretation of the measured data (or transformation) requires an analysis of the incoherence effects affecting the tomographic procedure. In particular, the measured map needs not be CP. If quantitative information is missing, our perturbation theory analysis of superoperator eigenvalue spectra can be used to determine an effective distribution of unitaries characterizing the process, provided a good mathematical model is available. While this requires a significant effort to measure a 3-qubit superoperator, it is certainly feasible within present experimental capabilities. Lastly, the knowledge of the distribution of control parameters should finally allow us to design more efficient control sequences aimed at counteracting these deleterious effects.

## V Acknowledgements

This work was supported by ARO, DARPA, NSF and the Cambridge-MIT institute. Correspondence and requests for materials should be addressed to D. G. Cory (e-mail: dcory@mit.edu). We thank Marcos Saraceno for valuable discussions.


format: text

text:
# New constraints on the mass composition of cosmic rays above 10¹⁷ eV from Volcano Ranch measurements

## I Introduction

The measurement of the mass composition of cosmic rays above $`10^{17}`$ eV is a challenging problem. This information is as important as the energy spectrum and the anisotropy in determining cosmic ray origin. One must know the likely mass range of a particular data set before one can interpret anisotropy information confidently, given the influence of galactic and intergalactic magnetic fields. Our knowledge of the mass composition of cosmic rays above $`10^{17}`$ eV remains very limited. Recent re-interpretation of measurements of the lateral distribution of water-Čerenkov signals made at Haverah Park  suggests a composition of 34% protons and 66% iron in the range $`2\times 10^{17}`$ \- $`10^{18}`$ eV. This contrasts with earlier claims, from observations made using Fly’s Eye, that the composition changes from a heavy mix around $`3\times 10^{17}`$ eV to a proton dominated flux around $`10^{19}`$ eV . At Yakutsk, both the inferred values of the depth of shower maximum ($`X_{max}`$) and the muon density favor a composition change from a mixture of heavy and light components to light composition over the same energy region . From HiRes/MIA data , there are claims that there is a rapid change from a heavy to a light composition between 0.1 and 1.0 EeV. A recent analysis by the HiRes collaboration of data collected in the energy range between $`10^{18}`$ \- $`10^{19.4}`$ eV  is consistent with a nearly constant purely protonic composition. The fraction of protons, however, decreases when they interpret their data using sibyll2.1. On the other hand, the AGASA group have argued for a “mixed” unchanging composition from 1 - 10 EeV  (using mocca for the simulations). A recent analysis of the muon component in air showers with aires/qgsjet98 around $`10^{19}`$ eV by the AGASA collaboration indicates a relatively light average composition .

The source of the discrepancy between different experiments is not understood and it is important to resolve the issue, because of its implications for cosmic ray models of origin, acceleration and propagation. Volcano Ranch data may provide a path for further understanding.

Following the successful re-examination of the Haverah Park data  with modern shower models, we report here a similar analysis using the Volcano Ranch data, collected by Linsley  to determine the shape of the lateral distribution of air showers. This is the first attempt to examine the Volcano Ranch data with the results of Monte Carlo calculations, using Monte Carlo tools that were unavailable when the data were recorded in 1970. It is timely as the situation on mass composition above $`10^{17}`$ eV remains confused and the steepness of the lateral distribution is sensitive to the depth of maximum of the shower, and therefore to the primary composition and to the character of the initial hadronic interactions.

To simulate the development of the air showers, we have used the aires  code (version 2.4.0), with the hadronic interaction generator qgsjet98 . The results of the simulated showers were convolved with a simulation of the detector response made using geant4 . A comparison of two hadronic generators (qgsjet98 and sibyll2.1 ) was presented in . Both give satisfactory descriptions of the data, but we have preferred to use qgsjet98 because this model has been shown to be consistent with experimental data at energies up to 10 PeV and beyond

## II The Volcano Ranch Array

The pioneering Volcano Ranch instrument consisted of an array of scintillation counters. The array was operated in three configurations from 1959-1976 at the MIT Volcano Ranch station located near Albuquerque, New Mexico (atmospheric depth 834 g $`\mathrm{cm}^2`$). One of its many distinctions was the detection of the first cosmic ray with an energy estimated at $`10^{20}`$eV . The final configuration, of relevance here, comprised 80 detectors of surface area 0.815 $`\mathrm{m}^2`$, scintillator thickness of 9.032 g $`\mathrm{cm}^2`$ laid out on a hexagonal grid with a separation of 147 m (Fig1). This configuration allowed precise measurement of the lateral distribution of the detector signals. The steepness of the lateral distribution, and its fluctuations, can be used to explore the primary mass composition as in . Fortunately, in his various writings, Linsley has left unusually detailed descriptions of his equipment, together with examples of events and a description of his data reduction methods.

## III Lateral distribution function

A generalized version of the Nishimura-Kamata-Greisen (NKG) formula was used to describe the lateral distribution of particles at ground in minimum ionizing particles per square metre ($`\mathrm{mips}\mathrm{m}^2`$) for Volcano Ranch data . This lateral distribution function is given as

$$S_{\mathrm{VR}}(r)=\frac{N_{fit}}{r_m^2}C(\alpha ,\eta )\left(\frac{r}{r_m}\right)^\alpha \left(1+\frac{r}{r_m}\right)^{(\eta \alpha )}$$
(1)
normalized to shower size $`N_{fit}`$ with

$$C=\frac{\mathrm{\Gamma }(\eta \alpha )}{2\pi \mathrm{\Gamma }(2\alpha )\mathrm{\Gamma }(\eta 2)}.$$
(2)
Here $`r_m`$ is the Molière radius, which is $``$ 100 m for the Volcano Ranch elevation. $`\eta `$ and $`\alpha `$ are parameters that describe the logarithmic slope of this function.

From a subset of 366 showers detected with the array, the form of $`\eta `$ as a function of zenith angle $`\theta `$ and shower size $`N_{fit}`$ was found to be :

$$\eta (\theta ,N_{fit})=a+b(\mathrm{sec}\theta 1)+clog_{10}(\frac{N_{fit}}{10^8})$$
(3)
with $`a=3.88\pm 0.054`$, $`b=0.64\pm 0.07`$, and $`c=0.07\pm 0.03`$ where a fixed value of $`\alpha =1`$ was adopted.

## IV Simulation of the detector response of the Volcano Ranch array

The aires code provides a realistic air shower simulation system, which includes electromagnetic algorithms  and links to different hadronic interactions models. As mentioned above, we have used the qgsjet98 model for nuclear fragmentation and inelastic collisions. For the highest energy showers, the number of secondaries becomes so large that it is prohibitive in computing time and disk space to follow and store all of them. Hillas  introduced a non-uniform statistical sampling mechanism which allows reconstruction of the whole extensive air shower from a small but representative fraction of secondaries that are fully tracked. Statistical weights are assigned to the sampled particles to account for the energy of the discarded particles. This technique is known as “statistical thinning”. The aires code includes an extended thinning algorithm, which has been explained in detail . The present work has been carried out using, in most cases, an effective thinning level $`ϵ_{th}=E_{th}/E_{prim}=10^7`$ which is sufficient to avoid the generation of spurious fluctuations and to provide a statistically reliable sample of particles far from the shower core. All shower particles with energies above the following thresholds were tracked: 90 keV for photons, 90 keV for electrons and positrons, 10 MeV for muons, 60 MeV for mesons and 120 MeV for nucleons and nuclei.

We have generated a total of 1735 proton and iron showers with zenith angles in the range $`\mathrm{sec}\theta `$ = 1.0 - 1.5 and primary energies between $`10^{17}`$ eV and $`10^{19}`$ eV, to match the Volcano Ranch data. To simulate the response of the detectors of the array to the ground particles, we utilized the general-purpose simulation toolkit geant4. Our procedure follows the prescription in , where the detector response to electrons, gamma, and muons is simulated in the energy range 0.1 to $`10^5`$ MeV and for five bins per decade of energy. The results of air shower simulations are convolved with the detector response to obtain the scintillator yield expressed in $`\mathrm{mips}\mathrm{m}^2`$. The computed lateral distributions of particles and the corresponding signal from the scintillators for photons, electrons and muons are displayed in Figure 2.

A comparison between the lateral distribution measurements and proton/iron showers simulated with aires/qgsjet98 including the scintillator response was presented previously . Each simulated shower was thrown a maximum of 100 times on to the simulated Volcano Ranch array with random core positions in the range 0 - 500 m from the array center.

With the thinning method used in Monte Carlo shower generators, when particles reach the thinning energy just one of them is followed and multiplied by the corresponding weight at the end. Thus, to simulate the response of the detectors correctly, it is necessary to perform a smoothing of the densities of the ground particles around the position of each detector. All particles in a “sampling zone” around a given detector are selected and the statistical weight, as obtained from aires, is multiplied by the “sampling ratio” $`A_{detector}/A_{sampling}`$ where $`A_{sampling}`$ is the area of the “sampling zone” and $`A_{detector}`$ is the corresponding detector area. This is equivalent to sampling particles on a larger area to get a realistic density around the detector position. As the densities depend mainly on the distance to the shower axis, the sampling area over which simulated particles are gathered is such that this ratio varies from about 0.1 at 100 m to about 0.001 at 1 km. As a first check on the validity of our approach, the data of a single large event were compared with calculations . Further checks between data and Monte Carlo were performed as the one in Figure 3. In this plot we present a comparison between lateral distribution measurements  and a $`10^{19.1}`$ eV proton (left) and iron (right) shower simulated with aires/qgsjet98, including the scintillator response of the detectors in the Volcano Ranch array configuration. It has been shown  that the fluctuation of the density of shower particles far from the core is quite small and that the density at 600 m, S(600m), depends only on primary energy. We normalize the showers to the value of S(600m) in order to decouple the normalization factor from the parameters related to the shape of the lateral distribution which change with primary mass but only slightly with shower size. The agreement between data and Monte Carlo is good and gives confidence in the procedures used.

## V Derivation of the primary mass composition

The nature of the primaries that initiate air showers is difficult to establish from the average properties of the data. For example, an average property can be explained with a mass composition of a single species (A) or by an appropriate mixture of species. However, with the Volcano Ranch array, accurate measurements of $`\eta `$ were made on a shower-by-shower basis for fixed bins of zenith angle separated by 80 g cm<sup>-2</sup> . Thus the fluctuations of $`\eta `$ can be used to break this degeneracy. Linsley determined the precision of each measurement of $`\eta `$ and reported the average value of this quantity for each zenith angle bin.

The average error in $`\eta `$ from the fit made to the simulated lateral distributions ($`\sigma _{sim}=0.029`$), is smaller than the one reported by Linsley($`\sigma =0.072`$). To include within the simulation the effect of data reconstruction, we smeared each value of $`\eta `$ calculated by Monte Carlo using a Gaussian with a width chosen so that Linsley’s overall uncertainties in $`\eta `$ are reproduced i.e. $`\sigma _{smear}=0.0662`$ is found by quadrative subtraction of the average values of $`\sigma _{sim}`$ and $`\sigma `$.

This is a minor correction since the measurement accuracy is so much smaller than the intrinsic shower-to-shower spread (r.m.s.= 0.19). Thus, for each value of $`\eta `$ found from the Monte Carlo calculation, we have a corresponding and realistic estimate of its “experimental” uncertainty. We are thus able to make comparisons of Volcano Ranch data with our calculations.

As a further check, we have calculated the variation of $`\eta `$ with shower size and zenith angle with Monte Carlo and made comparisons with the Volcano Ranch data. The number of particles at ground level ($`N_{fit}`$) is obtained from a fit to the lateral distribution function (with $`\alpha =1`$) for fixed bins of zenith angle. The variation of $`\eta `$ with $`N_{fit}`$ from the calculation has been compared with the average functional form of $`\eta `$ given by Equation 3. The results of this comparison for $`\mathrm{sec}\theta =`$ 1.0 - 1.1 and 1.3 - 1.4 can be seen for proton and iron showers in Figure 4. The error bars indicate the r.m.s. spread of data which is very much greater than the r.m.s. spread of the mean. The shaded band represents the fit to Volcano Ranch data, including the errors given in Equation(3) for $`a`$ and $`b`$.

The variation of $`\eta `$ with zenith angle is shown in Figure 5 for events with shower size in the range $`\mathrm{log}N_{fit}=7.68.6`$ (left) and $`\mathrm{log}N_{fit}=8.69.6`$ (right). One can see that the average form of $`\eta `$ over a realistic range of mass composition, from proton to iron, is well represented by the simulations. The error bars represent the r.m.s. spread as before.

### A Fitting VR data mass composition using finite Monte Carlo samples of different primaries.

We can estimate the primary mass that describes the Volcano Ranch data, assuming a bi-modal composition, using a maximum likelihood fit for the best linear combination of pure iron and pure proton samples to match the data sample. The available data are in bins of $`\eta `$ : the number of data points in several bins is small, so a $`\chi ^2`$ minimization is inappropriate. A maximum likelihood technique assuming Poisson statistics was adopted. The probability of observing a particular number of events $`d_i`$ in a particular bin is given by $`\mathrm{exp}^{f_i}f_i^{d_i}/d_i!`$ where $`f_i`$ is the predicted value for the number of events in this particular bin. If we assume a bi-modal composition of proton and iron with fractions $`P_{Fe}`$ and $`P_p`$ then $`f_i=C(P_{Fe}+P_p)`$ where $`C`$ is the overall normalization factor between numbers of data and Monte Carlo events. Estimates of the fractions $`P_j`$ are found by maximizing $`\mathrm{ln}(L)=\mathrm{\Sigma }d_i\mathrm{ln}(f_i)\mathrm{ln}(d_i!)f_i`$. Our Monte Carlo samples are at least ten times larger than the data sample to avoid effects of finite Monte Carlo data size.

In Figure 6 we compare the Monte Carlo results with the Volcano Ranch data points for near vertical showers. As can be seen the tail at large $`\eta `$ in the comparison with iron indicates that a lighter component is required to fit the experimental data. The best fit gives a mixture with (89 $`\pm `$ 5) % of iron, with a corresponding percentage of protons, and this distribution of $`\eta `$ is shown in Figure 7. One detail that Linsley did not describe is the distribution of shower sizes that comprise the data set. What is recorded is that the median energy was $`10^{18}`$ eV and that the shower sizes are between $`4\times 10^7`$ and $`6\times 10^9`$. This corresponds to an energy range approximately between $`10^{17}`$ eV and $`10^{19}`$ eV. To reproduce the data set, a differential energy spectrum of slope -3.5 was chosen for the primary energy spectrum. Simulations show that the whole area of the array was active above $`10^{18}`$ eV. Using this spectrum and the condition that the median energy be $`10^{18}`$ eV, the calculated threshold energy is approximately $`10^{17.7}`$ eV. Below $`10^{18}`$ eV the energy distribution of the data set was smoothed (always under the condition of median energy to be $`10^{18}`$ eV) assuming that the effective area of the detector increases with energy from the threshold up to $`10^{18}`$ eV. The conclusions about mass composition that we presently draw, are constrained by uncertainties in the details of the energy distribution of events recorded at VR. Additionally they are constrained by the hadronic model assumed.

The systematic error arising from our lack of knowledge of the energy distribution of the events has been estimated by repeating the fitting procedure with different energy spectra. From this analysis we estimate a systematic error of 12%. An additional source of systematic error is related to uncertainties in the hadronic interaction model: following the discussion in  which relates to the use of qgsjet98 rather than qgsjet01, the systematic shift in the fraction of iron is 14%. Showers that are calculated using qgsjet01 are found to develop higher in the atmosphere so that the fraction of Fe estimated is reduced from (89 $`\pm `$ 5)% to (75 $`\pm `$ 5)%.

## VI Comparison with other data

In Figure 8 we present the fraction of iron reported by Fly’s Eye, Agasa A100 and A1 (using sibyll as the hadronic model)  and the results for Haverah Park  using qgsjet98 and qgsjet01. Also shown is the mean composition with the corresponding error we get for the Volcano Ranch energy range using qgsjet98 and an estimation of what we would get using qgsjet01 as it was done in Haverah Park .

The inconsistency that exists between several experiments which span different energy ranges, use different techniques is doubtless enhanced as different hadronic models are used for the interpretation of the raw data. The application of a consistent hadronic model  brings the results of AGASA into better agreement with the Fly’s Eye conclusion, while their original analysis stated that there is no indication of changing composition. In Figure 8 the iron fraction for Fly’s Eye and AGASA corresponds to the re-interpretation under the sibyll hadronic model, including triggering efficiency effects; the original analysis from AGASA is not presented here.

While Haverah Park, Volcano Ranch and Akeno-AGASA infer $`X_{max}`$, and hence the overall composition, from properties of secondary particle distribution at ground, Fly’s Eye and HiRes experiments observe a image of the longitudinal shower profile and derive $`X_{max}`$ directly. Nonetheless the estimates can be biased due to the poor knowledge of atmospheric properties as recent studies of atmospheric profiles have suggested .

## VII Conclusions

Measurements of the steepness of the lateral distribution $`\eta `$ were made at Volcano Ranch on a shower-to-shower basis for fixed bins of zenith angle. We have compared the measured distribution of $`\eta `$ to our Monte Carlo results for proton and iron primaries using qgsjet98 including the scintillator response of the detectors in the Volcano Ranch array. Our ability to reproduce Volcano Ranch lateral distribution measurements give us confidence that our analysis procedure is correct.

The cosmic ray mass composition, deduced from Volcano Ranch data, is compatible with mean fraction (89 $`\pm `$ 5(stat) $`\pm `$ 12(sys)) % of iron in a bi-modal proton and iron mix, in the whole energy range $`10^{17.7}`$ eV to $`10^{19}`$ eV, mean energy $`10^{18}`$ eV. Following the discussion in , we estimate that this fraction would be reduced to 75 %, with the same qgsjet01 model adopted.

As mentioned above, some experiments find indications that the composition is becoming lighter with energy. However, in a recent work using Haverah Park data above $`10^{19}`$ eV, it was claimed that the observed time structure of the shower front can be best understood if iron primaries are dominant here . The differences between measurements on mass composition needs to be addressed further if more solid conclusions on the origin, acceleration or propagation of cosmic rays are to be reached. The rate of change with energy and the average mass are still under debate.

###### Acknowledgements.

We warmly thank the late John Linsley for the discussions that we had about his data during the 27th ICRC in Hamburg, and for his permission to proceed with an analysis along the lines described above. However, the interpretation given is that of the authors alone, as it was not possible to discuss the outcome of the work with Linsley before his death. Work at Universidad de La Plata is supported by ANPCyT (Argentina), at Northeastern University by the U.S. National Science Foundation and at the University of Leeds by PPARC, UK (PPA/Y/S/1999/00276). TPM gratefully acknowledges the support of the American Astronomical Society in the form of a travel grant to attend the 28th ICRC, Japan. MTD thanks the John Simon Guggenheim Foundation for a fellowship.


text:
# ”Qualms” from Lavenda, cond-mat/0311270

*

Abstract: The smoke (German: Qualm) produced by Lavenda obscures the simple and fundamental physics of non-extensive systems as proposed in .
The paper  addresses phase-transitions in systems with a fixed, finite, volume, here the Potts-($`q=3`$)-lattice gas on a $`50\times 50`$ lattice, details in the cited papers . Therefore, the volume is irrelevant for the transition. We have a system with the two control-parameters $`E,N`$. For this case Lavenda makes several wrong statements:

1. ”The author discards the fundamental properties of concavity”. Yes, this is the essence of my new extended formulation of thermo-statistics: At phase-separation the entropy \[$`S(E)`$ not $`S(E)/N`$(!)\] of any system is convex (bimodal). Therefore, in conventional thermo-statistics one gets in the thermodynamic limit Yang-Lee singularities. These are certainly not helpful to give a detailed insight into this most interesting physics: The original objects of Thermodynamics, steam engines, work exactly at phase-separation. Sailing on the surface of the sea would, presumably be impossible if there were no phase-separation. In hot nuclei this is the region where all the interesting fragmentation phenomena occur, from nucleon evaporation, to fission, to multifragmentation, and finally to vaporization.
2. ”The curvature eigenvalues $`\lambda _1\lambda _2`$ are devoid of any physical meaning.” On the contrary: for a totally concave entropy surface $`S(E,N)`$ the condition $`\lambda _1<0`$ guarantees that Lavenda’s equation (4) for any $`u,v`$ at fixed volume $`V`$ is fulfilled:

$$S_{EE}\times u^2+2S_{EN}\times uv+S_{NN}\times v^20$$
(1)
This can also be seen by starting from Lavenda’s equation for the eigen-curvatures :

$$\lambda _{1,2}=\frac{1}{2}(S_{NN}+S_{EE})\pm \frac{1}{2}\sqrt{(S_{NN}S_{EE})^2+4S_{EN}^2}$$
(2)
If both eigenvalues are $`<0`$ then evidently both $`S_{NN}`$ and $`S_{EE}`$ must be negative. Therefore, the single condition $`\lambda _1<0`$ suffices to guarantee relation (1), and of course also Lavenda’s conditions (2+3).
3. ”The inverse temperature $`\beta `$ has no meaning in the microcanonical ensemble”. Yes, of course, however, what was discussed in my paper  is the weight of a given energy in the Laplace transform $`ET`$:

$$Z(T,\mu ,V)=_0^{\mathrm{}}𝑑Ee^{S(E)\beta E},$$
(3)
with the external parameter $`\beta =1/T`$. The weight $`e^{S(E)\beta E}`$ has a single maximum in mono-phases, and at phase-separation it is bimodal, where it has at least two sharp maxima separated by the latent heat with a minimum in between (convexity).
4. ”Extensivity of the entropy is already incorporated in Boltzmann–Planck’s principle $`S\mathbf{=}k\mathrm{𝐥𝐧}W`$. First, extensivity and additivity are two different things. For two independent and non-interacting systems the microcanonical Boltzmann-Planck entropy is of course additive, $`S=S_1+S_2`$. Extensivity means something else: $`S(N)`$ scales as $`S(N)N`$. This is usually demanded in the thermodynamic limit but is, of course, invalid in non-extensive situations.


text:
# Power law for ensembles of stock prices

## 1 Introduction

Recent empirical findings regarding financial time series have suggested that financial markets can be considered complex systems in which a large number of agents participate in trade and interplay with each other . A goal of econophysics is to model the stochastic processes of financial markets as complex systems and to reproduce the statistical properties observed in the time evolution of financial prices by using tools of theoretical physics such as the master equation and the Fokker-Plank equation. Although most empirical studies so far have investigated the statistical properties of the time series of a single stock or index from theoretical point of view, it is more important to examine the statistical properties of ensembles of stocks. Lillo and Mantegna \[2-4\] studied the ensemble return distribution, and showed some statistical properties of the ensembles of stocks returns during the market’s crash and rally days. In this paper we focus our attention on the statistical properties of ensembles of stock prices. Using $`1200`$ stocks traded in the Tokyo Stock Exchange, we formed ensembles of daily stock prices.

## 2 Power-law for ensembles of stock prices

Figure 1(a) shows a so-called Zipf-plot for an ensemble of stock prices of 1200 companies on May 1, 1989, on a double logarithmic scale, where $`S(n)`$ denotes the stock prices in descending order ($`S(1)`$ being the stock with the highest price, $`S(2)`$ the stock with the second highest price, and so on). The Zipf-plot for stock prices higher than 800 yen shows a straight line and is well described by $`S(n)n^\beta `$ with $`\beta =0.46`$. Figure 1(b) shows the complementary cumulative distribution of the stock price ensemble on May 1, 1989. The tail of the distribution can be well described by a power-law decay, $`P(S>x)x^\alpha `$ with $`\alpha =2.17`$. The power-law exponent obtained by power-law regression fits the complementary cumulative distribution function, where the fit is for all $`x`$ larger than $`1000`$111Note that the Zipf plot and cumulative plot are equivalent. The relationship, $`\beta =1/\alpha `$ holds theoretically. Thus the power-law exponent calculated from $`\beta `$ is equal to $`2.175`$. This value is very close to the power-law exponent estimated from the stock price ensemble’s complementary cumulative distribution.. To confirm the robustness of the above analysis, we repeated this analysis for each of the trading days in a 5-year period from January 4, 1988, to December 30, 1992, which corresponds to the periods of bubbles and crashes in the Japan stock market. They are all consistent with a power-law asymptotic behavior, and the power-law exponent $`\alpha `$ moved in the range of $`1.7<\alpha <2.2`$. Figure 1(c) and Figure 1(d) show the movement of the power-law exponent $`\alpha `$, and the movement of the mean value of the stock prices in the period from January 4, 1989, to December 30, 1992. If follows from the figures that the power-law exponent $`\alpha `$ is positively correlated with the mean value of stock prices $`\overline{S}`$. The correlation coefficient is equal to $`0.88`$.

## 3 Concluding remarks

These findings suggest that the financial markets are itinerating over a series of the non-equilibrium stationary states characterized by a power-law distribution. We found here that the power-law exponent $`\alpha `$ is positively correlated with the stock price’s mean value. This finding suggests that the stock market bubbles in the period from 1988 to 1990 are defined as the generally extraordinary rise of stock prices. In our recent work on the Japanese land market , we quantitatively investigated the statistical properties of an ensemble of land prices in Japan in the period from 1981 to 2002, corresponding to a period of bubbles and crashes. We found that the tail of the complementary cumulative distribution function of the ensemble in the high price range is well described by a power-law distribution, $`P(S>x)x^\alpha `$, and the power-law exponent $`\alpha `$ is negatively correlated with the mean value of land prices. Therefore, the real estate bubble is not defined as the generally extraordinary rise of land prices, but rather as an extraordinary enlargement of the inequality of land prices. Why did we obtain the empirically contradictory results regardless of the bubbles that occurred in the same period in Japan? The answer could be as follows. At the peak of the Japanese stock market in December 1989, Japanese stocks had a total market value of about 4 trillion dollars, almost 1.5 times the value of all U.S. equities and close to 45 percent of the world’s equity market capitalization . Thus, the Japanese stock market bubbles in the late of 1980s may be defined as the extraordinary rise of market capitalization of the Japanese stock market among the world’s stock markets. This conjecture remains to be tested.

## 4 Acknowledgement

Financial support by the Japan Society for the Promotion of Science under the Grant-in-Aid, No. 06632 is gratefully acknowledged. All remaining errors, of course, are ours.


text:
# Linear Perturbations in Brane Gas Cosmology

## I Introduction

Understanding the behavior of strings in a time dependent background has been a subject of much interest and has been pursued in a number of differing ways. One scenario, known as Brane Gas Cosmology (BGC) is devoted to understanding the effect that string and brane gases could have on a dilaton-gravity background in the early Universe bv ; vafa ; bgc ; isotropization ; stable ; extended . In bv , it was suggested that the energy associated with the winding of strings around the compact dimensions would produce a confining potential for the scale factor and halt the cosmological expansion111This was later shown quantitatively in vafa .. The analysis of BGC was initially performed at the level of the equations for homogeneous and isotropic cosmology. The results were recently extended to the case of anisotropic cosmology in isotropization . There, it was shown that string gases can result in three dimensions growing large and isotropic due to string annihilation while the other six dimensions remain confined. In stable  it was shown that by considering both momentum and winding modes of strings, the six confined dimensions can be stabilized at the self-dual radius, where the energy of the string gas is minimal. This result demonstrated that, in BGC, the volume moduli of the extra dimensions can be stabilized in a natural and intuitive way.

To date, however, all analyses have been performed at the level of the equations for homogeneous cosmology. The string sources are usually represented by a perfect fluid with homogeneous energy and pressure densities given by the mass spectrum of the strings (see e.g. bgc ; stable ). One may worry that inhomogeneities of string sources (in particular strings winding around the confined dimensions) as a function of the unconfined spatial directions could lead to serious instabilities which would ruin one of the main successes of BGC, namely the prediction that three directions become large, leaving the other six confined uniformly as a function of the coordinates of the large spatial sections.

As a first step towards addressing this concern, we here study the evolution of linear fluctuations about a BGC background. We introduce linear matter and metric fluctuations, derive the linear perturbation equations, and study their stability on long wavelength (wavelengths larger than the Hubble radius) 222On shorter wavelengths, we expect that the motion of the strings will smear out potential instabilities in a way analogous to how the motion of light particles (“free-streaming”) leads to a decay of short wavelength fluctuations in standard cosmology (see e.g. Peebles  for a review). Note that since the dilaton plays a vital role in the stabilization process, we must also consider the effect of inhomogeneities in the dilaton field. Our main result is that the long wavelength fluctuations do not grow in time - there are no instabilities at the level of a linearized analysis.

In Section 2, we present the equations of BGC and review the background solution that leads to $`3+1`$ dimensions growing large and $`6`$ dimensions being stabilized at the self-dual radius. In Section 3, we introduce the fluctuations and derive the equations for linearized perturbations about the time dependent homogeneous background of BGC. The long wavelength solutions are given in Section 4, leaving a detailed description for the Appendix. We conclude with a discussion of the results in Section 5.

## II Background Solution

Our starting point is the low energy effective action for the bulk space-time with string matter sources vafa ,

$$S=\frac{1}{4\pi \alpha ^{}}d^Dx\sqrt{g}e^{2\phi }\left(R+4(\phi )^2\frac{1}{12}H^2\right)+S_\text{m},$$
(1)
where $`R`$ denotes the Ricci scalar, $`g`$ is the determinant of the background metric, $`\phi `$ is the dilaton field, and $`H`$ is the field strength of an antisymmetric tensor field. The action of the matter sources is denoted by $`S_\text{m}`$. For example, with $`D=10`$ this is the low energy effective action of type II-A superstring theory. For the purposes of this paper we will ignore the effects of branes, since it will be the winding and momentum modes of the string that ultimately determine the dimensionality and stability of space-time bgc . Here, we will ignore the effects of fluxes 333See Campos:2003ip  for inclusion of fluxes in the scenario., i.e. we set $`H=0`$.

This action yields the following equations of motion,

$`R_\mu ^\nu +2_\mu ^\nu \phi =8\pi M_p^2e^{2\phi }T_\mu ^\nu ,`$
$`R+4_\kappa ^\kappa \phi 4_\kappa \phi ^\kappa \phi =0,`$ (2)
where $``$ is the covariant derivative.

We will work in the conformal frame with a homogeneous metric of the form

$$ds^2=e^{2\lambda (\eta )}\left(d\eta ^2\delta _{ij}dx^idx^j\right)e^{2\nu (\eta )}\delta _{mn}dx^mdx^n,$$
(3)
where $`(\eta ,x^i)`$ are the coordinates of $`3+1`$ space-time and $`x^m`$ are the coordinates of the other six dimensions, both of which can be taken to be isotropic isotropization . The scale factors $`a(\eta )`$ and $`b(\eta )`$ are given by $`\lambda \mathrm{ln}(a)`$ and $`\nu \mathrm{ln}(b)`$.

We consider the effect of the strings on the background through their stress energy tensor

$$T_\mu ^\nu diag(\rho ,p_i,p_m),$$
(4)
where $`\rho `$ is the energy density of the strings, $`p_i`$ ($`i=1\mathrm{}3`$) is the pressure in the expanding dimensions and $`p_m`$ ($`m=4\mathrm{}9`$) is the pressure in the small dimensions (because of our assuption of isotropy of each subspace, there is only one independent $`p_i`$ and one independent $`p_m`$).

Strings contain winding modes, momentum modes and oscillatory modes. However, since the energies of the oscillatory modes are independent of the size of the dimensions, and since the winding modes and momentum modes dominate the thermodynamic partition function at very small and very large radii of the spatial dimensions, we shall here neglect the oscillatory modes. In the absence of string interactions, the contributions to the stress tensor coming from the string winding modes and momentum modes ($`T_{\mu \nu }^\text{w}`$ and $`T_{\mu \nu }^\text{m}`$ respectively) are separately conserved,

$`T_{\mu \nu }=T_{\mu \nu }^\text{w}+T_{\mu \nu }^\text{m}`$
$`^\mu T_{\mu \nu }^\text{w}=0^\mu T_{\mu \nu }^\text{m}=0.`$ (5)
The conservation equations take the form

$$\rho ^{\text{w},m}+\underset{i=1}{\overset{9}{}}\lambda _i^{}(\rho ^{\text{w},m}p_i^{\text{w},m})=0,$$
(6)
where the derivatives are with respect to the conformal time $`\eta `$, and where for the moment we consider 9 independent scale factors.

Expressing (II) in terms of the metric (3) and the stress tensor (4), we find the following system of equations,

$`3\lambda ^{\prime \prime }6\nu ^{\prime \prime }+6\lambda ^{}\nu ^{}6\nu ^2+\phi ^{\prime \prime }\lambda ^{}\phi ^{}=8\pi M_p^2e^{\phi +2\lambda }\rho `$ (7)
$`\lambda ^{\prime \prime }+2\lambda ^2+6\lambda ^{}\nu ^{}+\lambda ^{}\phi ^{}=8\pi M_p^2e^{\phi +2\lambda }p_i`$ (8)
$`\nu ^{\prime \prime }+6\nu ^2+2\lambda ^{}\nu ^{}+\phi ^{}\nu ^{}=8\pi M_p^2e^{\phi +2\lambda }p_m`$ (9)
$`6\lambda ^{\prime \prime }12\nu ^{\prime \prime }24\lambda ^{}\nu ^{}42\nu ^26\lambda ^2\phi ^2+2\phi ^{\prime \prime }+8\lambda ^{}\phi ^{}+12\phi ^{}\nu ^{}=0.`$ (10)
The explicit forms of the energy density and pressure were given in stable 444The equations here are related to Eq. (18) in stable  by the volume factor $`V=e^{3\lambda +6\nu }`$, e.g. $`\rho =\frac{E}{V}`$,

$`\rho =3\mu N^{(3)}e^{2\lambda 6\nu }+3\mu M^{(3)}e^{4\lambda 6\nu }+6\mu N^{(6)}e^{3\lambda 5\nu }+6\mu M^{(6)}e^{3\lambda 7\nu },`$ (11)
$`p_i=\mu N^{(3)}e^{2\lambda 6\nu }+\mu M^{(3)}e^{4\lambda 6\nu },`$ (12)
$`p_m=\mu N^{(6)}e^{3\lambda 5\nu }+\mu M^{(6)}e^{3\lambda 7\nu },`$ (13)
where $`\mu `$ is a constant, $`N^{(3)}`$ and $`M^{(3)}`$ are the numbers of winding and momentum modes in the large directions, and $`N^{(6)}`$ and $`M^{(6)}`$ in the six small directions.

We are interested in solutions that stabilize the internal dimensions, while allowing the three large dimensions to expand. Such solutions were discussed in stable , where it was shown that the winding and momentum modes of the strings lead naturally to stable compactifications of the internal dimensions at the self dual radius, while the other three dimensions grow large since the string modes winding them have annihilated, allowing for expansion. Thus, we will set $`N^{(3)}=0`$. At the self dual radius, the number of winding modes is equal to the number of momentum modes (i.e. $`N^{(6)}=M^{(6)}`$) and the pressure vanishes ($`p_m=0`$).

In Ref. stable , the solutions subject to the above conditions on the winding and momentum numbers were found numerically. In this paper, we wish to study the stability of these solutions towards linear perturbations in the time interval when the internal dimensions have stabilized and the large dimensions give power law expansion. In the following section, we will derive the equations for the linear fluctuations. The coefficients in these equations depend on the background solution. We will use analytical expressions which approximate the numerically obtained solutions of stable . We restrict our initial conditions so that the evolution preserves the low energy and small string coupling assumptions ($`g_se^{2\phi }1`$).

A approximation to a typical solution of the equations (7-10) is found to be of the form

$`\lambda (\eta )=k_1\mathrm{ln}(\eta )+\lambda _0\text{o}ra(\eta )=a_0\eta ^{k_1},`$
$`\phi (\eta )=k_2\mathrm{ln}(\eta )+\phi _0,`$ (14)
where the constants $`k_1`$, $`k_2`$, $`\lambda _0`$ and $`\phi _0`$ depend on the choice of initial conditions. We have made use of $`\nu =\nu ^{}=\nu ^{\prime \prime }=0`$, $`N^{(3)}=0`$, $`N^{(6)}=M^{(6)}`$, $`p_m=0`$. Note that in this limit (9) is trivially satisfied. An example of a solution yielding stabilized dimensions and three dimensions growing large corresponds to $`k_1=\frac{1}{9}`$ and $`k_2=\frac{9}{7}`$. The numerical solution of stable  and the analytical approximation used in this paper are compared in Fig. 1, for the above values of the constants $`k_1`$ and $`k_2`$.

## III Scalar Metric Perturbations

In this section we consider the growth of scalar metric perturbations (see e.g. MFB  for a comprehensive review of the theory of cosmological perturbations) due to the presence of string inhomogeneities. We are interested in the case where the fluctuations depend only on the external coordinates and conformal time, not on the coordinates of the internal dimensions. For simplicity we work in the generalized longitudinal gauge in which the metric perturbations are only in the diagonal metric elements 555As discussed e.g. in Dorca , for scalar perturbations depending on all spatial coordinates it would be inconsistent to choose the perturbed metric completely diagonal, and one would have to add a metric coefficient to the $`dtdx^m`$ terms, where $`x^m`$ are the coordinates of the internal dimensions. However, as discussed in GG96 , if the fluctuations are independent of the coordinates $`x^m`$, as in our case, the coefficient can be chosen to vanish, and thus the perturbed metric is completely diagonal.. Thus, the metric including linear fluctuations is given by

$$ds^2=e^{2\lambda (\eta )}\left((1+2\varphi )d\eta ^2(12\psi )\delta _{ij}dx^idx^j\right)e^{2\nu (\eta )}(12\xi )\delta _{mn}dx^mdx^n.$$
(15)
The dilaton $`\phi `$ also fluctuates about its background value $`\phi _0`$. The dilaton fluctuation $`\chi `$ is determined by

$$\phi =\phi _0+\delta \phi \chi \delta \phi .$$
(16)
In the above, the fluctuating fields $`\chi ,\varphi ,\psi `$ and $`\xi `$ are functions of the external coordinates $`x^i`$ and time, i.e.

$$\chi =\chi (\eta ,x^i),\varphi =\varphi (\eta ,x^i),\psi =\psi (\eta ,x^i),\xi =\xi (\eta ,x^i).$$
(17)
The perturbations of the matter energy momentum tensor result from over-densities and under-densities in the number of strings. From (11)-(13) and noting that we are interested in the case when $`N^{(3)}=0`$ and $`M^{(6)}=N^{(6)}`$ we find,

$`\delta \rho =\delta \rho _\text{w}+\delta \rho _\text{m},`$ (18)
$`\delta \rho _\text{w}=30\mu N\xi e^{3\lambda }+18\mu N\psi e^{3\lambda }+6\delta N^{(6)}e^{3\lambda },`$ (19)
$`\delta \rho _\text{m}=42\mu N\xi e^{3\lambda }+18\mu M\xi e^{4\lambda }+18\mu N\psi e^{3\lambda }+12\mu M\psi e^{4\lambda }+6\mu \delta M^{(6)}e^{3\lambda }+`$ (20)
$`+3\mu \delta Me^{4\lambda },`$ (21)
$`\delta p_\lambda =6\mu M\xi e^{4\lambda }+4\mu M\psi e^{4\lambda }+\mu \delta Me^{4\lambda },`$ (22)
$`\delta p_\nu =2\mu N\xi e^{3\lambda }\mu \delta N^{(6)}e^{3\lambda }+\mu \delta M^{(6)}e^{3\lambda },`$ (23)
where we define666Notice that we must be careful to distinguish between the perturbed quantities $`\delta N^{(6)}`$ and $`\delta M^{(6)}`$. $`NN^{(6)}=M^{(6)}`$ and $`MM^{(3)}`$. The fluctuations $`\delta N^{(6)}`$, $`\delta M^{(6)}`$, and $`\delta M`$ are taken as functions of both conformal time and the external space, e.g. $`\delta N^{(6)}=\delta N(\eta ,x^i)`$.

It follows from (6) that the perturbed sources obey modified conservation equations for both the winding and momentum modes,

$$\delta \rho ^{\text{w},m}+\underset{i=1}{\overset{9}{}}\lambda _i^{}(\delta \rho ^{\text{w},m}\delta p_i^{\text{w},m})+\underset{i=1}{\overset{9}{}}\delta \lambda _i^{}(\rho ^{\text{w},m}p_i^{\text{w},m})=0,$$
(24)
where $`\delta \lambda =a^1\delta a=\psi `$ and $`\delta \nu =b^1\delta b=\xi `$ are spatial variations.

We rewrite (II) to take the more familiar form of the Einstein and dilaton equations, namely

$`R_\mu ^\nu {\displaystyle \frac{1}{2}}\delta _\mu ^\nu R=e^{2\phi }T_\mu ^\nu 2g^{\alpha \nu }_\mu _\alpha \phi +2\delta _\mu ^\nu \left(g^{\mu \nu }_\mu _\nu \phi g^{\mu \nu }_\mu \phi _\nu \phi \right),`$
$`g^{\mu \nu }_\mu \phi _\nu \phi {\displaystyle \frac{1}{2}}g^{\mu \nu }_\mu _\nu \phi ={\displaystyle \frac{1}{4}}e^{2\phi }T_\mu ^\mu ,`$ (25)
where we invoke Planckian units (i.e. $`8\pi M_p^2=1`$). Plugging the perturbed metric (15) and dilaton into these equations, making use of the background equations of motion, and linearizing the equations about the background (i.e. keeping only terms linear in the fluctuations) yields the following set of equations:

$`\stackrel{}{}^2\psi +3\stackrel{}{}^2\xi 9\xi ^{}3\psi ^{}3^2\varphi ={\displaystyle \frac{1}{2}}e^{2\phi +2\lambda }\left(2\chi T_0^0+\delta T_0^0\right)6\varphi \phi ^{}`$
$`3\psi ^{}\phi ^{}6\xi ^{}\phi ^{}\stackrel{}{}^2\chi +3\chi ^{}+2\varphi \phi ^22\chi ^{}\phi ^{},`$ (26)
$$_i\psi ^{}+3_i\xi ^{}+_i\varphi 3_i\xi =\frac{1}{2}e^{2\phi +2\lambda }\delta T_0^i+_i\varphi \phi ^{}_i\chi ^{}+_i\chi ,$$
(27)
$$_i_j\left(\varphi \psi 6\xi 2\chi \right)=0ij,$$
(28)
$`\left(_i^2\stackrel{}{}^2\right)\left(\varphi \psi 6\xi \right)2\psi ^{\prime \prime }6\xi ^{\prime \prime }4\psi ^{}6\xi ^{}2{}_{}{}^{2}\varphi 4^{}\varphi 2\varphi ^{}`$
$`=e^{2\phi +2\lambda }\left(2\chi T_i^i+\delta T_i^i\right)+2_i^2\chi 4\varphi \phi ^{\prime \prime }2\varphi ^{}\phi ^{}4\varphi \phi ^{}4\psi ^{}\varphi ^{}12\xi ^{}\phi ^{}`$
$`+2\chi ^{\prime \prime }2\stackrel{}{}^2\chi +2\chi ^{}+4\varphi \phi ^24\chi ^{}\phi ^{},`$ (29)
$`\stackrel{}{}^2\varphi +5\stackrel{}{}^2\xi 5\xi ^{\prime \prime }+2\stackrel{}{}^2\psi 3\psi ^{\prime \prime }10\xi ^{}3\varphi ^{}9\psi ^{}6^2\varphi 6^{}\varphi `$
$`=e^{2\phi +2\lambda }\left(2\chi T_m^m+\delta T_m^m\right)4\varphi \phi ^{\prime \prime }2\varphi ^{}\phi ^{}8\varphi \phi ^{}6\psi ^{}\phi ^{}10\xi ^{}\phi ^{}+2\chi ^{\prime \prime }`$
$`2\stackrel{}{}^2\chi +4\chi ^{}+4\varphi \phi ^24\chi ^{}\varphi ^{},`$ (30)
$`2\varphi \phi ^2+2\phi ^{}\chi ^{}+\varphi \phi ^{\prime \prime }6\psi \phi ^{}+{\displaystyle \frac{1}{2}}\varphi ^{}\phi ^{}2\varphi \phi ^{}{\displaystyle \frac{3}{2}}\psi ^{}\phi ^{}3\xi ^{}\phi ^{}{\displaystyle \frac{1}{2}}\chi ^{\prime \prime }`$
$`+{\displaystyle \frac{1}{2}}\stackrel{}{}^2\chi \chi ^{}={\displaystyle \frac{1}{4}}e^{2\phi +2\lambda }\left(2\chi T+\delta T\right),`$ (31)
where $`TT_\mu ^\mu `$ is the trace. The modified conservation equations (24) take the form

$`{\displaystyle \frac{d}{d\eta }}\left(\delta N^{(6)}\right)=7N\xi ^{},`$ (32)
$`42\mu N\xi ^{}72\mu M\xi \lambda ^{}e^\lambda 48\mu M\psi \lambda ^{}e^\lambda +12\mu M\psi ^{}e^\lambda +6\mu {\displaystyle \frac{d}{d\eta }}\left(\delta M^{(6)}\right)`$
$`12\mu \delta M\lambda ^{}e^\lambda +3\mu {\displaystyle \frac{d}{d\eta }}\left(\delta M\right)e^\lambda =0.`$ (33)
These equations give us the evolution of the metric perturbations $`\varphi `$, $`\psi `$, and $`\xi `$ in terms of the matter perturbations $`\chi `$, $`\delta \rho `$, and $`\delta p_i`$. At first glance, it may appear that the above system is over-determined since we have eight equations for seven unknowns. However, as is the case in standard cosmology, the conservation equations are not independent of the Einstein equations. Thus, we can choose to keep only one of the modified conservation equations and our system will be consistent.

## IV Evolution of Fluctuations

In this section we focus on the evolution of long wavelength perturbations, since instabilities on these scales would be the most dangerous for the success of BGC at providing a mechanism which allows exactly three spatial dimensions to become large and stabilizes the radius of the other dimensions at a microscopic value. We focus on the time interval $`\eta 1`$ when the hierarchy in scales between the large and the small spatial dimensions has already developed, and where we can use the analytical approximations to the background dynamics given in Section II.

Since we are focusing on long wavelength fluctuations, we can neglect all terms with spatial gradients in Equations (III)-(III). Making use of (28) we can eliminate one of the scalar metric perturbations, $`\varphi `$, from the rest of the Equations (III)-(III). We take Equation (III) as a constraint on the initial data, noticing that it only contains first derivatives of the perturbation variables. By using the background solution (14) we find that the remaining equations represent a coupled system of harmonic oscillators with time dependent coefficients. With these approximations, the equations become 777The full equations can be found in the Appendix.

$`2\chi ^{\prime \prime }2\psi ^{\prime \prime }6\xi ^{\prime \prime }{\displaystyle \frac{230}{7\eta }}\xi ^{}{\displaystyle \frac{176}{21\eta }}\psi ^{}{\displaystyle \frac{230}{21\eta }}\chi ^{}{\displaystyle \frac{12868}{3969\eta ^2}}\chi {\displaystyle \frac{6434}{3969\eta ^2}}\psi {\displaystyle \frac{12868}{1323\eta ^2}}\xi =0,`$ (34)
$`3\psi ^{\prime \prime }5\xi ^{\prime \prime }2\chi ^{\prime \prime }{\displaystyle \frac{244}{21\eta }}\psi ^{}{\displaystyle \frac{1978}{63\eta }}\xi ^{}{\displaystyle \frac{718}{63\eta }}\chi ^{}{\displaystyle \frac{2672}{1323\eta ^2}}\psi {\displaystyle \frac{5344}{441\eta ^2}}\xi {\displaystyle \frac{5344}{1323\eta ^2}}\chi =0,`$ (35)
$`{\displaystyle \frac{1}{2}}\chi ^{\prime \prime }{\displaystyle \frac{250}{63\eta }}\chi ^{}+{\displaystyle \frac{9}{7\eta }}\psi ^{}{\displaystyle \frac{170}{49\eta ^2}}\chi {\displaystyle \frac{43}{49\eta ^2}}\psi {\displaystyle \frac{510}{49\eta ^2}}\xi =0.`$ (36)
Notice that the string source terms, $`\delta N`$ and $`\delta M`$, have disappeared from the equations. This is due to the fact that, for long wavelength fluctuations, the terms due to purely gravitational dynamics (e.g. the damping terms in the evolution of the radius of the internal dimensions as a consequence of the expansion of the large dimensions) are more important than the matter sources. This is not too surprising based on the results of the theory of linear cosmological perturbations MFB  which show that on super-Hubble scale matter is pulled along by gravity, but that it is the purely gravitational dynamics which determines the growth rate of the fluctuations. In analogy, we find here that the terms representing the string matter sources are sub-leading in the equations (see Appendix). This realization is the key physical reason which leads to our ultimate conclusion that also in BGC there are no instabilities of long wavelength fluctuations.

The above set of coupled equations can be solved exactly, and we find the leading terms to be

$`\psi {\displaystyle \frac{c_1}{\eta ^{0.06}}}+{\displaystyle \frac{c_2}{\eta ^{0.36}}},`$ (37)
$`\xi c_3+{\displaystyle \frac{c_4}{\eta ^{0.06}}},`$ (38)
$`\chi c_5+{\displaystyle \frac{c_6}{\eta ^{0.06}}}.`$ (39)
This is our main result. First, we have found that the perturbation of the internal geometry $`\xi `$ is constant to leading order. Thus, the stabilization mechanism resulting from considering winding and momentum modes remains viable in the presence of string inhomogeneities and dilaton fluctuations. Secondly, the geometry of the $`3+1`$ dimensions growing large does not suffer from exponential instabilities. In fact, we have power-law decay for $`\psi `$ and the dilaton perturbation $`\chi `$ decays to a constant.

The crucial point that led to this behavior was the fact that the string source terms were sub-leading in the equations above. Physically this means that the number of winding and momentum modes moving in or out of a region is negligible compared to the expansion rate and the evolution of the perturbations. Therefore, there is an averaging of the long wavelength perturbations that leads to a smoothing effect and the string inhomogeneities are rendered harmless.

## V Conclusions

We have considered the evolution of linearized string inhomogeneities and dilaton fluctuations in the background of Brane Gas Cosmology (BGC). We considered fluctuations which are independent of the spatial coordinates of the internal dimensions. We have derived the perturbation equations for BGC in the longitudinal gauge and solved them in the long wavelength approximation. We have found that for long wavelengths, the effects of the string sources is subleading compared to the purely gravitational and dilaton terms in the fluctuation equations. This generalizes the result of the conventional theory of cosmological perturbations in Einstein gravity in four space-time dimensions which states that on scales larger than the Hubble radius, self-gravitational effects dominate the dynamics, and that the matter sources are simply dragged along by the metric. In particular, there are no instabilities of the background solution of BGC towards such long wavelength fluctuations. Thus, we find that the predictions of BGC are robust towards the effects of matter and dilaton fluctuations. In particular, the stabilization mechanism for the extra dimensions (volume modulus) remains operative in the presence of these inhomogeneities, while the $`3+1`$ dimensional space-time continues to grow without instabilities.

Although this is a promising result for BGC, there are still many questions to be answered. Our analysis has been limited to the linear regime, while non-linear effects may be crucial. Our linearized analysis corresponds to taking the localized string sources and smearing them out over a length scale larger than the Hubble radius to obtain a string gas description. It would be interesting to see if these results hold when considering individual strings and their interactions and in what limits the string gas approximation holds. Also, BGC offers a dynamical way to stabilize the volume modulus, but we have not addressed the issue of stabilizing the dilaton and we have also ignored the role of fluxes. We leave these questions for further investigations.

###### Acknowledgements.

SW would like to thank Alan Guth and Steve Gubser for useful discussions. RB was supported in part by the U.S. Department of Energy under Contract DE-FG02-91ER40688, TASK A. SW was supported in part by the NASA Graduate Student Research Program.

## Appendix A Perturbation equations

In this Appendix we present a detailed derivation of the solutions presented in section 4. In the long wavelength limit we can neglect gradient terms and the equations (III)-(III) reduce to

$`9\xi ^{}3\psi ^{}3^2\varphi ={\displaystyle \frac{1}{2}}e^{2\lambda +2\phi }\left(2\chi T_0^0+\delta T_0^0\right)6\varphi \phi ^{}3\psi ^{}\phi ^{}6\xi ^{}\phi ^{}+`$
$`+3\chi ^{}+2\varphi \phi ^22\chi ^{}\phi ^{},`$ (40)
$$_i_j\left(\varphi \psi 6\xi 2\chi \right)=0ij,$$
(41)
$`2\psi ^{\prime \prime }6\xi ^{\prime \prime }4\psi ^{}6\xi ^{}2^2\varphi 4^{}\varphi 2\varphi ^{}=e^{2\lambda +2\phi }\left(2\chi T_i^i+\delta T_i^i\right)+`$
$`4\varphi \phi ^{\prime \prime }2\varphi ^{}\phi ^{}4\varphi \phi ^{}4\psi ^{}\phi ^{}12\xi ^{}\phi ^{}+2\chi ^{\prime \prime }+2\chi ^{}+4\varphi \phi ^24\chi ^{}\phi ^{},`$ (42)
$`5\xi ^{\prime \prime }3\psi ^{\prime \prime }10\xi ^{}3\varphi ^{}9\psi ^{}6^2\varphi 6H^{}\varphi =e^{2\lambda +2\phi }\left(2\chi T_m^m+\delta T_m^m\right),`$
$`4\varphi \phi ^{\prime \prime }2\varphi ^{}\phi ^{}8\varphi \phi ^{}6\psi ^{}\phi ^{}10\xi ^{}\phi ^{}+2\chi ^{\prime \prime }+4\chi ^{}+4\varphi \phi ^24\chi ^{}\phi ^{},`$ (43)
$`2\varphi \phi ^2+2\phi ^{}\chi ^{}+\varphi \phi ^{\prime \prime }6\psi \phi ^{}+{\displaystyle \frac{1}{2}}\phi ^{}\varphi ^{}2\varphi \phi ^{}{\displaystyle \frac{3}{2}}\psi ^{}\phi ^{}3\xi ^{}\phi ^{}{\displaystyle \frac{1}{2}}\chi ^{\prime \prime }\chi ^{}`$
$`={\displaystyle \frac{1}{4}}e^{2\lambda +2\phi }\left(2\chi \mu T_\mu ^\mu +\delta T_\mu ^\mu \right).`$ (44)
Making use of the background solution (14) and (41) to eliminate the scalar metric perturbation $`\varphi `$, we find the following set of equations describing the evolution;

$`{\displaystyle \frac{61}{21\eta }}\chi ^{}{\displaystyle \frac{61}{7\eta }}\xi ^{}{\displaystyle \frac{88}{21\eta }}\psi ^{}+\left({\displaystyle \frac{11114}{1323\eta ^2}}{\displaystyle \frac{3000}{\eta ^{\frac{176}{63}}}}{\displaystyle \frac{12000}{\eta ^{\frac{169}{63}}}}\right)\chi +\left({\displaystyle \frac{36000}{\eta ^{\frac{169}{63}}}}{\displaystyle \frac{9000}{\eta ^{\frac{176}{63}}}}{\displaystyle \frac{11114}{441\eta ^2}}\right)\xi +`$
$`+\left({\displaystyle \frac{6000}{\eta ^{\frac{176}{63}}}}{\displaystyle \frac{5557}{1323\eta ^2}}{\displaystyle \frac{18000}{\eta ^{\frac{169}{63}}}}\right)\psi {\displaystyle \frac{3}{\eta ^{\frac{169}{63}}}}\delta M^{(6)}{\displaystyle \frac{3}{\eta ^{\frac{169}{63}}}}\delta N^{(6)}{\displaystyle \frac{3}{2\eta ^{\frac{176}{63}}}}\delta M^{(3)}=0,`$ (45)
$`2\chi ^{\prime \prime }2\psi ^{\prime \prime }6\xi ^{\prime \prime }{\displaystyle \frac{230}{7\eta }}\xi ^{}{\displaystyle \frac{176}{21\eta }}\psi ^{}{\displaystyle \frac{230}{21\eta }}\chi ^{}+\left({\displaystyle \frac{2000}{\eta ^{\frac{176}{63}}}}{\displaystyle \frac{12868}{3969\eta ^2}}\right)\chi +\left({\displaystyle \frac{4000}{\eta ^{\frac{176}{63}}}}{\displaystyle \frac{6434}{3969\eta ^2}}\right)\psi +`$
$`+\left({\displaystyle \frac{6000}{\eta ^{\frac{176}{63}}}}{\displaystyle \frac{12868}{1323\eta ^2}}\right)\xi {\displaystyle \frac{1}{\eta ^{\frac{176}{63}}}}\delta M^{(3)}=0,`$ (46)
$`3\psi ^{\prime \prime }5\xi ^{\prime \prime }2\chi ^{\prime \prime }{\displaystyle \frac{244}{21\eta }}\psi ^{}{\displaystyle \frac{1978}{63\eta }}\xi ^{}{\displaystyle \frac{718}{63\eta }}\chi ^{}{\displaystyle \frac{2672}{1323\eta ^2}}\psi +\left({\displaystyle \frac{2000}{\eta ^{\frac{169}{63}}}}{\displaystyle \frac{5344}{441\eta ^2}}\right)\xi +`$
$`{\displaystyle \frac{5344}{1323\eta ^2}}\chi +{\displaystyle \frac{1}{\eta ^{\frac{169}{63}}}}\delta N^{(6)}{\displaystyle \frac{1}{\eta ^{\frac{169}{63}}}}\delta M^{(6)}=0,`$ (47)
$`{\displaystyle \frac{1}{2}}\chi ^{\prime \prime }{\displaystyle \frac{250}{63\eta }}\chi ^{}+{\displaystyle \frac{9}{7\eta }}\psi ^{}+\left({\displaystyle \frac{170}{49\eta ^2}}{\displaystyle \frac{6000}{\eta ^{\frac{169}{63}}}}\right)\chi +\left({\displaystyle \frac{9000}{\eta ^{\frac{169}{63}}}}{\displaystyle \frac{43}{49\eta ^2}}\right)\psi +\left({\displaystyle \frac{510}{49\eta ^2}}{\displaystyle \frac{15000}{\eta ^{\frac{169}{63}}}}\right)\xi +`$
$`{\displaystyle \frac{3}{\eta ^{\frac{169}{63}}}}\delta N^{(6)}=0.`$ (48)
The first equation, containing only first derivatives, is taken as a constraint on the initial data. We are interested in the late time behavior, i.e. $`\eta 1`$. Keeping only the leading order terms the remaining equations can be approximated as,

$`2\chi ^{\prime \prime }2\psi ^{\prime \prime }6\xi ^{\prime \prime }{\displaystyle \frac{230}{7\eta }}\xi ^{}{\displaystyle \frac{176}{21\eta }}\psi ^{}{\displaystyle \frac{230}{21\eta }}\chi ^{}{\displaystyle \frac{12868}{3969\eta ^2}}\chi {\displaystyle \frac{6434}{3969\eta ^2}}\psi {\displaystyle \frac{12868}{1323\eta ^2}}\xi =0,`$ (49)
$`3\psi ^{\prime \prime }5\xi ^{\prime \prime }2\chi ^{\prime \prime }{\displaystyle \frac{244}{21\eta }}\psi ^{}{\displaystyle \frac{1978}{63\eta }}\xi ^{}{\displaystyle \frac{718}{63\eta }}\chi ^{}{\displaystyle \frac{2672}{1323\eta ^2}}\psi {\displaystyle \frac{5344}{441\eta ^2}}\xi {\displaystyle \frac{5344}{1323\eta ^2}}\chi =0,`$ (50)
$`{\displaystyle \frac{1}{2}}\chi ^{\prime \prime }{\displaystyle \frac{250}{63\eta }}\chi ^{}+{\displaystyle \frac{9}{7\eta }}\psi ^{}{\displaystyle \frac{170}{49\eta ^2}}\chi {\displaystyle \frac{43}{49\eta ^2}}\psi {\displaystyle \frac{510}{49\eta ^2}}\xi =0.`$ (51)
This set of coupled equations can be solved exactly and we find the leading terms to be

$`\psi {\displaystyle \frac{c_1}{\eta ^{0.06}}}+{\displaystyle \frac{c_2}{\eta ^{0.36}}},`$ (52)
$`\xi c_3+{\displaystyle \frac{c_4}{\eta ^{0.06}}},`$ (53)
$`\chi c_5+{\displaystyle \frac{c_6}{\eta ^{0.06}}}.`$ (54)
As a consistency check, we can plug this result back into the original equations (A)-(48). In addition we must consider the evolution of $`\delta N^{(6)}`$, $`\delta M^{(6)}`$ and $`\delta M`$ given by (32) and (33). We find that it was consistent to neglect the time evolution of the string sources (i.e. $`\delta N`$ and $`\delta M`$) compared with the expansion of the background and the evolution of the perturbations. In addition, we find that at late times this holds as an exact solution of the perturbation equations. In this way, we find a self consistent solution for the behavior of the long wavelength linear perturbations.