text:
# Quantum initial value representations using approximate Bohmian trajectories

## I Introduction

Over the past few years there has been significant effort in the development of the so-called hydrodynamic or quantum trajectory based approach to quantum mechanics. rabitz98a ; qtm1 ; mayor ; bw1 ; wb1 ; erb ; maddox ; maddox3 ; irene1 ; irene2 ; irene3 ; irene4 ; hw ; tw1 ; bwreview ; maddox2 ; makri  This formulation is based upon Bohm’s “hidden variable” representation bohm ; debroglie ; holland  is initiated by substituting the amplitude-phase decomposition of the time-dependent quantum wavefunction, $`\psi (y,t)=R(y,t)\mathrm{exp}(iS(y,t)/\mathrm{})`$, into the time-dependent Schrödinger equation,

$`i\mathrm{}{\displaystyle \frac{\psi }{t}}=\left({\displaystyle \frac{\mathrm{}^2}{2m}}^2+V\right)\psi `$ (1)
Substituting Eq. 1 into Eq. 2 and separating the real and imaginary components yields equations of motion for $`\rho `$ and $`S`$,

$`{\displaystyle \frac{\rho }{t}}`$ $`=`$ $`(\rho S)/m`$ (2)
$`{\displaystyle \frac{S}{t}}`$ $`=`$ $`{\displaystyle \frac{SS}{2m}}+V{\displaystyle \frac{\mathrm{}^2}{2m}}{\displaystyle \frac{1}{\sqrt{\rho }}}^2\sqrt{\rho }=0`$ (3)
where $`\rho =|R|^2`$ is the probability density and $`S`$ is the phase. These last two equations can be easily identified as the continuity equation (Eq. 2) and the quantum Hamilton-Jacobi equation (Eq. 3). In Eq. 3 we can identify three contributions to the action. The first two are the kinetic and potential energies and the third is the quantum potential $`Q`$ bohm ; holland . This term is best described as a shape kinetic energy since it is determined by the local curvature of the quantum density. In numerical applications, computation of the quantum potential is frequently rendered more accurate if derivatives are evaluated using the log-amplitude $`C=\mathrm{log}(R)`$ ($`C`$ is also referred to as the $`C`$-amplitude). In terms of derivatives of this amplitude, the quantum potential is $`Q=\mathrm{}^2(^2C+(C)^2)/2m`$.

Eq. 2 and 3 are written in terms of partial derivatives in time. In order to compute the time-evolution of $`\rho `$ (or equivalently $`C`$) and $`S`$ along a give path, $`x(t)`$ we need to move to a moving Lagrangian frame by transforming $`d_t=_t+\dot{x}(t)_x`$. The Lagrangian form of the hydrodynamic equations of motion resulting from this analysis are given by:

$`{\displaystyle \frac{d\rho }{dt}}`$ $`=`$ $`\rho v`$ (4)
$`{\displaystyle \frac{dv}{dt}}`$ $`=`$ $`{\displaystyle \frac{1}{m}}(V+Q)`$ (5)
$`{\displaystyle \frac{dS}{dt}}`$ $`=`$ $`{\displaystyle \frac{1}{2m}}(S)^2(V+Q)`$ (6)
in which the derivative on the left side is appropriate for calculating the rate of change in a function along a fluid trajectory. Eq. 6 is a Newtonian-type equation in which the flow acceleration is produced by the sum of the classical force, $`f_c=V`$, and the quantum force is $`f_q=Q`$

The quantum trajectories obey two important noncrossing rules: (1) they cannot cross nodal surfaces (along which $`\psi =0`$ ); (2) they cannot cross each other. (In practice, because of numerical inaccuracies, these conditions may be violated.) Because of these two conditions, the quantum trajectories are very different from classical paths and represent a geometric optical rendering of the evolution of the quantum wave function.

The primary difficulty in implementing a quantum trajectory based approach is in constructing the quantum potential. In spite of 3-4 years of intense effort by a number of groupsrabitz98a ; qtm1 ; mayor ; bw1 ; wb1 ; erb ; maddox ; hw ; tw1 ; bwreview ; maddox2  , including our own, there has yet to emerge a satisfactory and robust way to compute $`Q`$ and the quantum force $`\stackrel{}{}Q`$ directly from the quantum trajectories in multidimensional systems. This severely limits the applicability of a quantum trajectory based approach for computing exact quantum dynamics for all but the most trivial systems.

In this paper, we show that quantum trajectories can be used to construct exact quantum initial value representations (IVR) of correlation functions. The basic equations are extremely straightforward to derive and follow directly from a Bohm trajectory based description. The heart and sole of our approach, however, is the use of an cummulant/derivative propagation scheme tw2  that generates local equations of motion for a given quantum trajectory and all the derivatives of $`C`$ and $`S`$ along the path, allowing us to compute the correlation function on a trajectory by trajectory basis. In our numerical implementation of this approach, we examine a highly anharmonic one dimensional system (an inverted gaussian well) and find that the Bohm-IVR approach is surprisingly accurate in computing correlation functions even though at long-times the quantum trajectories themselves show crossings.

## II The Bohm Initial Value Representation

While in principle, the quantum wavefunction or the density matrix tells us everything we can know about a given physical process, often what we are interested in is a correlation function or an expectation value. This can be equivalently cast as either a time-correlation or a spectral function via Fourier transform:

$$\widehat{A}(\omega )=[\widehat{A}(t)\widehat{A}].$$

Where

$`\widehat{A}(t)\widehat{A}`$ $`=`$ $`\psi _o|e^{+iHt/\mathrm{}}\widehat{A}e^{iHt/\mathrm{}}\widehat{A}|\psi _o`$ (7)
is the time-correlation function for some quantum operator $`\widehat{A}`$ following initial state $`\psi _o`$ evolving under Hamiltonian, $`H`$.

Within the Bohmian framework, the synthesis of a quantum wavefunction along a quantum trajectory, $`x(t)`$, is given by

$`\psi (x,t)`$ $`=`$ $`{\displaystyle 𝑑x(t)\delta (xx(t))J^{1/2}(x(t),x(0))}`$ (8)
$`\times `$ $`e^{iS(x(0),x(t))/\mathrm{}}\psi _o(x(0)).`$
Here,

$`J(x(t),x(0))=e^{+_0^tv𝑑s}`$ (9)
is the Jacobian for the transformation of a volume element $`dx(0)`$ to $`dx(t)`$ along path $`x(t)`$, $`dx(t)=J(x(0),x(t))dx(0)`$. Its time evolution follows directly from the continuity equation with initial condition $`J(x(0),x(0))=1`$. Likewise, $`S`$ is the action integral integrated over the quantum trajectory connecting $`x(0)`$ and $`x(t)`$. Using this, we can rewrite Eq. 7 as an integral over starting points for Bohmian trajectories:

$`\widehat{A}(t)\widehat{A}={\displaystyle 𝑑x(0)𝑑y(0)A(x(t),y(t))J_x^{1/2}J_y^{1/2}e^{iS(x(0),x(t))/\mathrm{}}\psi ^{}(x(0))e^{iS_A(y(0),y(t))/\mathrm{}}\psi _A(y(0))}`$ (10)
where $`|\psi _A=\widehat{A}|\psi `$ and we assume that $`\widehat{A}`$ is a scalar quantum operator with matrix elements $`x|\widehat{A}|y=A(x,y)`$. In this expression, the forward-going quantum path along $`y(0)`$ to $`y(t)`$ is subject to a quantum potential derived from the evolution of $`\psi _A(t)`$ (with action integral $`S_A(y(0),y(t))`$ and Jacobian $`J_y`$), where as the reverse path along $`x(0)`$ to $`x(t)`$ has a quantum potential derived from the evolution of $`\psi (t)`$ (with action integral $`S(x(0),x(t))`$ and Jacobian $`J_x`$.) Thus, even if $`x(0)=y(0)`$ the forward path ($`y(s)`$) will be very different than the reverse path. We can interpret this in the following way. At time $`t=0`$, we act on the initial wave function with operator $`\widehat{A}`$. If $`\psi `$ is not an eigenstate of $`\widehat{A}`$, then we have a new wave function $`\widehat{A}|\psi =|\psi _A`$ which evolves out to time $`t`$ where we act again with $`\widehat{A}`$ and $`y(t)`$ is instantly scattered to $`x(t)`$ where it evolves back in time under the reversed evolution of $`\psi (t)`$. The IVR correlation function then measures the probability amplitude that the final forward-backward state at $`x(0)`$ was a member of the initial state at $`y(0)`$. This is illustrated schematically in Fig. 1.

So far, no approximations have been made. The simplest correlation function is that the autocorrelation function of the wavefunction which monitors the overlap between the initial wavefunction $`\psi _o`$ and the time-evolved wavefunction.

$`c(t)`$ $`=`$ $`\psi (0)|\psi (t)`$ (11)
$`=`$ $`{\displaystyle }dx(0)\psi _o^{}(x(t)\psi _o(x(0))J^{1/2}e^{iS/\mathrm{}}`$
Here, $`\psi _o(x)`$ is the initial wave function. This is directly analogous to the semi-classical IVR developed extensively over the past number of years. makri ; gr1 ; mchild ; makri02 ; miller03  In the SC version, one estimates the quantum propagator via the well known Van Vleck expression vvleck .

$`c(t)`$ $``$ $`{\displaystyle 𝑑x(0)\psi _o^{}(x(t))K_{SC}(x(t),x(0);t)\psi _o(x(0))}.`$ (12)
Here, $`K_{SC}(x(t),x(0);t)`$ is the semi-classical propagator

$`K_{SC}={\displaystyle \underset{cl.paths}{}}{\displaystyle \frac{C^{1/2}}{(2\pi i\mathrm{})^{D/2}}}e^{iS_{cl}i\kappa \pi /2}`$ (13)
where the summation is over all classical paths connecting $`x(0)`$ to $`x(t)`$ with classical action $`S_{cl}`$,

$`S_{cl}={\displaystyle _0^t}{\displaystyle \frac{m}{2}}\dot{x}^2(s)V(x(s))ds,`$ (14)
$`C`$ is the matrix of the negative second variations of the classical action with the end points,

$`C`$ $`=`$ $`{\displaystyle \frac{^2S_{cl}}{x(0)x(t)}}`$ (15)
and $`\kappa `$ gives the Morse index obtained by counting the number of conjugate points along each trajectory. What is appealing in the Bohm-IVR expression above is that one does not need to compute $`C`$ or $`\kappa `$ along each trajectory. This is not surprising since we are dealing with a purely quantum mechanical object and $`C`$ and $`\kappa `$ are semi-classical quantities.

A closer connection between Eq. 12 and Eq. 11 can be drawn by writing $`C`$ as

$`C`$ $`=`$ $`{\displaystyle \frac{p(t)}{x(0)}}=J{\displaystyle \frac{^2S_{cl}}{x(t)^2}}`$ (16)
Thus, the SC-IVR expression can be written as

$`c_{sc}(t)`$ $`=`$ $`{\displaystyle \frac{1}{(2\pi i\mathrm{})^{D/2}}}{\displaystyle 𝑑x(0)\underset{cl.paths}{}J^{1/2}\left|\frac{^2S_{cl}(t)}{x^2(t)}\right|^{1/2}e^{iS_{cl}/\mathrm{}i\kappa \pi /2}\psi _o^{}(x(t))\psi _o(x(0))}.`$ (17)
Comparing Eq. 17 to Eq. 11 we first note the appearance of $`J^{1/2}`$ in both expressions. Here, as in the case above, $`J`$ is the Jacobian for the transformation of a volume element $`dx(0)`$ to $`dx(t)`$ along the *classical path* $`x_{cl}(t)`$. What is most striking, however, is that in the exact Bohm IVR expression one has a *single trajectory connecting each initial point to each final point* where as in the semi-classical case, we have to consider every possible trajectory that starts at $`x(0)`$. In the semi-classical case, this can take the form of a phase-space integration, as in the widely used Herman-Kluk expressionhk ; hkd ; mh86 ; gr1 ; mchild ; batista ; makri02 ; gara-light .

As mentioned above, one of the primary difficulties encountered in developing numerical methods based upon the Bohm trajectory approach is in computing the quantum potential, $`Q`$, and its derivatives on a generally unstructured mesh of points evolving according to the Bohmian equations of motion. Moreover, the trajectories themselves can be highly erratic and kinky requiring that one use very accurate numerical integration methods to achieve an accurate solution of the equations of motion. makri  Consequently, obtaining an accurate estimate of the time-evolving wavefunction (via the synthesis equation above) is quite difficult. On the other hand, one may be able to get by with a less robust and less exact calculation of the quantum trajectories if one is only interested in computing time-correlation functions. We next turn our attention towards a hierarchy approach for generating approximate quantum trajectories.

## III Derivative Propagation

Recently, Trahan, Hughes, and Wyatt tw2  introduced an interesting approach for computing quantum trajectories by simultaneously evolving the spatial derivatives of both $`C`$ and $`S`$ along a given trajectory. This notion is similar in spirit to the hydrodynamic moment expansions of the density matrixmaddox ; maddox3 ; irene1 ; irene2 ; irene3 ; irene4  and follows directly by taking derivatives of $`C/t`$ and $`S/t`$ with respect to $`x`$ and transforming to the moving Lagrangian representation. Following this notion, the time-evolution of the $`n`$-th spatial derivatives of $`C`$ and $`S`$ at $`x_k(t)`$ are given by

$`{\displaystyle \frac{dC_k^{(n)}}{dt}}`$ $`=`$ $`{\displaystyle \frac{1}{2m}}\left(S_k^{(n+2)}+2{\displaystyle \underset{j=0}{\overset{n}{}}}\left(\begin{array}{cc}n& \\ j& \end{array}\right)S_k^{(n+1j)}C_k^{(j+1)}\right)+S_k^{(1)}C_k^{(n+1)}/m`$ (20)
$`{\displaystyle \frac{dS_k^{(n)}}{dt}}`$ $`=`$ $`{\displaystyle \frac{1}{2m}}\left({\displaystyle \underset{j=0}{\overset{n}{}}}\left(\begin{array}{cc}n& \\ j& \end{array}\right)S_k^{(j+1)}S_k^{(n+1j)}\right)`$ (23)
$`+`$ $`{\displaystyle \frac{\mathrm{}^2}{2m}}\left(C_k^{(n+2)}+{\displaystyle \underset{j=0}{\overset{n}{}}}\left(\begin{array}{cc}n& \\ j& \end{array}\right)C_k^{(j+1)}C_k^{(n+1j)}\right)V_k^{(n)}+S_k^{(1)}S_k^{(n+1)}/m`$ (26)
Our notation is that $`f_k^{(n)}=^nf`$ evaluated at $`x_k(t)`$ and $`(_j^n)`$ is the binomial coefficient. These equations, along with $`dx_j/dt=S_j^{(1)}/m`$ are an infinite hierarchy of equations for both $`C`$ and $`S`$ and their derivatives evaluated along the path $`x_j(t)`$. In this hierarchy, the equation for the $`n^{th}`$ order derivative of $`C`$ involves up to the $`n+2`$ and order derivative of $`S`$ and up to the $`n+1`$ order derivative of $`C`$. Likewise, the equation of motion for the $`n^{th}`$ order derivative of $`S`$ involves up to 2 higher order terms of $`C`$ and one higher order term of $`S`$. The $`C_k^{(n)}`$ and $`S_k^{(n)}`$ terms are in essence $`n^{th}`$-order cumulants of the wavefunction evaluated at $`x_k(t)`$. Note also that the hierarchy equations are *l*ocal in that they do not couple different trajectories. Each trajectory can be computed as an independent quantity and the task of computing multiple quantum trajectories can be easily distributed over multiple machines.

Imbedded in these equations of motion is the mechanism by which complex structure emerges within the wave function since even simple wave functions with few cumulants will rapidly evolve into wave functions with high order cumulants. Thus, in order to close the equations of motion, we need to introduce an artificial boundary condition which limits the extent of the hierarchy. The most straight-forward scheme is to simply require any term, $`C_j^{(n)}`$ or $`S_j^{(n)}`$ with $`n>n_{order}`$ to be set to zero in the initial equations of motion. One can also introduce “absorbing boundary conditions” by damping higher-order terms, e.g. $`dS_j^{(n)}/dt=\lambda S_j^{(n)}`$.

To compute correlation functions, we first start off by sampling $`\psi _o`$ on a uniform grid of DVR points such that $`\psi _{on}=\psi _o(x_n(0))w_n^{1/2}`$ where $`w_n^{1/2}`$ is the gaussian quadrature weight associated with the $`n`$-th grid point. Thus, the correlation function can be evaluated by gaussian quadrature over a grid defined by the initial positions of the Bohm trajectories. For example, for the auto-correlation function of the wavefunction:

$`c(t)={\displaystyle \underset{n}{}}w_n\psi _o^{}(x_n(t))\psi _o(x_n(0))J_n^{1/2}e^{iS_n/\mathrm{}}`$ (27)
where $`J_n`$ and $`S_n`$ are the Jacobian and quantum action along $`x_n(s)`$ which starts at $`x_n(0)`$ and ends at $`x_n(t)`$.

Written in terms of the derivative terms above, the wave function correlation becomes simply

$`c(t)={\displaystyle \underset{j=1}{\overset{N}{}}}w_je^{2C_j^{(0)}(t)}e^{iS_j^{(0)}(t)/\mathrm{}}\psi _o^{}(x_j(t))\psi _o(x_j(0))`$ (28)
where $`C_j^{(n)}`$ and $`S_j^{(n)}`$are computed using the above hierarchy equations along a finite set of trajectories. As each trajectory is independent, the sum can evaluated by distributing the trajectory calculations amongst several nodes, each running a single quantum trajectory. This last expression along with its implementation comprise the central results of this paper. In the following section we examine a series of one dimensional test cases.

## IV Sample Calculations

As discussed above, quantum interference in the wave function leads to very complex structure in the quantum potential. Even in the simple double-slit experiment example, the quantum potential due to the interference between the two gaussian components of the wave function varies rapidly as the underlying quantum trajectories avoid one another. holland  Consequently, in numerical applications of the quantum trajectory approach, nodes and interferences in the wave function have present considerable hurdles in the development of generalized multi-dimensional applications of this method.

Here, we take three related test cases with increasing anharmonicity and test the ability of the Bohm-IVR approach with derivative propagation to obtain accurate correlation functions. In each case we consider an electron ($`m=1`$) in one-dimensional potentials, chosen such that the curvature at the bottom is identical in each case.

$`V_{HO}`$ $`=`$ $`x^2/21`$ (29)
$`V_{QO}`$ $`=`$ $`x^2/2+0.01x^41`$ (30)
$`V_G`$ $`=`$ $`e^{x^2/2}`$ (31)
$`V_{HO}`$ is the harmonic well, $`V_{QO}`$ is harmonic well perturbed by a small quartic term, and $`V_G`$ is an inverted gaussian well. In each case, the initial state $`\psi (x)=\mathrm{exp}(c+is/\mathrm{})`$ is gaussian with $`c`$ and $`s`$ given by

$`c(x)={\displaystyle \frac{1}{4}}\mathrm{log}(\beta /\pi )\beta (xx_o)^2`$ (32)
$`s(x)=0`$ (33)
with $`x_o=1`$ and $`\beta =1`$. (Note, atomic units used throughout.) A plot of each potential along with the location of the lowest energy bound states is show in Fig. 2.

### IV.1 Harmonic Oscillator

Harmonic systems are particularly useful as test cases since the an initial gaussian wave function retains its gaussian form for all time. Furthermore, any quantity of interest can be computed analytically thereby providing a useful benchmark for the method. In Fig. 3 we show the quantum trajectories computed via the derivative propagation scheme outlined above. The trajectories display the coherent oscillation one expects for Bohm trajectories in this system. This is not surprising since a low order truncation of the derivative propagation equations is perfectly valid for this case. These agree perfectly with the analytical results for this system and we can integration the derivative propagation equations out to very long times without significant loss of accuracy. The initial points were chosen as evenly spaced gaussian quadrature points.

The solid curve in Fig. 4 shows the time correlation function of the wavefunction and its Fourier transform. Here, we performed the calculation out to $`t=20\pi `$ corresponding to 10 complete oscillation cycles. Both $`C(t)`$ and its spectral representation $`\stackrel{~}{C}(\omega )`$ agree perfectly with analytical results with $`\stackrel{~}{C}(\omega )`$ accurately mapping out the harmonic progression of the energy levels in this well $`\omega _n=\omega (n+1/2)`$ as well as the expansion coefficients for the projection of the initial wavefunction on to the eigenstates of the well, $`c_n=|n|\psi _o|^2`$. For the case at hand, the first 4 coefficients are $`c=\{0.606531,0.303265,0.0758163,0.0126361\}`$ and agree nicely with the roughly 8:4:1 ratio of the first three peaks in $`C(\omega )`$.

### IV.2 Quartic Perturbed Harmonic

In this case, we apply a slight anharmonic squeeze to the harmonic well. In an anharmonic system, an initial gaussian wave function retains its gaussian form for short period of time before bifurcating into a multitude of components due to interferences. Such features are most likely to occur where the underlying classical trajectory crosses a conjugate or focal point. Such crossings are strictly forbidden for quantum Bohmian trajectories and lead to rapidly varying quantum forces over a short spatial region. Consequently, interference effects in anharmonic systems has proven to be somewhat of a bugbear for implementing quantum trajectory approaches based upon fitting the quantum potential (and its derivatives) to a simple polynomial over a set of quantum trajectories.

The derivative propagation trajectories also retain their coherences for a short period of time. However, after only one or two oscillations, crossings become more frequent and the ensemble of trajectories degrades into an incoherent mess. Had we performed this with the moving-weighted least squares (MWLS) or other methods which extracts the quantum potential from an ensemble of trajectories, the results would have been cataclysmic. Our experience has indicated that even small errors in constructing the quantum force from the ensemble can lead to a rapid amplification of error in the propagation. Remarkably, even though the ensemble has apparently lost coherence amongst its members, the individual trajectories exhibit non-classical deflections, and quasi-periodic behavior one expects from truly quantum trajectories.

Turning towards the autocorrelation and its Fourier transform. $`C(t)`$ shows some recurrences at $`t2\pi n`$ as in the harmonic case, but these rapidly die out after 5-6 beats. This limits our ability to resolve the eigenenergies in the Fourier transform. Nonetheless, at least 4 peaks can be distinguished.

### IV.3 Inverted Gaussian Potential

In this final example, $`V_G`$ is an inverted gaussian and supports only two eigenstates (at $`E=0.59386`$ and $`E=0.0356576`$). In this example we used 30 trajectories since the majority of the trajectories near the leading edges of the distribution leave the potential region almost immediately. In Fig. 5 we show 10 trajectories that started in the center of the initial wavepacket. As time progresses, many of these eventually escape the potential well leaving only one of the original 30 trajectories in the well.

The long-time accuracy of the ensemble is not so important since we are ultimately interested in computing $`C(t)`$. This we show in Fig. 4 along with its Fourier transform. One can clearly identify 6 recursion peaks spaced roughly every 0.35 ps. These correspond to the nearly harmonic motion of the trajectories closest to the center of the ensemble.

## V Discussion

For low dimensional systems, there are arsenals of robust and computationally efficient methodologies for computing exact quantum dynamics for both time dependent and time-independent systems. However, these methods are severely limited by dimensionality since the computational effort and storage required to perform an exact quantum calculation increases almost exponentially with dimensionality. Thus, approximate methods, such as those introduced here, provide a viable path way towards higher dimensions with increasing complexity. While we have focused entirely upon one-dimensional examples, the Bohm-IVR and derivative propagation equations are exact and can be easily implemented in a multi-dimensional setting.

What do the individual trajectories mean?  Are the trajectories computed via the derivative propagation methods really Bohm trajectories? The answer, I believe is a definite yes and no. Yes, they are Bohm trajectories, but not for a initial single gaussian wavepacket evolving in a given potential. They correspond to trajectories arising from a diverging set of wave functions where each wave function component is constrained to have a particular functional form dictated by truncating the hierarchy. In essence, by truncating the hierarchy, one introduces dephasing and decoherence into the dynamics of a quantum system. This is evident in the time-correlation function for the particle in the harmonic + quartic potential. The time correlation function for a fully coherent wavepacket should show stronger recurrences and eventually be able to resolve the eigenenergies of the system upon Fourier transform. One evidence for this is that for bound systems other than the harmonic oscillator, the autocorrelations we compute using the derivative propagation decay to zero after a few recurrences, and never recover coherence.

For exact quantum mechanics, this implicit loss of coherence at long-times clearly undesirable. However, for condensed phase systems, especially in cases where we have a low dimensional quantum subsystem in contact with a high-dimensional “bath” of oscillators, reduction of the density matrix of the subsystem to a gaussian form in the coherences occurs spontaneously due to decoherence. maddox ; irene3 ; irene4 . Consequently, so long as the “physical” decoherence time is smaller than the “truncation induced ” decoherence time, time-correlation functions computed using the Bohm-IVR approach we describe herein are expected to quite accurate.

###### Acknowledgements.

This work was supported by the National Science Foundation and the Robert A. Welch Foundation. The author wishes to acknowledge discussions with Jeremy Maddox regarding this work and Prof. Bob Wyatt for providing us with a preprint of Ref.tw2 .


format: text

text:
# Ising-like Spin Anisotropy and Competing Antiferromagnetic - Ferromagnetic Orders in GdBaCo2O5.5 Single Crystals

On leave from \]Institute of Semiconductor Physics, Siberian Branch of RAS, 630090 Novosibirsk, Russia.

## Abstract

In $`R`$BaCo2O5+x compounds (R is rare earth), a ferromagnetic-antiferromagnetic competition is accompanied by a giant magnetoresistance. We study the magnetization of detwinned GdBaCo2O5.5 single crystals, and find a remarkable uniaxial anisotropy of Co3+ spins which is tightly linked with the chain oxygen ordering in GdO0.5 planes. Reflecting the underlying oxygen order, CoO2 planes also develop a spin-state order consisting of Co3+ ions in alternating rows of $`S=1`$ and $`S=0`$ states. The magnetic structure appears to be composed of weakly coupled ferromagnetic ladders with Ising-like moments, which gives a simple picture for magnetotransport phenomena.

Transition-metal oxides exhibit a complex interplay between charge, spin, orbital and lattice degrees of freedom, which is at the heart of many fascinating phenomena such as colossal magnetoresistance (MR) in manganites. Phenomenologically, the colossal MR originates from a magnetic-field induced transition between delicately balanced and competing phases that dramatically differ in resistivity CMR\_rev . However, the microscopic origin of these competing phases still remains far from being clear, owing to the complexity of the manganites; for example, magnetic fields (which induce a ferromagnetic (FM) spin alignment in manganites) have been found to melt the charge order melting , change the orbital order with accompanying Jahn-Teller distortions Nojiri , and modify the scale and topology of domains in microscopically heterogeneous phases CMR\_rev ; inhomo .

In order to clarify the MR mechanisms operating in transition-metal oxides, one may look for compounds that also exhibit competing phases and giant MR, but possess fewer degrees of freedom. Recently, an intriguing antiferromagnetic (AF) - ferromagnetic competition accompanied with a giant MR has been found in cobalt-oxide compounds RBaCo2O5+x (where R is rare earth) Martin ; Troy ; Respaud ; Akahoshi . While cobaltites with a variable spin state of Co-ions Martin ; Troy ; Respaud ; Akahoshi ; Vogt ; Moritomo ; Korotin  are in general no simpler than manganites, RBaCo2O5+x with $`x0.5`$ possesses a plenty of attractive features that may turn it into a model system. For $`x0.5`$, oxygen ions in ROx layers order into alternating filled and empty rows running along the $`a`$ axis Martin ; Moritomo ; Akahoshi , thus allowing one to deal with almost perfect crystal structure. No considerable structural distortions are detected at the AF-FM transition as well Moritomo ; Fauth . Moreover, at $`x0.5`$ the valence state of Co-ions approaches 3+, so that charge ordering and electronic phase separation typical for mixed-valence systems CMR\_rev  become irrelevant for this composition.

In this Letter, we present a first study of the static magnetization in detwinned GdBaCo2O5.5 single crystals. We find a strong Ising-like spin anisotropy which is clearly manifested in the FM phase, where Co spins are aligned along the oxygen-chain direction with a net magnetic moment of $`1`$ $`\mu _\text{B}/`$Co. This value indicates a 1:1 ratio of Co3+ ions in the $`S=0`$ and $`S=1`$ spin states; based on this result and a recent structural study Frontera , we discuss that the oxygen chain ordering causes the $`S=1`$ Co3+ ions to be arranged into ferromagnetic 2-leg ladders that are separated by non-magnetic layers. The AF-FM competition and giant MR in GdBaCo2O5.5 appear to reflect a relative magnetic order among weakly coupled FM ladders, which can be rather easily altered by temperature or magnetic fields.

Using the floating-zone technique, we have succeeded in growing high-quality GdBaCo2O5+x single crystals suitable for magnetization and transport measurements. The crystals were annealed in a flow of oxygen at 470C, the temperature being carefully tuned unpubl  to provide the oxygen stoichiometry $`x=0.50`$. Parallelepiped samples with all faces adjusted to the crystallographic planes with a 1-accuracy were prepared by cutting and polishing under Laue x-ray diffraction control. In GdBaCo2O5.5, the oxygen ordering in GdO0.5 planes induces a tetragonal-to-orthorhombic transition, which is accompanied by heavy twinning of crystals that mixes the $`a`$ and $`b`$ orthorhombic axes. To detwin crystals, we slowly cooled them under a uniaxial pressure of $`0.15`$ GPa from $`260^{}`$C, using an optical microscope to control the twin removing twins ; note that these crystals are fragile and only few survive the detwinning procedure. According to X-ray measurements, the remaining fraction of misoriented domains (which characterizes the quality of detwinning) was 4-5%. Magnetization measurements were carried out using a SQUID magnetometer at fields up to 7 T applied along $`a`$, $`b`$ or $`c`$ axis. Throughout this paper, the magnetization coming from Co ions is determined by subtracting the contribution of Gd-ions, assuming their ideal paramagnetic (PM) behavior with total spin $`S=7/2`$ \[see inset in Fig. 1(a)\]; the latter is a good approximation since no ordering of Gd3+ moments is detected down to 1.7 K. The magnetoresistance was measured by a four-probe method on twinned crystals.

Figure 1(a) shows the magnetization of untwinned GdBaCo2O5.5 single crystals measured along $`a`$, $`b`$, and $`c`$ axes in the field-cooling (FC) process. Below 300 K, a net ferromagnetic component appears in the Co sublattice and suddenly vanishes at $`260`$ K, indicating successive PM-FM-AF transitions Martin ; Troy ; Moritomo ; Respaud ; Akahoshi . A remarkable feature of the FM state is not only the narrow temperature window where it shows up, but also a very strong anisotropy: We have found that the net FM moment appears only along the $`a`$ axis, and even the 7-T field is not enough to turn over the magnetic moment along the $`b`$ or $`c`$ axis. Note that a small magnetization along the $`b`$ axis in Fig. 1(a) comes mostly from residual misoriented domains. This behavior suggests the spin system in GdBaCo2O5.5 to be Ising-like, which drastically simplifies the understanding of magnetic ordering.

The balance of FM and AF ordering in GdBaCo2O5+x turns out to be quite delicate, and can be easily affected by temperature, magnetic fields, or even subtle variation in stoichiometry. Magnetic fields applied along the spin-easy $`a`$ axis stabilize the FM state and shift the FM-AF transition to lower temperatures \[Fig. 1(b)\]. Nevertheless, it would not be correct to consider the FM and AF orders as simply competing ones. If this switching were originating from equally strong and competing AF and FM exchange interactions, spin fluctuations would inevitably be enhanced in the vicinity of the FM-AF phase boundary, smearing the transition. Isothermal magnetization curves in Fig. 2(a) demonstrate, however, that although the FM-AF balance is subtle, and fairly weak fields ($`1.5`$ T at 240 K) are capable of recovering the FM order, the AF-FM transition remains always sharp, showing that thermal fluctuations are irrelevant here. This behavior clearly indicates that the observed AF-FM switch occurs within the ordered spin state and is governed by reorientation of spins in one of weakly-coupled sublattices.

The fact that the AF-FM transition remains sharp, even when it is induced by a weak field $`\mu _\text{B}HkT`$ \[Fig. 2(a)\], unambiguously points to a strong hierarchy of spin interactions in GdBaCo2O5.5. A strong ferromagnetic interaction aligns spins within each sublattice, while a much weaker antiferromagnetic (at $`T<260`$ K) interaction provides a subtle coupling between the sublattices, which can be broken by an applied field. Figure 2(c) shows that a magnetic field $`Ha`$ required to overcome the AF coupling grows roughly linearly upon cooling, from zero at $`T260`$ K up to $`20`$ T at $`T=0`$ (Fig. 2(c)). Whatever the temperature, however, the 7-T field $`Hb`$ or $`Hc`$, cannot compete with the spin anisotropy and causes just a partial tilting of spins (in each sublattice) from their easy axis, thus giving linear $`M(H)`$ curves in both FM and AF regions (Fig. 2(b)).

Qualitatively, the overall magnetic behavior is quite reminiscent of that in canted antiferromagnets, such as La2CuO4 LCO , where slight spin canting brings about a weak ferromagnetism. However, the magnitude of the field-induced FM moment, reaching $`0.6`$ $`\mu _\text{B}/`$Co at $`T=205`$ K \[Fig. 2(a)\], is apparently inconsistent with such weak spin canting, and points to a different origin of the ferromagnetism. A rough extrapolation to $`T=0`$ suggests a saturated magnetic moment of $`1`$ $`\mu _\text{B}/`$Co \[Fig. 2(d)\], which corresponds to a 1:1 mixture of low-spin (LS: $`t_{2g}^6`$, $`e_g^0`$; $`S=0`$) and intermediate-spin (IS: $`t_{2g}^5`$, $`e_g^1`$; $`S=1`$) states of Co3+ ions, if a simple FM spin order is realized. A similar conclusion, that Co3+ ions exhibit a 1:1 ratio of LS and IS states below the metal-insulator transition at $`T360`$ K Respaud ; Frontera ; Martin , can be reached based on the Curie-Weiss fitting of the PM susceptibility in the temperature range 300-360 K \[Fig. 1(b)\]. It is worth noting that polycrystalline RBaCo2O5.5 samples demonstrate smaller FM moments Troy ; Respaud ; Akahoshi , seemingly inconsistent with the expected Co spin states, which has been one of the mysteries of the FM state. The discovered Ising anisotropy, that prevents moments from being seen along the $`b`$ and $`c`$ axes, readily resolves this discrepancy.

To understand why Co3+ ions in GdBaCo2O5.5 even at zero temperature exhibit two different spin states, one should consider the oxygen ordering in GdO0.5 planes, sketched in Fig. 3(a). The alternating filled and empty oxygen chains create two types of structural environment, octahedral and pyramidal, for Co ions; one of these environments stabilizes the LS ground state ($`t_{2g}^6`$, $`e_g^0`$; $`S=0`$), while the other one makes the IS state ($`t_{2g}^5`$, $`e_g^1`$; $`S=1`$) preferable. While our experiments can hardly distinguish whether the IS ground state is realized in pyramidal or octahedral positions, a structural study by Frontera et al. Frontera  suggests the IS state for Co ions in pyramidal sites. In this case, magnetic Co ions form 2-leg ladders extended along the $`a`$ axis \[black spheres in Fig. 3(a)\] and separated by non-magnetic layers.

The magnetic ordering in insulators is known to be predominantly caused by the superexchange (SE) interaction, whose sign for each pair of ions can be estimated using Goodenough-Kanamori rules Anderson . The FM spin order in ladders can be explained by the orbital ordering among IS Co-ions, an example of which is shown in Fig. 3(a). Owing to the reduced dimensionality (quasi-1D/2D) of the ladders, the FM order develops quite gradually upon cooling, being subject to strong fluctuations, as is evident in the $`M(H)`$ curves in Fig. 2(a). Eventually, FM-ordered ladders with spins aligned along the $`a`$ axis are formed; however, whether a macroscopic magnetic moment will emerge or not depends on the relative orientation of moments between these ladders, which can be only ferro- or antiferromagnetic due to the Ising nature of the spins. Experimentally, the interaction between different FM ladders, mediated by spinless-Co-O layers \[Fig. 3(b)\], turns out to be antiferromagnetic, bringing about the AF ground state. Of course, our magnetization data cannot tell along exactly which axis, $`b`$ or $`c`$, or both, the ladder stacking is AF, but this does not matter for the qualitative picture.

Quite naturally, the inter-ladder coupling is weak, and switching from the ground-state AF order to the FM one can be induced by magnetic fields or temperature. To understand the role of temperature, one should consider the thermally-excited states: upon heating, a certain amount of LS Co3+ ions become Co2+, Co4+, or change their spin state. Whatever the case, each excited ion acquires a non-zero spin, and replaces the weak SE interaction between two IS Co-ions on the neighboring ladders with usual SE interactions, providing a strong bridge between spin-ordered ladders (Fig. 3(c)). Regardless of whether the excited spin couples ferro- or antiferromagnetically with the spins in ladders, the symmetry of the bridge makes it certain that the additional coupling between the ladders is FM. Therefore, thermally-excited spins should inevitably induce an AF-FM transition at some temperature, if fluctuations will not kill the FM order in the ladders first. It is worth noting also that the AF ordering can be suppressed not only by magnetic fields or thermal excitations, but by changing the oxygen stoichiometry as well: any deviation from $`x=0.5`$ also introduces Co2+ or Co4+ ions and should shift the AF-FM phase to lower temperature, which we have indeed observed experimentally unpubl .

To complete the picture of magnetic ordering in GdBaCo2O5.5, it is useful to estimate the scale of different magnetic interactions. The strongest one is definitely the FM superexchange between IS Co3+-ions within ladders, which can be estimated from the Curie temperature $`T_c`$: The molecular-field theory for 3D systems would give $`J/k_\text{B}150`$ K for $`T_c=300`$ K; however, for the quasi-1D/2D magnetic ordering, $`J`$ should be roughly twice as large Fisher . The weak AF coupling between adjacent FM ladders, $`J_w`$, can be directly evaluated from the critical field of the AF-FM transition: $`J_w/k_\text{B}`$ reaches $`15`$ K at $`T=0`$, where $`H_c20`$ T (Fig. 2(c)). In fact, it is this weak AF coupling that makes possible the competition and easy switching between the AF and FM orders in GdBaCo2O5.5. Lastly, based on the magnetization anisotropy in the FM or AF state, we can estimate the spin anisotropy energy, which appears to be unusually large: The rotation of a spin from the spin-easy $`a`$ axis to the $`b`$ axis requires $`\mathrm{\Delta }E_{ab}/k_\text{B}1520`$ K at $`T=260`$ K, which grows up to $`80100`$ K at $`T=2`$ K; the energy necessary to turn a spin along the $`c`$ axis seems to be at least two times larger. Thus, the spin anisotropy energy appears to be about ten times larger than $`J_w`$, which explains the observed Ising-like behavior: the magnetic fields capable of inducing the AF-FM transition are still too weak to rotate the spins away from the $`a`$ axis.

One might wonder how the reorientation of weakly coupled ladders affects the charge transport, bringing about a giant magnetoresistance. The MR mechanism may be quite simple: GdBaCo2O5.5 appears to be a narrow-gap insulator, where the carrier generation goes through formation of Co2+ \- Co4+ pairs, and the excitation energy for these states may well depend on the magnetic order. Indeed, in a low-spin Co-O layer ($`ac`$), each thermally-excited state strongly couples with two adjacent spin-ordered ladders, providing them with a FM bridge; correspondingly, if the ladders’ moments are AF oriented, the resulting frustration should increase the energy of the excited state by $`2J`$. Thus, the relative ordering of adjacent FM ladders is capable of significantly changing the insulating-gap size. Applied magnetic fields align the FM ladders and reduce the insulating gap, which results in step-like increase in the number of carriers and decrease in resistivity \[Fig. 4(a)\]. Figure 4(b) demonstrates that the activation energy is actually diminished by the magnetic field, and thus the MR grows roughly exponentially upon decreasing temperature.

Two salient points should be emphasized in the MR behavior. The first one is the cooperative nature of the MR, whereby an apparently small energy of the magnetic field $`g\mu _\text{B}H/k_\text{B}`$, $``$ several K, is capable of changing the carriers’ activation energy by several hundreds K. The second point is the extremely large MR anisotropy: the 14-T field $`Hc`$ can do nothing comparable to the MR caused by $`Ha`$ \[Fig. 4(c)\], which confirms the remarkable Ising-like spin anisotropy, which could be overcome only by magnetic fields in the 100-T range.

The revealed Ising-like spin behavior together with the absence of significant structural disorder turn RBaCo2O5.50 into a model system for studying the competing magnetic interactions and accompanying MR phenomena. Moreover, RBaCo2O5.50 compounds may represent an intriguing realization of a new MR scheme – a kind of “magnetic field-effect transistor” – where the charge-carrier injection into a 2D semiconducting channel is controlled by a magnetic state of neighboring “ligands”.

###### Acknowledgements.

We thank S. Komiya and K. Segawa for technical assistance and I. Tsukada for fruitful discussions. A.A.T. gratefully acknowledges support from JSPS.


text:
# Gauge-Invariant Initial Conditions and Early Time Perturbations in Quintessence Universes

## I Introduction

The advent of high precision data Spergel:2003cb  of the cosmic microwave background (CMB) anisotropies permits detailed tests of the composition and shape of the primordial density fluctuations. The most popular models of inflationary cosmology predict adiabatic fluctuations Mukhanov:xt ; Bardeen:qw ; Liddle:cg . More elaborate models lead to an admixture of adiabatic and isocurvature fluctuations Lyth:2002my ; Bartolo:2002vf . The time evolution of adiabatic and non-adiabatic fluctuations is well understood for a universe composed of radiation, baryons, cold dark matter (CDM) and neutrinos Ma:1995ey . In the context of quintessence Wetterich:fm ; Ratra:1987rm ; Caldwell:1997ii , the behaviour of the field fluctuation has been studied in several works Viana:1997mt ; dave\_phd ; Malquarti:2002iu ; Abramo:2001mv ; Kawasaki:2001nx . Initial conditions have been proposed in Perrotta:1998vf  for the case of negligible quintessence contribution in the early universe. We present here a systematic treatment of initial conditions for quintessence models which differs from that of Perrotta:1998vf  in approach and interpretation.

Our basic setting assumes that small deviations from homogeneity are generated during a very early stage of the big bang, typically an inflationary epoch. During the following radiation dominated period the wavelength of the relevant fluctuations is far outside the horizon. Apart from this, we will not use any further constraint on the primordial fluctuations. Only the spectra of a certain number of “dominant” modes can possibly influence events such as emission of the CMB and its anisotropies since the other modes decay. The information about these dominant modes therefore constitutes the initial conditions for practical purposes. Primordial information beyond the dominant modes is effectively lost and not observable. The detailed time of specification of the initial conditions is therefore irrelevant as long as it is much shorter than the time of matter-radiation equality.

During the period relevant for the discussion of the initial conditions the universe is radiation dominated. However, our approach allows for the presence of scalar fields which evolve like radiation at early times or are subdominant. Consequently, our results hold for a wide class of quintessence models, including those with non-negligible $`\mathrm{\Omega }_q`$ at early times Caldwell:2003vp . In fact, we only use a “tracking” property Steinhardt:nw  for the background of homogenous quintessence, namely that its equation of state $`w_q=p_q/\rho _q`$ is almost constant and determined only by the energy densities of the radiation and matter components. The parameters $`w_q`$ and $`\mathrm{\Omega }_q=1\mathrm{\Omega }_m\mathrm{\Omega }_\nu \mathrm{\Omega }_\gamma `$ will therefore be the only parameters of the quintessence model that influence the early time evolution of small fluctuations. This makes our analysis model independent to a large extent.

We will formulate the evolution equations for the perturbation variables as a first order differential matrix equation:

$$\frac{\text{d}}{\text{d}\mathrm{ln}x}𝑼=A(x)𝑼,$$
(1)
where the vector $`𝑼`$ contains all perturbation variables and the matrix $`A(x)`$ encodes the evolution equations. In doing so, we relate the problem of finding initial conditions and dominant modes to the familiar language of eigenvalues and eigenvectors. This formulation makes “mode-accounting” transparent by counting the degeneracy of the largest eigenvalue. We find four dominant modes that remain regular at early times. For physical reasons, we choose a basis using adiabatic, CDM isocurvature, baryon isocurvature and neutrino isocurvature initial conditions. As we will show, adiabaticity between CDM, baryons and photons implies adiabaticity of quintessence. There is therefore no pure quintessence isocurvature mode. In addition, using the matrix formulation reveals facets of the modes that otherwise remain obscured.

In order to avoid the appearance of gauge modes, we will use the gauge-invariant formalism Bardeen:kt ; Kodama:bj ; Mukhanov:1990me ; Durrer:2001gq . In contrast to earlier work, we find it more appropriate to specify the initial conditions and time evolution of the quintessence field in terms of the gauge-invariant density contrast and velocity, thus unifying the language for all species. As anticipated, the quintessence density perturbation remains constant at super-horizon scales for adiabatic initial conditions. In contrast to this, the field fluctuation follows a simple power law in conformal time that only depends on the quintessence equation of state.

We will proceed as follows: in section II we give the gauge-invariant perturbation equations for a radiation-dominated universe containing radiation, cold dark matter, neutrinos, baryons in the tight coupling limit and tracking quintessence. We express the evolution in matrix form in II.2. In section III.1, we classify the modes and determine them in sections III.2, III.3 and III.4. To illustrate the effect of non-adiabatic contributions to the CMB spectrum, we plot a few spectra for different initial conditions in section IV. A summary of our findings is given in section V. In Appendix A, we derive the perturbation equations used in detail, while Appendices B and C discuss supplementary issues.

## II The Perturbation Equations

In the following we adopt the gauge-invariant approach as devised by Bardeen Bardeen:kt . It is not difficult to obtain the initial conditions in any gauge from the corresponding gauge-invariant quantities given here. In Appendix A, we summarize the definitions of the perturbation variables and sketch the derivation of the evolution equations. It turns out that the evolution is best described as a function of $`xk\tau `$, where $`\tau `$ is the conformal time and $`k`$ the comoving wavenumber of the mode. We assume that at early times, the universe expands as if radiation dominated. This assumption is well justified for small $`\mathrm{\Omega }_q`$ at early times, as well as for potentials that are essentially exponentials at the time of interest, regardless of $`\mathrm{\Omega }_q`$. The assumption is certainly not justified for models in which quintessence is dominating the universe at early times with equation of state $`w_q1/3`$. For such (slightly exotic) models, the following steps would need to be modified.

### II.1 Full Set of Equations

Assuming tracking quintessence we obtain the following set of equations (for a derivation, see Appendix A):

$`\mathrm{\Delta }_c^{}`$ $`=`$ $`x^2\stackrel{~}{V}_c,`$ (2)
$`\stackrel{~}{V}_c^{}`$ $`=`$ $`2\stackrel{~}{V}_c+\mathrm{\Psi },`$ (3)
$`\mathrm{\Delta }_\gamma ^{}`$ $`=`$ $`{\displaystyle \frac{4}{3}}x^2\stackrel{~}{V}_\gamma ,`$ (4)
$`\stackrel{~}{V}_\gamma ^{}`$ $`=`$ $`{\displaystyle \frac{1}{4}}\mathrm{\Delta }_\gamma \stackrel{~}{V}_\gamma +\mathrm{\Omega }_\nu \stackrel{~}{\mathrm{\Pi }}_\nu +2\mathrm{\Psi },`$ (5)
$`\mathrm{\Delta }_b^{}`$ $`=`$ $`x^2\stackrel{~}{V}_\gamma ,`$ (6)
$`\mathrm{\Delta }_\nu ^{}`$ $`=`$ $`{\displaystyle \frac{4}{3}}x^2\stackrel{~}{V}_\nu ,`$ (7)
$`\stackrel{~}{V}_\nu ^{}`$ $`=`$ $`{\displaystyle \frac{1}{4}}\mathrm{\Delta }_\nu \stackrel{~}{V}_\nu {\displaystyle \frac{1}{6}}x^2\stackrel{~}{\mathrm{\Pi }}_\nu +\mathrm{\Omega }_\nu \stackrel{~}{\mathrm{\Pi }}_\nu +2\mathrm{\Psi },`$ (8)
$`\stackrel{~}{\mathrm{\Pi }}_\nu ^{}`$ $`=`$ $`{\displaystyle \frac{8}{5}}\stackrel{~}{V}_\nu 2\stackrel{~}{\mathrm{\Pi }}_\nu ,`$ (9)
$`\mathrm{\Delta }_q^{}`$ $`=`$ $`3(w_q1)[\mathrm{\Delta }_q+3(1+w_q)\{\mathrm{\Psi }+\mathrm{\Omega }_\nu \stackrel{~}{\mathrm{\Pi }}_\nu \}`$ (10)
$`+\{3{\displaystyle \frac{x^2}{3(w_q1)}}\}(1+w_q)\stackrel{~}{V}_q],`$
$`\stackrel{~}{V}_q^{}`$ $`=`$ $`3\mathrm{\Omega }_\nu \stackrel{~}{\mathrm{\Pi }}_\nu +{\displaystyle \frac{\mathrm{\Delta }_q}{1+w_q}}+\stackrel{~}{V}_q+4\mathrm{\Psi },`$ (11)
with the gauge-invariant Newtonian potential $`\mathrm{\Psi }`$ given by

$$\mathrm{\Psi }=\frac{\underset{\alpha =c,b,\gamma ,\nu ,q}{}\mathrm{\Omega }_\alpha (\mathrm{\Delta }_\alpha +3(1+w_\alpha )\stackrel{~}{V}_\alpha )}{\underset{\alpha =c,b,\gamma ,\nu ,q}{}3(1+w_\alpha )\mathrm{\Omega }_\alpha +\frac{2x^2}{3}}\mathrm{\Omega }_\nu \stackrel{~}{\mathrm{\Pi }}_\nu .$$
(12)
We denote the derivative $`\text{d}/\text{d}\mathrm{ln}x`$ with a prime. The gauge-invariant energy density contrasts $`\mathrm{\Delta }_\alpha `$, the velocities $`\stackrel{~}{V}_\alpha `$ and the shear $`\stackrel{~}{\mathrm{\Pi }}_\nu `$ are the ones found in the literature Bardeen:kt ; Kodama:bj ; Durrer:2001gq , except that we factor out powers of $`x`$ from the velocity and shear defining $`\stackrel{~}{V}V/x`$ and $`\stackrel{~}{\mathrm{\Pi }}_\nu x^2\mathrm{\Pi }_\nu `$. This factoring out leads to the particularly simple form of the system of equations for $`x1`$ (see also Appendix A). It does, however, exclude modes with diverging $`\mathrm{\Psi }`$ at early times such as a neutrino velocity mode Rebhan:1994zw . The index $`\alpha `$ runs over the five species in our equations, namely cold dark matter, baryons, photons, neutrinos and quintessence, denoted with the subscript $`q`$. We assume tight coupling between photons and baryons. The equation of state $`w=\overline{p}/\overline{\rho }`$ takes on the values $`w_c=w_b=0`$, $`w_\gamma =w_\nu =1/3`$ and $`w_q`$ is left as a free parameter. Equations (2), (4), (6) and (7) can be regarded as continuity relations between the density fluctuations and the velocity. We obtain equations (10) and (11) from the perturbed Klein-Gordon equation of the quintessence scalar field expressed in terms of $`\mathrm{\Delta }_q`$ and $`V_q`$, the energy density and velocity perturbations as defined in Appendix A.

### II.2 Matrix Formulation and Dominant Modes

Conceptually, it is convenient to note that the above set of equations can be concisely written in matrix form according to Equation (1) where the perturbation vector is defined as

$$𝑼^T(\mathrm{\Delta }_c,\stackrel{~}{V}_c,\mathrm{\Delta }_\gamma ,\stackrel{~}{V}_\gamma ,\mathrm{\Delta }_b,\mathrm{\Delta }_\nu ,\stackrel{~}{V}_\nu ,\stackrel{~}{\mathrm{\Pi }}_\nu ,\mathrm{\Delta }_q,\stackrel{~}{V}_q).$$
(13)
The matrix $`A(x)`$ can easily be read off from equations (2)-(11). This enables us to discuss the problem of specifying initial conditions in a systematic way.

The initial conditions are specified for modes well outside the horizon, i.e. $`x1`$. In this case, the r.h.s. of equations (2), (4), (6) and (7) can be neglected, provided $`\stackrel{~}{V}_\alpha `$ does not diverge $`x^2`$ or faster for $`x^20`$. The evolution matrix $`A(x)`$ loses any explicit $`x`$ dependence for $`x^20`$. Yet, it still depends on $`x`$ via terms involving $`\mathrm{\Omega }_c,\mathrm{\Omega }_b`$ and $`\mathrm{\Omega }_q`$. By our assumptions on quintessence, the term involving $`\mathrm{\Omega }_q`$ is either a constant (for $`w_q=1/3`$) or negligible (yet, in Appendix C, we extend the treatment to include models with considerable $`\mathrm{\Omega }_q`$ and $`w_q1/3`$.) In both cases $`\mathrm{\Omega }_q`$ can be approximated by a constant ($`\mathrm{\Omega }_q=0`$ for $`w_q<1/3`$) and $`\mathrm{\Omega }_c`$, $`\mathrm{\Omega }_b`$ vanish $`x`$. In leading order, the matrix $`A`$ becomes therefore $`x`$-independent for very early times. In fact, the general solution to Equation (1) in the (ideal) case of a truly constant $`A`$ would be

$$𝑼(x)=\underset{i}{}c_i\left(\frac{x}{x_0}\right)^{\lambda _i}𝑼^{(i)},$$
(14)
where $`𝑼^{(i)}`$ are the eigenvectors of $`A`$ with eigenvalue $`\lambda _i`$ and the time independent coefficients $`c_i`$ specify the initial contribution of $`𝑼^{(i)}`$ towards a general perturbation $`𝑼`$. As time progresses, components corresponding to the largest eigenvalues $`\lambda _i`$ will dominate. Compared to these “dominant” modes, initial contributions in the direction of eigenvectors $`𝑼^{(i)}`$ with smaller $`\mathrm{Re}(\lambda _i)`$ decay. It therefore suffices to specify the initial contribution $`c_i`$ for the dominant modes, if one is not interested in very early time behaviour shortly after inflation. In our case, the characteristic polynomial of $`A(x)`$ indeed has a fourfold degenerate eigenvalue $`\lambda =0`$ in the limit $`x^20`$, independent of $`\mathrm{\Omega }_c`$, $`\mathrm{\Omega }_b`$ and $`\mathrm{\Omega }_q`$.111For $`w_q=1`$ we find another eigenvalue with $`\lambda =0`$. We will ignore this special case in what follows. While it is not feasible to obtain the remaining six eigenvalues by analytic means, we have checked numerically for a wide range of $`\mathrm{\Omega }_\gamma `$, $`\mathrm{\Omega }_\nu `$, $`\mathrm{\Omega }_b`$, $`\mathrm{\Omega }_c`$, $`\mathrm{\Omega }_q`$ and $`w_q`$ that the remaining eigenvalues have indeed negative real parts and contributions from the corresponding eigenvectors towards a general perturbation $`𝑼`$ will therefore decay according to Equation (14). We can improve the analytic description of the dominant modes by taking corrections $`x`$ into account.

As $`\mathrm{\Omega }_c\mathrm{\Omega }_bx`$, it is appropriate to split $`A(x)`$ according to the scaling with $`x`$,

$$A=A_0+xA_1,$$
(15)
where $`A_0`$ and $`A_1`$ are constant and $`xA_1`$ contains the small, time-dependent corrections from terms involving $`\mathrm{\Omega }_c`$ and $`\mathrm{\Omega }_b`$. We may also write222This form is not an ansatz, but dictated by Equation (1), once the dependence of $`A(x)`$ on $`x`$ is given. the eigenvectors as a series in $`x`$,

$$𝑼=𝑼_0+x𝑼_1.$$
(16)
Inserting Equations (15)-(16) in Equation (1), we get

$$A_0𝑼_0=0,$$
(17)
and

$$𝑼_1=(A_0𝟙)^\mathrm{𝟙}𝔸_\mathrm{𝟙}𝕌_\mathrm{𝟘}.$$
(18)
Equation (18) is easy to solve, once $`𝑼_0`$ has been determined (we discuss the possibility of a vanishing $`𝑼_0`$ in Appendix C). We see from eqn. (17) that to constant order the solutions of eqn. (1) are indeed given by eigenvectors to the eigenvalue $`\lambda =0`$. We should emphasize that the vectors $`𝐔_0`$ do not evolve in time if their corresponding eigenvalues are $`\lambda =0`$. Thus, the perturbations remain constant in the super-horizon regime during radiation domination in this approximation. If we include the next-to-leading order contribution to $`𝐔`$, the eigenvectors do evolve and we can no longer apply eq. (14). These corrections are, however, small as long as we are deep in the radiation dominated era due to the small contributions of baryons, radiation and quintessence during this era. Given a set of initial conditions in the form of coefficients for the four dominating modes at $`z_{initial}`$ we can find the perturbations at some later time (provided the modes are still super-horizon sized and we have radiation domination). In leading order, the coefficients will remain the same while in next-to-leading order we can use the evolution of $`𝐔`$ to compute the coefficients for $`z<z_{initial}`$. If initial conditions are specified with accuracy of next-to-leading order one therefore has to specify $`z_{initial}`$ as well. In leading order this is unecessary for $`z`$ in a wide range long before last scattering.

### II.3 Constraint Equations to Leading Order

Equation (17) is equivalent to setting the l.h.s. of Equations (2)-(11) equal to zero and using $`\mathrm{\Omega }_c=\mathrm{\Omega }_b=x^2=0`$. Then Equations (2), (4), (6) and (7) are automatically satisfied (provided $`\stackrel{~}{V}_\alpha `$ does not diverge $`x^2`$ or faster), and Equations (3),(5),(8)-(11) yield non-trivial constraints for the components of $`𝑼_0`$:

$`2\stackrel{~}{V}_c\mathrm{\Psi }`$ $`=`$ $`0,`$ (19)
$`1/4\mathrm{\Delta }_\gamma \stackrel{~}{V}_\gamma +\mathrm{\Omega }_\nu \stackrel{~}{\mathrm{\Pi }}_\nu +2\mathrm{\Psi }`$ $`=`$ $`0,`$ (20)
$`1/4\mathrm{\Delta }_\nu \stackrel{~}{V}_\nu +\mathrm{\Omega }_\nu \stackrel{~}{\mathrm{\Pi }}_\nu +2\mathrm{\Psi }`$ $`=`$ $`0,`$ (21)
$`8/5\stackrel{~}{V}_\nu 2\stackrel{~}{\mathrm{\Pi }}_\nu `$ $`=`$ $`0,`$ (22)
$`3\mathrm{\Omega }_\nu \stackrel{~}{\mathrm{\Pi }}_\nu +\mathrm{\Delta }_q/(1+w_q)+3\stackrel{~}{V}_q+3\mathrm{\Psi }`$ $`=`$ $`0,`$ (23)
$`3\mathrm{\Omega }_\nu \stackrel{~}{\mathrm{\Pi }}_\nu +\mathrm{\Delta }_q/(1+w_q)+\stackrel{~}{V}_q+4\mathrm{\Psi }`$ $`=`$ $`0.`$ (24)
In the above, all quantities are considered only to constant order. (we have omitted the subscript ’$`0`$’ for notational convenience.) In particular, there is no contribution of CDM and baryons to $`\mathrm{\Psi }`$ at constant order. Note that, apart from $`w_q`$, no model-specific parameters occur in any of these equations so the modes will be independent of the type of quintessence as long as the scalar field is in a regime with approximately constant $`w_q`$. We note that for $`w_q`$ substantially smaller than $`1/3`$ the quintessence fraction $`\mathrm{\Omega }_q`$ changes with time. By the assumption that the universe expands as if radiation dominated, the quintessence contribution would however be small in this case and its contribution to $`\mathrm{\Psi }`$ can be neglected (see Appendix C for an extended discussion).

We mention that for $`w_q=1/3`$, quintessence evolves the same way as radiation, therefore $`\mathrm{\Omega }_q`$ does not change in this case. If $`w_q=1/3`$, quintessence has the same influence on the scale factor $`a`$ as a curvature term in an open universe. However, the geometry is still flat and one can distinguish an open universe from this quintessence model by measuring the position of the first acoustic peak in the CMB.

## III The modes in detail

### III.1 Classifying the modes

While any basis for the subspace spanned by the eigenvectors with eigenvalue $`\lambda =0`$ can be used to specify the initial conditions, it is still worthwhile to use a basis that is physically meaningful. Following the existing literature, we use the gauge-invariant entropy perturbation Kodama:bj

$$S_{\alpha :\beta }=\frac{\mathrm{\Delta }_\alpha }{1+w_\alpha }\frac{\mathrm{\Delta }_\beta }{1+w_\beta },$$
(25)
between two species $`\alpha `$ and $`\beta `$, as well as the gauge-invariant curvature perturbation on hyper-surfaces of uniform energy density of species $`\alpha `$ Bardeen:qw ; bardeen\_isocurvature\_book ; Lyth:2002my ; Wands:2000dp

$$\zeta _\alpha =\left(H_L+\frac{1}{3}H_T\right)+\frac{\delta \rho _\alpha }{3(1+w_\alpha )\overline{\rho }_\alpha },$$
(26)
in order to classify the physical modes. On slices of uniform total energy density, the curvature perturbation is correspondingly

$$\zeta _{tot}=\left(H_L+\frac{1}{3}H_T\right)+\frac{_\alpha \delta \rho _\alpha }{_\alpha 3(1+w_\alpha )\overline{\rho }_\alpha }.$$
(27)
In our variables, these expressions take on the manifestly gauge-invariant form

$$\zeta _\alpha =\frac{\mathrm{\Delta }_\alpha }{3(1+w_\alpha )},\zeta _{tot}=\frac{_\alpha \mathrm{\Delta }_\alpha \mathrm{\Omega }_\alpha }{_\alpha 3(1+w_\alpha )\mathrm{\Omega }_\alpha }.$$
(28)
If $`\zeta _{tot}=0`$, energy density perturbations do not generate curvature. It is therefore clear that such a perturbation is a perturbation in the local equation of state. One should note that the definition of $`\zeta _{tot}`$ is different from that of Mukhanov:1990me :

$$\zeta _{MFB}=\frac{2}{3}\frac{^1\dot{\mathrm{\Psi }}+\mathrm{\Psi }}{(1+w)}+\mathrm{\Psi }.$$
(29)
However, one may verify that this quantity coincides with $`\zeta _{tot}`$ in the super-horizon limit for a flat universe martin .

### III.2 The Adiabatic Mode

The first (rather intuitive) perturbations one would try to find are adiabatic perturbations, which are specified by the adiabaticity conditions $`S_{\alpha :\beta }=0`$ for all pairs of components. In our case, this results in eleven constraints333Without requiring quintessence to be adiabatic, we have six constraints from equations (19)-(24) plus three constraints from eq. (30) plus one constraint from the overall normalization, which is fixed by choosing a specific value for $`\mathrm{\Delta }_\gamma `$. for the ten components of $`𝑼_0`$. It is a priori not clear that this has a solution so we will not include quintessence in the adiabaticity requirement. Requiring adiabaticity between CDM, baryons, neutrinos and radiation,

$$\mathrm{\Delta }_\nu =\mathrm{\Delta }_\gamma =\frac{4}{3}\mathrm{\Delta }_c=\frac{4}{3}\mathrm{\Delta }_b,$$
(30)
and using the six constraint Equations (19)-(24), we obtain

$$\left(\begin{array}{c}\mathrm{\Delta }_c\\ \stackrel{~}{V}_c\\ \mathrm{\Delta }_\gamma \\ \stackrel{~}{V}_\gamma \\ \mathrm{\Delta }_b\\ \mathrm{\Delta }_\nu \\ \stackrel{~}{V}_\nu \\ \stackrel{~}{\mathrm{\Pi }}_\nu \\ \mathrm{\Delta }_q\\ \stackrel{~}{V}_q\end{array}\right)_{\text{adiabatic}}=C\left(\begin{array}{c}3/4\\ (5/4)𝒫\\ 1\\ (5/4)𝒫\\ 3/4\\ 1\\ (5/4)𝒫\\ 𝒫\\ 3(1+w_q)/4\\ (5/4)𝒫\end{array}\right),$$
(31)
where $`𝒫=(15+4\mathrm{\Omega }_\nu )^1`$ and $`C`$ is an arbitrary constant. From $`\mathrm{\Delta }_q/\mathrm{\Delta }_\gamma =3(1+w_q)/4`$ we conclude that quintessence is automatically adiabatic if CDM, baryons, neutrinos and radiation are adiabatic, independent of the quintessence model for as long as we are in the tracking regime. As all components are non-vanishing, we do not quote the next to leading order contributions from $`x𝑼_1`$.

### III.3 Neutrino Isocurvature

Having found the adiabatic vector, one could specify three additional linearly independent vectors satisfying the constraint Equations (19)-(24). This would complete the basis. It is, however, appropriate to choose modes that may be generated by physical processes. These modes are in general not orthogonal but span the eigenspace of $`\lambda =0`$. Modes that may be generated by physical processes are isocurvature modes. A given mode is an isocurvature mode, if the gauge-invariant curvature perturbation $`\zeta _{tot}`$ vanishes, i.e. $`\zeta _{tot}=0`$. In order to distinguish different isocurvature modes from one another, we require that the other species are adiabatic with respect to each other, i. e. $`S_{\alpha :\beta }=0`$ except for quintessence and one species $`\sigma `$, which has non-vanishing $`S_{\sigma :\gamma }`$.

Let us first consider the neutrino isocurvature mode. For this, we require that CDM, baryons and radiation are adiabatic, while $`S_{\nu :\gamma }0`$ and that the gauge-invariant curvature perturbation vanishes:

$$\zeta _{tot}=0,\mathrm{\Delta }_c=\mathrm{\Delta }_b=\frac{3}{4}\mathrm{\Delta }_\gamma .$$
(32)
Using this and Equations (19)-(24) leads to

$$\left(\begin{array}{c}\mathrm{\Delta }_c\\ \stackrel{~}{V}_c\\ \mathrm{\Delta }_\gamma \\ \stackrel{~}{V}_\gamma \\ \mathrm{\Delta }_b\\ \mathrm{\Delta }_\nu \\ \stackrel{~}{V}_\nu \\ \stackrel{~}{\mathrm{\Pi }}_\nu \\ \mathrm{\Delta }_q\\ \stackrel{~}{V}_q\end{array}\right)_{\text{neutrino iso.}}=C\left(\begin{array}{c}3/4\\ \mathrm{\Omega }_\gamma 𝒫\\ 1\\ \left(\mathrm{\Omega }_\gamma +\mathrm{\Omega }_\nu +\frac{15}{4}\right)𝒫\\ 3/4\\ \mathrm{\Omega }_\gamma /\mathrm{\Omega }_\nu \\ \frac{15}{4}𝒫\mathrm{\Omega }_\gamma /\mathrm{\Omega }_\nu \\ 3𝒫\mathrm{\Omega }_\gamma /\mathrm{\Omega }_\nu \\ 0\\ \mathrm{\Omega }_\gamma 𝒫\end{array}\right).$$
(33)
It is important to note that we did not require quintessence to be adiabatic. One can see from the neutrino isocurvature vector that $`\mathrm{\Delta }_q=0`$, and as a consequence quintessence is not adiabatic with respect to either neutrinos, radiation, baryons or CDM. Hence, we could just as well have labeled this vector “quintessence isocurvature”. We cannot require adiabaticity between neutrinos, CDM, baryons and radiation and hope to obtain a “pure” quintessence isocurvature vector since, as we have seen in the discussion of the adiabatic mode, these requirements lead to quintessence being adiabatic as well.

### III.4 CDM isocurvature and baryon isocurvature

The CDM isocurvature mode is characterized by $`S_{c:\gamma }0`$, $`\zeta _{tot}=0`$ and adiabaticity between photons, neutrinos and baryons:

$$\zeta _{tot}=0,\mathrm{\Delta }_\gamma =\mathrm{\Delta }_\nu =\frac{4}{3}\mathrm{\Delta }_b.$$
(34)
Using this and Equations (19)-(24) yields

$$𝑼_0^T(\text{CDM iso.})=(1,0,0,0,0,0,0,0,0,0).$$
(35)
This vector fulfills $`\zeta _{tot}=0+𝒪(\mathrm{\Omega }_c)`$, which is in line with our approximation since $`\mathrm{\Omega }_c1`$. Similarly, for the baryon isocurvature mode, we require $`S_{b:\gamma }0`$, $`\zeta _{tot}=0`$ and adiabicity between photons, neutrinos and baryons. The resulting vector reads

$$𝑼_0^T(\text{baryon iso.})=(0,0,0,0,1,0,0,0,0,0).$$
(36)
As all but one of the components of $`𝑼_0`$ are vanishing for CDM isocurvature and baryon isocurvature, we use Equation (18) to obtain the next to constant order solution for CDM isocurvature

$$\left(\begin{array}{c}\mathrm{\Delta }_c\\ \stackrel{~}{V}_c\\ \mathrm{\Delta }_\gamma \\ \stackrel{~}{V}_\gamma \\ \mathrm{\Delta }_b\\ \mathrm{\Delta }_\nu \\ \stackrel{~}{V}_\nu \\ \stackrel{~}{\mathrm{\Pi }}_\nu \\ \mathrm{\Delta }_q\\ \stackrel{~}{V}_q\end{array}\right)_{\text{CDM iso.}}=C\left(\begin{array}{c}1\\ \mathrm{\Omega }_c(4\mathrm{\Omega }_\nu 15)𝒰/12\\ 0\\ (15/4)\mathrm{\Omega }_c𝒰\\ 0\\ 0\\ (15/4)\mathrm{\Omega }_c𝒰\\ 2\mathrm{\Omega }_c𝒰\\ \mathrm{\Omega }_c(15+2\mathrm{\Omega }_\nu )(1+w_q)𝒰\\ \mathrm{\Omega }_c𝒰𝒱\end{array}\right),$$
(37)
where $`𝒰=(30+4\mathrm{\Omega }_\nu )^1`$ and $`𝒱=[10545w_q+4\mathrm{\Omega }_\nu (3w_q1)]/[36(w_q1)]`$. Similarly, we find for baryon isocurvature

$$\left(\begin{array}{c}\mathrm{\Delta }_c\\ \stackrel{~}{V}_c\\ \mathrm{\Delta }_\gamma \\ \stackrel{~}{V}_\gamma \\ \mathrm{\Delta }_b\\ \mathrm{\Delta }_\nu \\ \stackrel{~}{V}_\nu \\ \stackrel{~}{\mathrm{\Pi }}_\nu \\ \mathrm{\Delta }_q\\ \stackrel{~}{V}_q\end{array}\right)_{\text{baryon iso.}}=C\left(\begin{array}{c}0\\ \mathrm{\Omega }_b(4\mathrm{\Omega }_\nu 15)𝒰/12\\ 0\\ (15/4)\mathrm{\Omega }_b𝒰\\ 1\\ 0\\ (15/4)\mathrm{\Omega }_b𝒰\\ 2\mathrm{\Omega }_b𝒰\\ \mathrm{\Omega }_b(15+2\mathrm{\Omega }_\nu )(1+w_q)𝒰\\ \mathrm{\Omega }_b𝒰𝒱\end{array}\right).$$
(38)
Note that these vectors are not constant since $`\mathrm{\Omega }_b`$ and $`\mathrm{\Omega }_c`$ both evolve in time. We observe that the corrections to $`U`$ are indeed proportional to $`\mathrm{\Omega }_c`$ or $`\mathrm{\Omega }_b`$ as expected. This result holds for all tracking quintessence models with $`w_q=1/3`$ or $`w_q0`$ during the radiation dominated period. For intermediate values $`0<w_q<1/3`$ the deviation from the leading behaviour scales $`x^\alpha `$, $`\alpha <1`$. Obviously, the adiabatic, CDM isocurvature, baryon isocurvature and neutrino isocurvature vectors $`𝑼_0`$ are linearly independent. We have therefore identified four modes corresponding to the fourfold degenerate eigenvalue zero of $`A(x)`$. These four vectors span the subspace of dominant modes in the super-horizon limit, and there are no more linearly independent vectors that satisfy the constraints (19) \- (24). Arbitrary initial perturbations may therefore be represented by projecting a perturbation vector $`𝑼`$ at initial time into the subspace spanned by the four aforementioned vectors, as this is the part of the initial perturbations which will dominate as time progresses.

Figure 1 demonstrates that the early time behaviour is well described by our analytic formulae. The analytic results agree very well with the simulation for early times, when the mode is outside the horizon. In the lower graph, we plot the equation of state $`w_q`$. The quintessence model used is parameterized by an equation of state $`w_q(a)=0.95+0.75(1a)`$, leading to $`w_q(\mathrm{early})=0.2`$ and according to (98), $`X\tau ^{0.8}`$. This differs from reference Perrotta:1998vf . 444In Perrotta:1998vf  it is stated that the quintessence fluctuation in Newtonian gauge scales $`\tau ^2`$ for adiabatic initial conditions. This does not agree with our results in Appendix B. Actually, equation (101) of Perrotta:1998vf  includes a factor $`\phi _{t0}`$, which, interpreted as a dynamical quantity $`\text{d}\phi /\text{d}t`$ (and not fixed at some initial time $`t_0`$), leads to a power law in $`\tau `$ which is then consistent with our result of Appendix B.

We see that including quintessence does not add a new dominant mode. The two additional modes added by the fluctuations of the scalar field are both subdominant and decay with negative eigenvalue $`\lambda _i`$. This is due to the fact that none of the perturbation equations for quintessence equate to zero in the superhorizon limit. This holds for non-tracking quintessence models as well. Let us investigate this in detail. For all the other fluid components, $`\mathrm{\Delta }_a^{}=0`$ in the super-horizon limit, but for quintessence we get from eq. (68) that $`\mathrm{\Delta }_q^{}=3(c_{s(q)}^2w_q)\mathrm{\Delta }_q3w_q\mathrm{\Gamma }_q`$. For tracking quintessence, we obtain from equation (85) that $`c_{s(q)}^2=w_q`$ and we find

$$\mathrm{\Delta }_q^{}=3w_q\mathrm{\Gamma }_q$$
(39)
Since $`\mathrm{\Gamma }_q`$ does not vanish except for $`w_q=1`$ (see eq. (84)), this does not equate to zero. 555Note that $`w_q=0`$ does not lead to $`\mathrm{\Delta }_q^{}=0`$. Hence, due to the non-vanishing entropy perturbation of quintessence there is no additional dominant mode. 666We have not yet investigated the relationship between decaying quintessence modes and the background evolution.

## IV Isocurvature initial conditions and the CMB

We illustrate the influence of different initial conditions on the CMB with an example. For an analysis of experimental data and a possible isocurvature contribution to the CMB we refer the reader to Enqvist:2000hp ; Trotta:2002iz ; Gordon:2002gv . Here, we merely whish to show the qualitative features of the different modes. We use a modified version of cmbeasy Doran:2003sy ; Seljak:1996is  to compute CMB spectra corresponding to different initial conditions for an early quintessence cosmology with parameters as in model A of Caldwell:2003vp . We set the spectral index of the isocurvature modes identical to the spectral index of the pure adiabatic mode, $`n_s=0.99`$. The resulting spectra are plotted in Fig. 2. The spectrum of the pure CDM isocurvature mode decays quickly when going to small scales as has been found in previous works Bucher:2000hy ; Langlois:ar ; Amendola:2001ni . The neutrino isocurvature mode shows prominent peaks at higher multipoles than the adiabatic mode with different peak ratios. For the mixed initial conditions with only small isocurvature contribution, the shape of the curve remains more or less the same. A small admixture of isocurvature fluctuations leads to a decrease of power at larger multipoles if the overall normalization is fixed at $`l=10`$. Comparison with the WMAP data in the same figure shows that non-adiabatic initial perturbations are strongly constrained. Clearly, pure isocurvature initial conditions are inconsistent with CMB observations.

## V Conclusion

We have investigated perturbations in a radiation-dominated universe containing quintessence, CDM, neutrinos, radiation and baryons in the tight coupling limit. The perturbation evolution has been expressed as a differential equation involving a matrix acting on a vector comprised of the perturbation variables. This formulation leads to a systematic determination of the initial conditions. In particular, we find that due to the presence of tracking scalar quintessence no additional dominant mode is introduced. This fact is beautifully transparent in the matrix language. Indeed, contributions of higher order in $`xk\tau `$ towards a perturbation vector $`𝑼`$ can easily be determined by solving a simple matrix equation once the constant part of $`𝑼`$ has been determined.

In total, we find four dominant modes and choose them as adiabatic, CDM isocurvature, baryon isocurvature and neutrino isocurvature. For the neutrino isocurvature mode, quintessence automatically is forced to non-adiabaticity. Hence, we could have as well labeled the neutrino isocurvature mode as quintessence isocurvature. To demonstrate the influence on the cosmic microwave background anisotropy spectrum, we have calculated spectra for all modes. Clearly, non-adiabatic contributions are severely constrained by the data. A detailed study may provide ways to put additional constraints on quintessence models or tell us more about the initial perturbations after inflation.

###### Acknowledgements.

We thank Robert R. Caldwell, Pier Stefano Corasaniti, Karim A. Malik and Roberto Trotta for helpful discussions. M. Doran was supported by NSF grant PHY-0099543 at Dartmouth and PHY-9907949 at the KITP. C.M. Müller and G. Schäfer are supported by GRK grant 216/3-02.

## Appendix A Gauge-invariant perturbation equations

In this Appendix we will explain the derivation of equations (2)-(12) in detail.

### A.1 The general story

First, we briefly summarize the gauge-invariant approach of Bardeen:kt ; Kodama:bj ; Durrer:2001gq . Perturbing a homogenous Friedman universe, one classifies fluctuations according to their transformation properties with respect to the rotation group. In flat spacetime, we may expand the perturbation variables in terms of harmonic functions york . With $`Q_{,i}=Q/x^i`$ one defines

$$Q_i(𝒌,𝒙)k^1Q(𝒌,𝒙)_{,i}$$
(40)
and

$$Q_{ij}(𝒌,𝒙)k^2Q(𝒌,𝒙)_{,ij}+\frac{1}{3}\delta _{ij}Q(𝒌,𝒙),$$
(41)
where the $`Q(𝒌,𝒙)`$ are eigenfunctions of the Laplace-operator, $`^2Q_k(𝒙)=k^2Q_k(𝒙)`$ and in spatially flat universes $`Q=\mathrm{exp}(i𝒌𝒙)`$. As modes with different $`𝒌`$ decouple in linear theory, we will not display the $`𝒌`$-dependence of $`Q`$ in the following. The scalar parts of vector and tensor fields can then be written as

$$B_i=BQ_i$$
(42)
and

$$H_{ij}=H_LQ\delta _{ij}+H_TQ_{ij}.$$
(43)
respectively.

In this work, we are only interested in scalar fluctuations because scalar quintessence will not influence vector or tensor modes. The most general line element for a perturbed Robertson-Walker metric may be written as

$`\text{d}s^2=a(\tau )^2[(1+2A)\text{d}\tau ^22B_i\text{d}\tau \text{d}x^i+`$
$`(\delta _{ij}+2H_{ij})\text{d}x^i\text{d}x^j],`$ (44)
where in the scalar case $`B_i`$ and $`H_{ij}`$ are given by equations (42) and (43). The gauge transformation of a tensor $`T`$ is given by Bardeen:kt ; Kodama:bj ; Mukhanov:1990me ; Durrer:2001gq ; Doran:2003sy

$$\stackrel{~}{T}(x)=T(x)L_ϵ\overline{T},$$
(45)
where $`L_ϵ`$ is the Lie derivative. The transformation vector $`ϵ`$ can be decomposed as

$$\stackrel{~}{\tau }=\tau +T(\tau )Q(𝒙),$$
(46)
$$\stackrel{~}{x}^i=x^i+L(\tau )Q^i(𝒙),$$
(47)
where $`L`$ and $`T`$ are arbitrary functions of $`\tau `$. The transformation properties of the metric perturbations are given by Bardeen:kt ; Doran:2003sy

$`\stackrel{~}{A}`$ $`=`$ $`AT\dot{T},`$ (48)
$`\stackrel{~}{B}`$ $`=`$ $`B+\dot{L}+kT,`$ (49)
$`\stackrel{~}{H}_L`$ $`=`$ $`H_LT{\displaystyle \frac{k}{3}}L,`$ (50)
$`\stackrel{~}{H}_T`$ $`=`$ $`H_T+kL,`$ (51)
where a dot denotes the derivative with respect to conformal time $`\tau `$ and $`\dot{a}(\tau )/a(\tau )`$. The functions $`L`$ and $`T`$ can be used to eliminate two of the metric perturbations. Popular choices are $`A=B=0`$ for the synchronous gauge and $`B=H_T=0`$ for the longitudinal gauge.

From equations (48)-(51) one can construct the gauge-invariant Bardeen potentials Bardeen:kt

$`\mathrm{\Psi }=Ak^1\sigma k^1\dot{\sigma },`$ (52)
$`\mathrm{\Phi }=H_L+{\displaystyle \frac{1}{3}}H_Tk^1\sigma ,`$ (53)
with $`\sigma k^1\dot{H}_TB`$. It is worthwhile to note that in longitudinal gauge, for which $`B=H_T=\sigma =0`$, the perturbed metric takes on the simple form

$$\begin{array}{c}ds_{(long.)}^2=a(\tau )^2[(1+2\mathrm{\Psi }Q)d\tau ^2\hfill \\ \hfill +(1+2\mathrm{\Phi }Q)\delta _{ij}dx^idx^j].\end{array}$$
(54)
With $`\mathrm{M}_{\overline{P}}(8\pi G)^{1/2}`$ denoting the reduced Planck mass, Einstein’s equation reads

$$T_\nu ^\mu =\mathrm{M}_{\overline{P}}^2\left(R_\nu ^\mu \frac{1}{2}\delta _\nu ^\mu R\right),$$
(55)
where the energy momentum tensor of a perfect fluid is given by

$$T_\nu ^\mu =p\delta _\nu ^\mu +(\rho +p)u^\mu u_\nu +\pi _\nu ^\mu .$$
(56)
The covariant 4-velocity is $`u_i=a(v(\tau )B)Q_i`$. We define the energy density contrast $`\delta `$ by $`\rho =\overline{\rho }(1+\delta (\tau ))Q`$, the spatial trace by $`p\delta _j^i=\overline{p}(\tau )(1+\pi _L(\tau )Q)\delta _j^i`$ and the traceless part by $`\pi _j^i=\overline{p}\mathrm{\Pi }Q_j^i`$. Therefore the components of the energy momentum tensor are

$`T_{\mathrm{\hspace{0.33em}\hspace{0.33em}0}}^0`$ $`=`$ $`\overline{\rho }(1+\delta Q),`$ (57)
$`T_{\mathrm{\hspace{0.33em}\hspace{0.33em}0}}^i`$ $`=`$ $`\overline{\rho }(1+w)vQ^i,`$ (58)
$`T_i^0`$ $`=`$ $`\overline{\rho }(1+w)(vB)Q_i,`$ (59)
$`T_j^i`$ $`=`$ $`\overline{p}\left[\left(1+\pi _LQ\right)\delta _j^i+\mathrm{\Pi }Q_j^i\right].`$ (60)
Given the gauge-transformation properties of $`\delta `$, $`v`$ and $`\pi _L`$ Bardeen:kt ; Kodama:bj ; Mukhanov:1990me ; Durrer:2001gq ; Doran:2003sy , one can construct gauge-invariant quantities for the energy density contrast $`\mathrm{\Delta }`$, the velocity $`V`$ and the entropy perturbation $`\mathrm{\Gamma }`$. These are given by

$`\mathrm{\Delta }`$ $`=`$ $`\delta +3(1+w)\left(H_L+{\displaystyle \frac{1}{3}}H_T\right),`$ (61)
$`V`$ $`=`$ $`vk^1\dot{H}_T,`$ (62)
$`\mathrm{\Gamma }`$ $`=`$ $`\pi _L{\displaystyle \frac{c_s^2}{w}}\delta .`$ (63)
Here, $`c_s^2\overline{p}/\overline{\rho }`$ is the adiabatic sound speed. From the conservation of the zero component of the energy momentum tensor $`_\mu \overline{T}_{\mathrm{\hspace{0.33em}\hspace{0.33em}0}}^\mu =0`$ we obtain

$$\frac{\dot{\overline{\rho _\alpha }}}{\overline{\rho _\alpha }}=3(1+w_\alpha ),$$
(64)
where $`w=\overline{p}/\overline{\rho }`$ is the equation of state of the particular species. The perturbed Einstein equations in gauge-invariant variables are Bardeen:kt ; Kodama:bj ; Mukhanov:1990me ; Durrer:2001gq ; Doran:2003sy :

$`a^2\overline{\rho }\mathrm{\Delta }`$ $`=`$ $`2\mathrm{M}_{\overline{P}}^2k^2\mathrm{\Phi }`$ (65)
$`3a^2\overline{\rho }(1+w)(k^1V\mathrm{\Phi }),`$
$`a^2(\overline{\rho }+\overline{p})V`$ $`=`$ $`2\mathrm{M}_{\overline{P}}^2k\left(\mathrm{\Psi }\dot{\mathrm{\Phi }}\right),`$ (66)
$`a^2\overline{p}\mathrm{\Pi }`$ $`=`$ $`\mathrm{M}_{\overline{P}}^2k^2(\mathrm{\Phi }+\mathrm{\Psi }),`$ (67)
In the above, it is understood that the quantities $`\mathrm{\Delta },V`$ and $`\mathrm{\Pi }`$ are the sum of the contributions of all species $`\alpha `$. Using $`\dot{w}=(c_s^2w)\dot{\overline{\rho }}/\overline{\rho }`$ and (64) we get from $`T_{\mathrm{\hspace{0.33em}\hspace{0.33em}0};\mu }^\mu =0`$ that

$$\dot{\mathrm{\Delta }}+3(c_s^2w)\mathrm{\Delta }+kV(1+w)+3w\mathrm{\Gamma }=0,$$
(68)
and from $`T_{i;\mu }^\mu =0`$

$`\dot{V}=(3c_s^21)V+k[\mathrm{\Psi }3c_s^2\mathrm{\Phi }]`$
$`+{\displaystyle \frac{c_s^2k}{1+w}}\mathrm{\Delta }+{\displaystyle \frac{wk}{1+w}}\left[\mathrm{\Gamma }{\displaystyle \frac{2}{3}}\mathrm{\Pi }\right].`$ (69)
### A.2 Gauge-invariant quintessence perturbations

The scalar quintessence field is decomposed into a background and fluctuation part according to $`\phi (\tau ,𝒙)=\overline{\phi }(\tau )+\chi (\tau ,𝒙)`$. The fluctuation can be promoted to a gauge-invariant quantity by defining the gauge-invariant quintessence field fluctuation $`X\chi \dot{\overline{\phi }}k^1\sigma `$. The field dynamics is governed by the Klein-Gordon equation. For the background, it reads:

$$\ddot{\overline{\phi }}=2\dot{\overline{\phi }}a^2V^{}(\phi ),$$
(70)
while the perturbation obeys the equation of motion

$$\begin{array}{c}\ddot{X}=\dot{\overline{\phi }}(\dot{\mathrm{\Psi }}3\dot{\mathrm{\Phi }})2a^2V^{}(\phi )\mathrm{\Psi }\hfill \\ \hfill (a^2V^{\prime \prime }(\phi )+k^2)X2\dot{X}.\end{array}$$
(71)
From the energy momentum tensor for the quintessence field

$$T_\nu ^\mu =\phi ^{,\mu }\phi _{,\nu }\delta _\nu ^\mu \left(\frac{1}{2}\phi ^{,\alpha }\phi _{,\alpha }+V(\phi )\right),$$
(72)
using $`\phi =\overline{\phi }+X`$ and the longitudinal gauge metric, one gets

$`\delta T_{0}^{0(lon.)}`$ $`=`$ $`\left[a^2\left(\dot{\overline{\phi }}^2\mathrm{\Phi }\dot{X}\dot{\overline{\phi }}\right)V^{}(\phi )X\right]Q`$ (73)
$`\delta T_{0}^{i(lon.)}`$ $`=`$ $`a^2k\dot{\overline{\phi }}XQ^i.`$ (74)
Using the definition of $`\mathrm{\Delta }`$, equation (61) in longitudinal gauge and $`\overline{\rho }_q+\overline{p}_q=a^2\dot{\overline{\phi }}^2`$ one can read off from equation (73) the gauge-invariant expression

$$\mathrm{\Delta }_q=(1+w_q)\left[3\mathrm{\Phi }\mathrm{\Psi }+\dot{X}\dot{\overline{\phi }}^1\right]+XV^{}(\phi )\overline{\rho }_q^{(1)}.$$
(75)
In the same manner, one gets from equation (74) and the fact that $`v^{(long.)}=V`$ the relation

$$V_q=k\dot{\overline{\phi }}^1X.$$
(76)
Taking the time derivative of equations (75) and (76) and using the equation of motion (71), one obtains the evolution equations

$`\dot{\mathrm{\Delta }}_q=(1`$ $`+`$ $`w_q)[{\displaystyle \frac{2a^2V^{}(\phi )}{\dot{\overline{\phi }}}}({\displaystyle \frac{\mathrm{\Delta }_q}{1+w_q}}3\mathrm{\Phi })`$ (77)
$`+`$ $`({\displaystyle \frac{6a\dot{a}V^{}(\phi )}{k\dot{\overline{\phi }}}}k)V_q]+{\displaystyle \frac{\dot{w_q}\mathrm{\Delta }_q}{1+w_q}}`$
and

$$\dot{V}_q=k\left[\frac{\mathrm{\Delta }_q}{1+w_q}3\mathrm{\Phi }+\mathrm{\Psi }\right]+2V_q.$$
(78)
Equation (77) depends on the specific quintessence model through $`V^{}`$ and $`\dot{\overline{\phi }}`$. We can however make progress in the case of nearly constant $`w_q`$: Many quintessence models have solutions for which $`\phi `$ approaches an attractor solution irrespectively of its initial value. For these tracking quintessence models Wetterich:fm ; Ratra:1987rm ; Steinhardt:nw , the equation of state of the quintessence field $`w_q`$ is nearly constant during radiation domination. We will use this vanishing of $`\dot{w}_q`$ in the following to derive relations to simplify equation (77). Considering $`a^2\dot{\overline{\phi }}^2=(1+w_q)\rho _\phi `$ it follows using the Friedman equation $`3a^2\mathrm{M}_{\overline{P}}^2^2=\rho `$ that

$$\dot{\overline{\phi }}=[3(1+w_q)\mathrm{\Omega }_q]^{\frac{1}{2}}\mathrm{M}_{\overline{P}},$$
(79)
and hence

$$\frac{\ddot{\overline{\phi }}}{\dot{\overline{\phi }}}=\frac{\text{d}}{\text{d}\tau }\mathrm{ln}\dot{\overline{\phi }}=\frac{1}{2}\frac{\dot{\mathrm{\Omega }}_q}{\mathrm{\Omega }_q}+\frac{\dot{}}{},$$
(80)
where we have neglected a term involving $`\dot{w}_q`$. We will in the following assume that at early times, the universe expands as if radiation dominated. In this case, $`=\tau ^1`$ and inserting the above equation (80) into the equation of motion (70), one finds

$$\frac{a^2V^{}}{\dot{\overline{\phi }}}=\frac{3(1w_q)}{2\tau }.$$
(81)
Using this relation (81), the evolution equation for $`\mathrm{\Delta }_q`$ becomes

$$\begin{array}{c}\dot{\mathrm{\Delta }}_q=3(w_q1)\frac{k}{x}[\mathrm{\Delta }_q3(1+w_q)\mathrm{\Phi }\hfill \\ \hfill +\{3\frac{x^2}{3(w_q1)}\}(1+w_q)\stackrel{~}{V}_q],\end{array}$$
(82)
whereas the one for the velocity remains almost unaltered while we move to the reduced velocity $`\stackrel{~}{V}_q`$:

$$\dot{\stackrel{~}{V}}_q=\frac{k}{x}\left[\frac{\mathrm{\Delta }_q}{1+w_q}3\mathrm{\Phi }+\mathrm{\Psi }\right]+\tau ^1\stackrel{~}{V}_q.$$
(83)
Note that $`\mathrm{\Gamma }_q`$ does not usually vanish. Instead, we obtain

$$w_q\mathrm{\Gamma }_q=(1c_{s(q)}^2)\left[\mathrm{\Delta }_q3(1+w_q)\mathrm{\Phi }+3\frac{\dot{a}}{a}(1+w_q)\frac{V_q}{k}\right]$$
(84)
with the sound speed of quintessence given by

$$c_{s(q)}^2=\dot{\overline{p}}_q/\dot{\overline{\rho }}_q=w_q\frac{1}{3}\frac{a}{\dot{a}}\frac{\dot{w}_q}{1+w_q}$$
(85)
### A.3 Matter and Radiation

Setting $`w=c_s^2=\mathrm{\Gamma }=0`$ in equations (68) and (A.1), we obtain the cold dark matter evolution equations

$`\dot{\mathrm{\Delta }}_\mathrm{c}`$ $`=`$ $`kx\stackrel{~}{V}_\mathrm{c},`$ (86)
$`\dot{\stackrel{~}{V}}_\mathrm{c}`$ $`=`$ $`{\displaystyle \frac{k}{x}}(\stackrel{~}{V}_\mathrm{c}+\mathrm{\Psi }).`$ (87)
The multipole expansion of the neutrino distribution function Ma:1995ey ; Durrer::fund  can be truncated beyond the quadrupole at early times. In terms of density, velocity and shear, it is given by Durrer::fund ; Doran:2003sy

$`\dot{\mathrm{\Delta }}_\nu `$ $`=`$ $`{\displaystyle \frac{4}{3}}kx\stackrel{~}{V}_\nu ,`$ (88)
$`\dot{\stackrel{~}{V}}_\nu `$ $`=`$ $`{\displaystyle \frac{k}{x}}\left({\displaystyle \frac{1}{4}}\mathrm{\Delta }_\nu \stackrel{~}{V}_\nu {\displaystyle \frac{1}{6}}x^2\stackrel{~}{\mathrm{\Pi }}_\nu +\mathrm{\Psi }\mathrm{\Phi }\right),`$ (89)
$`\dot{\stackrel{~}{\mathrm{\Pi }}}_\nu `$ $`=`$ $`{\displaystyle \frac{k}{x}}\left({\displaystyle \frac{8}{5}}\stackrel{~}{V}_\nu 2\stackrel{~}{\mathrm{\Pi }}_\nu \right).`$ (90)
Deep in the radiation dominated era, for which the initial conditions here are derived, Compton scattering tightly couples photons and baryons Peebles:ag ; Kodama:bj . The coupling leads to $`V_\mathrm{b}=V_\gamma `$ and the evolution equations become Kodama:bj

$`\dot{\mathrm{\Delta }}_\gamma `$ $`=`$ $`{\displaystyle \frac{4}{3}}kx\stackrel{~}{V}_\gamma ,`$ (91)
$`\dot{\stackrel{~}{V}}_\gamma `$ $`=`$ $`{\displaystyle \frac{k}{x}}\left({\displaystyle \frac{1}{4}}\mathrm{\Delta }_\gamma \stackrel{~}{V}_\gamma +\mathrm{\Psi }\mathrm{\Phi }\right),`$ (92)
$`\dot{\mathrm{\Delta }}_\mathrm{b}`$ $`=`$ $`kx\stackrel{~}{V}_\gamma .`$ (93)
As the photon quadrupole and all higher photon multipoles are suppressed during tight coupling, it follows that $`\mathrm{\Phi }`$ is given from Einstein’s equation by

$$\mathrm{\Phi }=\mathrm{\Psi }\mathrm{\Omega }_\nu \stackrel{~}{\mathrm{\Pi }}_\nu ,$$
(94)
where we have used the Friedmann equation. Finally, the Poisson equation (65) in terms of the various species is

$$\mathrm{\Psi }=\frac{\underset{\alpha =c,b,\gamma ,\nu ,q}{}\mathrm{\Omega }_\alpha (\mathrm{\Delta }_\alpha +3(1+w_\alpha )\stackrel{~}{V}_\alpha )}{\underset{\alpha =c,b,\gamma ,\nu ,q}{}3(1+w_\alpha )\mathrm{\Omega }_\alpha +\frac{2x^2}{3}}\mathrm{\Omega }_\nu \stackrel{~}{\mathrm{\Pi }}_\nu ,$$
(95)
where the index $`\alpha `$ runs over all species. Rewriting the evolution equations (86) - (93) in terms of $`\text{d}/\text{d}\mathrm{ln}x`$ and replacing $`\mathrm{\Phi }`$ by means of (94), one arrives at (2)-(11).

## Appendix B Early time quintessence field fluctuations

While throughout this work, we describe quintessence perturbations by the variables $`\{\mathrm{\Delta }_q,V_q\}`$, one could instead use the field fluctuation and its time derivative $`\{X,\dot{X}\}`$. In this section, we will give analytic expressions for $`X`$ and $`\dot{X}`$ in the case of tracking quintessence for super-horizon modes. We will do so assuming that $`\mathrm{\Psi }`$ and $`\mathrm{\Phi }`$ are at least almost constant. As this is not the case for CDM isocurvature and baryon isocurvature, the following steps do not apply in these modes. Furthermore, we will assume that the universe expands as if radiation dominated during the time of interest. In this case, $`=\tau ^1`$, $`\mathrm{\Omega }_q\tau ^{13w_q}`$ and hence by means of equation (79) $`\dot{\overline{\phi }}\tau ^{\frac{1}{2}(1+3w_q)}`$. Using this, we infer from equation (81) that $`V^{}\tau ^{\frac{1}{2}(7+3w_q)}`$. In addition, a straightforward calculation using (80) and (81) yields

$$a^2\tau ^2V^{\prime \prime }=a^2\tau ^2\frac{\text{d}V^{}}{\text{d}\tau }\frac{\text{d}\tau }{\text{d}\phi }=\frac{3}{4}(1w_q)(7+3w_q).$$
(96)
The equation of motion for $`X`$ (71) contains a term $`\dot{\overline{\phi }}\left(\dot{\mathrm{\Psi }}3\dot{\mathrm{\Phi }}\right)`$, which by assumption we may drop. In addition, we see from equation (96), that for super-horizon modes, $`a^2V^{\prime \prime }k^2`$ (except for $`w_q`$ very close to $`1`$), and hence the equation of motion reduces to

$$\ddot{X}=2a^2V^{}\mathrm{\Psi }a^2V^{\prime \prime }X2\frac{\dot{a}}{a}\dot{X}.$$
(97)
Using the power law behaviour in $`\tau `$ of $`V^{},V^{\prime \prime }`$ and $`a`$, as well as equations (81) (96), one finds the particular solution

$$X(\tau )=\frac{\tau }{2}\mathrm{\Psi }\dot{\overline{\phi }},$$
(98)
as well as two complementary solutions that may be added to obtain the general solution

$$\begin{array}{c}X(\tau )=\frac{\tau }{2}\mathrm{\Psi }\dot{\overline{\phi }}+c_1\tau ^{\frac{1}{2}\left(1\sqrt{14a^2\tau ^2V^{\prime \prime }}\right)}\hfill \\ \hfill +c_2\tau ^{\frac{1}{2}\left(1+\sqrt{14a^2\tau ^2V^{\prime \prime }}\right)}.\end{array}$$
(99)
The mode proportional to $`c_2`$ is at least as rapidly decaying as the one proportional to $`c_1`$. Using the explicit form of $`4a^2\tau ^2V^{\prime \prime }`$, equation (96), we find that $`\sqrt{14a^2\tau ^2V^{\prime \prime }}`$ is imaginary if $`w_q[\frac{2}{3}(1+\sqrt{6}),\frac{2}{3}(1\sqrt{6})]`$, which holds for all scalar quintessence models of current interest. Hence, the complementary modes decay $`1/\sqrt{\tau }`$ in an oscillating manner.

Coming back to the dominating particular solution (98), Figure 1 shows that the accuracy of this analytic result is indeed high at early times, when compared to numerical simulations.

Inserting the solution (98) and its time derivative into equation (75), we find the simple expression

$$\mathrm{\Delta }_q=3(1+w_q)\left(\mathrm{\Phi }\frac{1}{2}\mathrm{\Psi }\right),$$
(100)
which is just a restatement of eqn. (23) and (24). Hence, the energy density contrast in tracking quintessence models remains constant on super horizon scales, provided the gravitational potentials are constant to good approximation.

## Appendix C Extended Matrix Formulation

For simplicity, we have limited the discussion in section (II.2) to cases where either $`w_q=1/3`$, or quintessence contributions to $`A(x)`$ are neglected. Here, we will discuss cases for which $`w_q<1/3`$, while the background expands radiation dominated. In this case, $`\mathrm{\Omega }_q\tau ^{(13w_q)}`$ and we can split the matrix in three parts according to their scaling with $`x`$:

$$A(x)=A_0+xA_1+x^{(13w_q)}A_q.$$
(101)
Again, Equation (1) will lead to a solution vector of the form

$$𝑼(x)=𝑼_0+x𝑼_1+x^{(13w_q)}𝑼_q.$$
(102)
Substituting this into Equation (1) and keeping only leading orders in $`x`$, we get

$`A_0𝑼_0`$ $`=`$ $`0,`$ (103)
$`A_1𝑼_0+A_0𝑼_1`$ $`=`$ $`𝑼_1,`$ (104)
$`A_q𝑼_0+A_0𝑼_q`$ $`=`$ $`(13w_q)𝑼_q.`$ (105)
While the conclusion regarding $`𝑼_0`$ and $`𝑼_1`$ are still the same as in section (II.2), we see that quintessence may introduce a correction

$$𝑼_q=\left[A_0(13w_q)𝟙\right]^1A_q𝑼_0.$$
(106)
This contribution $`x^{(13w_q)}𝑼_q`$ could in principle dominate over $`x𝑼_1`$ for $`w_q>0,\mathrm{\Omega }_q>\mathrm{\Omega }_c`$. However, the contribution is only of interest for the CDM isocurvature and baryon isocurvature modes, as it is otherwise negligible compared to the constant order. Yet for CDM isocurvature and baryon isocurvature, $`A_q𝑼_0=0`$. Therefore, the discussion below applies, leading to $`𝑼_q=0`$ for CDM isocurvature and baryon isocurvature modes. One order higher in $`x`$, there may be a contribution. Yet this is in any case a higher order contribution, which we may neglect.

Finally, we briefly discuss the case of vanishing $`𝑼_0`$. This only concerns possible subdominant modes. Equation (104) then yields $`A_0𝑼_1=𝑼_1`$, i.e. $`𝑼_1`$ is an eigenvector of $`A_0`$ with eigenvalue $`\lambda =1`$. As $`A_0`$ does not have such an eigenvector, we are led to conclude that Equation (1) does not have a regular solution involving $`𝑼_1`$, if $`𝑼_0=0`$. Turning to Equation (105), we similarly conclude that $`𝑼_q`$ needs to be a eigenvector of $`A_0`$ with $`\lambda =(13w_q)`$ for vanishing $`𝑼_0`$. For $`w_q<1/3`$ this is once again excluded and for $`w_q=1/3`$, we just regain the results of section II.2.