{'fraction_non_alphanumeric': 0.030320182638063433, 'fraction_numerical': 0.010037546461387522, 'mean_word_length': 4.554018652415383, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Adiabatic perturbations in coupled scalar field cosmologies

We present a comprehensive and gauge invariant treatment of perturbations around cosmological scaling solutions for two canonical scalar fields coupled through a common potential in the early universe, in the presence of neutrinos, photons and baryons, but excluding cold dark matter. This setup is relevant for analyzing cosmic perturbations in scalar field models of dark matter with a coupling to a quintessence field. We put strong restrictions on the shape of the common potential and adopt a matrix-eigensystem approach to determine the dominant perturbations modes in such models. Similar to recent results in scenarios where standard cold dark matter couples to quintessence, we show that the stability of the adiabatic perturbation mode can be an issue for this class of scalar field dark matter models, but only for specific choices of the common potential. For an exponential coupling potential, a rather common shape arising naturally in many instances, this problem can be avoided. We explicitly calculate the dominant perturbation modes in such scenarios.


I. INTRODUCTION
Scalar fields can appear in modern cosmology in many forms. They are typically associated with phenomena which can not be attributed to ordinary forms of matter, most commonly with an accelerated expansion of the universe. Thus they are particularly important in the very early universe for modeling the era of inflation [1][2][3][4][5][6], and the very late universe, where scalar quintessence models are one of the most promising candidates for explaining dark energy [7][8][9][10][11][12][13].
For the long time between these eras, scalar fields typically play only a supporting role in the cosmic evolution. From reheating onwards, throughout the entire radiation dominated era, their energy densities only contribute a small percentage to the total energy density of the universe [14][15][16][17][18]. The same holds true for the matterdominated era up to redshifts of about 5, with the exception of scalar field models of dark matter [19][20][21][22][23][24][25]. Still, much effort has ben put into analyzing scalar field dynamics in these eras. The reason for this lies in the so called coincidence problem, which adresses the question why the energy densities of dark energy and matter are roughly of the same order of magnitude only in the recent past. In the ΛCDM model this seems like a very precise finetuning of the cosmological constant. Scalar field models of dark energy provide a possible solution to this problem, since they can exhibit a behavior often called tracking or scaling [26][27][28][29], where the field equations have solutions for which the dynamics follow a specific trajectory largely irrespective of the initial conditions chosen, thereby explaining why dark energy became dominant only recently.
In this paper we are interested in such scaling solutions involving multiple scalar fields which couple through their common potential and the evolution of perturbations in such scenarios. Similar models have been investigated before in the context of inflation [30,31], i.e. without any additional matter content. We consider the dynamics of multiple scalar fields in the presence of ordinary matter, i.e baryons, neutrinos and photons, but exclude cold dark matter or a cosmological constant. This setting is relevant for models where the entire dark sec-tor, i.e. both dark energy and dark matter is described by scalar fields.
If one looks beyond the homogeneous and isotropic background cosmology, one quickly finds that the the inhomogeneities of scalar fields in such models are of great importance for the study of cosmological perturbations in general. This starts with the fact that all structures seen in the universe today are now widely believed to have been sourced by quantum fluctuations of scalar fields during the era of inflation, but does not end there. Scalar quintessence models can influence late-time structure growth [32,33], but even during the era of radiation domination, where scalar fields are subdominant, their presence can impact how perturbations in the baryonic or dark matter sector evolve. While standard minimally coupled one-field scaling quintessence models seem to be unproblematic in this respect [34], recent studies have shown that a coupling between a scalar field and dark matter can result in an instability of the adiabatic perturbation mode [35,36], potentially rendering such scenarios difficult to reconcile with observational bounds, which clearly demand largely adiabatic perturbations [37][38][39][40][41][42][43][44][45][46]. It is therefore interesting to investigate whether such problems can also arise in models where multiple scalar fields couple to each other, e.g. in coupled scalar field dark matter models, and we do so in our analysis below. This paper is organized as follows: In section II we derive the generic shape of the potential necessary for multiple canonical scalar fields to exhibit scaling solutions in the presence of a background fluid. In section III we study the evolution of linear perturbations in the superhorizon regime in such models, in the presence of photons, baryons and neutrinos. We restrict ourselves to the case of two scalar fields here for simplicity. In section IV, we apply the results to one particular choice of scalar potential with an exponential coupling. This is the relevant potential for the coupled scalar dark matter model investigated in an accompanying paper [47]. Finally, in section V we investigate alternative potentials and construct an explicit example for which the adiabatic perturbation mode is unstable. We present our conclusions in section VI.

II. COUPLED CANONICAL SCALAR FIELDS AND SCALING SOLUTIONS
In this section we investigate exact scaling solutions in FLRW-cosmologies involving multiple coupled canonical scalar fields and a background matter fluid, defined by an energy density ρ m and a constant equation of state ω m . In the literature the definition of a scaling solution is not completely unambiguous. In the case of a single scalar field it can either mean a solution for which the scalar energy density scales like the scale factor to some constant power, or, and this can be considered a stricter definition, a solution for which the scalar energy density scales exactly like the background fluid. In this work we adopt the stricter definition and consider an exact scaling solution to be a scenario in which all scalar energy densities scale like the background fluid.
It is well known that the only potential providing an exact scaling solution for a single canonical scalar field is an exponential potential. This even holds true in the presence of a coupling to the background fluid [48][49][50][51]. Furthermore, even in the presence of non-canonical kinetic terms -however still restricted to those yielding second order field equations -the form of possible Lagrangians yielding scaling solutions can be strongly restricted [51].
We will now show that in the case of multiple canonical coupled scalar fields similar restrictions for the common potential can be found. To do this we employ an approach already used in the earlier papers cited above, slightly adjusted to fit our scenario. We start by considering the following action where X ϕi = 1 2 D µ ϕ i D µ ϕ i , n denotes the number of scalar fields ϕ i present and S m denotes the matter action. The scalar field equations for this action in the FLRW-background read where a prime denotes a derivative with respect to conformal time and h is the conformal hubble parameter h = a /a. The common potential V can in principal be identified with any one of the scalar field energy densities or even split in some arbitrary fashion, but we choose to assign it to ϕ 1 for simplicity and define scalar energy densities and pressure densities as follows: ρ ϕj = p ϕj = X ϕj , j = 2, ..., n .
As we assume no direct coupling between the background fluid and the scalar sector, the equations governing the background-evolution can be written as where N = ln(a) and conservation of the total energymomentum tensor implies that the scalar couplings have to satisfy The definition of the couplings q i in this manner is of course a matter of convention, here we remain consistent with the general formulas given in appendix A. In an exact scaling scenario we demand Note that the scalar equations of state can differ from the background equation of state even in an exact scaling scenario if a coupling between the scalar fields is present, only the combined scalar equation of state has to fulfill if the scalar density parameters are non-zero. From equations (5) and (8) we can directly conclude that and Since X ϕi ∝ (dϕ i /dN ) 2 ρ tot , equations (10) and (11) directly give Furthermore all couplings have to be constant by virtue ofequations (6) and (10): Now we can easily derive a differential equation for the potential V by again employing equations (6) and (10): which directly gives where ξ i = ϕ i − ci c1 ϕ 1 and f is some arbitrary (smooth) function. One can quickly check that this solution is not in conflict with the the demand of constant couplings simply by noting that by virtue of equation (12) all ξ i have to be constant for a scaling solution and by convention (3) we therefore have (16) for j > 2 and dynamics corresponding to a scaling solution.
One comment might be in order here: In the derivation of equation (15) we made some implicit assumptions. The most obvious one is that ω ϕ1 = 1, i.e. the scalar potential energy density should be non-vanishing. As we will see in section IV, purely kinetic fixed points with ω ϕ1 = 1 do exist, and we expect them to be present for a large class of potentials, including many not of the shape derived here. The only restriction here is that the common potential should become zero somewhere, at least asymptotically. Furthermore one might suggest that we have excluded solutions for which cosmological expansion is purely scalar field dominated. This is however not the case, as one can simply replace all ω m appearing in the above equations with a constant effective equation of state ω eff and reach the same conclusions concerning the potential shape.
To clarify this result: If you want to construct a potential allowing for exact scaling solution with non-vanishing potential energy density involving n multiple canonical scalar fields coupled through their potential, you first have to check if your potential fulfills equation (15) for some set of constants c 1 , ..., c n . If so, a scaling solution can exist.

III. EARLY UNIVERSE PERTURBATIONS
We will now turn our attention to the evolution of perturbations in the early universe in models containing multiple non-minimally coupled canonical scalar fields. In order to keep things simple, we restrict our attention to the case of two fields and call them ϕ and χ, deviating from the conventions used in section II. We also split the potential differently, this time assigning a ϕ-dependent part to ρ ϕ and a common part to ρ χ , i.e.
By "early universe" we mean an era during which the cosmic expansion proceeds as if dominated by radiation, i.e. h(η) = η −1 , with η being conformal time. We will consider a background-evolution during which both fields follow an exact scaling bahaviour, in particular the equations of state are assumed to be constant. Note that this does not imply a strongly subdominant role of the scalar fields or equations of state of exactly 1/3. The coupling between the two fields enables a wide range of possible ω ϕ 's and ω χ 's while still allowing for a radiation-like expansion (a(η) ∝ η) and non-negligible scalar field contributions to the energy-densities already in very simple models, as we will see below.
In addition to the two coupled scalar fields we include neutrinos, photons and baryons into our model, but exclude a possible cold dark matter contribution. This sets our study apart from former systematic studies of early universe perturbation modes, which have usually been limited to cosmologies containing either no scalar fields [52], only scalar fields (usually in the context of inflation) [30] or a single quintessence field [31,34,36,53,54]. Leaving out a possible cold dark matter contribution is realistic, if at least one of the two scalar fields changes its dynamics away from a scaling solution and acts like a dark matter component during the later stages of its evolution. If this is the case, the coupling arising from the common potential will result in a coupled quintessence model in the late universe [11,55,56], more specifically a coupled quintessence-scalar-dark-matter model [25]. It is precisely these kinds of models that require the analysis presented here to draw the initial conditions for numerical evolutions of cosmological perturbations, a study we performed in an accompanying work [47].
Before moving on to our analysis, let us recap some relevant results from the previous studies. As was shown in [31], even in generic non-minimally coupled models (which of course includes ours), purely adiabatic perturbations remain purely adiabatic on superhorizon scales, no matter what the interaction is. Here a mode is called adiabatic if all entropy perturbations vanish. This includes relative entropy perturbations between different components of the cosmic fluid as well as internal entropy perturbations for a single component -these can exist for non-perfect fluids (see section 4 in appendix A for more details). However, as was noted in ref. [34], the existence of such an adiabatic mode is non-trivial, since demanding all entropy perturbations to vanish typically introduces more constraints than we have degrees of freedom in the equations, in particular if the number of species considered is large. One should understand that this is not in conflict with the theorem found by Weinberg in [57,58], which states that an adiabatic superhorizon mode always exists, no matter the content of the universe. The apparent contradiction stems from Weinberg's definition of an adiabatic superhorizon-mode, which he takes to mean a mode for which the total curvature perturbation ζ remains constant, which is the case if the total entropy perturbation Γ tot = 0 (see appendix A). This is a weaker condition than demanding all entropy perturbations to vanish, which is what we take adiabatic to mean in this work.
Furthermore this result is not sufficient to explain a suppression non-adiabatic modes, since any, presumably initially small, admixture of such modes could in principle eventually outgrow the adiabatic perturbations. Such a growth can happen for example during the 'adjustment' phases in minimally coupled tracking quintessence models when the quintessence field does not yet follow the tracker trajectory. However, a sufficiently long tracking regime typically erases these modes in this case [53].
The situation in non-minimally coupled models of quintessence appears to be more complicated. We are aware of only one study which systematically analyzes the evolution of superhorizon-perturbations in a coupled quintessence model [36]. It appears that in such scenarios there are certain sections of parameter space for which the adiabatic mode becomes unstable, i.e. other, faster growing perturbation modes exist. Since this analysis includes a conventional cold dark matter component and uses a specific commonly used form of the coupling given by Q ϕ = −βρ cdm , the question if similar issues will arise in coupled scalar field models of dark matter remains open, and we will address it below.

A. Basic Formalism
The basic idea we employ in our study of early universe superhorizon perturbations is one already used in earlier works [34,53,54]. We write the differential equations governing the evolution of linear perturbations in a convenient matrix form, i.e.
Here we introduced a convenient new time variable x ≡ k/h and combined all relevant perturbative quantities into a single perturbation vector U (x). Which quantities these are depends on the approximations one choses to use in the neutrino-, photon-and baryonic sectors. Here we employ the simplest (lowest order) version of the tight-coupling approximation for photons and baryons, which yields a common velocity potential for both components and excludes all higher momenta of the Boltzmann hierarchy for photons. In the case of neutrinos we truncate the Boltzmann-expansion after the quadrupole, leaving an additional anisotropic stress contribution. More details and precise definitions of the perturbative quantities used in this chapter can be found in appendix A. Since the scalar fields only have two degrees of freedom, the components of the perturbation vector can be chosen to be the energy density contrasts ∆ α and velocity potentials V α for all components of the cosmic fluid as well as anisotropic stress for neutrinos, i.e.
The entries of the perturbation matrix A can now be read off from the differential equations governing the evolution of these variables. For the photon-, baryon-and neutrino-sector they read These equations can be derived from the generic equations given in appendix A by setting ω eff = 1/3. Note that the equation for the velocity potential V γb is different from the one usually employed, e.g. in refs. [34,36]. This version is the correct one if one makes no assumption about the ratio Ω b /Ω γ . Both version agree in the limit Ω b /Ω γ → 0. The gravitational potentials are then given by The generic equations of motion for the scalar fields are quite complicated (see appendix A), but they can be simplified considerably in an exact scaling scenario. The background quantities appearing in the scalar sector are the equations of state ω ϕ and ω χ , the adiabatic sound speeds c 2 a,ϕ and c 2 a,χ , the couplings q ϕ and q χ , their time-derivatives and the second derivative of the common potential V 2,ϕϕ .
In an exact scaling scenario we assume constant equations of state, i.e. we set ω ϕ = ω χ = 0. This implies that the adiabatic sound speeds are given by c 2 a,ϕ = ω ϕ and c 2 a,χ = ω χ . Furthermore the couplings q ϕ and q χ have to be constant by equation (13) and are of course related by equation (7). Finally the scalar equations of state are related by equation (9) and we only have to deal with the second derivative V 2,ϕϕ of the common potential. For a generic scalar field evolution this can of course take any, in principle time-dependent, form, but in a scaling scenario it merely introduces another constant into the equations, which can be seen as follows.
We start by investigating the generic potential for exact scaling solutions given by equation (15). In order to split our potential according to equation (17) we have to rewrite f as where ξ = χ − c 2 /c 1 ϕ and µ 1 , µ 2 are some suitable constants. Note that this already shows that V 1 has to be an exponential potential. We can take the second derivative of this potential with respect to ϕ and obtain (for (31) Since ξ is constant during a scaling solution, it is obvious that V 2,ϕϕ ∝ V ∝ ρ tot in such scenarios and we can thus define with some constant r ϕ . Note that r ϕ is in principle a new independent constant depending on the functional form of g and can generally generally not be related to q ϕ , as it depends on the second derivative g , whereas only the first derivative g enters into q ϕ . Now we can finally write down the equation for the scalar sector in a form where the only remaining constant parameters are ω ϕ , r ϕ and the density parameters. They The coefficients of the matrix A(x) can be read off from equations, (21) - (26) and (33) - (36) after replacing the gravitational potentials. Before analyzing the solutions for the perturbations we have to investigate the background quantities appearing in A(x). Those are generally x-dependent, but in this analysis we are interested only in the superhorizon-limit defined by x 1 and we can therefore work with Taylorexpansions in x. To see that, first note that we assume a radiation-like expansion with a scalar scaling solution in the early universe, a scenario only slightly disturbed by the baryonic density parameter, which is small, but grows linearly in x: This implies that all the other density parameters are constant to leading order, but decrease slowly at linear order in x . Evaluating Friedmanns equations for a flat universe in the form Ω ν + Ω γ + Ω ϕ + Ω χ + Ω b = 1 order by order in x then yields to leading order Using this, it becomes obvious that we can expand the matrix A(x) as a type of Taylor-series in x and recover a constant matrix at leading order: With this approximation we can immediately write down the leading order solution for equation (19), which reads where U (i) 0 are eigenvectors of the matrix A 0 and λ i the corresponding eigenvalues. With this as a starting point, one can now easily determine higher order corrections to the found solution by expanding the eigenvectors as follows: where the first order correction is then given by and the other corrections can be calculated in a similar fashion, but we will not need them here. Of course, the general solution is given by A short comment concerning the different approximations used here might be in order. As explained in appendix A we use the simplest version of the tight coupling approximation for photons and baryons, i.e. we only include coupling terms to leading order (that is zeroth order) of hτ c 1, where τ c = (an e σ T ) −1 is the Thomson interaction timescale. One might question whether it is consistent to go to next to leading order in the x-expansion of the matrix A(x), but not include the next order in the TCA, i.e. terms suppressed by hτ c . As it turns out, for the relevant wavenumbers we have x = k/h hτ c during radiation domination. This breaks down for very small k, but for realistic cosmologies the boundary lies at roughly 10 −4 Mpc −1 , which is already far below what is currently observable.

Eigenvalues
In order to determine the dominant perturbation modes in the early universe in coupled two scalar field models one simply has to determine the eigenvalues λ i of the matrix A 0 with the largest real parts and the corresponding eigenvectors. In principle these depend on the 4 parameters determining the scaling solution, which can be taken to be ω ϕ , r ϕ , Ω ϕ and Ω χ . The remaining density parameters for neutrinos and photons are fixed by the fact that we assume a FLRW metric without curvature, i.e. Ω ν + Ω γ + Ω ϕ + Ω χ = 1 and the well known relation Ω ν = 21/8 × (4/11) 4/3 Ω γ , which is valid after electronpositron annihilation [59]. However, four of the eigenvalues are always given by {λ 1 , λ 2 , λ 3 , λ 4 } = {0, 0, 0, −1}, irrespective of the values of these parameters. The remaining 6 eigenvalues are very complicated functions of the 4 free parameters and quoting them is not very enlightening. We will however present numerical results below, obtained simply by running through a grid of the four independent parameters {ω ϕ , r ϕ , Ω ϕ , Ω χ }. For the wide range of parameters we investigated, 4 of the 6 remaining eigenvalues had real parts which were bound from above by −0.5. Only two eigenvalues can have a real part bigger than zero, a case will will call strongly growing from now on. As we will see below, the nullspace of A 0 always contains an adiabatic mode. It is precisely the two potentially strongly growing modes which can render the adiabatic mode unstable, which would make the corresponding scenario difficult to reconcile with observations [37][38][39][40][41][42][43][44][45][46]. On could of course conceivably come up with reasons why these modes are initially strongly suppressed compared to the adiabatic mode, so that it does not become dominant during the radiation dominated era, but such a fine-tuning of initial conditions is often not very natural and generally undesirable.
To visualize the effect of the coupling between the two scalar fields we have plotted the real part of one of the two critical eigenvalues in Figures 1 and 2. As can be seen, the values become bigger than zero for large deviations of ω ϕ from the uncoupled value of 1/3, i.e. for big couplings. Furthermore, as one would expect, this effect becomes stronger for larger scalar field density parameters. In the uncoupled case, even the critical eigenvalues have real parts bound by −0.5 from above and the corresponding modes are thus subdominant. This shows that a strong coupling can be relevant for early universe perturbations, because it can potentially change the dominant perturbation mode, rendering the adiabatic mode unstable.
One should note however, that just because we can find a set of parameters which yield a strongly growing eigenvalue in this numerical treatment does not necessarily mean that we can construct a model actually yielding a scaling solution with such a feature. The problem is of course, that for any given scaling solution some of the parameters treated here as free might well be related, or scaling solutions with a radiation like expansion might not even exist or be confined to specific regions of parameter-space. Later in this work we will examine one specific potential where the r ϕ − ω ϕ relation is fixed to the solid red line in Figures 1 and 2, thus excluding the region of strongly growing eigenvalues, and a second one represented by the dashed red line, for which the adiabatic mode becomes unstable for strong couplings.

Eigenvectors
Since the generic eigenvalues have very complicated forms in general, so do the corresponding eigenvectors. This also holds for the eigenvectors belonging to the two critical eigenvalues, representing potentially dominant modes. Things look a lot better for the nullspace of A 0 , which of course includes all dominant perturbations modes in the absence of strongly growing modes. One is free to choose any suitable basis of this three-dimensional space, we choose to categorize the modes according following ref. [34], i.e. by means of the total curvature and relative entropy perturbations.
In this categorization a mode is called an isocurvature mode if the total curvature perturbation vanishes, i.e.
and adiabatic if all relative and internal entropy perturbations, i.e.
There are some deviating definitions of the total curvature pertubation in the literature [60], but they agree in the superhorizon limit in a flat universe. This categorization splits the nullspace into a twodimensional isocurvature space and one adiabatic mode. The adiabatic eigenvector is given by where R = (15 + 4Ω ν ). This is very similar to the adiabatic mode in minimally coupled quintessence scenarios [34].
The isocurvature-subspace can be further decomposed into one baryon-isocurvature and one neutrino/scalarisocurvature mode.
The baryon isocurvature mode is characterized by the demand that all species except for baryons should be adiabatic, i.e. S b,α = 0 for all α, but S α,β = 0 for α, β = b. For our set of scaling scenarios this necessarily implies that the scalar internal entropy perturbations also vanish, i.e. Γ ϕ = Γ χ = 0. In the corresponding eigenvector all entries except for ∆ b vanish at leading order, we therefore quote the result to subleading order in x here: where we defined the quantities P = (1 + 2Ω γ + 2Ω ν ) and X = (−2(1 − Ω χ ) + (1 − 3ω ϕ )Ω ϕ /4). This vector differs considerably from the ones found in previous studies, which is due not only to the presence of two coupled scalar fields, but also to the improved treatment of the tight coupling approximation, which changes the subleading order contributions to the perturbation matrix A(x). Finally, the neutrino/scalar isocurvature mode is characterized by which S γ,b = 0, as this is the only relative entropy perturbation not involving neutrinos or scalar fields. Enforcing this necessarily implies S ϕ,χ = Γ ϕ = Γ χ = 0, with all other relative entropy perturbations non-vanishing. The corresponding eigenvector reads at leading order: Note that due to the convenient choice of perturbation variables, the dependence on the coupling only appears at subleading order in the baryon isocurvature mode through ω ϕ . At leading order, the nullspace of the leading order perturbation matrix A 0 is completely independent of r ϕ or q ϕ .

IV. COUPLED EXPONENTIAL POTENTIALS
In this section we investigate one specific example of two coupled canonical scalar fields given by We chose this particular shape for two reasons. First, it arises in the cosmon-bolon model for quintessence and dark matter [25], which gets investigated more deeply in an accompanying paper [47]. Furthermore, as we will show in the next section, this simple potential is already the most generic case, in the sense that all exact scaling solutions existing for other potentials are already present for this exponential form.
One can easily check that this potential fulfills equation (15) with c 1 = 3(1 + ω eff )/α, c 2 = 4(2β − α)/αλ and For the subsequent analysis we split the potential as follows with and assign energy-and pressure-densities accordingly:

A. Exact scaling solutions
We start by finding all exact scaling solutions for the case of two coupled canonical scalar fields ϕ and χ in the presence of a radiation fluid, which is the relevant case for studying physics in the early universe. The equations of motion governing this system can be directly read off from equations (2) and (5). As is common when trying to find exact scaling solutions, we now employ the dynamical systems approach [26,66] and introduce the following new variables: where these variables are subject to the constraint since are assuming a flat universe. The equations of motion can now be rewritten as and a fixed point is characterized by dx/dN = dy/dN = du/dN = dv/dN = 0. When solving the resulting algebraic equations we can reduce the number of fixed points by eliminating some redundancies. First we assume that α > 0 and λ > 0, both can be easily achieved by a suitable sign-change in the fields ϕ and χ. Furthermore we have y > 0 and v > 0 by definition. Restricting ourselves to these ranges, the complete set of all fixed points is given in Table I. These fixed points can be split up into the purely scalar field dominated fixed points K 1 to S 3 and radiation-like fixed points R 1 to R 4 . We expect the kinetically dominated fixed points K 1 and K 2 to always exist, irrespective of the exact shape of the scalar potential, as long as it vanishes asymptotically. Since the kinetic energies are non-zero whereas the potential energies vanish, this requires an 'extended region of a zero potential' where the fields can roll freely, a scenario reached asymptotically in the limit ϕ → ∞ for our potential. Furthermore the point R1 corresponds to both scalar fields sitting at a root of the potential, in the case of ϕ at infinity. This point is present even more generically, all that is required of the potential is a (possibly asymptotic) root. However, these points are always unstable, as we will see below.
The remaining points S 1 to S 3 and R 2 to R 4 are the ones only existing for potentials of the form derived in section II and among these, only the ones exhibiting a radiation-like expansion, i.e. R 2 to R 4 , can result in a realistic early cosmology.
The stability of the fixed points can be analyzed by investigating the Jacobian matrix for x, y, u and v for each fixed point. The typical classification distinguishes between nodes, spirals and saddle-points as well as stable and unstable points [26,66]. Here we are not interested in the precise details and simply call points for which the Jacobian matrix has only positive (or zero) eigenvalues unstable, and points for which all eigenvalues are negative stable. The conditions for existence and stability of the fixed points are given in Table II.
As is common in scenarios such as ours, the stable fixed points split up the parameter space into disjunct sections. For a given set of parameters α, λ and β we therefore have a unique attractive scaling solution towards which the cosmological evolution will adjust itself relatively quickly in the early universe. To avoid this fixed point for the prolonged period of radiation domination, one would necessarily have to start with initial conditions very far away from the fixed point configuration. The splitting of the parameter space can be seen in

B. Dominant perturbation modes
The dominant perturbation modes in this potential can now be obtained in the same fashion as the general results in section III, but this time with the simplifications and The generic results for the eigenvalues do of course still hold, but as is shown in Figures 1 and 2, the particular shape of the potential leads to a form of the coupling term r ϕ that excludes the region of strongly growing eigenvalues. This holds for all stable fixed points, with the exception of R3, where Ω χ = 0 and which therefore requires some additional treatment. However, a vanishing Ω χ clearly requires ω ϕ = 1/3 by virtue of equation (9), effectively reducing the scenario to an uncoupled one, where all eigenvalues are non-positive. Thus for the coupled exponential potential there is a huge region of the parameter space which is not only viable as an early cosmology at the background level, but for which no strongly growing perturbation modes exist and the adiabatic mode is therefore stable. With the results presented here one can therefore easily write down the early universe background evolution for any set of parameters α, β and λ and also the dominant perturbation modes, given by the nullspace of A 0 , with the density parameters replaced by the ones given in Table I. Figure 3: Parameter space of stable fixed points for the coupled exponential potential with λ = 1, the structure remains valid for all λ < 2. The shaded region shows the points which allow for a realistic early universe cosmology, i.e. a radiationlike expansion.

V. NON EXPONENTIAL MODELS
We now move on to provide short arguments why all exact scaling solutions for coupled two scalar field cosmologies are effectively the ones found for the coupled exponential potential in section IV. We start by recalling the generic shape of the common scalar potential required for the existence of scaling solutions given in equation (15) together with the splitting (17), which implies via equation (30) that the common potential is given by where ξ = χ−σϕ with σ = c 2 /c 1 and α = 3(1+ω m )M/c 1 . Our aim is to rewrite the field equations in the same form used in section IV, in order to recover a set of algebraic equations whose fixed points determine the existing scaling solutions. For non-exponential potentials we can of course not express the potential derivatives in terms of the potential, and we therefore need a new set of variables, defined by The evolution equations for these variables are very simple and read ds n dN = (s n+1 − s 1 s n ) √ 6 (u − σx) .
(68) Figure 4: Parameter space of stable fixed points for the exponential coupling potential with λ = 3, the structure remains valid for all λ > 2. The shaded region shows the points which allow for a realistic early universe cosmology, i.e. a radiationlike expansion.
The field equations can now be rewritten as in section IV. They read: An algebraic solution for the scaling solutions can now be found by demanding dx/dN = dy/dN = du/dN = dv/dN = ds n /dN = 0. Luckily, the only s n -term relevant for the x, y, u and v-equation is s 1 , which has to be a constant. Thus, equations (69) to (72) give exactly the same system of algebraic equations as equations (60) to (63), with the replacements Let us quickly recap what this tells us. Any exact scaling solution for two coupled canonical scalar fields is effectively given by one of the points found above for the exponential cross-coupling potential. One should be careful not to take this conclusion too far. We emphasize that the conditions that all s n have to be static are not necessarily met by all the points above for all potentials of the shape given in equation (66). However, equations (68) and (73) tell us that all points for which u = x(2β − α)/λ, i.e. S3 and R4, do exist for all such models. Furthermore, as was discussed in section IV, the points K1, K2 exist if V 2 vanishes asymptotically and R1 does, if V 2 exhibits a real root. The stability of scaling solutions may also be very different for different potentials. While the stability-analysis in the x, y, u and v-directions proceeds precisely as in section IV, the stability of the s n provides new constraints, which depend on the functional form of the potential.
When it comes to perturbations around the possible fixed points, in particular the question of the stability of the adiabatic perturbation mode, we can draw no conclusions in general. The reason for this is, as we have seen above, a second order potential derivative, denoted by r ϕ enters the equations, which does not appear at the background level. In the language used here, this corresponds to knowledge about the term s 2 , which is not relevant for finding the scaling solutions.
To sum up: We have shown that all scaling solutions for two coupled canonical scalar fields are included in the set of fixed points found for the exponential potentials in section IV, and that two of these points exist for all potentials of the allowed shape. However, the stability conditions for the fixed points can change for different potentials, and so can the growth of linear perturbations modes around them.

A. 'Power law couplings'
As an example of a model involving strongly growing perturbation modes, we now investigate potentials with a power law coupling. By this we mean models for which g(ξ) = λξ m /M m , where obviously m has to be even in order for the potential to be bound from below. In this case the system of algebraic equations is closed after adding the equation for s 1 , since and of course s n = 0 for n ≥ m. In this case the evolution equations for s 1 reads where we already used β as defined in equation (73). As we have just seen, the point R4 exists for all choices of g, and a stability analysis reveals that stability for the power-law coupling requires As we have restricted ourselves to α > 0, it is immediately clear that stability is not possible for m = 2, but already for m = 4, the parameter range of stability is quite big, but requires a positive coupling.
In order to see how perturbations around such a scaling solutions behave, we need to calculate r ϕ . A relatively straightforward calculation shows that Here we used the superscripts pc for the 'power law coupling' models investigated in this section, and exp for the exponential coupling model investigated in section IV. Clearly a power law coupling lowers the r ϕ (ω ϕ ) for positive values of r ϕ curve compared to the exponential coupling model, which pushes it into the regime where strongly growing modes exist. We have visualized this in Figures 1 and 2  We have thus shown that it is possible to construct coupled scalar field models where stable scaling solutions with strongly growing perturbation modes exist, i.e. the adiabatic mode becomes unstable.

VI. CONCLUSION
In this paper we have analyzed exact scaling solutions in cosmologies containing multiple canonical scalar fields coupled through their potential at both the background level and at the level of linear perturbations. At the background level we have formulated the restrictions the demand that realistic scaling solutions should exist puts on the shape of the common potential and analyzed the existing scaling solutions for one specific choice of potential, and exponential coupling, in detail. As it turns out, the resulting set of scaling solutions is complete in the sense that all possible alternative potentials can not add new solutions to the list. However, the behavior of the scaling solutions concerning stability and growth of perturbations around such solutions can be different for alternative potentials.
At the level of linear perturbations we have analyzed the growth behavior of perturbations around scaling solutions in some detail. In particular we have shown that a strong coupling can dramatically change the growth of perturbation modes. For the specific case of an exponentially coupled potential however we have seen that there are no fast growing modes, i.e. no modes growing faster than the adiabatic mode demanded to be dominant by current observations. In this sense, the adiabatic mode is stable in this model, any initially small admixture of nonadiabatic modes will remain small. This is not a trivial fact, as we have demonstrated by constructing a scenario where a fast growing perturbation mode around a stable scaling solution exists. In this case one would have to find a mechanism to suppress the initial contribution of the fast growing isocurvature modes very strongly in order to remain in agreement with observations. The analysis presented here is of particular relevance for coupled scalar field models of dark matter, for which it is conceivable that they follow such a scaling solution during the early universe. In an accompanying paper we have investigated such models further [47], building on the results presented here.
When studying linear perturbations in a cosmological setting, the number of different gauges and conventions for perturbative quantities one can employ is quite large. In order to avoid confusion, we will use this appendix to set up the notation we are using throughout this paper. We will work within a manifestly gauge-invariant approach and use the first part of this section to specify our definitions, setting up Einstein's equations and the equations of energy-and momentum-conservation for a generic cosmology containing several non-minimally coupled fluids. In a second part we will turn to the perturbations for two canonical scalar fields with a common potential and show how they can be mapped onto the fluid description.

a. Metric perturbations
The most generic linearly perturbed metric around a FLRW-background contains scalar, tensor and vector modes (see e.g. refs. [60,61]). Since we are investigating a theory of scalar fields, we will restrict our attention to the scalar perturbations. With this constraint the line element of the linearly perturbed FLRW-metric reads where η denotes conformal time, defined by η(t) ≡ t t0 dt a(t ) with t 0 being some suitable initial cosmic time. Under a gauge-transformation η → η + ξ 0 (η, x) , where denotes a derivative with respect to conformal time. Note that we have ignored the vector component of gauge transformation acting on the spatial coordinates, since it does not affect scalar perturbations. In the subsequent calculations we will always resort to the commonly used gauge-invariant Bardeen potentials [62], defined as where σ is given by σ ≡ −G + F , i.e. σ → σ − ξ 0 under the above gauge transformation.

b. Fluid perturbations
For any (non-minimally coupled) component of the cosmic fluid the most generic form of the linearly perturbed energy momentum tensor reads Here ρ α , p α and u α have the usual meanings of the total (i.e. background and perturbed) energy density, pressure density and four-velocity, respectively. π µ ν is the anisotropic stress tensor, which can be decomposed into trace free scalar part, a vector part and a tensor part. Again ignoring the vector and tensor contributions we end up with the following expression for the spatial components while all other components vanish. Writing the irrotational part (again neglecting vector modes) of the 4velocity u µ as we obtain the scalar part of the stress energy tensor to linear order: Here we have used a bar (likeρ) to mark background quantities. The anisotropic stress tensor π i α j is already gauge-invariant, whereas scalar quantities (like δρ α or δp α ) transform like under a generic gauge transformation and the velocity potential obeys We will therefore define a gauge-invariant version for scalar quantities (i.e. δρ α and δp α ) via and also define a gauge-invariant momentumperturbation via These definitions are of course not unique, there is a huge number of gauge invariant metric and fluid quantities one could employ (see e.g. ref. [61] for a quite comprehensive overview). We choose to work with the definitions that correspond to the standard quantities used when working in the longitudinal gauge (defined by F = σ = 0). In the following sections the superscript "(gi)" will be dropped since we will be using gauge-invariant quantities exclusively. We will also stop employing the bar to denote the background quantities, because the total densities will not appear in any of the equations anymore.

Equations of motion a. Einsteins Equations
At the background level Einstein's equations M 2 G µν = T µν give rise to the usual Friedmann equations, which read At linear order around a FLRW-background one obtains the usual equations for the Fourier modes of the two Bardeen potentials: In accordance with common practice we do not explicitly mark the Fourier modes in the perturbative quantities. The two Bardeen potential are related by Here M denotes the reduced Planck mass.

b. Energy-and momentum (non-)conservation
At the background level the equations of energymomentum conservation D µ T µν gives rise to total energy conservation, which reads ρ tot + 3h(ρ tot + p tot ) = 0 .
(A21) However, in a multicomponent cosmology each component obeys its own energy (non-)conservation equation, which reads where the sum of all coupling terms is of course subject to the constraint from equation (A21), which gives At the background level the assumption of spatial homogeneity and isotropy require that and we define a dimensionless version of the coupling via Evaluating equation (A22) then simply gives In order to write down the generic perturbed equations of energy-and momentum-conservation, we first have to specify the generic forms for energy-and momentum transfer at the linear level. This is done by perturbing the coupling four-vector Q, which reads Note that we have again dropped a possible vector contribution to the spatial part of the perturbed coupling fourvector. While f α is already gauge-invariant, α transforms as α → α − Q α Qα ξ 0 . We define a gauge-invariant momentum transfer via The equation of energy conservation then reads at the linear level where we have used the dimensionless coupling defined in equation (A25). The equation of momentum conservation gives 3. New variables a. Generic equations Following ref. [63] we now introduce new variables for all perturbative fluid quantities, which are particularly well suited to investigate the perturbation modes in the early universe in our model. We start with a redefined density contrast and velocity potential given by The equations of motion (A30) and (A31) can now be rewritten in terms of these new variables. This gives: where we have introduced the adiabatic sound speed and the internal entropy perturbation (see next section for more details), given by as well as a dimensionless quantity for the anisotropic stressΠ Einstein's equations (A17) -(A20) can of course easily be adapted as well:

b. Equations for neutrinos, photons and baryons
At this point we can already write down the equations relevant for photons, neutrinos and baryons in the early universe, but after electron-positron annihilation. Strictly speaking, the evolution of neutrinos and photons is described in terms of a multipole-expansion of the respective phase-space distribution functions. Following ref. [34] we truncate the neutrino expansion after the quadrupole,which leaves an additional anisotropic stress contributionΠ ν in addition to ∆ ν and V ν . Since scalar fields, even when coupled, do not develop any anisotropic stress and the corresponding quantity for the baryonphoton plasma can be ignored for the era under consideration, we haveΠ tot =Π ν and the evolution of this quantity is governed by the following equation [52]: Here we closed the equation by setting all higher order moments in the Boltzmann expansion to zero. Other ways of closing the equations are more accurate for later times, but for the very early times considered here this is accurate enough [64]. The corresponding equations of energy-and momentum conservation can be read off from the generic equations (A34) and (A35), but with the substitutions ω ν = c 2 a,ν = 1/3 and Γ ν = q ν = 0 to get The situation is slightly more complicated for photons and baryons, which form a strongly coupled plasma before the era of decoupling. The full equations for a fluid description (i.e. after truncation of the photon-Boltzmann expansion after the dipole) can be obtained by using the following replacement in the generic equations: where τ c h −1 is the very short Thomson interaction timescale (τ c = (an e σ τ ) −1 ) and c 2 s,b is the (nonadiabatic) baryonic sound speed. As is common practice, we employ the tight coupling approximation (TCA), which can be found in many forms in the literature (e.g. [64,65]) to simplify the interaction terms. The basic feature of this approximation stems from the fact that the momentum exchange between photons and baryons involves the very small interaction timescale τ c , satisfying hτ c 1 in the early universe. The TCA to any given order then expresses the equations of energy-and momentum-conservation for photons and baryons to the corresponding order in hτ c . We use only the lowest order approximation here, which yields V γb ≡ V γ = V b and Π γ = 0. Redoing the calculations presented in ref. [64] to zeroth order for our variables gives the following equations for the remaining fluid-variables: where R = 4Ω γ /3Ω b 1. Note that the last equation is different from the corresponding equation presented in ref. [34], where the usual uncoupled photon equation is used. Both equations agree to leading order in an 1/R-expansion, but in the main text we will need the corresponding expansion to subleading order, and this is the correct equation to use.
A further simplification can be made by setting the baryonic sound-speed to zero, an approximation valid in the early universe [52,64].
is what is relevant for us. Interestingly, it is easy to show that (see e.g. ref. [63]) on the superhorizon scales i.e. the total curvature perturbation is constant for adiabatic modes. 5. Canonical scalar fields with a common potential a. Basic equations Let us move on to the case of two canonical scalar fields ϕ and χ with a common potential V (ϕ, χ). Variation of the scalar part of the action (A62) with respect to the scalar field yields the Klein-Gordon equations For the spatially homogeneous and isotropic background we can insert the FLRW-metric and drop all spatial derivatives to obtain the scalar field equations equations as given in equation (2). At the level of linear perturbations we can split the scalar fields into background parts and perturbations: ϕ = ϕ 0 + δϕ and χ = χ 0 + δχ. In order to remain manifestly gauge-invariant we will work with the redefined field perturbations X = δϕ − ϕ σ and Y = δχ − χ σ . (A65) Plugging this and the perturbed FLRW-metric (equation (A1)) into the Klein-Gordon equation yields the following (manifestly gauge invariant) linearized field equations: Following section III (but keeping things slightly more general) we split the potential into two parts V (ϕ, χ) = V 1 (ϕ, χ) + V 2 (ϕ, χ), where the exact definitions of V 1 and V 2 are in principle arbitrary. We define a corresponding splitting of the energy-momentum tensor by T µν ϕ = D µ ϕD µ ϕ − g µν 1 2 D τ ϕD τ ϕ + V 1 (ϕ, χ) , (A68) T µν χ = D µ χD µ χ − g µν 1 2 D τ χD τ χ + V 2 (ϕ, χ) . (A69) An evaluation at the background level (inserting the FLRW-metric and dropping spatial derivatives) yields a perfect fluid-form of this tensor with the identifications (54) and (55). At the linear level we again recover the fluid form, this time with the following assignments: δρ ϕ = 1 a 2 (ϕ X −Φϕ 2 + a 2 V 1,ϕ X + a 2 V 1,χ Y , (A70) In order to find the correct expressions for energy-and momentum-transfer we have to evaluate the equations of energy-and momentum conservation D µ T µν α = Q ν α . At the background level (and using the definition (A25)) we recover At the linear level we can use the linearized field equations (A66) and (A67) to finally obtain rather complicated expressions for the energy-and momentum transfer. Since we assume that no additional couplings are present for the scalar fields (except of course to gravity) we know that D µ T µν ϕ + D µ T µν χ = 0. From this it is easy to see that energy-and momentum transfer of the two fields are related by and we can restrict ourselves to quoting only the expressions for the perturbed cosmon-coupling here: