\section{Introduction}

Observations, such as SNe Ia \cite{sn}, WMAP \cite{cmb}, SDSS
\cite{sdss}, etc., indicate that
the expansion of the universe is accelerating at late time \cite{Frieman:08, Chen:09}.
Within the framework of general relativity, this acceleration of the universe may be explained by
supposing that the energy density of the universe is dominated by an unknown form of energy,
whose pressure is sufficiently negative, known as dark  energy.
It has been suggested that dark energy may be in the form of cosmological constant, exotic fluid, evolving scalar fields, etc.
In general, scalar field models of dark energy, such as  quintessence\cite{Wetterich:88}- \cite{Caldwell:98} or k-essence \cite{Chiba:00,Picon:00},
may be motivated by high energy physics theory.
For example, quintessence may be a scalar field with the exponential potential
$V(\phi)=V_0 e^{-\lambda \phi}$, that arises
in superstring and Kaluza-Klein theories \cite{Kanti:99}.
Here, $\lambda$ is a constant parameter and the reduced Planck mass is set to be unity.

For a single quintessence field with an exponential potential \cite{Ferreira:97,Copeland:98},
  the evolution of the field in the early epoch will converge to the scaling solution for a wide range of initial conditions
if $\lambda > \sqrt{3(1+w_f)}$, where $w_f$ is the equation of state parameter of the dominant fluid.
When the scaling solution is reached, the energy density of the field will evolve in the same way as that of the dominant fluid \cite{Ferreira:97,Copeland:98}.
Hence, although the scaling behaviour may provide a solution to the cosmic coincidence problem,
it cannot lead to an accelerated expansion of the universe.
The accelerated expansion of the universe can be achieved if the field evolves towards the field dominated solution with $\lambda < \sqrt 2$.
Unfortunately, the range $\lambda > \sqrt{3(1+w_f)}$ does not overlap with $\lambda < \sqrt 2$, so that
the accelerating solution cannot be reached from the scaling solution.

In general, there is no reason why the present accelerated expansion of the universe should be driven by a single scalar field.
It has been shown that, for multiple scalar fields with exponential potentials, although the potential of each field
is too steep to drive an accelerated expansion, the cumulative dynamics of these fields can lead to an accelerated expansion of the universe \cite{Liddle:98}.
Based on this idea, multiple scalar fields with exponential potentials, with the slope of each potential  in agreement with the prediction from particle 
physics, can drive an accelerated expansion of the universe.
This idea has been originally proposed in the assisted inflationary scenario, and has also been used in dark energy model.
This model of dark energy is assisted dark energy \cite{Kim:05,Ohashi:09},
and in this model, the scaling regime can be followed by accelerated epoch.

Based on the motivation from high energy physics, it is possible that there exists an internal structure  in the field space of multiple scalar fields, e.g., in DBI models \cite{Easson:09}.
 Some features of this internal structure, which may be represented by Kinetic interactions between  fields, have been studied using a 
simple phenomenological form of kinetic interaction in several dark energy models \cite{Chimento:08} - \cite{Saridakis:10}.
It has been shown that allowing kinetic interaction in assisted dark energy can 
have only minor effects on its background dynamics \cite{Bruck:09}.

In addition to the background dynamics, the features of dark energy also depend on the dynamics of its perturbations.
In some cases, the perturbations in dark energy can leave interesting imprints on CMB anisotropies.
For example, in the case of single scalar field dark energy, the entropy (isocurvature) perturbations in dark energy can modify the ISW effect
such that the CMB power spectrum is suppressed at low multipoles,
if the  entropy perturbations in dark energy and  the curvature perturbations are initially anti-correlated
and dark energy field is nearly frozen through out the evolution of the universe \cite{GorHu:04} - \cite{kk:07}.
In this case, the dark energy field is required to be nearly frozen otherwise the entropy perturbation in dark energy damps away
before it can affect the ISW effect.
However, we will show that, for two-field  dark energy models, it is possible for the entropy perturbations 
to survive until the present epoch for a generic evolution of both fields in some case.

The dynamics of perturbations in multiple scalar field system have been studied by many authors \cite{Malik:05} - \cite{Langlois:086}.
In many cases, the dynamics of the field perturbations are conveniently described by decomposing
the field perturbations into adiabatic and entropy modes \cite{Gordon:00} - \cite{Langlois:086}.

In this work, we are interested in the behaviour of the perturbations, in particular the entropy perturbations, in the two-field assisted dark energy model.
We use the formulation in \cite{Gordon:00, Langlois:08} to decompose the perturbations in dark energy into adiabatic and entropy modes.
    We also allow  kinetic interaction between the fields in our consideration.

This paper is organized as follows.
In section 2, we present the phenomenological Lagrangian of assisted dark energy including the kinetic interaction between fields.
In section 3, we first present the evolution equations for perturbations in general two-field dark energy model,
and then apply them to the assisted dark energy model with mixed kinetic terms.
We end with the conclusion in section 4.

\section{Assisted dark energy with mixed kinetic terms}

The general action for multiple scalar field dark energy with noncanonical kinetic terms can be written as \cite{Langlois:08}
\begin{equation}
S =  \int {\rm d}^4 x \sqrt{-g}\(\frac{1}{2} R  +  P(X,\varphi^I)\)\,,
\label{act1}
\end{equation}
where index $I = 1,2, ... N$ denotes the I-th scalar field, the reduced Planck mass is set to unity and
\begin{equation}
X=-\frac{1}{2} G_{IJ} \nabla_\mu \varphi^I \nabla^\mu \varphi^J\,.
\label{def-x}
\end{equation}
where  $G_{IJ}(\varphi^K)$ is the metric in field space and is used to raise or lower the field index $I$.
The summation convention on both field and space-time indices is implicit.
The inverse metric $G^{IJ}$ can be computed from the relation $G_{IK}G^{KJ} = \delta_I^J$.
In the following, we will simplify our consideration by supposing that $G_{IJ}$ is constant.
From (\ref{act1}), one can derive the energy-momentum tensor for dark energy \cite{Langlois:08},
\begin{equation}
T^{\mu \nu}=P\,  g^{\mu \nu}+P_{,X} \, G_{IJ}\, \nabla^{\mu} \varphi^I \nabla^{\nu} \varphi^J\,,
\label{tmn}
\end{equation}
where subscript $,_X$ denotes  derivative with respect to $X$, and the evolution equation for the scalar fields,
\begin{equation}
G_{IJ} \nabla_{\mu}(P_{,X} \nabla^{\mu} \varphi^J) + P_{,I}=0\,,
\label{kg1}
\end{equation}
where the subscript $,_I$ denotes  derivative with respect to $\phi^I$.
In the following, we will use subscript $,_A$ to denote  derivative with respect to any variable $A$,
and use subscript $,_I$ where $I$ is the field index  to denote  derivative with respect to field $\varphi^I$.

In a spatially flat FLRW universe, the evolution equation for scalar fields takes the form
\begin{equation}
\ddot \phi^I + \(3H+\frac{\dot P_{,X}}{P_{,X}}\)\dot \phi^I - \frac{1}{P_{,X}}G^{IJ} P_{,J}=0\,,
\label{kg2}
\end{equation}
where a dot denotes  derivative with respect to time $t$, $\phi^I$ is the homogeneous part of $\varphi^I$ and the Hubble parameter $H$ is given by
\begin{equation}
H^2 = \(\frac{\dot a}{a}\)^2 = \frac 13 \sum_s \rho_s\,,
\label{hub}
\end{equation}
where the index $s = r, m, d$ represents radiation, matter and dark energy respectively.
From (\ref{tmn}), the energy density of dark energy is given by
\begin{equation}
\rho_d = 2 X P_{,X} - P\,.
\label{rhod}
\end{equation}
Multiplying (\ref{kg2}) by $\dot\phi_I$, we obtain
\begin{equation}
\dot X + 6c_e^2 H X + \frac{\rho_{d,I} \dot\phi^I}{\rho_{d,X}} = 0\,,
\label{dotx}
\end{equation}
where $c_e^2 = P_{,X} / \rho_{d,X}$ is the effective sound speed of dark energy.

The important ingredient in  assisted dark energy model is that the
dark energy scales as the dominant fluid, i.e., undergoes the scaling regime, in the early epoch.
This implies that the equation of state parameter $w_d = P / \rho_d$ of dark energy is constant, and therefore $P \propto \rho_d$
in the scaling regime. From (\ref{hub}), one sees that $\rho_d \propto P \propto H^2 \propto t^{-2}$.
Hence, in the simplest case where $P_{,X}$ is constant during the scaling regime, it follows from (\ref{rhod}) that $X \propto t^{-2}$,
and consequently  $\dot\phi^I \propto 1/t$ and $\phi^I \propto \ln t$,
where $\phi^I$ are the I-th fields that contribute to the dynamics of dark energy during the scaling regime.
If we further suppose that the sound speed $c_e$ is also constant during the scaling regime,
(\ref{dotx}) shows that $\rho_{d, I} \dot\phi^I \propto t^{-3}$,  which is possible when 
the $\phi^I$-dependent parts of $P$ depend on
\begin{equation}
Q = \(\sum_J c_J{\rm e}^{\pm \lambda_J \phi^J}\)\,,
\label{p_i}
\end{equation}
where $c_J$ and $\lambda_J$ are constant.
During the scaling regime, $P_{,X}$ is supposed to be constant, so that 
we can write $P \propto X g(X, Q)$,
where $g(X,Q)$ and $X g_{,X}(X,Q)$ are constant during scaling regime.
A possible form of $g(X,Q)$ is $g(X,Q) = g(Y)$,
where $g(Y)$ is a polynomial function of $Y$, and $Y$ is the combination of $X$ and (\ref{p_i})
such that it is constant during scaling regime.
From the above consideration, one sees that dark energy can scale as the dominant fluid in the early epoch if its Lagrangian takes the form
\begin{equation}
P(X, \phi^I) = X g(Y)\,,
\label{p-assi}
\end{equation}
where, in general, $Y$ is the sum of the following expressions
\begin{eqnarray}
Y_a &=& X\(\sum_J c_J{\rm e}^{\lambda_J\phi^J}\),\qquad 
Y_b = \frac{X}{\sum_Jc_J{\rm e}^{-\lambda_J \phi^J}}, 
\nonumber\\
Y_c &=& \frac 1X \(\sum_J c_J{\rm e}^{-\lambda_J\phi^J}\), \qquad 
Y_d = \frac 1{X\(\sum_Jc_J{\rm e}^{\lambda_J \phi^J}\)}\,,
\label{y}
\end{eqnarray}
where $\lambda_J\phi^J >0$.

Let us now compute $\Omega_d$ during the radiation and matter eras,
and $w_d$ during the dark energy era.
Substituting (\ref{p-assi}) into (\ref{dotx}), we get
\begin{equation}
\hspace*{-2.5cm}
- 6 H P_{,X} X = \dot X\rho_{d,X} + \rho_{d, I} \dot\phi^I
= \dot X\(2X^2g_{,XX} + 5X g_{,X} + g\) + X\(2X g_{, X I} + g_{,I}\)\dot\phi^I\,,
\label{dotx-sol}
\end{equation}
where the index $I$ runs over the fields that contribute to the dynamics of dark energy.
It can be seen that, during radiation and matter domination, the above equation can be satisfied if $\phi^I = (2/ \lambda_I)\ln t +$ constant
and hence $\dot\phi^I = (2 / \lambda_I t)$.
Using this expression for $\dot\phi^I$, (\ref{p-assi}) and (\ref{y}), one can show that
\begin{equation}
X g_{,X} = \frac t2 g_{,I}\dot\phi^I,
\qquad {\rm and}\qquad 
X g_{,XX} = \frac t2 g_{, X I}\dot\phi^I - g_{,X}\,.
\label{gx-gi}
\end{equation}
Substituting the above equations into (\ref{dotx-sol}), we get
\begin{equation}
\Omega_d = \frac{\rho_d}{3H^2} = 3(1+w_f)\frac{P_{,X}}{\lambda_{\rm eff}^2}\,,
\label{omd}
\end{equation}
where
\begin{equation}
\frac 1{\lambda_{\rm eff}^2} = G_{I I'} \frac 1{\lambda_I \lambda_{I'}}\,.
\label{leff}
\end{equation}
In the above equation, the index $I'$ also runs over the fields that contribute to the dynamics of dark energy.
Similarly, one can show that $w_d$ during dark energy domination is given by
\begin{equation}
w_d = -1 + \frac{\lambda_{\rm eff}^2}{3P_{,X}}\,.
\label{wdf}
\end{equation}
It follows from (\ref{leff}) and (\ref{wdf}) that for suitable choices of $G_{I I'}$, 
$\lambda_{\rm eff}$ decreases and hence $w_d$ gets closer to $-1$ when more fields contribute to the dynamics of dark energy.
This shows that the multi-field dark energy described by (\ref{p-assi}) has assisted behaviour.
Note that the form  of Lagrangian in (\ref{p-assi}) does not coincide with that given in \cite{Ohashi:09}.
This is because  the internal structure in field space as in (\ref{def-x}) can give rise to the coupling terms
between the I-th field and the other fields such as $\dot\phi^I\dot\phi_{I'}$, $\(\dot\phi^{I'}\)^2 \e^{\lambda_I\phi^I}$ ($I'\neq I$) etc..
These terms are absent  in the Lagrangian in \cite{Ohashi:09}.
However, according to the above consideration, the scaling solution and assisted behaviour exist even when such terms appear in the Lagrangian.
In general,  the form of Lagrangian (\ref{p-assi}) cannot be reduced to that given in \cite{Ohashi:09} by choices in
the forms of $G_{IJ}$, $g(Y)$ or $Y$.
Nevertheless, in some cases, the both Lagrangians can become identical (see e.g. (\ref{quint})).

In the scenario for assisted dark energy, $\lambda_{\rm eff}$ must be significantly larger than unity in the early epoch because
 during radiation era $\Omega_d$ is required to be smaller than 0.045 according to BBN constraints \cite{Bean:01}.
For accelerated expansion of the universe, $\lambda_{\rm eff}$ has to decreases at late time.
This is possible when some of the subdominant fields, that had no contributions to the dynamics of dark energy
in the early epoch, can start to contribute to the dynamics of dark energy at late time.

We now show that although (\ref{p-assi}) has assisted behaviour, 
not all the possible forms of $Y$ presented above can realize the scenario for assisted dark energy.
It can be shown that (\ref{kg2}) can be written as
\begin{eqnarray}
\ddot\phi^I &=& - 3c_e^2 H \dot\phi^I
- \frac 12 \frac{\dot\phi^I}{X}\frac{c_e^2}{P_{,X}}\dot\phi^K\rho_{d, K}
- \frac 12 \frac{\dot\phi^I}{XP_{,X}}\dot\phi^K P_{,K}
+ \frac{P_{, K} G^{IK}}{P_{,X}}\,,
\label{phid1}\\
\ddot\phi^J &=& - 3c_e^2 H \dot\phi^J
- \frac 12 \frac{\dot\phi^J}{X}\frac{c_e^2}{P_{,X}}\dot\phi^K\rho_{d, K}
- \frac 12 \frac{\dot\phi^J}{XP_{,X}}\dot\phi^K P_{,K}
+ \frac{P_{, K} G^{JK}}{P_{,X}}\,.
\label{phid2}
\end{eqnarray}
Here, $\phi^I$ are the I-th dominant  fields that contribute to the dynamics of dark energy in the early epoch,
and $\phi^J$ are the J-th subdominant fields in the early epoch.
For simplicity, we set $G^{IJ} = \delta^{IJ}$,
Substituting (\ref{p-assi}) into (\ref{phid1}) and (\ref{phid2}),
one can see that, for all possible form of $Y$, all terms on the RHS of (\ref{phid1}) are proportional to $t^{-2}$ if $\phi^I$ are in the scaling regime.
Nevertheless, for some forms of $g(Y)$ and $Y$, e.g., $g(Y) = -1 + Y = -1 + Y_a$,
the fourth term on the RHS of (\ref{phid2}) decreases with time much faster the the remaining terms
and the fourth term on the RHS of (\ref{phid1}).
We will consider this case below.
Since $\phi^J$ are subdominant, $|\dot\phi^J| \ll |\dot\phi^I|$.
According to the scenario for assisted dark energy, the expansion of the universe will start to accelerate when some of $\phi^J$ begin to be dominant.
Hence, in order to drive a present accelerated expansion of the universe, $|\dot\phi^J|$ should be comparable with $|\dot\phi^I|$ about the present epoch.
It follows from (\ref{phid1}) and (\ref{phid2}) that the ratio $|\dot\phi^J| / |\dot\phi^I|$ can increase
if $\ddot\phi^J/\dot\phi^J$ are less negative than $\ddot\phi^I/\dot\phi^I$.
For convenience, let the forth term on the RHS of (\ref{phid1}) and (\ref{phid2}) be denoted by $F_I$ and $F_J$ respectively.
In the case where $\dot\phi^J > 0$ initially, $F_I$ and $F_J$ are possitive while the remaining terms in their equations are negative.
This implies that $\ddot\phi^J/\dot\phi^J$ will be less negative than $\ddot\phi^I/\dot\phi^I$
if the initial conditions for $\phi^J$ are chosen such that $F_J / \dot\phi^J > F_I /\dot\phi^I$.
Since $F_J / \dot\phi^J$ decrease with time much faster than $F_I / \dot\phi^I$
and $F_J$ decrease faster than the remaining terms in their equation,
the ratio $|\dot\phi^J| / |\dot\phi^I|$ will increase only for a short period of time
if the initial values of $F_J/\dot\phi^J$ are not significantly larger than $F_I/\dot\phi^I$.
If the initial values of $F_J/\dot\phi^J$ increase, the  rate of growth of $|\dot\phi^J| / |\dot\phi^I|$ will  increase.
Nevertheless, the period of increasing $|\dot\phi^J| / |\dot\phi^I|$ does not expand so much.
As a result, $\phi^J$ will start to be dominant in the early time when the initial value of $F_J /\dot\phi^J$ increase.
Using  similar consideration, one will see that $\phi^J$ also cannot be dominant at the right time for the case where $\dot\phi^J <0$ initially.

In the following, we restrict ourselves to the case of two fields.
In the early epoch, the evolution of entropy perturbations in this dark energy model can be
studied analytically using (\ref{p-assi}).
However,  numerical integration is required to study the effects of entropy perturbations on CMB power spectrum,
so that we have to specify the form of $P(X,\phi^I)$.
We consider both the cases of quintessence and k-essence whose forms of $P(X,\phi^I)$ are presented in the following sections.

\subsection{The quintessence with exponential potential}

If we set $g(Y) = 1 - Y = 1 - Y_c$, (\ref{p-assi}) becomes
\begin{equation}
P = X - \e^{-\lambda_1 \phi^1} - \e^{-\lambda_2\phi^2}\,.
\label{quint}
\end{equation}
Here, we have set $c_1 = c_2 =1$ for convenience.
If we choose the form of $G_{IJ}$ as
\begin{equation}
G_{IJ} = \delta_{IJ} + \(\delta_{I1}\delta_{J2} + \delta_{I2}\delta_{J1}\) \alpha\,,
\label{gij}
\end{equation}
where $|\alpha | < 1$ is constant,
the Lagrangian (\ref{quint}) becomes that considered in \cite{Bruck:09} and \cite{Ohashi:09} when $\alpha =0$.
We will use this form of $G_{IJ}$ in our numerical integration for both quintessence and k-essence models.

\subsection{The k-essence model}

We now set $g(Y) = - 1 + Y = - 1 + Y_a - M^4 Y_c$, 
where $M$ is constant, so that (\ref{p-assi}) becomes
\begin{equation}
P = - X + X^2 \(\e^{\lambda_1 \phi^1} + c_2\e^{\lambda_2\phi^2}\) - M^4\(\e^{-\lambda_1 \phi^1} + c_2 \e^{-\lambda_2\phi^2}\)\,,
\label{kesse}
\end{equation}
where we have set $c_1$ in $Y_a$ and $Y_c$ equal to unity.
The above Lagrangian has the form of the Lagrangian for dilatonic ghost condensate
model studied in \cite{Piazza:04}.
According to the previous consideration, $g(Y) = -1 + Y_a$ cannot realize the scenario for assisted dark energy
because the fourth term on the RHS of (\ref{phid2}) decreases too rapidly compared with the others.
Nevertheless, due to the term $M^4Y_c$ in the expression for $g(Y)$ in this section,
the fourth term on the RHS of (\ref{phid2}) does not decrease much faster than the remaining terms.
As a result, $\phi^2$ can contribute to the dynamics of dark energy at late time, i.e., the assisted scenario can be realized.

The value of $c_2/ M^4$ that can realize the scenario for assisted dark energy,
depends on the value of $M^4$, the starting time  for the numerical integration
and the contribution from $\phi^2$ to the dynamics of dark energy at initial stage,
i.e., the initial value of $\dot\phi^2/\dot\phi^1$ and $\e^{\lambda_2\phi^2 - \lambda_1\phi^1}$.
In our consideration, we set $M^4 = 1$ and $c_2 = 10^{-40}$.

\section{The perturbations in assisted dark energy}

We will mainly perform our calculations in the longitudinal gauge \cite{Bardeen:80}, in which the line element is 
\begin{equation}
{\rm d} s^2 =- (1+2\Psi){\rm d} t^2
+a(t)^2 (1-2\Phi)\delta_{ij}{\rm d} x^i {\rm d} x^j
= g_{\mu\nu}{\rm d} x^\mu {\rm d} x^\nu\,.
\label{long}
\end{equation}
The evolution equation for the field perturbations $\delta\phi^I = \varphi^I(t,{\rm x}) - \phi^I(t)$ can be obtained by expanding (\ref{kg1})
around the FLRW background.
From (\ref{kg1}) we write
\begin{equation}
G_{IJ} g^{\mu\nu} \nabla_{\mu}(P_{,X} \nabla_{\nu} \phi^J) + P_{,I}=0\,.
\end{equation}
Expanding the field variable around its homogeneous background and substituting $g^{\mu\nu}$ from (\ref{long}) into the above equation,
we obtain the  following linearized equation for $\delta\phi^I$ in  Fourier space:
\begin{eqnarray}
P_{,X} && \[\delta\ddot\phi^I + \(3H + \frac{\dotP_{,X}}{P_{,X}}\)\delta\dot\phi^I + \frac{k^2}{a^2}\delta\phi^I\]
+ \frac 1{a^3} \(a^3\dot\phi^I\delta P_{,X}\)^{\displaystyle{\cdot}}
\nonumber\\
&&= P_{,X}\dot\phi^I\(3\dot\Phi + \dot\Psi\) + 2 G^{IJ}P_{,J}\Psi + G^{IJ}\delta P_{,J}\,,
\label{pkg}
\end{eqnarray}
where $\delta P_{,X}$ and $\delta P_{,J}$ are the perturbations in $P_{,X}$ and $P_{,J}$ respectively.
In the single field case this equation is in agreement with \cite{Christopherson:09, Felice:10}.
The evolution of the metric perturbations is governed by the perturbed Einstein equation \cite{Doran:03},
\begin{eqnarray}
\fl
3H\(\dot\Phi + H \Psi\) + \frac{k^2}{a^2}\Phi &=& - \frac 12 \sum_s\delta\rho_s\,,
\label{ein1}\\
\dot\Phi + H \Psi &=& \frac 12 a k^{-1} \sum_s\(1 + w_s\)\rho_s v_s,\,,
\label{ein2}\\
\Phi - \Psi &=& a^2 k^{-2}\sum_s \(\rho_s\Pi_s\)\,,
\label{ein3}
\end{eqnarray}
where $\delta\rho_s$, $v_s$ and $\Pi_s$ are the density perturbations, velocity perturbations and anisotropic perturbations in specie $s$ respectively.
For dark energy, $\Pi_d = 0$ and
\begin{eqnarray}
\fl
\delta\rho_d = 2\delta\(XP_{,X}\) - \delta P 
&=& P_{,X}\delta X + 2X\deP_{,X} - P_{, J}\delta\phi^J
= \rho_{d,X}\delta X + \rho_{d, J}\delta\phi^J\,,
\label{del-r}\\
(1+w_d)\rho_d v_d &=& \rho_d u_d = P_{,X} ka^{-1} \dot\phi^J\delta\phi_J\,,
\label{del-v}
\end{eqnarray}
where
\begin{eqnarray}
\delta X &=& \dot\phi^J\dot{\delta\phi}_J - 2X\Psi\,,
\label{del-x}\\
\delta P &=& P_{,X}\delta X + P_{,J}\delta\phi^J\,,
\label{del-p}\\
\deP_{,X} &=& P_{,XX}\delta X + P_{,XJ}\delta\phi^J\,.
\end{eqnarray}
Using (\ref{pkg}), (\ref{del-r}) and (\ref{del-v}), we obtain the evolution equation for $\delta_d = \delta\rho_d / \rho_d$ and $u_d$,
\begin{eqnarray}
\dot\delta_d &=& - 3H\(\frac{\delta P}{\delta\rho_d} - w_d\)\delta_d + 3\(1+w_d\)\dot\Phi - \frac ka u_d\,,
\label{dot-del}\\
\dot u_d &=& H\(3w_d - 1\) u_d + \frac ka \frac{\delta P}{\rho_d} + (1+w_d) \frac ka \Psi\,.
\label{dot-u}
\end{eqnarray}
As expected,  these equations are the evolution equations for the perturbations in a perfect fluid whose pressure perturbations are $\delta P$.
These equations can be closed if the relation between $\delta P$ and $\delta\rho_d$ is specified.
For the the Lagrangian of the form $P(X, \varphi^I)$, it is convenient to compute this relation by
decomposing the field perturbations into adiabatic and entropy  modes.

\subsection{Decomposition into adiabatic and entropy modes}

The field perturbations can be decomposed into adiabatic and entropy modes using
a set of suitable orthonormal bases in the field space.
According to \cite{Gordon:00,Langlois:08}, we choose the first basis as
\begin{equation}
\e^I = \frac{\dot\phi^I}{\sqrt{2 X}}\,.
\label{ei}
\end{equation}
The other basis can be specified by solving the equations
 \begin{equation}
{\rm s}^I G_{IJ} {\rm s}^J = {\rm s}_I {\rm s}^I = 1,
\qquad {\rm and} \qquad
\e^I G_{IJ} {\rm s}^J = \e^I {\rm s}_I = 0\,.
\label{si}
\end{equation}
Using these bases, the field perturbations can be decomposed as
\begin{equation}
\delta\phi^I = \delta\sigma\e^I + \delta S{\rm s}^I\,.
\end{equation}
 From the orthogonality of these bases, the adiabatic mode $\delta\sigma$ and entropy mode $\delta S$ can be written in terms of the field perturbations as
\begin{equation}
\delta\sigma = \e^I G_{IJ}\delta\phi^J = \e_J\delta\phi^J,,\qquad {\rm and} \qquad
\delta S = {\rm s}_J \delta\phi^J\,.
\end{equation}
If this decomposition is applied to the background field, the perturbations $\delta\phi^I$ should refer to a variation in time.
Hence, we can write
\begin{equation}
\dot\sigma =  \e_J\dot\phi^J = \sqrt{2X},
\qquad {\rm and} \qquad
\dot S = {\rm s}_J\dot\phi^J = 0\,.
\end{equation}
This shows that the background entropy field is constant, which is in agreement with the fact that the ``relative'' entropy
perturbations $\delta S$ is gauge invariant.

Multiplying (\ref{kg2}) by $\e_I$, we obtain the evolution equation for $\sigma$
\begin{equation}
\ddot\sigma = - \(3H + \frac{\dotP_{,X}}{P_{,X}}\)\dot\sigma + \frac{P_{,\sigma}}{P_{,X}}
= -3c_e^2 H\dot\sigma - \frac{\rho_{d, \sigma}}{\rho_{d, X}}\,,
\end{equation}
where $P_{,\sigma} = P_{,I}\e^I$.
On the other hand, if we multiply (\ref{kg2}) by ${\rm s}_I$, we will get a useful relation
\begin{equation}
\dot{\rm s}^I = - \frac{P_{,S}}{\dot\sigma P_{,X}}\e^I = - \dot\theta \e^I\,,
\end{equation}
where $P_{,S} = P_{,I}{\rm s}^I$.
Using the relation $\e_J{\rm s}^J = 0$, we obtain
\begin{equation}
\dot\e^I = \dot\theta{\rm s}^I\,.
\end{equation}
The above equations can be used to compute the relation between $\delta P$ and $\delta\rho_d$.
The pressure perturbations $\delta P$ in (\ref{del-p}) can be written in terms of $\delta\sigma$ and $\delta S$ as
\begin{equation}
\delta P = P_{,X}\delta X + P_{,\sigma}\delta\sigma + P_{,S}\delta S\,.
\end{equation}
If we also write $\delta\rho_d = \rho_{d, X}\delta X + \rho_{d,\sigma}\delta\sigma + \rho_{d, S}\delta S$ and use
the relation $P_{,X} = c_e^2\rho_{d, X}$, we will get
\begin{equation}
\delta P = c_e^2 \delta\rho_d + \(P_{,\sigma} - c_e^2\rho_{d, \sigma}\)\delta\sigma
+ \(P_{,S} - c_e^2\rho_{d, S}\)\delta S\,.
\label{dp0}
\end{equation}
Since $\dot S = 0$, $P_{,\sigma} - c_e^2\rho_{d, \sigma} = \(\dot P - c_e^2\dot\rho_d\)/\dot\sigma$.
Using the definition of the adiabatic sound speed $c_a^2 = \dot P / \dot\rho_d$, and 
the conservation equation for $\rho_d$, $\dot\rho_d = - 3H\(1+w_d\)\rho_d$,
the second term on the RHS of the above equation becomes
\begin{equation}
3H\(c_e^2 - c_a^2\)P_{,X}\dot\sigma\delta\sigma = 3 H \frac ak\(c_e^2 - c_a^2\)\rho_d u_d\,.
\end{equation}
Here, we have written $u_d$ in terms of the adiabatic field as
\begin{equation}
\rho_d u_d = P_{,X} ka^{-1} \dot\phi^J\delta\phi_J =  P_{,X} ka^{-1} \dot\sigma\delta\sigma\,.
\end{equation}
From (\ref{rhod}), we can write the third term on the RHS of (\ref{dp0}) in terms of the pressure only.
Finally, we have
\begin{equation}
\delta P = c_e^2\delta\rho_d + 3 \frac ak H \(c_e^2 - c_a^2\)\rho_d u_d + \[P_{,S}\(1 + c_e^2\) - 2c_e^2 X P_{,XS}\]\delta S\,.
\label{dp-ds}
\end{equation}
If $\delta S = 0$, the form of this relation will be the same as the one for single scalar field derived in \cite{GorHu:04}.
This relation implies that (\ref{dot-del}) and (\ref{dot-u}) can be closed if the evolution of $\delta S$ is known.

In order to derive the evolution equation for $\delta S$, the following relations \cite{Langlois:08} are required,
\begin{eqnarray}
\dot P_{,S} &=& -\frac{P_{,S} P_{,\sigma}}{\dot\sigma P_{,X}}+P_{,\sigma S}\dot\sigma+P_{,XS}\dot\sigma\ddot\sigma\,,
\\
P_{,\sigma S} &=&  P_{,IJ}\e^I{\rm s}^J, \qquad {\rm and} \qquad P_{,XS} = P_{,X J}{\rm s}^j\,.
\end{eqnarray}
The evolution equation for $\delta S$ is computed by multiplying (\ref{pkg}) with ${\rm s}_I$.
The result is
\begin{equation}
\ddot{\delta S} + \(3H + \frac{\dot{P_{,X}}}{P_{,X}}\)\dot{\delta S} + \(\frac{k^2}{a^2} + \mu_S^2 + \frac{\Xi^2}{c_e^2}\)\delta S
= - \frac{\Xi}{\dot\sigma P_{,X}}\delta\rho_{dv}\,,
\label{dot-s}
\end{equation}
where
\begin{eqnarray}
\mu_S^2 &\equiv & -\frac{P_{,SS}}{P_{,X}} - \frac{1}{2 c_e^2 X}\frac{P_{,S}^2}{P_{,X}^2}+2\frac{P_{,XS} P_{,S}}{P_{,X}^2}\,,
\\
\Xi &\equiv & \frac{1}{\dot \sigmaP_{,X}}\[(1+c_e^2)P_{,S} - c_e^2 P_{,XS}\dot\sigma^2\]\,,
\end{eqnarray}
and $\delta\rho_{dv}$ is the density perturbations of dark energy in the dark energy rest frame, given by
\begin{equation}
\hspace*{-2.5cm}
\delta\rho_{dv} = \delta\rho_d + 3 H \frac ak \rho_d u_d
= \rho_{d, X}\(\dot\sigma\dot\delta\sigma - 2X\Psi - \ddot\sigma\delta\sigma\) - \[P_{,S} \(\frac 1{c_e^2} + 1\) - 2XP_{,XS}\]\delta S\,.
\end{equation}
In the second equality, the expression of $\delta\rho_d$ in terms of the adiabatic and entropy fields,
\begin{eqnarray*}
\hspace*{-2.5cm}
\delta\rho_d = \rho_{d,X}\delta X + \rho_{d, J}\delta\phi^J
= \rho_{d, X}\(\dot\sigma\dot\delta\sigma - 2X\Psi \) + \rho_{d, \sigma}\delta\sigma - \[P_{,S} \(\frac 1{c_e^2} + 1\) - 2XP_{,XS}\]\delta S\,,
\end{eqnarray*}
has been used.
For the early epoch, it is possible to solve
(\ref{dot-del}), (\ref{dot-u}) and (\ref{dot-s}), which completely describe the dynamics of perturbations in dark energy, analytically 
Nevertheless,  numerical  integration is required to study the behavior of perturbations in dark energy about the present epoch.
In order to solve these equations numerically, the initial conditions for $\delta_d$, $u_d$, $\delta S$ and $\dot{\delta S}$
in the early epoch must be specified.
Since the contribution from entropy perturbations $\delta S$ is small in the early era,
it is possible to use the following adiabatic initial conditions for $\delta_d$ and $u_d$ if $w_d$ is approximately constant during the initial stage.
\begin{equation}
\frac{\delta_d}{1 + w_d} = \frac{\delta_r}{1 + w_r},\qquad {\rm and}\qquad \frac{u_d}{1+ w_d} = v_r\,,
\label{adia-con}
\end{equation}
where $\delta_r = \delta\rho_r/\rho_r$ and $v_r$ is the density contrast and velocity perturbations of radiation respectively.
Supposing that dark energy is a fundamental  field which has already existed since inflationary era,
the initial magnitude of $\delta S$ may be estimated by studying the dynamics of perturbations in dark energy during inflation.
 
\subsection{The initial conditions from inflation}

Here, we would like to roughly estimate the initial magnitude of $\delta S$ by studying the evolution of perturbations in dark energy during inflation.
In the following calculation, we assume for simplicity that inflation is driven by a single canonical inflaton field, the energy density of dark energy is subdominant,
dark energy and inflaton are weakly coupled via  gravity only.
Since it is convenient to study the dynamics of perturbations during inflation using spatially flat gauge,
we write the evolution equation for perturbations (\ref{pkg}) 
in spatially flat gauge using the relations \cite{Malik:01, Riotto:02}
\begin{equation}
\Psi_{\rm flat} = \Psi_{\rm long} + \(\frac{\Phi_{\rm long}}{H}\)^{\displaystyle{\cdot}}\,,
\qquad
Q^I = \delta\phi_{\rm flat}^I = \delta\phi_{\rm long}^I + \dot\phi^I\frac{\Phi_{\rm long}}{H}\,.
\end{equation}
where subscripts ``flat'' and ``long'' denote spatially flat and longitudinal gauge respectively.
After performing a transformation , the metric perturbations $\Psi_{\rm flat}$ can be eliminated using the relation \cite{Gordon:00}
\begin{equation}
\Psi_{\rm flat} = \frac{1}{2H}\dot\psi Q_f\,,
\label{psiflat}
\end{equation}
where $\psi$ is the background inflaton field
and $Q_f$ is the perturbed inflaton field in spatially flat gauge.
It is not hard to show that (\ref{pkg}) in spatially flat gauge
is compatible with the action
\begin{eqnarray}
\hspace*{-2cm}
S^{(2)} &=& \frac{1}{2} \int {\rm d}t\, {\rm d}^3{\rm k} \,a^3
\[\(P_{,X} G_{IJ} + P_{,XX} \dot\phi_I \dot\phi_J \)\dot Q^I \dot Q^J
-\frac{P_{,X} k^2}{a^2}G_{IJ}Q^I Q^J 
\right.
\nonumber\\
\hspace*{-2cm}
&&  
  \left.
  + P_{,IJ}Q^I Q^J - C_K Q_f Q^K
- D_K\dot Q_f Q^K
+ 2 P_{,XJ} \dot\phi_I Q^J \dot Q^I\]\,,
\label{flat1}
\end{eqnarray}
whose form is similar to that in \cite{Langlois:08},
except for the third, fourth and fifth terms in the action.
This is because the metric perturbations $\Psi_{\rm flat}$ in (\ref{psiflat})
are sourced by the inflaton field not by $Q^I$ in our case.
Thus, to simplify our task, we follow the calculation in \cite{Langlois:08}.
Here, inflaton and dark energy are coupled via the perturbed metric, and the coupling coefficients are given by 
\begin{equation}
\hspace*{-2cm}
C_K = \frac{2X}{H} P_{,XK}\dot\psi
-\frac{1}{a^3}\[\frac{a^3}{H}P_{,X}\(1+\frac{1}{c_e^2}\)\dot\psi\dot\phi_K\]^{\displaystyle{\cdot}}\,
\qquad
D_K =\frac{P_{,X}}{H}\(1 - \frac 1{c_e^2}\)\dot\phi_K\dot\psi\,.
\label{c-inflaton}
\end{equation}
The action (\ref{flat1}) can be written in terms of the adiabatic and entropy fields as
\begin{eqnarray}
\hspace*{-2cm}
 S^{(2)} 
&=&   \frac{1}{2} \int  {\rm d}t\, {\rm d}^3{\rm k} \,a^3\
\[ \frac{P_{,X}}{c_e^2}\(\dotQ_{\sigma} - \dot\theta Q_{S}\)^2 +P_{,X} \(\dot Q_{S} + \dot\thetaQ_{\sigma}\)^2
- \frac{P_{,X} k^2}{a^2}Q_MQ^N
 \right.\nonumber\\
\hspace*{-2cm}
  &+&\left. 
P_{, MN} Q^MQ^N - C_M Q_f Q^M - D_M \dot{Q}_f Q^M
+2\dot\sigma \(\dotQ_{\sigma} - \dot\theta Q_{S}\)P_{,XM}Q^M\]\,,
\label{flat2}
\end{eqnarray}
where indices $M, N$ run over $\sigma, S$.
Using the new fields which are canonically normalized,
\begin{equation}
v_{\sigma}=\frac{a \sqrt{P_{,X}}}{c_e}\, Q_{\sigma} = \beta Q_{\sigma}\,,\qquad \,v_S = a\,\sqrt{P_{,X}}\, Q_{S} = \gammaQ_{S}\,,
\label{v}
\end{equation}
together with conformal time $\tau = \int d a /a$, the second-order action (\ref{flat2}) can be rewritten in the form
\begin{eqnarray}
\label{flat3}
\hspace*{-1.9cm}
S^{(2)}=\frac{1}{2}\int {\rm d}\tau {\rm d}^3{\rm k} 
&&\[ 
  v_\sigma^{\prime\, 2}+ v_S^{\prime\, 2} -2\xi v_\sigma^\prime v_S-k^2 c_e^2 v_\sigma^2 -k^2 v_S^2 
\right. \\
\hspace*{-1.9cm}
&& \left.
+\Omega_{\sigma\sig}v_\sigma^2+\Omega_{SS} v_S^2+2\Omega_{\sigma S}v_\sigma v_S
+ \bar{C}_\sigma Q_f v_\sigma + \bar{C}_S Q_f v_S
+ \bar{D}_\sigma Q'_f v_\sigma
\]\,,
\nonumber
\end{eqnarray}
where a prime denotes  derivative with respect to conformal time, and
\begin{eqnarray}
\hspace*{-2.7cm}
\xi&=& \frac{a}{c_e}\Xi=\frac{a}{\dot\sigma P_{,X} c_e}[(1+c_e^2)P_{,S} - c_e^2 \dot \sigma^2 P_{,XS}],
\nonumber\\
\hspace*{-2.7cm}
\Omega_{\sigma\sig} &=& \frac{\beta''}{\beta} + a^2c_e^2\dot\theta^2 + \frac{a^2c_e^2}{P_{,X}}\(P_{,\sigma\sig} - \(a^3 \dot\sigma P_{,X\sigma}\)^{\displaystyle{\cdot}}\)\,,
\quad  
\Omega_{\sigma S}= \(\frac{\beta'}{\beta} + a\frac{\ddot\sigma}{\dot\sigma}\)\xi\,, 
\quad
\Omega_{SS}=\frac{\gamma''}{\gamma}-a^2 \mu_S^2\,,
\nonumber\\
\hspace*{-2.7cm}
\bar{C}_\sigma &=&  \frac{a^3 c_e}{\sqrt{P_{,X}}}\[
\[\(1+\frac{1}{c_e^2}\)\frac{\dot\psi}{H}\]^{\displaystyle{\cdot}} P_{,\sigma}
- \frac{2X}{H} P_{,X\sigma}\dot\psi\],
\quad
\bar{C}_S =
\frac{a^2\dot\sigma \sqrt{P_{,X}}}{c_e}\frac{\dot\psi}{H}\xi,
\nonumber\\
\hspace*{-2.7cm}
\bar{D}_\sigma &=& \frac{a^2\dot\sigma}{P_{,X}}\(\frac{1}{c_e} - c_e\)\frac{\dot\psi}{H}\,.
\end{eqnarray}
The evolution equation for entropy mode can be derived from (\ref{flat3}) as
\begin{equation}
\hspace*{-1cm}
v_S''+\xi  v_\sigma'+\(k^2- \frac{\gamma''}{\gamma}+a^2\mu_S^2\) v_S - \(\frac{\beta'}{\beta} + a\frac{\ddot\sigma}{\dot\sigma}\) \xi v_\sigma = \frac{a^2\dot\sigma \sqrt{P_{,X}}}{2c_e}\frac{\dot\psi}{H}\xi Q_f\,.
\label{dvs}
\end{equation}
In order to solve this equation, the terms that are propportional to $\xi$ and $\mu_S^2$ are supposed to be negligible, i.e.,
$|\xi | \ll aH$ and $\mu_S^2 \ll H^2$, which is possible when one dark energy field is subdominant compared with the other.
We will discuss this issue in the next section.
According to the scenario for assisted dark energy, the field, which undergoes scaling regime during radiation era,
has larger contribution to the energy density as compared to the other before it enters the scaling regime.
Since, before the scaling regime starts, the subdominant field tends to be nearly frozen \cite{Ohashi:09}
while the dominant one can undergo a long period of kinetic regime (see figure (\ref{fig:1})),
the contribution of the dominant field to the energy density of dark energy does not tend to decrease when looking back in time.
Therefore, it is possible that the field that is dominant during radiation era is also dominant during inflation.
 If $|\xi | \ll a H$, $\mu_S^2\ll H^2$ and the dark energy fields slowly evolve, i.e., $P_{,X}$ and $c_e$
evolve slower than the scale factor, the approximate solution of (\ref{dvs}) is
\begin{equation}
v_S(k) \approx  \frac{1}{\sqrt{2k}}e^{-ik \tau }\(1-\frac{i}{ k \tau}\)\,.
\label{vsk}
\end{equation}
For a slow rolling inflaton, it can be shown that \cite{Riotto:02}
\begin{equation}
|Q_f| \approx \frac 1{a\sqrt{2 k}}\,,
\end{equation}
on large scales, so that
\begin{equation}
\hspace*{-2cm}
||\delta S || = ||Q_{S} || \simeq \frac 1{\sqrt{P_{,X}}}||Q_f || = \frac 1{\sqrt{P_{,X}}}\frac{\dot\psi}{H}||{\cal R}|| = \frac 1{\sqrt{ P_{,X}}} \sqrt{2\epsilon} ||{\cal R}||
\simeq \frac 1{\sqrt{ P_{,X}}} \sqrt{0.125 R} ||{\cal R}||\,,
\label{init-ds}
\end{equation}
where ${\cal R}$ is the curvature perturbations on the comoving hypersurfaces, $\epsilon$ is the slow roll parameter of inflaton,
$R$ is the relative amplitude  of the tensor to  scalar perturbations,
and $||\cdots ||$ denotes $\sqrt{\langle \(\cdots\)^2\rangle}$ and $|\cdots |$
when $Q_{S}$ and $Q_f$ are uncorrelated and fully correlated respectively.
Here, $\langle\cdots\rangle$ refers to the ensemble average.
In the second equality, we have used the fact that $\delta S$ is gauge invariant.
It follows from this equation that the magnitude of $\delta S$ can be large if $|P_{,X} |$ is small,
which is possible when dark energy slowly evolves in many k-essence models,
e.g., in the case of Lagrangian in (\ref{kesse}),
$P_{,X} \rightarrow 0$ when $w_d \rightarrow -1$ if $M^4 = 1$.
This result can be used to estimate the initial magnitude of $\delta S$ during the radiation era,
because, as we will see in the next section, $\delta S$ tends to be constant for generic evolution of dark energy as long as one of the scalar fields is subdominant.

Since we use the Lagrangian of the form $P(X, \varphi^I)$ for dark energy in the above calculation,
the adiabatic and entropy mode of the dark energy perturbations propagates with the speed of sound $c_e$ and speed of light respectively.
However, in more general case in which the Lagrangian takes the form $P(X^{IJ}, \varphi^K)$ \cite{Langlois:086}
where $X^{IJ}= -(1/2)\nabla_\mu \varphi^I \nabla^\mu \varphi^J$,
the propagation speed of adiabatic and entropy modes can be different from the ones considered here.
For example,  in the case of DBI Lagrangian \cite{Langlois:084, Langlois:086},
both adiabatic and entropy modes propagate with the speed of sound, so that
the gradient term in the second order perturbed action is multiplied by sound speed.
In this case, the amplitude of the entropy field $\delta S$ can be amplified compared with the adiabatic field if $|c_e| < 1$.
However, in our case, it follows from (\ref{init-ds}) that the magnitude of $\delta S$ can be enhanced compared with the curvature perturbation 
when $|P_{,X} | < 1$, because, in order to obtain the canonically normalized field variable,
we rescale the entropy field $\delta S$ by the scale factor and square root of is kinetic coefficient $\sqrt{P_{,X}}$,
and the curvature perturbation is generated by canonical inflaton field.

In general, it is not convenient to use the value of $Q_{\sigma}$ during inflation to define the initial conditions for $\delta\sigma$ in the radiation era,
because the evolution of $Q_{\sigma}$ depends on the evolution of dark energy after inflation.
Hence, we define the initial conditions for the $\delta\sigma$ and also $\dot{\delta\sigma}$ through 
the adiabatic conditions for the perturbations in the total energy density and velocity given in (\ref{adia-con}).
However, strictly speaking, (\ref{adia-con}) can be justified if the curvature perturbation and the perturbations in dark energy have correlation.
Nevertheless, in the above calculation, the perturbations in dark energy and curvature perturbation are uncorrelated
because we suppose that $|\xi| \ll a H$, $\mu_S^2\ll H^2$ and the inflaton slowly evolves.
Several models for generating correlated dark energy perturbations have been proposed \cite{GorHu:04, Moroi:04}.
To simplify our task, instead of considering these models in detail, we use the model in \cite{Moroi:04} to estimate
how (\ref{init-ds}) is modified when the perturbations in dark energy and curvature perturbation are correlated.
Based on the model in \cite{Moroi:04}, the correlated dark energy perturbations can be generated
if the curvaton field is coupled with the dark energy fields such as $g_I\varphi^I \zeta$
where $\zeta$ is the curvaton field and $g_I$ is the coupling coefficient.
In the case where $|\xi| \ll a H$, $\mu_S^2\ll H^2$ and $P_{,X}$ together with $c_e$ evolves slower than the scale factor,
the dynamics of $Q_{S}$ and $Q_{\sigma}$ should not be much different from the dynamics of perturbations in dark energy field in  \cite{Moroi:04}.
Hence, comparing the result in \cite{Moroi:04} with (\ref{init-ds}),
one expects that (\ref{init-ds}) are modified only by multiplication factor whose magnitude is smaller than unity,
such as $|\delta S| < \sqrt{0.125 R}|{\cal R}|/\sqrt{P_{,X}}$,
 when the entropy perturbations in dark energy and curvature perturbation are correlated.
Thus, in the following consideration, we assume that the perturbations in dark energy and curvature perturbation are fully correlated
and use the RHS of the last equality in (\ref{init-ds}) to give a rough estimate of the upper bound of $|\delta S|$.


\subsection{Evolution of the entropy perturbations	 during the radiation and matter epochs}

In order to study the evolution of the entropy perturbations,
the form of ${\rm s}^I$ should be known.
In the case of two dimensional field space, applying the constraint equations for basis ${\rm s}^I$ in (\ref{si}) to a symmetric $G_{IJ}$ yields
\begin{equation}
{\rm s}^1 = \mp \frac{G_{22}\e^2 + G_{12}\e^1}{\sqrt{\det G}}
\qquad {\rm and} \qquad 
{\rm s}^2 = \pm \frac{G_{11}\e^1 + G_{12}\e^2}{\sqrt{\det G}}\,,
\label{si-gen}
\end{equation}
where $\det G > 0$.
Obviously, there are two possible choices of ${\rm s}^I$ in the two dimensional field space.
These choices correspond to two ${\rm s}^I$ pointing in opposite directions.
Nevertheless, these different choices of ${\rm s}^I$ do not lead to  different dynamics of perturbations in dark energy.
This is because $\delta S$ and the derivative $\partial^n /\partial S^n$, where $n$ is an odd integer, change sign for these different choices of ${\rm s}^I$,
so that the evolution equation for $\delta S$ in (\ref{dot-s}), the relation between $\delta P$ and $\delta\rho_d$
in (\ref{dp-ds}) and the expression for $\delta\rho_d$ in terms of adiabatic and entropy modes are the same for these different choices.
However, it follows from (\ref{init-ds}) that
these different choices of ${\rm s}^I$ imply a change in the relative sign between $\delta S$ and the curvature perturbations.
Nevertheless, this relative sign cannot be computed explicitly from the choices of ${\rm s}^I$.
In our consideration, we choose ${\rm s}^1 = - (G_{22}\e^2 + G_{12}\e^1)/\sqrt{\det G}$
and ${\rm s}^2 = (G_{11}\e^1 + G_{12}\e^2)/\sqrt{\det G}$,
and set the ratio $\delta S / {\cal R}$ to be negative initially,
i.e., $\delta S$ and the curvature perturbation are anti-correlated.
We set $\delta S / {\cal R}< 0$ initially because it leads to
the suppression of the CMB power spectrum at low multipoles for our choice of Lagrangians (\ref{quint})
and (\ref{kesse}) in our numerical integration.
However, in general, initial $\delta S / {\cal R}< 0$  does not necessarily lead to the suppression of the 
CMB power spectrum at low multipoles.

In the radiation and matter epochs, when the field $\phi^2$ is subdominant compared with $\phi^1$,
the magnitude of ${\rm s}^I$ can be estimated by solving (\ref{kg2}).
Rewriting (\ref{kg2}) as
\begin{equation}
\frac{d}{d t}\(a^3P_{,X} \dot\phi^I\) = a^3 G^{IJ}P_{, J}\,.
\end{equation}
When $\phi^2$ is subdominant, $|P_{,1}| \gg |P_{,2}|$, and hence
\begin{equation}
\hspace*{-1.5cm}
\dot\phi^1 = \frac{G^{11}}{a^3P_{,X}}\int dt \(a^3 P_{,1}\) + \frac{c_1}{a^3P_{,X}}\,,
\qquad
\dot\phi^2 = \frac{G^{12}}{a^3P_{,X}}\int dt \(a^3 P_{,1}\) + \frac{c_2}{a^3P_{,X}}\,,
\end{equation}
where $c_1$ and $c_2$ are the constants of integration.
This result implies that
\begin{equation}
\dot\phi^2 = - \frac{G_{12}}{G_{22}}\dot\phi^1 + \frac 1{a^3P_{,X}}\(c_2 + \frac{G_{12}}{G_{22}}c_1\)\,.
\label{p2p1}
\end{equation}
Here, we have used
\begin{equation}
G^{12} =  - \frac{G_{12}}{\det G}
\quad
{\rm and}
\quad
G^{11} = \frac{G_{22}}{\det G}\,.
\end{equation} 
Substituting (\ref{p2p1}) into (\ref{def-x}), we get
\begin{equation}
X = - \frac 12 \(\dot\phi^1\)^2 \(G_{11} - \frac{G_{12}^2}{G_{22}}\) - \frac{G_{22}}2 \[\frac 1{a^3P_{,X}}\(c_2 + \frac{G_{12}}{G_{22}}\)\]^2\,.
\label{x-ear}
\end{equation}
In the case of our interest where $-1 < w_d <1$, the first term on the RHS of this equation will not decrease faster than the second term,
so that the second term will be negligible compared with the first term as the universe evolves.
Thus, we ignore this term in the following calculation.
Substituting (\ref{p2p1}) into (\ref{si-gen}), we obtain
\begin{equation}
|{\rm s}^1| \ll |{\rm s}^2| \simeq \left |\frac{\sqrt{\det G}}{G_{22}}\e^1\right |\,,
\qquad {\rm where}\qquad 
|\e^1| \simeq \left |\sqrt{\frac{G_{22}}{\det G}}\right | \gg |\e^2|\,.
\label{s1ll}
\end{equation}
Here, the expression for $|\e^1|$ is computed using (\ref{x-ear}) and (\ref{ei}).
In the following calculations, we suppose that $|{\rm s}^2| \sim |\e^1| \sim {\cal O}(1)$.

Let us now come to the case of assisted dark energy whose Lagrangian is given by (\ref{p-assi}).
Using this Lagrangian, we estimate the magnitudes of $\mu_S^2$ and $\Xi$ which
correspond to the amplitude of the ``effective mass'' of $\delta S$ and the coupling between $\delta S$ and $\delta\rho_d$.
From (\ref{gx-gi}),
one can derive the following relations in the case where $\phi^2$ is subdominant:
\begin{equation}
|\lambda_1| X g_{, X} \simeq g_{,1}\,,
\qquad {\rm and} \qquad
|\lambda_1 g_{,1}| \sim |g_{,11}| \,.
\label{xgx2g1}
\end{equation}
Using these relations, $|P_{,1}| \gg |P_{,2}|$ and (\ref{s1ll}), if $\lambda_1$ is not very much larger than unity, it can be seen that
\begin{eqnarray}
\hspace*{-2.2cm}
X |P_{,XS}| &=& X |\(g + X g_{,X}\)_{,I}{\rm s}^I| \simeq X \left |\(g + \frac{g_{,1}}{\lambda_1}\)_{,I}{\rm s}^I\right | \sim |P_{,S}|\,.
\label{pxs-ps}\\
\hspace*{-2.2cm}
|P_{,SS}| &=& |{\rm s}^1{\rm s}^1 P_{,11} + 2{\rm s}^1{\rm s}^2 P_{,12} + {\rm s}^2{\rm s}^2 P_{,22}| = |P_{,11}B_{SS}| \sim |\lambda_1 P_{,1}B_{SS}| \simeq |\lambda_1 P_{,\sigma}B_{SS}|,
\label{pss-psig}
\end{eqnarray}
where $B_{SS} \equiv {\rm s}^1{\rm s}^1 + 2{\rm s}^1{\rm s}^2P_{,12}/P_{,11} + {\rm s}^2{\rm s}^2P_{,22}/P_{,11}$
so that $|B_{SS}|\ll 1$ when $\phi^2$ is subdominant.
The relation (\ref{pxs-ps}) is valid when $P_{,XS} \neq 0$.
From $P_{,1}\gg P_{,2} and $(\ref{s1ll}), it is easy to see that
\begin{equation}
|P_{,S} | = |P_{,1}{\rm s}^1 + P_{,2}{\rm s}^2| = |P_{,1}B_S| \simeq |P_{,\sigma}B_S|\,,
\label{ps-psig}
\end{equation}
where $ B_S \equiv {\rm s}^1 + {\rm s}^2P_{,2}/P_{,1}$
so that $|B_S|\ll 1$ when $\phi^2$ is subdominant.
Using (\ref{xgx2g1}), (\ref{p-assi}) and $(1+w_d)\rho_d = 2XP_{,X}$, we get
\begin{equation}
\hspace*{-2cm}
\frac{P_{,\sigma}}{\rho_d} = \frac{X g_{,\sigma}}{2X^2 g_{,X} + X g} \simeq \frac{\lambda_1 X g_{,X}}{2 X g_{,X} + g} \simeq (1 - w_d)\frac{\lambda_1}2\,,
\quad {\rm and}\quad
\frac{P_{,\sigma}}{XP_{,X}} \simeq \lambda_1\frac{1 - w_d}{1 + w_d}\,.
\label{psig-rd}
\end{equation}
Hence, $|P_{,S} |/ \rho_d \simeq \lambda_1(1 - w_d)|B_S|/2 \ll 1$,
implying that the contribution of $\delta S$ to $\delta P$ is negligible in the early epoch.
In the early epoch, when the contribution from $\phi^2$ to $\rho_d$ is significantly small such that $|B_{SS}| \ll 1+w_d$
and $|B_S| \ll 1+w_d$, it follows from (\ref{pss-psig}), (\ref{ps-psig}) and (\ref{psig-rd}) that
\begin{eqnarray}
\left |\frac{P_{,S}}{\dot\sigma P_{,X} \sqrt{\rho_d}}\right | &\lesssim & \left |\frac{P_{,\sigma}B_S}{X P_{,X}}\right | 
\simeq (1 - w_d)\frac{|\lambda_1 B_S|}{1+w_d} \ll 1\,,
\label{ps-px}\\
\left |\frac{P_{,SS}}{\rho_dP_{,X}}\right | &\lesssim & \left |\frac{\lambda_1P_{,\sigma}B_{SS}}{XP_{,X}}\right |
 \simeq (1 - w_d)\frac{\lambda_1^2 |B_{SS}|}{1+w_d} \ll 1\,.
\label{pss-px}
\end{eqnarray}
In the above calculations, we have used $|\dot\sigma | \lesssim \sqrt{\rho_d}$ instead of $|\dot\sigma| \sim \sqrt{\rho_d}$,
because $|\dot\sigma|$ can be much smaller than $\sqrt{\rho_d}$ when $w_d\rightarrow -1$ in some cases, e.g.,
in the case of quintessence models where $P = X - V$.
Using (\ref{pxs-ps}), (\ref{ps-px}), (\ref{pss-px})
 and $\sqrt{\rho_d} < H$ together with $c_e^2 \sim {\cal O}(1)$,
it can be shown that $|\Xi | \ll H$ and $\mu_S^2 \ll H^2$.

When $|\Xi | \ll H$ and $\mu_S^2 \ll H^2$, the evolution equation for $\delta S$ (\ref{dot-s}) on super horizon scales becomes
\begin{equation}
\frac{d}{d t} \(a^3P_{,X}\dot{\delta S}\) \approx 0\,.
\end{equation}
The solution of this equation is
\begin{equation}
\dot{\delta S} = \frac{d_1}{a^3P_{,X}}\,,
\qquad
\delta S = \int \frac{d_1 dt}{a^3P_{,X}} + {\rm const.} \sim - \frac{\dot{\delta S}}{H} + {\rm const}\,,
\label{ds-long}
\end{equation}
where $d_1$ is the constant of integration.
To evaluate the integral in the above expression  for $\delta S$, $P_{,X}$ is assumed to be approximately constant.
Since, according to the above equations, $|\dot{\delta S} / H|$ decreases with time,
and $|\dot{\delta S} / H| < |\delta S|$ on large scales during inflation because the dark energy fields slowly evolve,
we get $\delta S \simeq$ constant on large scales after the inflationary era.

To check the validity of the above consideration, we solve (\ref{kg2}), (\ref{dot-del}), (\ref{dot-u}) and (\ref{dot-s})
numerically using the modified version of CMBEASY \cite{cmbeasy}.
In the following, the cosmological parameters for the background are chosen such that
$\Omega_b^{(0)} = 0.044$, $\Omega_{\rm CDM}^{(0)} = 0.216$, $\Omega_d^{(0)} =0.74$, $h =0.72$.
In our numerical integration, the form of $G_{IJ}$ is chosen according to (\ref{gij}), so that
\begin{equation}
{\rm s}^1 = - \frac{\e^2 + \alpha \e^1}{\sqrt{1 - \alpha^2}} \qquad {\rm and} \qquad 
{\rm s}^2 = \frac{\e^1 + \alpha \e^2}{\sqrt{1 - \alpha^2}}\,.
\label{s12}
\end{equation}

We first consider the case where $\alpha =0$ using the quintessence model whose Lagrangian is given by (\ref{quint}).
The evolutions of $w_d = P / \rho_d$ and $w_2$, which is the equation of state parameter of subdominant field $\phi^2$,
for various initial conditions are plotted in figure (\ref{fig:1}).
There are four types of the initial conditions in figure (\ref{fig:1}).
These four types correspond to the four possible combinations of the thick and thin lines in this figure,
e.g., the thick solid together with thin long dashed lines  correspond to the case where
dark energy is initially in the potential regime $w_d \approx -1$, while $\phi^2$ is in the kinetic regime $w_2 \approx 1$ initially.
All of these initial conditions, except the case where $w_d \approx -1$ while $w_2 \approx 1$ initially,
 lead to the same evolution of $\delta S$ in the early epoch, as shown by thin solid line in figure (\ref{fig:3}),
This evolution of $\delta S$ is also similar to that in the case of $\alpha \neq 0$ represented by thick solid line in the same figure.
In the exceptional case where $w_d \approx -1$ while $w_2 \approx 1$ initially,
the kinetic energy of $\phi^2$ can be larger than that of $\phi^1$ initially
although $\phi^2$ gives a negligible contribution to the energy density $\rho_d$, so that
$|\e^2 | > |\e^1|$ at the initial stage, and therefore the above analytical consideration is not applicable.
The evolution of $\delta S$ for this case is represented by thin long dashed line in figure (\ref{fig:3}).
This figure shows that $\delta S$ decreases until $w_2$ is close to $-1$, i.e., $|\e^1|$ becomes larger than $|\e^2|$.
Due to the source term on the RHS of (\ref{dot-s}), $\delta S$ does not decrease to zero.
Thus, this case does not imply the vanishing of $\delta S$.

We now consider the case where $\alpha \neq 0$ using the k-essence model whose Lagrangian is given by (\ref{kesse}).
The initial conditions for k-essence in figure (\ref{fig:2}) are chosen such that 
either $\dot\phi^2 = \pm \alpha \dot\phi^1$ or $\dot\phi^2 = 0$ initially,
while dark energy is either in potential regime $w_d \approx -1$ or approximately in the scaling regime $w_d\approx w_f$ at initial stage. 
In the case where dark energy starts in scaling regime,
figure (\ref{fig:2}) shows that $\dot\phi^2$ converges to $\dot\phi^2 = -\alpha\dot\phi^1$ rather quickly.
Although in the cases where $\dot\phi^2 = \alpha\dot\phi^1$ and $\dot\phi^2 = 0$ initially,
$\phi^2$ is subdominant and dark energy can be approximately in the scaling regime during the initial stage,
the relation $\dot\phi^2 \simeq -\alpha\dot\phi^1$ is not satified.
Hence, in these cases, $\delta S$ is not constant during the initial stage as shown by thick long dashed and dashed lines in figure (\ref{fig:3}).
However, after the initial transient, $\delta S$ becomes constant.
We also consider the case where $w_d \approx -1$ and $\dot\phi^2 = -\alpha \dot\phi^1$ initially.
In this case $\delta S$ is also constant in the early epoch.
The above results for $\alpha \neq 0$ also hold in the case of quintessence model.
\begin{figure}[ht]
\includegraphics[height=0.4\textwidth, width=0.9\textwidth,angle=0]{fig1.eps}
\caption{
The evolution of $w_d$ and the equation of state parameter $w_2$ of $\phi^2$
for various initial conditions in the case of quintessence model.
The tick lines represent $w_d$, while the thin lines represent $w_2$.
In this plot, $\lambda_1 = 10$, $\lambda_2 = 1$ and $\alpha = 0$.
}
\label{fig:1}
\end{figure}
\begin{figure}[ht]
\includegraphics[height=0.4\textwidth, width=0.9\textwidth,angle=0]{fig2.eps}
\caption{
The evolution of $w_d$ and $\dot\phi^2/\dot\phi^1$ for the case of k-essence model.
The thick line represents the evolution of $w_d$, while the thin lines correspond to $\dot\phi^2/\dot\phi^1$.
The tin solid, thin long dashed and thin dashed lines
correspond to $\dot\phi^2= -\alpha\dot\phi^1$, $\dot\phi^2= 0$ and $\dot\phi^2= \alpha\dot\phi^1$ at initial time respectively.
In this plot, $\lambda_1 = 50$, $\lambda_2 = 1$ and $\alpha = 0.5$.
}
\label{fig:2}
\end{figure}
\begin{figure}[ht]
\includegraphics[height=0.4\textwidth, width=0.9\textwidth,angle=0]{fig3.eps}
\caption{
The evolution of the long wavelength mode of $\delta S$ for models in figures (\ref{fig:1}) and (\ref{fig:2}).
The thin and thick lines correspond to the cases of $\alpha = 0$ and $\alpha = 0.5$ respectively.
The thin long dashed line is the case where $w_2 \approx 1$ while $w_d \approx -1$ initially.
The thick solid, thick long dashed and thick dashed lines
correspond to $\dot\phi^2= -\alpha\dot\phi^1$, $\dot\phi^2= 0$ and $\dot\phi^2= \alpha\dot\phi^1$ at initial time respectively.
}
\label{fig:3}
\end{figure}

From the above consideration, it can be concluded that in the early epoch
if $\phi^2$ is subdominant compared with $\phi^1$,
the entropy perturbations in dark energy $\delta S$ tend to be constant on large scales,
and consequently can survive until the present epoch for generic evolution of dark energy.
The constancy of the entropy perturbations is not due to the particular form of Lagrangian (\ref{p-assi}),
but is a consequence of $|\Xi | \ll H$ and $\mu_S^2 \ll H^2$,
which can generally occur in two-field system when one field is subdominant compared with another.

\subsection{The effects of entropy perturbations on the CMB power spectrum}

Having consider the evolution of $\delta S$ in the early epoch,let us now study the evolution and effects of $\delta S$ about the present era.
Since $\Xi$ and $\mu_S^2$ depend on $P_{,S}$ and $P_{,XS}$,
the evolution of $\delta S$ and the influence of $\delta S$ on $\delta_d$ are governed by these quantities.
However, since the contributions from $P_{,S}$ and $P_{,XS}$ are approximately the same, sometime we consider $P_{,S}$ only.
As have been shown, $|P_{,S} / \rho_d| \ll 1$ in the early era.
Using  similar approach one expects that $|P_{,S} / \rho_d| < 1$ in the present epoch.
This is because when $\phi^2$ is dominant compared with $\phi^1$ at present, we have $|\e^2| > |\e^1|$, $|{\rm s}^1| > |{\rm s}^2|$
and $|P_{,2}| > |P_{,1}|$ so that 
$|P_{,S} | \sim | P_{,1}{\rm s}^1 + P_{,2}{\rm s}^2| < |P_{,2}\e^2| \sim |P_{,\sigma}| \sim \rho_d$.
The ratio $|P_{,S} / \rho_d|$ can be of order of unity during the transition from $\phi^1$ - dominant to $\phi^2$ - dominant,
because $|{\rm s}^2|$ and $|P_{,2}|$ are of order of $|{\rm s}^1|$ and $|P_{,1}|$ respectively.
This rough estimate is confirmed by the numerical results in figure (\ref{fig:4}).
Hence, $\delta S$ can have a contribution to the dynamics of $\delta_d$ only during the transition stage,
because the contribution from the entropy perturbations to the pressure perturbations in (\ref{dp-ds}) depends on the magnitude of $P_{,S}$ and $P_{,XS}$.
This implies that in the general two-field dark energy models,
although the entropy perturbations in dark energy can be constant until the present epoch
because one field is subdominant compared with the other in the early era,
the entropy perturbations cannot significantly influence the dynamics of density perturbations in the universe
unless the subdominant field becomes more important at late time.

The evolution of $\delta S$ can influence the evolution of $\delta_d$ via the relation $\delta P/\rho_d \propto \tilde{P_{,S}}\delta S$ given in (\ref{dp-ds}).
Here, $\tilde{P_{,S}} = \[P_{,S}\(1 + c_e^2\) - 2c_e^2 X P_{,XS}\]/\rho_d$.
We follow CMBEASY code to normalize all perturbation variables such that ${\cal R} = -1$.
If the value of $\tilde{P_{,S}}\delta S$ is sufficiently large and possitive, $\delta S$ will give a possitive contribution to $\delta P$,
and also $\dot u_d$ in (\ref{dot-u}).
Consequently, it follows from (\ref{dot-del}) that $\delta_d$ can decrease more rapidly on large scales when $\tilde{P_{,S}}\delta S$ increases.
This qualitative consideration is independent of the form of $P(X,\phi^I)$ and
in agreement with  the numerical results in figure (\ref{fig:5}).
In contrast, if $\tilde{P_{,S}}\delta S$ is sufficiently large and negative, its contributions should increase $\dot{\delta_d}$ on large scales.

For illustration, in figures (\ref{fig:4}) and (\ref{fig:5}), we plot the evolutions of $P_{,S}/\rho_d$, $\tilde{P_{,S}}\delta S$ and $\delta_d$
from the numerical integration for k-essence model whose Lagrangian is given by (\ref{kesse}).
Since the evolutions of these quantities for our choice of quintessence and k-essence models have approximately the same features,
we consider the case of k-essence model only.

We set $\lambda_1 = 50$ and $\lambda_2 =1$ in the numerical integration.
The initial conditions for the background field are chosen such that dark energy is in the scaling regime
while $\dot\phi^2 = -\alpha\dot\phi^1$ initially.
For the perturbed quantities, the adiabatic initial conditions given in (\ref{adia-con}) are used for $\delta_d$ and $u_d$.
The initial value of $\delta S$ are chosen such that the CMB power spectrum at $\ell =2$ for $\alpha =0$ is maximally suppressed.
This initial value can be in agreement with the initial value from inflation if $|P_{,X} | < 1$ during inflation.
In our case, initial $\delta S / {\cal R} < 0$ corresponds to
$\tilde{P_{,S}}\delta S / {\cal R} < 0$ at the initial time, and leads to $\tilde{P_{,S}}\delta S / {\cal R} < 0$ during the transition stage.
In general, the relation between the sign of $\tilde{P_{,S}}\delta S$ during the transition stage and
the relative signs of $\delta S$ or $\tilde{P_{,S}}\delta S$ to ${\cal R}$ at the initial time depends on the form of  Lagrangian for dark energy.

Since the effects of entropy perturbations in dark energy are most important on large scales 
we consider the perturbation mode whose wavenumber is $k \simeq 2\times 10^{-4}$ Mpc${}^{-1}$.
This mode enters the Hubble radius about the present epoch.
\begin{figure}[ht]
\includegraphics[height=0.4\textwidth, width=0.9\textwidth,angle=0]{fig4.eps}
\caption{
The evolution of $P_{,S} / \rho_d$. 
for various $\alpha$.
The solid, long dashed and dashed lines represent $\alpha = 0, 0.5$ and 0.8 respectively.
}
\label{fig:4}
\end{figure}
\begin{figure}[ht]
\includegraphics[height=0.4\textwidth, width=0.9\textwidth,angle=0]{fig5.eps}
\caption{
The evolution of $\tilde{P_{,S}}\delta S$ (thick lines) and $\delta_d$ (thin lines).
the solid, long dashed and dashed lines correspond to $\alpha =0$, $0.5$ and $0.8$ respectively.
For comparison , the evolution of $\delta_d$ in the case where $\delta S = 0$ is plotted using dotted line.
In this plot, $k \simeq 2\times 10^{-4}$ Mpc$^{-1}$.
}
\label{fig:5}
\end{figure}

The effects of entropy perturbations on the density perturbations in the universe can be studied using the relation between
the metric perturbations and $\delta_d$ in (\ref{ein1}).
Since we are interested in the perturbations on large scales
and radiation contributes a negligible fraction of the total energy density during the transition stage,
we write (\ref{ein1}) as
\begin{equation}
\dot\Phi = - H\Phi
- \frac H2 \(\Omega_m\delta_m + \Omega_d\delta_d\)\,,
\label{dot-phi}
\end{equation}
where the subscript $m$ denotes matter.
The relation $\Phi \sim \Psi$ has been used in the above equation because radiation is negligible.
Hence, if the contributions from density perturbations is ignored for a while,
it can be read from (\ref{dot-phi}) that the amplitude of metric perturbations decreases with time due to the expansion of the universe.
The decreases of $|\Phi|$ and $|\Psi|$ which correspond to the ISW effect,
partially cancel the ordinary Sachs-Wolfe (OSW) effect.
It follows from (\ref{dot-phi}) that the inclusion of density perturbations can modify the evolution of $\Phi$,
such that if $\Phi / \delta_d > 0$, $|\Phi|$ will decrease faster, and $|\Phi|$ will decrease more slowly or increase if $\Phi /\delta_d < 0$
\cite{GorHu:04, Moroi:04, kk:07}.
When $|\Phi|$ decreases faster, the ISW contributions will be larger, leading to more cancellations between ISW and OSW effects.
In this case, the CMB power spectrum is suppressed  at low multipoles.
In contrast, OSW effect is less canceled by ISW effect when $|\Phi|$ decreases more slowly,
and OSW and ISW can add up when $|\Phi|$ increases.
Consequently, the CMB power spectrum is enhanced  at low multipoles in these cases.

Using the same model of k-essence and initial conditions as above,
we plot the evolutions of $\Omega_d\delta_d$ and $\dot\Phi$ in figure (\ref{fig:6})
and plot the CMB power spectrum in figure (\ref{fig:7}).
In our numerical integration, $\Phi$ is negative because all perturbation variables are normalized such that ${\cal R} = -1$.
Thus, in figure (\ref{fig:6}), $\dot\Phi$ increases when $\Omega_d\delta_d < 0$.
This means that due to the contribution from $\delta_d$, $|\phi|$ decreases more rapidly, i.e., the ISW effect is enhanced.
The enhancement of ISW effect can lead to the suppression of the CMB power spectrum
at low multipoles as shown in figure (\ref{fig:7}).
These are in agreement with the above qualitative consideration
\begin{figure}[ht]
\includegraphics[height=0.7\textwidth, width=0.9\textwidth,angle=0]{fig6.eps}
\caption{
The evolution of $\Omega_d\delta_d$ is shown in the upper panel, while
the evolution of $\dot\Phi$  is shown in the lower panel.
In the both panels, the solid, long dashed and dashed lines correspond to $\alpha =0$, $0.5$ and $0.8$ respectively.
The evolution of both quantities  in the case where $\delta S = 0$ is denoted by dotted line.
In this plot, $k \simeq 2\times 10^{-4}$ Mpc$^{-1}$.
}
\label{fig:6}
\end{figure}
\begin{figure}[ht]
\includegraphics[height=0.4\textwidth, width=0.9\textwidth,angle=0]{fig7.eps}
\caption{
The CMB power spectra for the k-essence model.
The solid, long dashed and dashed lines correspond to $\alpha = 0$, $0.5$ and $0.8$ respectively.
The case where $\delta S =0$ while $\alpha =0$ is denoted by dotted line.
}
\label{fig:7}
\end{figure}

From the above analysis and figures (\ref{fig:5}) - (\ref{fig:7}),
one can conclude that if the curvature perturbation ${\cal R}$ and the contribution from entropy perturbations
in dark energy, which is proportional to $\tilde{P_{,S}}\delta S$, have opposite signs after the matter era,
the ISW effect can be enhanced and consequently the CMB power spectrum can be suppressed at low multipoles.
In contrast, if ${\cal R}$ and $\tilde{P_{,S}}\delta S$ have the same signs after the matter era, the CMB power spectrum at low multipoles can be enhanced.
This is in agreement with literature \cite{GorHu:04, Moroi:04}.
However, for two-field dark energy models, the contribution from entropy perturbations in dark energy depends on $\tilde{P_{,S}}$,
i.e., the evolution of dark energy, during the transition from $\phi^1$ - dominant to $\phi^2$ - dominant.
Hence, in general, $\tilde{P_{,S}}\delta S / {\cal R} < 0$ at the initial time does not necessarily lead to
$\tilde{P_{,S}}\delta S / {\cal R} < 0$ after the matter era, and consequently the suppression of the CMB power spectrum at low multipoles.

In order to study the influence of $\alpha$ on the effect of entropy perturbations, we plot
the fractional difference in the CMB power spectra for k-essence and quintessence in figure (\ref{fig:8}).
The fractional difference in the CMB power spectra between  
the case of $\delta S = 0$ and the case of $\delta S \neq 0$ is expressed as the ratio
\begin{equation}
\Delta C_\ell = \frac{C_\ell - C_\ell^*}{C_\ell^*}\,,
\end{equation}
where $C_\ell^*$ is the power spectrum for the case of $\delta S = 0$ and
$C_\ell$ is the power spectrum for the case of $\delta S \neq 0$.

In this figure, the initial value for $\delta S$ in the case of quintessence
is also chosen such that the CMB power spectrum at $\ell =2$ for $\alpha =0$ is maximally suppressed.
This initial value of $\delta S$ is much larger than a possible initial value from inflation  given in (\ref{init-ds}),
because $P_{,X} =1$ for quintessence. If the initial value in (\ref{init-ds}) is used,
the entropy perturbations have small contribution to ISW effect and hence have no effect on CMB power spectrum.
This figure shows that the suppression of the CMB power spectrum depends on $\alpha$ in the case of k-essence,
but weakly depends on $\alpha$ in the case of quintessence.
The influence of $\alpha$ on the suppression of CMB power spectrum strongly depends on the form of $g(Y)$ in (\ref{p-assi}).
\begin{figure}[ht]
\includegraphics[height=0.4\textwidth, width=0.9\textwidth,angle=0]{fig8.eps}
\caption{
The fractional difference in the CMB power spectra for k-essence and quintessence models.
}
\label{fig:8}
\end{figure}

\section{Conclusion}

In this work, we study dynamics of entropy perturbations in two-field assisted dark energy model.
Based on the scenario for this dark energy model in which one field is subdominant compared with the other in the
early era, we show that the entropy perturbations in this two-field system tend to be constant on large scales in the early era,
and hence can survive until the present epoch, for a generic evolution of dark energy during the radiation and matter eras.
This behaviour of the entropy perturbations is preserved even when the fields are coupled via kinetic interaction.
Although in our case the constancy of the entropy perturbations is obtained using the Lagrangian which
has scaling behaviour, this result can also be obtained from more general form of Lagrangian.
This is because the constancy of the entropy perturbations is the consequence of the subdomination
of one field compared with the other in the early epoch.
However,  the entropy perturbations cannot significantly influence the dynamics of density perturbations in the universe
unless the subdominant field in the early epoch becomes more important at late time.

Since, for assisted dark energy, the subdominant field in the early epoch becomes dominant at late time,
the entropy perturbations can influence the dynamics of metric perturbations so that the ISW effect is modified.
The influences of the entropy perturbations on the ISW effect depend  on both the initial condition for $\delta S$ and the evolution of dark energy during 
the transition between $\phi^1$- domination and $\phi^2$- domination after the matter era.
Assuming correlations between the entropy and curvature perturbations,
the ISW effect will be enhanced and therefore the CMB power spectrum will be suppressed at low multipoles
if the contribution from entropy perturbations and the curvature perturbations have opposite signs,
i.e., $\tilde{P_{,S}}\delta S/{\cal R} < 0$, after the matter era.
For our choice of models in the numerical integration,
the negative  $\tilde{P_{,S}}\delta S/{\cal R}$ at initial time leads to negative $\tilde{P_{,S}}\delta S/{\cal R}$ after the matter era.
However, in general, initial $\tilde{P_{,S}}\delta S/{\cal R} < 0$ does not necessarily imply $\tilde{P_{,S}}\delta S/{\cal R} < 0$ after the matter era.

Assuming dark energy  to exist since inflationary epoch,
the initial value of $|\delta S |$ estimated during inflationary epoch cannot be large enough to influence the ISW effects in the case of quintessence.
This is because $P_{,X} =1$ for usual quintessence models.
The initial value of $|\delta S |$ can be large enough to influence the ISW effect if $|P_{,X} |$ is significantly smaller than unity
during inflationary era, which is possible in many k-essence models.

If dark energy fields are coupled via kinetic terms, the modifications of CMB power spectrum due to entropy perturbations
also depend on the coupling strength $\alpha$.
For our choice of models in the numerical integration,
this dependence is significant in the case of k-essence, but not in the case of quintessence.
The influence of $\alpha$ on the CMB power spectrum strongly depends on the form of $g(Y)$ in (\ref{p-assi}).

\section*{Acknowledgments}

The author would like to thank A. Ungkitchanukit and anonymous referee for comments on the manuscript.
This work is supported by Thailand Research Fund (TRF) through grant TRG5280036.

