Warm dark energy

Motivated by some of the recent swampland conjectures, we study the implementation for the late time acceleration of the Universe of a mechanism developed by Anber and Sorbo in the context of primordial inflation, in which an axion field can slowly roll in a steep potential due to additional friction provided by its coupling to some U(1) gauge field. We first study the realization of this mechanism in N = 2 supergravity models resulting from string compactifications on Calabi--Yau manifolds. We then study the transition between matter domination and the axion domination, and show that indeed the backreaction of the produced gauge field can sufficiently slow the motion of the axion, so to produce the present accelerated era. We finally study the transition from a pre-inflationary matter or radiation domination to primordial inflation. In the regime that we could explore numerically, the evolution is characterized by stages of faster axion roll (and consequent bursts of gauge field amplification) intermitted by stages of slower roll, with a pattern that"oscillates'' about the steady state Anber and Sorbo solution, but that does not appear to relax to it.

for some constant c of order one. In this relation prime denotes derivative with respect to the quintessence or inflaton field φ, while M p is the reduced Planck mass.
While this is a constraint that is generically violated in both inflationary and dark energy models, it has revived the quintessence scenario and called for a renewed and more careful analysis of its realizations.
Phenomenological constraints on quintessence work as upper bounds rather than lower bounds on |V /V | [11]. This creates a (small) tension between observation and theory. In this work we explore the possibility that a large slope of the quintessence potential can be compatible with the data in presence of a coupling of the quintessence with light fields.
We choose a coupling so that the motion of φ leads to an amplification of these fields. These fields are therefore produced at the expense of the kinetic energy of φ, resulting in an additional source of friction for its motion, beside the standard one due to the background expansion.
The idea of exploiting this extra friction to obtain an accelerated expansion in otherwise too steep potentials has been well explored in the context of primordial inflation. The original mechanism has been dubbed warm inflation [12] and for this reason we call our implementation "warm dark energy". Warm inflation assumes that the particles produced by the inflaton are in a thermal state. This assumption is dropped in the more recent implementations of this idea, as for instance trapped inflation [13], or the mechanism of gauge field amplification by an axion inflaton [14].
In this work, we explore this latter mechanism, developed by Anber and Sorbo in the context of primordial inflation [14], to the case of late time quintessence. The gauge fields are coupled to an axion-like particle, which couples to the gauge fields by means of the topological term 1 : The extra friction resulting from this interaction has also the effect of reducing the range of field values spanned by the inflaton, and can result in a sub-Planckian inflaton evolution 2 . For instance, for a typical axion potential V ∝ cos (φ/f ) (as in the so called model of natural inflation [22]), this coupling can result in an inflationary expansion also for an axion scale f that is a few orders of magnitude smaller than M p . On the other hand, a trans-Planckian axion scale of natural inflation, beside being in tension with the latest CMB results [23], can hardly be reconciled with quantum gravity [24][25][26][27][28][29][30][31][32][33][34].
Given the renewed interest in quintessence models because of the swampland programme, we revisit here axion models coupled to vector fields as a possible solution of the tension between swampland criteria and observations of the current phase of expansion of the Universe. We therefore investigate potentials satisfying (1.1), with c of order one, trying to make them compatible with observation thanks to the vector field couplings (1.2).
As in the inflationary application [14], the last term in (1.2) results in a significant gauge field amplification. This amplification occurs at the expense of the kinetic energy of the scalar field, and therefore it backreacts on the scalar field motion. Anber and Sorbo provided an analytic approximation for the backreaction term, that is valid only after a sufficiently prolonged period of accelerated expansion, with nearly constantφ. This is not the case at the onset of the present accelerated expansion. For this reason, we resort to a numerical simulation of this effect. Previous numerical studies either performed full lattice simulations [35][36][37][38], in which both the gauge and the scalar field are inhomogeneous, or adopted a discretization scheme [39][40][41][42], in which the scalar field is taken to be homogeneous, while a discretized set of modes is evolved for the gauge field. All these works focus on inflation or (p)reheating after inflation, with the exception of ref. [38], which studied the production of dark matter gauge bosons through this mechanism. We also follow a discretization scheme as this second class of works.
In this work we are also interested in the string theory uplift of warm dark energy models. The first step to embed them in string theory is to find an appropriate supergravity model that reproduces them at sufficiently low energies. We therefore identified a class of N=2 theories that produce the action (1.2). We focus on N=2 supergravity theories because they are a natural framework to describe the effective theories of Calabi-Yau compactifications. In fact their truncations capture the closed string sector also in the presence of supersymmetric configurations of branes and orientifolds. We will see that in a fairly general setup one can generalize the N=1 construction of [43] to N=2 in a way that also clarifies the higherdimensional origin of some of the ingredients. Although we do not offer here a full string theoretic construction, we see that extended supersymmetry gives us a better control on the origin of crucial couplings in (1.2), like the inflaton-vector coupling.
The paper is structured as follows. In Section 2 we construct N=2 supergravities that contain the model (1.2). In Section 3 we study the background evolution in this model, and show that the gauge field production can sufficiently slow down the dark energy field, so to lead to a late time acceleration compatible with data. While this part is mostly focused on the transition period between matter domination and dark energy domination, in the following Section 4 we study the later dynamics inside the dark energy dominated phase. We also study the onset of primordial inflation, choosing parameters as in the Anber-Sorbo paper, and considering whether their solution can be dynamically achieved from an earlier period of matter or of radiation domination. In Section 5 we then summarize and discuss our results. The paper is concluded by two appendices, in which, respectively, we provide an analytic approximation for the amplified gauge field and we present the rescaled evolution equations used in our numerical integrations.

Supergravity realizations
A generic N=1 supergravity model is defined in terms of a Kähler potential K, defining the scalar σ-model, a superpotential W , giving the scalar-self interactions and generating masses for the fields in the chiral multiplets, the gauge-kinetic functions f ΛΣ , describing the nonminimal couplings in the vector kinetic terms, and the gauge group G. We are interested in embedding our supergravity models of inflation/quintessence where the inflaton/quintessence field couples to a vector field. We therefore need at least one chiral multiplet containing our scalar field and a vector multiplet, containing our gauge field. The details of this minimal construction have been worked out in [43] for a rather general choice of scalar potentials and we therefore use it as a starting point for our discussion.
In this section we want to go beyond the minimal construction, asking ourselves what are the N=2 supergravity models that can give rise to the required couplings between the inflaton/quintessence field and vector fields. There are several reasons to look for such extension.
The first one is to see if and how extended supersymmetry constrains such scenarios. More interestingly, we would like to better understand the string theory origin of models of this type. While 4-dimensional supergravity models with minimal supersymmetry like those in [43] are fairly easy to construct, their stringy origin is usually very difficult to assess. In particular, quite often, fields that have an axionic coupling to the vector fields do not have canonical kinetic terms and therefore most of the studies on inflation/quintessence scenarios with vector field couplings should be re-evaluated. On the other hand, generic string compactifications scenarios use Calabi-Yau manifolds as internal spaces, supplemented by various types of branes, fluxes and orientifolds. This means that the closed string sector is described by an N=2 supergravity theory in 4 dimensions, which is then truncated to N=1 by the presence of orientifolds. Analyzing N=2 supergravity models we can therefore address if the inflaton/quintessence field can come from the closed-string sector and in case of positive answer, identify the relevant compactification scenario.
To recover the N=1 models of [43] as truncations of some N=2 supergravity, we have to imagine that our N=1 supermultiplets derive from either N=2 vector or hypermultiplets. As a first step, in the following we focus on scenarios that involve only N=2 vector multiplets. We therefore analyze N=2 supergravities that contain the gravity multiplet (g µν , ψ A µ , A µ ), coupled to n V vector multiplets (A i µ , λ i , z i ), i = 1, . . . , n V . Note that in N=2 supergravity the gravity multiplet, in addition to the graviton g µν and to the two gravitini ψ A µ , A = 1, 2, contains also a vector field, the graviphoton, which generically mixes with the other vector fields. This imposes several restrictions to the scalar self-couplings and to the couplings of the scalars to the vector fields. The σ-model described by the scalars in the vector multiplets is a Kähler manifold, but of a restricted type. In detail, its Kähler potential depends on the scalar fields only via projective coordinates X Λ (z). The reason of this restriction is the fact that complex scalar fields z i appear in vector multiplets, but there are n V + 1 vectors because of the graviphoton and they all mix under duality transformations. This means that the correct way to describe vectors is by means of a duality-covariant set of fields A Λ µ , where Λ = 0, 1, . . . , n V , so as to include the graviphoton. At the same time, since duality transformations mix up all the vector fields, the scalar fields associated to them should transform accordingly, hence the need of new projective coordinates X Λ , holomorphic functions of n V complex coordinates z i . Actually, duality transformations can mix vector fields with their duals A µΛ and therefore one has to introduce also "dual" holomorphic coordinates F Λ (z i ), which, in specific frames, can be deduced from a prepotential function F (X), by taking its derivatives with respect to the projective coordinates F Λ = ∂F (X) ∂X Λ . The prepotential function and its derivatives (denoted by additional pedices involving capital greek indices) determine all the couplings of the vector multiplet fields. The outcome is that the Kähler potential must have the following restricted form 3 [44,45] Once the number of vector multiplets has been fixed, any arbitrary holomorphic function F (X), which scales quadratically F → λ 2 F when X → λX, is a legitimate prepotential from which one can derive the various couplings in the action. However, for simplicity and because it comes as the natural outcome of Calabi-Yau compactifications, we will restrict the prepotential function to the case of very special Kähler manifolds, i.e. those with prepotentials of the form where d ijk is a constant totally symmetric tensor that describes the various intersection numbers of the cycles on the Calabi-Yau manifold. The Kähler potential in terms of the physical fields can be recovered by gauge-fixing the projective coordinates (for instance X 0 = 1, X i = z i ): As discussed previously, the prepotential fixes all the couplings of the vector multiplets and hence it also fixes the non-minimal gauge couplings to the scalar fields withF µν = 1 2 µνρσ F ρσ . These couplings are dictated by the N=2 gauge-kinetic matrix N = R + iI, with I negative-definite, whose explicit form can be deduced from the second derivatives of the prepotential F ΛΣ = ∂ 2 F (X) ∂X Λ ∂X Σ . The explicit expression is the following [45]: where repeated indices are summed over. In order to identify the relevant N=2 models, we need to consider which fields survive the truncation process to N=1. For instance, the truncation to N=1 will get rid of the graviphoton A 0 µ (because the N=1 gravity multiplet does not contain any vector field) and therefore one needs at least two additional vector multiplets in the model, one from which we recover a chiral multiplet containing the inflaton/quintessence field and one from which we recover the surviving gauge field. Luckily, the rules for N=2 supergravity truncations to N=1 have been worked out in detail in [46]. We therefore refer to that paper for the derivation of such rules, while we simply apply them in the following.
As mentioned above, we focus our analysis to the N=2 models with a cubic prepotential of the form (2.2). These depend on a totally symmetric tensor d ijk , whose possible form has been studied and classified in a series of papers by Van Proeyen, de Wit and collaborators [47][48][49]. In particular one can check the form of the gauge-kinetic function (2.5) and of the Kähler potential (2.3) and see which models have fields with canonical kinetic terms coupling linearly to the vector fields through the matrix R. After a straightfoward analysis one can check that the simplest model one can construct that leads to the wanted truncation requires n V = 3, which implies that one has an extra chiral multiplet in the resulting N=1 truncation, whose scalar is going to play the role of a stabilizer field. In detail, one can use a homogeneous manifold with and keep in the truncation the scalar fields z 1 and z 3 , while truncating µ . We verified that the rules for a consistent truncation given in [46] are satisfied. In fact, for X 2 (z) = z 2 = 0, the Kähler covariant derivatives of the holomorphic sections with respect to the truncated scalar identically vanish The scalar σ-model metric factorizes and the terms involving z 2 also identically vanish Finally, the gauge kinetic functions that mix the vector fields we truncated with the surviving one also vanish These conditions fully decouple the residual vector multiplet from the chiral multiplets in the N=1 supersymmetry rules. The result of the truncation gives a N=1 model with a gauge-kinetic coupling that is linear in the surviving chiral multiplet scalar field as in [43] (2.10) As mentioned above, we actually want more, because we need it to be linear in the quintessence scalar field, which should be canonically normalized. In order to see this, we split the z 1 field into its real and imaginary part, while we also introduce for later convenience a reparameterization of the other scalar field in terms of the stabilizer S: We will shortly see that the field φ has a canonical kinetic term if α = 0 and it couples linearly to the topological term of the surviving vector field once (2.10) is taken into account. The peculiar parameterization of the N=2 scalars z i in terms of the N=1 fields as in (2.11) has been chosen in a way that allows to obtain a model such that S acts as a stabilizer field to set both α = S = 0 consistently in the equations of motion. This can be achieved as follows. The Kähler potential deriving from truncating (2.1) is (2.12) The last term can be removed by a Kähler transformation In the new frame we can then fix the Kähler potential to The scalar potential follows from the superpotential in the usual way (2.14) In order to achieve a consistent truncation to α = S = 0, we can follow a procedure similar to the one outlined in [50] and propose a superpotential that has a linear dependence on S, so that where F is an arbitrary holomorphic function of Φ, which we subject to conditions coming from the equations of motion and Λ is a mass-dimension 1 parameter. The scalar field sector of our theory is described by the σ-model metric with scalar potential V following from inserting (2.13) and (2.15) in (2.14). The special choice of superpotential guarantees that S = 0 is a critical point and that where prime means derivative of the argument. This can be further simplified positing where u is a holomorphic function of Φ and f is a mass-dimension 1 parameter. The critical point at α = 0 is assured if u ū − uū = 0 (2.19) and the truncated potential at S = α = 0 is which can be used to embed almost any potential in our construction. For instance, the potential follows from the choice u = sin Φ 2f . While having a consistent truncation, it is not obvious that setting S = α = 0 gives a good effective theory. Although we do not have a general behaviour of the masses as functions of F that guarantees us that our consistent truncation is also a good effective theory, we can work out the specifics of the models of interest to see that it becomes true with some tuning. For the cosine potential, for f = 10 −4 M p , at φ = π we have that φ is tachyonic, α is clearly stabilized, but S is tachyonic, too. This is usually corrected by higher-order corrections in the Kähler potential that generate an additional term of order |S| 4 /M 4 inside the second log, giving S a mass of order M (see for instance [51,52]). However, as one can see from Figure 2, the ratio of normalized masses is such that the scalar fields will immediately roll in the φ direction, while S develops a huge mass. In fact all the plots in Figure 2 are in terms of the dimensionless parameter and the scale of the plot shows that very rapidly S acquires a huge mass. When the field φ reaches π/2 it becomes massless. At that point the other fields are much heavier. The situation is a bit more complicated for the exponential potential but it can also be solved with some additional tuning. In this case it is indeed difficult to achieve a consistent truncation and at the same time preserving a hierarchy in the masses, because of the nature of the exponential potential. If one naively sets F = e − i w 2 Φ Mp , the fields α and φ acquire masses of the same order. However, if one posits with w developing a small imaginary part w = 1 − i g f , with g f /10, then S and α have the mass ratio (2.22) of order 10 10 , while φ is almost massless. In this case we do not have a consistent truncation at α = 0, but the fields are stabilized by the huge masses they acquired.
Since we discussed Kähler potentials of restricted form because of their origin from string theory compactifications on Calabi-Yau manifolds, we can ask ourselves if the models we just presented can actually be obtained in such string theory reductions. Effective theories for Calabi-Yau compactifications with branes and orientifolds have been discussed in detail in [53][54][55]. While it is not difficult to check that there are Calabi-Yau manifolds with the right geometric structure to give rise to the prepotential (2.2) with (2.6), one should be careful about the natural units of the coupling between the inflaton/quintessence scalar φ and the vector fields. Although in a dual setup, the analysis of [56] shows that such couplings are always of order 1/(2π) 2 for all gauge fields in the closed string sector and can be enhanced by a factor N when considering gauge fields on N branes, provided we identify them. One therefore would need a rather large number of branes to make the couplings work, but this in turn would require to take into account the backreaction due to their presence, thus questioning the supergravity approximation. We therefore conclude that the full uplift of these models is still an open question, though we provided here some of the necessary ingredients, at least in the low-energy effective theory.

Dynamical evolution and accelerated expansion from dissipation
Let us now discuss the cosmological consequences of the models studied in the previous Section. We consider only the bosonic part of the action, including only the dynamically relevant fields. In this view, the study presented in this and in the next section can be applied also beyond the supergravity constructions that we have studied so far, and it relies on the starting action (1.2). We want to study whether this model can account for the present accelerated expansion of the Universe, for the cases in which the potential of a scalar field φ in this model respects the condition in (1.1), and it is therefore too steep to lead to accelerated expansion, in absence of particle production. We show that a typical axionic coupling between the scalar field in this model φ and a U(1) gauge field can allow for sufficient dissipation, that slows down the scalar field to a sufficiently level so to behave as an effective (slowly varying) cosmological constant.
This section is divided in five parts. In Subsection 3.1 we present the Anber-Sorbo mechanism for gauge field amplification and backreaction. In Subsection 3.2 we present and discuss the concrete potential that we study in this work. The following three subsections are devoted to study the dynamics of this model, from the early matter dominated phase to the onset of the late time acceleration.

Gauge field amplification from a rolling axion, and backreaction
We consider the action where φ is an axionic field that enjoys a shift symmetry φ → φ + C, which is respected by its interaction with a U(1) gauge field, and which is broken by the potential term V . In this action, F µν is the gauge field strength, andF µν ≡ µναβ 2 √ −g F αβ is its dual (where the tensor µναβ is totally anti-symmetric, and it is normalized to 0123 = 1). We assume a flat isotropic and homogeneous (FLRW) geometry, with line element where a is the so called scale factor of the Universe, and τ is conformal time. We assume that φ is a cosmologically relevant field. Therefore, to respect the FLRW geometry, φ only depends on time, φ = φ (τ ) up to small fluctuations, that we disregard in this work. The motion of φ then modifies the dispersion relations for the gauge field: the gauge field mode function obeys the equation where A + and A − are, respectively, the right-handed and the left-handed gauge field polarizations, and k is the comoving momentum of the gauge field mode. Finally, H ≡ ∂τ a a 2 denotes the Hubble rate.
As seen from eq. (3.2), due to the motion φ (t) one gauge field polarization becomes tachyonic for momenta smaller than the threshold momentum where we have assumed that ∂ τ φ > 0, so that the tachyonic mode is 4 A + . For notational convenience, from now on we disregard the A − mode, and we relabel A + as A. For sufficiently large coupling (sufficiently small f ) this leads to a strong gauge field amplification, which in turn can backreact on the background dynamics. The produced gauge field modifies the evolution equation for φ (following from the extremization of (1.2) with respect to φ) through the typical axionic coupling [14] It also modifies the Friedman equation (namely, the 00 component of the Einstein equation) for the evolution of the scale factor through its energy density [14] 1 a ∂a ∂τ where we have also included the contribution of the energy density of matter, which scales as ρ ∝ a −3 , and the suffix 'in' refers to a moment in the matter dominated era. In these expressions, we have used the standard electromagnetic notation, with (3.6) simply for notational convenience, without necessarily requiring that the U(1) field is the Standard Model photon (or the hypercharge gauge boson). As we show in the next subsection below in a concrete example, the main backreaction effect is on the dynamics of φ, via eq. (3.4). The backreaction term is controlled by the parameter ξ, which is in turn related to the speed of φ (for a constant φ, the last term in (1.2) is a topological term, which does not contribute to the equations of motion). This generates a friction term analogous to the second term in eq. (3.4) (the so called "Hubble friction"). As we show below, this additional friction can slow down the motion of φ and allow for accelerated expansion, in a potential that would be too steep to lead to acceleration in absence of this effect.

A concrete example
Let us now study how the mechanism presented above works for a specific example. For definiteness, we consider an exponential potential We consider a sufficiently early time, when the energy density is dominated by that of matter, and the scalar φ is nearly frozen, due to Hubble friction, with a negligible energy. Eventually, the matter energy density is diluted by the expansion of the Universe, and the scalar field becomes dynamically relevant. For λ < √ 3, the scalar field eventually comes to dominate over matter, and acquires an equation of state (pressure over energy density) [57] This value is obtained only after the scalar field comes to dominate. Moreover, we are for the moment disregarding the coupling between φ and the gauge field. Accelerated expansion requires w < − 1 3 . This is however not enough to agree with the data. For this potential, comparison with data provides an upper bound on λ. Marginalizing over the scale of the potential, one obtains λ < ∼ 0.49, 0.80, and 1.02 at 68%, 95%, and 99.7%, respectively [11]. This is incompatible with the generic expectation coming from string theory, where λ √ 6 [7,9]. In our numerical study, we fix λ = 1, leading to w φ = −2/3 when the scalar field dominates. This is compatible with the data only at 3σ, and, we therefore study whether the additional friction from this production can allow for a more phenomenologically acceptable expansion law.
As can be verified a-posteriori, due to the friction from the Hubble expansion and the gauge field amplification, the evolution of the scalar field only spans a small region of the potential, and we could perform the present study for a linearization V = α + βφ of a generic potential. However, the current presentation is focused on the potential (3.7) for definiteness. We solve the background evolution equations for the scalar field, for the scale factor, and a set of N equations (3.2) for a grid of N different modes of the gauge field (each one characterized by a distinct comoving momentum k). We solve these equations in three ways, progressively including the various physical contributions. Firstly, we disregard the backreaction of the produced gauge fields, and we include only the energy density of matter in the evolution for the scale factor. This allows to obtain an analytic solution for the scalar field and scale factor evolution, that is accurate only at very early times, when the φ, A sector is completely subdominant 5 . This allows us to obtain the initial conditions for the numerical evolutions that we consider next. Secondly, we include the energy density of the scalar field in the evolution equations, but we still disregard the backreaction of the gauge fields. We now solve also the equations for the gauge fields in this approximation, and we evaluate the backreaction terms as a diagnostic (we compute it, but we do not include it in the evolution equation). This allows us to find the moment at which the approximation breaks down and, more importantly, the range of momenta of gauge field modes that dominate the backreaction term at this moment. Thirdly, we evolve the full set of equations, making sure that the sampled range of momenta for the gauge fields covers the dynamically relevant regions, and that the set of probed modes is sufficiently dense.
We discuss these three solutions in the next three subsections.

Analytic solutions for ρ m V and negligible backreaction
We consider an initial time τ in deep inside the matter dominated regime. With no loss of generality we can always set at this initial time. The first condition is allowed by the fact that the normalization of the scale factor is arbitrary (and unphysical) in a flat Universe. The second condition is allowed by the fact that any constant initial value for the scalar field can be reabsorbed in the value of V 0 in the exponential potential (3.7). At these early times, the motion of the scalar field is extremely small due to Hubble friction, and φ essentially contributes to the expansion of the Universe as a cosmological constant term V 0 . Therefore, the ratio between the energy density of matter and is the most immediate physical quantity that can be used to parametrize the initial time.
Keeping only the energy density of matter at the r.h.s. of eq. (3.5), one obtains the solution where an integration constant has been fixed to set the first of (3.9). To obtain the early time solution for the scalar field, we take equation (3.4), we insert in it the solution (3.11), we disregard the backreaction term, and we approximate the potential derivative as ∂V Mp as it is appropriate at the early times, when φ M p . The resulting equation is then solved by where one integration constant has been used to remove a decreasing mode (which, if included at times τ in , would have become negligible at τ in ), and another one to set the second of (3.9).

Solutions with negligible backreaction
We now improve over the solutions presented in the previous subsection by consistently solving eqs. (3.4) and (3.5) in absence of the backreaction terms (these are therefore exact solutions in the f → ∞ limit). The parameters V 0 , ρ m,in and M p can be then removed from the evolved system if we rescale φ = M p ϕ and if we use the dimensionless combinatioñ as time variable, the evolution equations become: (3.14) The initial time (3.11) corresponds to the rescaled initial timeτ in =ρ −1/6 m . To estimate when the transition between matter domination and scalar field domination occurs, we equate the last two terms in the second of (3.14), corresponding, respectively, to the potential energy of the scalar field and to the energy density of matter. As we are interested in a parametric estimate, we set ϕ = 0 in this comparison, and obtain where the early time analytic solution (3.11) has been used in the second estimate. This final estimate justifies the rescaling done in (3.13). Namely, we chose a dimensionless time variable that evaluates to 1 at the equality between matter and the scalar field. We solved the system of equations (3.14) numerically, fixingρ m = 10 12 and, as we wrote above, λ = 1. The conditions (3.9) provide the initial conditions for this evolution. In the left panel of Figure 2 we show the evolution of the scalar field. We see that φ M p in the dynamical range of interest, justifying the statement that we could have obtained analogous results by simply considering a linearization V = α+βφ rather than the exponential potential (the backreaction due to the gauge field production that we include in the next subsection further decreases the range spanned by φ). In the central panel of the figure we show the evolution of the fractional energy of the scalar field We recall that w < − 1 3 leads to accelerated expansion, and we note that the equation of state indeed starts from the early matter dominated value w = 0 and it then evolves to the value −2/3 at late times, in agreement with eq. (3.8). In the central and in the right panel, the evolutions as shown are a function of the number of e-folds of expansion, defined as

Growth of the backreaction term
We study here the backreaction effects of the gauge field production, encoded by the last two terms in eqs. (3.4) and (3.5). We evaluate them using the numerical solutions for φ and a derived in absence of this backreaction. Therefore, the solution for the backreaction terms that we obtain is valid only as long as it is negligible. This computation allows us to study the growth of the backreaction terms, and to find the value of the axion decay constant f for which the backreaction becomes important close to the transition point between the matter field and the scalar field domination. Moreover, it also allows us to understand what is the relevant range for the gauge field modes that we need to cover in the full numerical evolutions that we study in the next section (namely, what are the modes that dominate the backreaction terms when they become dynamically relevant). We start from eq. (3.2) for the gauge field modes, that we rewrite as We can insert in this expression the time derivative dφ dτ from the numerical system evolved in the previous subsection, or the early time analytic solution (3.12), valid during the matter dominated regime. In the latter case, rescaling the time as in eq. (3.13), we obtain We see that at sufficiently early times the first term dominates the parenthesis, and the mode is stable. The second term becomes dominant at sufficiently largeτ , leading to the amplification of the gauge mode. At any fixed time, the maximum momentum amplified is We normalize the comoving momentum to the value of k max at the estimated transition timeτ = 1,k In these variables, the equation for the gauge field modes becomes where the axion decay constant has also been rescaled according tof ≡ f /M p . In Appendix A we present an analytic estimate for the solution of this equation, and for the corresponding backreaction on the evolution of the scalar field. As a measure of the backreaction, we take the ratio between the fourth and third term in eq. (3.4), which is smaller than one in the regime of negligible backreaction, and of O (1) when backreaction is important. Thanks to the above rescalings, we know that we need to study the backreaction atτ = O (1) (the moment of approximate equality between the matter and the scalar field) and for momentak up to O (1) (the highest momentum excited at this time). In Appendix A we provide an analytic estimate for the backreaction term (see eqs. (A.4) and (A.6)). In figure 3, we plot the integrand B EB obtained from that estimate, against the one obtained from the numerical evolution of the system (3.14). As discussed above, the integrand is normalized so that the backreaction becomes relevant when the area underneath the curve approaches one. In all cases shown the integral under the solid line is much smaller than one. The vertical line shown in the plots separates the unstable (on the left) from the stable (on the right) modes, according to the analytic approximation (3.21). As explained in the appendix, our analytic estimate holds for the unstable regime. We note that, in this region, the approximated line matches the exact one well at early times, while the two lines start to diverge asτ approaches one. This is due to the fact that the analytic approximation assumes a matter dominated background, which stops to be true atτ 1. In the stable region to the right of the line, the modes are still close to the vacuum state, and we should also add the contribution of the A − mode, that would cancel against the one shown here in the high momentum vacuum regime. We do not implement this cancellation (namely, we do not include the backreaction of the modes A − in our computations), because this region always contributes a negligible amount of backreaction.
We see from the figure that, even if the analytic approximation does not reproduce the correct amplitude of the backreaction atτ 1, it still continues to correctly provide the momentum range of the unstable modes. Therefore, we can use these runs to understand the dynamically relevant range of momenta that must be included in the full numerical simulations that we discuss in the next subsection.

Solutions with backreaction included
The final step in our analysis is to solve the full system of equations including the backreaction term in the equation of motion of the scalar field. We want to investigate what values of the coupling constant f produce a sufficiently large backreaction that affects the motion of the scalar field before the present moment. The additional friction changes the equation of state parameter of the scalar field w φ ≡ p φ ρ φ with respect to what is obtained in the absence of backreaction (formally, in the f → ∞ limit).
The system of equations to be evolved as well as all the necessary re-scaling of the variables are presented in Appendix B. In order to compute the backreaction term in the equation of motion of the scalar field, we discretize the momenta in equally spaced intervals in log k space. The contributions of the various momenta are then added up using the trapezoid rule. We varied the number of discreet k modes between O (10) and O (1000) and determined that there is negligible difference in the evolution for anything larger than approximately 500 modes.
In Figure 4 a number of quantities of interest are displayed for a sample run of our numerical code. The showed run terminates about half an e-fold after the "present moment", which is defined as the moment in which the ratio of dark energy to matter is the observed one. The "present moment" is denoted in the figures with a vertical red line. We have reserved a more in-depth discussion of the evolution deep into the accelerated regime for the next section 4.
The main important feature that emerges from the numerical results is that an acceptable equation of state for w φ can indeed be obtained from this model due to the backreaction. Specifically, we notice from the bottom left panel that, before the backreaction becomes relevant (namely, for N < ∼ 8.5), w φ is climbing toward the late time attractor (3.8). Backreaction adds a strong friction to the motion of φ, that brings w φ back toward −1. We notice from the plots that backreaction manifests itself very suddenly (we verified that this is not an artifact  Again the moment the backreaction becomes dominant is visible by the sudden change in value that occurs at N ∼ 8.6. The red line again denotes the present moment. Bottom right panel : The particle production parameter defined in (3.2) of a poor sampling of the modes in the numerical evolution), rather than in a smooth way, due to its exponential growth. The scalar field reacts to this by performing some oscillations about an overall growing solution. These oscillations are too rapid to be seen clearly in figure  4 but it is possible to see a thickening of the blue line immediately after backreaction has become dominant in the equation of motion for the scalar field. The oscillations are more clearly visible in Figure 5, which is a zoomed version of the previous figure, centered around the moment in which backreaction becomes important.

Comparison with the data
The above results indicate that, even if one starts from a scalar field potential that does not provide a sufficiently fast acceleration (a large value of λ in the specific case (3.7) that we are studying), a sufficiently strong amplification of the vector field can result in an equation of state w φ sufficiently close to −1, as required by data.This suggests that, for any choice of λ, there can be a threshold value f thr (λ) below which the mechanism can result in an acceptable phenomenology. Finding the precise threshold value however requires a dedicated data comparison between this model and the data, which is beyond the scope of this paper. In particular, an input for the comparison with the supernovae data is the luminosity distance d L = a 2 0 a a 0 a da a 2 H(a ) . The distance can be derived from the expansion law a (N ), which in turn is one to one related to the evolution of the scalar field equation of state w φ (N ). In the present case, several choices of f results in significant oscillations of w φ , so that the simplest parametrizations existing in the literature cannot be used. This is for instance the case for the example shown in Figure 5, in whichf = 8.7 · 10 −4 was chosen.
Even without a dedicated data comparison, we can obtain an interval in which we expect to find the threshold value. Specifically, we can obtain a value f 1 for which we can be highly confident that the mechanisms agrees with data, and a greater value f 2 for which the mechanism is not in agreement with the data (we recall that increasing the value of f results in a smaller backreaction). This will allow us to indicate that the threshold value is between f 1 and f 2 . The value of f 2 can be obtained with relative ease. Specifically, by numerically evolving the system as described in the previous subsection, we find that forf = 9.25 · 10 −4 , or higher, the backreaction becomes noticeable only at the present time (we recall that by 'present time' we mean the moment in the evolution in which the ratio between the energies of the dark energy and of matter is the observed one). This is phenomenologically indistinguishable from the case of no gauge production, which, as we discussed above, is compatible with the data only at 3 σ (we recall that we are assuming λ = 1). Thereforef 2 = 9.25 · 10 −4 .
To find f 1 , we consider a sufficiently strong production so that w φ reaches a smooth evolution by the present time, with very suppressed oscillations. In this case, the late time evolution of w φ is well described by the parametrization shows the allowed region in the parameter space (w 0 , w a ) for different choice of data. We choose the dataset that results in the most stringent constraints. The the central panel is forf =f 1 = 9.25 · 10 −4 , leading to a significant evolution of w φ , so that the parametrization (3.25) cannot be employed for this choice off (we know that this choice is ruled out by data, since it leads to an appreciable backreaction only after the present time). The right panel is forf =f 1 = 2.2 · 10 −5 . This leads to a much smaller change in w φ in the last e-fold (see also Figure 6), and to an evolution inside the 1σ allowed region. Fig. 6 shows the evolution of the scalar field and of its equation of state for a coupling for which the oscillations in w φ took place well before the present time, and for which we argue that the late time evolution can be sufficiently well described by (3.25) for the purpose of data comparison. We support this claim through Figure 7. In the left panel of the figure we reproduce the 1 and 2σ contours from Figure 30, of [58], focusing on the w 0 ≥ −1 region. In the central panel we superimpose to these contours a line that encodes the evolution of w φ and of its derivative wrt the scale factor, for the choicef =f 2 = 9.25 · 10 −4 . We recall that for this choice, backreaction is important only after the present time, and the model is incompatible with data. The evolution encoded in the line lasts for one e-fold before the present time, with earlier times on the left of the line, and later times on the right (in the color version of the figure, the color of the line changes from blue to red as time progresses). Due to the strong evolution of w φ and of its derivative the parametrization (3.25) cannot be employed, and the main purpose of this panel is to provide an immediate comparison with the right panel. In the right panel we show an analogous evolution of w φ for the choicẽ f = 2.2 · 10 −5 . We note the very different scale shown now in the figure. As can also be seen from Figure 6, the equation of state is extremely close to −1 during the last e-folds, and its derivative has fast oscillations, that remain within the allowed 1σ region. Although, strictly speaking, the parametrization (3.25) requires an exactly linear evolution of w φ (a) we see from the horizontal scale that the fast and small oscillations of dw da around zero nearly average out on the net evolution of w φ (a). Due to this, and due to the fact that evolution is always within the 1σ allowed region, we conclude that this choice is compatible with the data. Therefore we setf 1 = 2.2 · 10 −5 .
To summarize, a rigorous data comparison for this mechanism requires a dedicated analysis, due to the nontrivial and oscillatory nature of w φ (a). We know that the mechanism will be compatible with the data for f smaller than a threshold value (which depends on the curvature of the potential, which in our case is controlled by the parameter λ). From our numerical study, performed with λ = 1, that have argued that the threshold value should be between 2.2·10 −5 M p and 9.25·10 −4 M p . A more refined computation of the precise threshold value f thr can be obtained with a dedicated analysis of the data for the present model.

Application to primordial inflation and the road to the Anber-Sorbo solution
In addition to the case studied in the previous section, we are interested in applying our numerical code to the transition between a pre-inflationary era and inflation. In the original Anber and Sorbo paper [14], the model is studied under the assumption that the backreaction term is always the main source of friction of the motion, and that the inflaton has a nearly constant speed, determined by the interplay between this term and the derivative of the inflaton potential. This results in a slow roll inflationary solution, with nearly constant inflaton speed. The backreaction term is computed in [14] by disregarding the departure from de Sitter expansion, and by providing a convenient approximation for the mode functions of the amplified gauge fields. One obtains [14] where the last term is the one due to backreaction, and the numerical coefficient is I ≡ 7!/(2 21 π) 2.4 × 10 −4 . Equating the backreaction term with the derivative of the inflaton potential leads to [14] ξ 1 2π which is the Anber-Sorbo (AS) solution in the regime of slow roll, sustained by the gauge field amplification. We want to understand whether the Anber-Sorbo inflationary solution can be dynamically attained if one starts from a pre-inflationary period of matter domination, during which the gauge field amplification, and its back reaction is negligible. Loosely speaking, this analysis extends to later times the one performed in the previous sections, since so far we have been interested to at most one e-fold of accelerated expansion (corresponding to the present dark energy domination) while now we want to discuss the evolution well inside the inflationary regime.

Numerical considerations
The system of equations to be evolved as well as all the necessary re-scaling of the variables are presented in Appendix B. Evolving our numerical code deep into the inflationary era poses a numerical challenge. The reason for this is that the more e-folds of inflation we cover, the greater is the momentum range that corresponds to the unstable modes that have to be included in the momentum integral in the equation of motion of the inflaton. The momentum threshold between stable and unstable modes is given by If, as a first approximation, we ignore the variation ofφ during inflation, we can see that during the course of inflation the threshold is increasing. Assuming that we evolve N 15 e-folds of inflation, the value k th will be a = e 15 3 × 10 6 times greater at the end compared to the initial value of the threshold. Evolving gauge modes of really high momenta slows the code down significantly since we are evolving gauge modes which spend most of the evolution in the stable regime, in which the evolution of the mode is highly oscillatory.
An additional problem associated to the high momenta is that the integral giving the backreaction is UV-divergent (as we did not renormalize the modes, since this was not needed for the computations performed in the previous section 6 ). Therefore, if we include a too large range of high momentum mode, instead of integrating over the instability associated with particle production, one integrates over a contribution associated with modes which are still in their vacuum configuration. This problem is also apparent from the first two panels of figure 3. We see that, at early times the integrand has a UV contribution that is dominant compared to the instability associated with particle production. For inflationary energy scales, and for very large momenta (which are necessary for consistent evolution deep in the inflationary regime) it is possible that the UV contribution might become dominant at early times. In principle, one way to overcome this problem would be to put a cut-off at the interface between the stable and unstable modes. We will see that this is not possible in our case because, due to a non monotonic evolution of dφ dt , there are gauge modes that become stable after they have been enhanced; therefore, their contribution would be incorrectly neglected if we included in the backreaction only modes that are unstable at any given time.
In order to overcome the problems mentioned above we will perform the evolution by choosing a convenient value of the potential 7 V 0 = 10 −20 M 4 p . As can be seen by eq (3.24) in order for the backreaction to become dominant the integrand has to overcome a suppression of the order of V 0 M 4 p . The smallness of the potential allows us to choose a higher value of the highest momentum, without risking the UV spurious divergence being dynamically relevant. Simultaneously it is beneficial to have a not too small value of the potential because one does not have to evolve too many e-folds for the gauge modes to overcome the suppression and thus all the various effects arise sooner. Our choice of V 0 is a compromise between these two opposite requirements.
Another way to improve the performance of the code is to keep the really high momenta "frozen" in their vacuum configuration until late in the evolution when they come close to the instability region. We realize this by separating the evolution into two parts. In the first part we evolve the numerical system from deep inside the matter dominated period until the early inflationary era, including only gauge modes that are relevant during this time period. We use the final values of the fields obtained in this run as the initial condition for a second run, that extends more deeply into the inflationary regime. In this second run we consider a much greater range of high momentum modes, that we expect to become unstable during this second evolution (these additional modes are initialized from their vacuum state; we verified that this is a good assumption at the start of the second evolution).

Numerical results for an exponential potential
We discuss here the results obtained for the run discussed in the previous subsection, for the same potential (3.7) that we have considered so far in this work (choosing λ = 1, and V 0 = 10 −20 M 4 p , as described above). For completeness we performed the analysis assuming an equation of state of matter (w = 0) as well as radiation (w = 1 3 ) for the initially dominating fluid. The two cases provide qualitatively similar results, with the main (trivial) difference that, starting from the same ratio between the energy density of the initial fluid and that of the scalar field, the transition from fluid domination to inflation takes place sooner in the case of initial radiation domination. We show in figures 8 and 9 some of the key quantities obtained in the evolution for the first case (initially, w tot = 0).
An interesting feature that characterizes the later inflationary stage, is the appearance of different stages in which the scalar field evolves more slowly, separated by quicker stages of faster evolution. This is visible in the top-left panel of figure 8, which shows the evolution of the scalar field, as well as in the bottom-right panel, which shows the evolution of the gauge production parameter ξ ≡ ∂τ φ 2af H . The evolution shows signs of self-similarity, and the different spikes in ξ (N ) (each spike corresponding to a stage of enhancedφ) have a nearly identical amplitude. The dashed line in the bottom right panels of both figures denotes the value for the particle production parameter obtained from the Anber-Sorbo analytical computation reproduced in our eq. (4.2). Qualitatively speaking, we see that the Anber-Sorbo predicted value is achieved in an "averaged" way, but that the evolution does not reach a stage in which eq. (4.2) is valid at all times. We have also experimented with various values of parameters in the code, and, particularly, with the amplitude of the potential varying between V 0 = 10 −20 M 4 p and V 0 = 10 −120 M 4 p . This patterns repeated in all cases that we studied, with a frequency that decreases with decreasing value of V 0 .
The Anber-Sorbo solution denotes the case in which the effect of gauge particle production and the backreaction to the equation of motion of the scalar field are in perfect equilibrium with each other. In our case, the transition to inflation is achieved far away from that equilibrium, as testified by the presence of the spikes in the gauge production parameter ξ. Another way to describe the observed effect is the following. If, at some point of the evolution of the scalar field the gauge field production parameter is below the Anber-Sorbo solution, that implies that there is less friction in the equation of motion of the scalar field and therefore the scalar field accelerates. Since the speed of the field increases, the gauge field production also increases. The production keeps growing until it overshoots the equilibrium level defined by the AS solution. This leads to a quick increase of the backreaction, that then significantly slows down the scalar field and the gauge production. This leads to a quick Bottom left panel : Evolution of the state parameter of the scalar field. Bottom right panel : Evolution of the particle production parameter defined in (3.2). The dashed line that is superimposed denotes the value of ξ as predicted by the AS solution (4.2). One can observe a series of spikes after N 16 that appear to be self similar. decrease of the gauge field amplification. The gauge modes produced up to that moment are then subject to redshift, which decreases the backreaction, and it allows the scalar field to speed up again. This process repeats itself and it does not show any sign of relaxation to an equilibrium (constant or nearly constantφ state) during the evolution that we could cover in our numerical integration.
The six panels of figure 9 demonstrate how the integrand defined in (3.24) varies with time. The first four panels are chosen to be around the transition time N 7.8 whereas the last two show later stages of the evolution. The vertical dashed line is plotted at the exact threshold between the stable and unstable modes at the time shown. In the top left panel the backreaction is still subdominant, and one observes a similar shape for the particle production as displayed in figure 3. In the top middle panel the backreaction has just become dominant and the speed of the inflaton has just decreased sharply. As a consequence the threshold between stability and instability decreased as well to a value on the left of the range shown in the panel 8 . The threshold between stability and instability then proceeds to  fluctuate wildly for the next ∆N ∼ 1.5 e-folds 9 with a positive average value that is slightly greater than 0.002. For this reason, the contribution to the backreaction from the modes atk 0.002 continues to grow with respect to what shown in the first three panels. This growth is visible in the fourth and fifth panel. On the contrary, the modes corresponding to the second positive peak in the third panel are stable during the evolution between the fourth and the fifth panel, and their contribution to the backreaction decreases (due to cosmological redshift). In fact, we observe that all the contributions with momentum greater thank 0.002 have been redshifted away as a −4 and they are too small to see in a linear scale 10 .
Finally, we can make a number of observations on the last panel shown in the figure: firstly, the bump aroundk 0.002 has decreased due to cosmological redshift. Secondly, we notice an enhancement present at aroundk 0.010. These are the highest modes that were enhanced immediately before the backreaction became important. These modes then became stable (and therefore their amplitude oscillated) while simultaneously redshifting away. Due to the redshift, the contribution of these modes does not appear to differ from zero in the fifth panel. Nonetheless, their amplitude is greater than that of the surrounding modes, field moves slowly due to Hubble friction, and backreaction effects are subdominant. Due to this one might be tempted to includee in the numerical evolution only the gauge modes that are unstable at any one given time. The sharp decrease of the threshold that takes place between the first two panels shows that this would be a mistake. The third panel shows modes that were previously amplified, and that are now from some time in the stability region. Their contribution to the backreaction is highly significant at this stage. 9 This is a consequence of the rapid oscillations ofφ. The same qualitative behavior can be observed in the last panel of Figure 4). 10 Modes withk 0.002 are formally unstable; namely, they have a tachyonic effective frequency in eq. (3.19). However, due to the smallness of k, this instability is extremely mild, and it is dominated by the dilution due to the redshift. and it acts as a seed for the amplification that leads to the highest bump seen in the last panel. Thirdly, we observe a new emerging bump at aroundk 0.2. This enhancement corresponds to newly unstable modes that are enhanced from their vacuum configuration only in the latest stages of the evolution shown. This is the bump that grows the fastest at the times immediately after those shown in the figure, and it will eventually become the single dominant bump, much like as in the first panel of the figure. This generates a new cycle in the evolution, in which the features shown in these panels repeat themselves for the higher,k = O 10 1 , comoving momenta.

Numerical results for a sinusoidal potential
In this subsection we want to examine whether the (approximate) self-similar inflationary evolution shown in Figure 8 is a consequence of the assumed exponential potential. We therefore replace the potential with a sinusoidal one, which is more typical for axion fields, and that was also chosen by Anber ad Sorbo [14]. Since f now denotes the scale of the potential, we follow [14] and we parametrize the coupling in the topological term with an additional positive parameter α, To study the proposed Anber-Sorbo solution, we assume the presence of N Abelian gauge fields, all coupled to φ with an identical coupling (4.5). Quite nontrivially, the presence of N gauge fields provides an increased friction term also in the evolution equation for the perturbation of φ, leading to a nearly scale invariant (as long as the Anber-Sorbo solution is achieved) spectrum of the primordial density perturbations with amplitude P ζ 0.05 N ξ 2 [14]. If this mechanism is realized in nature, the combination N ξ 2 therefore needs to be fixed so to produce the observed value P ζ 2.1 · 10 −9 [23]. In the presence of N gauge fields, equation (4.2) should also present the rescaling α → N α, so that we need to enforce We fix ξ 20 as the benchmark value chosen in the discussion by Anber and Sorbo. This then requires N = 6 · 10 4 identically coupled gauge fields. Eq. (4.6) can be then turned into a relation for the scale of the potential: where ϕ needs to be evaluated when the CMB modes left the horizon. We choose, ϕ f , and, after parametrically identifying the value of the potential with the energy scale of inflation, we find that this corresponds to "low-scale" inflationary scale V 0 = O (10 TeV). We modified our numerical code to the choices of parameters just described, and performed a numerical evolution, starting for definiteness from φ in = π 4 f . In Figure 10 we display the results for the case of transition from radiation domination to inflation. The code was separated in two runs like in the case of the exponential potential. The first run was performed from N = 0 until N = 13 e-folds and in a range of momenta between 2 · 10 −6 ≤k ≤ 1 · 10 −3 . The second run lasted from N = 13 until N = 25 and the range of momenta was 2 · 10 −6 ≤k ≤ 50. Bottom right panel : The particle production parameter displays similar features to the exponential potential. The transient effects are different but the "spikes" of particle production are still present. Again we can see that the AS-Solution (dashed line) is qualitatively achieved in an average way.
Unfortunately, running the code deeper into inflation would require including higher value of the momenta, which slows the numerical evolution to a prohibitive level. However, the amount of inflation that we could cover shows that the bursts of particle production are also present for the case of the sinusoidal potential, and that the AS inflationary solution is reached only in a (qualitatively) average fashion. The spikes in this case appear with a slower frequency than in the case of the exponential potential. This is a consequence of the smallness of the potential rather than the shape of the potential. We also experimented with various values of V 0 between ∼ 10 −20 M 4 p and ∼ 10 −120 M 4 p . These increased values are inconsistent with the CMB normalization, but they allow to obtain more bursts of gauge field amplification in the numerical integration. We could verify that also in this case sufficiently large values of V 0 lead to a qualitatively self-similar spikes that are analogous to the ones that appear in the exponential potential example (although clearly, in this case, the self similarity will disappear when φ approaches the minimum of the potential).

Conclusions
Cosmological data suggest that the Universe underwent two periods of accelerated expansion: one responsible for setting the initial conditions (background and perturbations) of the "standard" big-bang era, and the current one.
The observed cosmological perturbations require that the primordial expansion is not an exact de Sitter one. This is successfully explained by models of slow-roll scalar field inflation. On the contrary, the cosmological data are consistent with the assumption that the present Universe is dominated by a vacuum energy density, and strongly constrain any departure from this simplest possibility [23]. Despite of this, our inability to naturally obtain such a small vacuum energy as the one required by the Universe has motivated the construction of quintessence models where the vacuum energy is set to zero (typically, by some unspecified symmetry), and the current expansion is due to some additional source, denoted by "dark energy". Borrowing from the success of scalar field inflation, the simplest models of dark energy make use of scalar fields rolling in sufficiently flat potentials. Of particular interest are models of tracking quintessence, in which the evolution of the scalar field tracks that of the dominant background source (matter and radiation) for most of the cosmic history, and emerges only at late times [60]. The tracking mechanism can allow for an explanation of the coincidence problem, namely why the current value of the dark energy is comparable to that of dark matter.
Also in these interesting models, however, the mass of the quintessence field (namely, the curvature of the potential) needs to be set to an extremely small value (as compared to typical scales that arise in particle physics) set by the current Hubble parameter H 0 = O 10 −42 M p . In addition, as already discussed in the introduction, several theoretical considerations appeared in the recent literature pointing to the requirement that the slope of the scalar potential might be too steep to support accelerated expansion. This observation provides one of the main motivations for the present analysis. Also in this case we borrow an idea that, despite not being as popular as standard slow-roll inflation, is well known in the inflationary community, namely that of warm inflation [12]. In this framework, the scalar field potential is too steep to lead to accelerated expansion, if one accounts only for the Hubble friction. This problem is overcome by the fact that the scalar field is coupled to some additional field. The coupling, together with the motion of the scalar field, leads to the production of this additional field, at the expenses of the kinetic energy of the rolling scalar. This can significantly slow down the rolling scalar, allowing for an effective equation of state of the scalar that can provide accelerated expansion. The more recent literature has implemented this idea in several relatively simple QFT constructions, as for instance the Anber-Sorbo mechanism [14], trapped infaltion [13] (see also [61] for a reanalysis of cosmological perturbations in this model) and chromo-natural inflation [17]. Our main goal is to implement the mechanism of [14] in the context of the present expansion.
In this mechanism, the rolling field is a pseudo-scalar φ, coupled to a U(1) field via the topological term φFF . In the first part of this paper we present a novel implementation of this mechanism in the context of extended supergravity models. Most phenomenologically viable models coming from String Theory assume a compactification to 4 dimensions by means of a Calabi-Yau (CY) manifold as internal space, supplemented by the presence of branes and orientifolds. The final low-energy effective theory is given by N=1 supergravity, but of a restricted type. In fact the closed string sector modes experience the underlying CY geometry and their interactions are therefore better described by an N=2 model, possibly further restricted by the geometric details of the chosen CY. We show that warm dark energy scenarios can be constructed in this N=2 setup by choosing a CY with specific intersection numbers that allow for the appropriate couplings between the vector field and the quintessence scalar (see eqs. (2.2) and (2.6)). By introducing an orientifold projection we can also consistently embed in this scenario almost any scalar potential via an appropriate choice of superpotential, though some tuning may be required. While this does not provide a fully consistent uplift of warm dark energy models within string theory, we have clearly identified most of its necessary ingredients.
We account for the backreaction of the produced gauge field by solving the system numerically. As we disregard the perturbations of the scalar field in this paper, we can write linear exact equations (in the limit of homogeneous φ) for the U(1) vector field mode functions. We solve these equations for a dense grid of modes in momentum space. Focusing for definiteness on an exponential potential, we verified that a sufficiently strong φFF interaction can indeed successfully account for the current acceleration of the Universe.
The current accelerated stage has been taking place for only about one e-fold of expansion. Therefore, to perform this study it is enough to reach only the beginning of the accelerated stage. In our numerical simulations the expansion is initially dominated by the energy density of matter. This strongly suppresses the motion of the scalar field (just due to Hubble friction) and the amplification of the gauge fields. 11 The scalar field comes eventually to dominate the energy density of the Universe, as the energy of matter decreases due to the expansion. Compatibility with data in a steep potential requires that the backreaction from the produced gauge fields becomes important before the transition from matter to scalar field domination, allowing for at least an e-fold of sufficiently accelerated expansion.
While not necessarily for this quintessence application, we nonetheless exploited our numerical code to also study the evolution well inside the accelerated regime. As it is well known, simulations of inhomogeneous fields during inflation on a fixed grid are extremely expensive since, due to the rapid growth of the scale factor, they must cover a progressively larger and larger dynamical range as the inflationary expansion goes on. Therefore we could only simulate a limited number of e-folds of accelerated expansion. The simplest expectation from this study is that, in the strong backreaction regime, the background evolution should approach that of the Anber-Sorbo solution [14] (in a very different energy regime, in the quintessence application, or in the original regime, if we change the parameters of the code to simulate the primordial inflationary stage). The Anber-Sorbo solution is characterized by an adiabatically evolving inflaton speed, and by a nearly identical evolution of the different gauge modes. Specifically the evolution of a gauge mode of comoving momentum k 1 at the time t 1 is nearly identical to that of a mode of comoving momentum k 2 at the time t 2 when we plot it in terms of the physical momentum of the two modes, namely when k 1 a(t 1 ) = k 2 a(t 2 ) ). Therefore, in this solution, the system evolves in a nearly steady state fashion, which is the typical slow-roll inflationary expectation, characterized by only a very small breaking of time invariance (this is what generates nearly scale invariant signals).
On the contrary, the inflationary stages that emerged from our simulations are characterized by a highly nontrivial evolution of φ (t). As we approach the transition between matter and scalar field domination, the scalar field starts accelerating. The backreaction from the produced gauge fields, while still subdominant, grows extremely rapidly with time, and it becomes dominant all of a sudden, giving rise to a large burst of gauge field amplification. This suddenly slows down the gauge field, and actually makes it oscillate back and forth for some time. In this phase, the energy density of the produced gauge field redshifts away due to the expansion of the Universe, and no new fields are efficiently produced, due to the decreased speed of φ (the gauge field amplification is exponentially sensitive to the speed of φ, so even a small decrease of the speed can result in a very strong decrease of the gauge field production). Due to this, the backreaction becomes rapidly subdominant again, and the scalar field starts accelerating once more. The cycle we just described repeats itself over and over in our simulations, with no sign of convergence toward the Anber-Sorbo steady state solution. A useful diagnostic for this study is the evolution of the parameter ξ (t) ≡φ (t) 2f H(t) , which controls the gauge field amplification. While this quantity is nearly constant in the Anber-Sorbo solution, in our numerical evolutions ξ (t) "oscillates" about the Anber-Sorbo value, without showing signs of convergence. This evolution has some similarities with the transition between overshooting and undershooting evolution in tracking quintessence [60], with the big difference that in this second case the scalar field eventually reaches the steady state tracking solution, in which the Hubble friction balances against the slope of the potential. As a consequence, the tracking solution is an attractor for these models. The same is true for the inflationary solution in standard slow roll inflation. Our numerical evolutions do show sign of this convergence for the Anber-Sorbo solutions. Analogous "oscillations" between fast and slow evolution were also observed in refs. [39], that studied the transition during inflation between Hubble dominated friction and gauge field production dominated friction, and in refs. [41], that assumed that the scalar inflaton field had zero initial speed. These studies were limited to the inflationary case. One could have hoped that in our study, the extra friction due to the initially dominating matter field could have "gently" accompanied the system toward the steady state Anber-Sorbo evolution. However, we could not see this, within the limits of what we could numerically simulate. We believe that this is worth further study.
where we have normalized the gauge filed mode function so to have magnitude 1 at early (rescaled) times, y → 0. This is the equation of an oscillator with the time dependent frequency ω (y) ≡ 1 − y 5 . At the initial time, the frequency varies adiabatically, dω dy ω 2 , and the early time evolution ofÂ is well described by the adiabatic solution A adiabatic = e −i y dy ω = exp −i 2y 1 − y 5 7 + 5 7 where 2 F 1 is an ordinary hypergeometric function (according to the usual vacuum, prescription, we keep only the positive frequency term). The adiabaticity condition is violated at y −1/5 , but it is then valid again at later times. At these late times we therefore expect the solution to be a linear combination of the growing mode e + y dy √ −ω 2 and of the decreasing mode e − y dy √ −ω 2 . To obtain an accurate late time solution one should find the two coefficients of this linear combination, by studying the solution in the non adiabatic stage.
Our present goal is to obtain a rough estimate of the backreaction term, to estimate the range of parameters needed for the full numerical evolution (which, eventually, is the only computation that we want to perform accurately). Therefore, we simply use the growing term as an estimate, which is real (as it can be inspected directly). The integrand in the E · B backreaction term is proportional to the time derivative of the magnitude square of the mode function, namely to d dy |Â| 2 , up to constant rescalings.
We computed this quantity for the (A.3) estimate, and we compared it with the one obtained from a numerical integration of eq. (A.1). We found that a much better agreement is obtained when d dy |Â| 2 is divided by In Figure 11 we compare the approximation (A.4) against the numerical integration of eq. (A.1). As we are interested in λ,k = O (1), the parameter is of O f 5 M 5 p , which, as we should in the main text, could be as small as 10 −20 in the cases that we consider. As τ ranges from 0 to a few, the combination y 5 also ranges in this interval. We see that the approximation (A.4) is very accurate, and clearly adequate for our estimates.
In particular, we use this result to estimate the backreaction term in the evolution equation (3.4) of the scalar field, that in rescaled variables, reads We define as measure of backreaction the ratio between the fourth and third term in this equation, that we write as in eq. We plot this quantity in Figure 3 in the main text.

B System of equations solved
In this Appendix we present the system of equations that we use in our numerical integrations, using the number of e-folds of expansion N as "time variable", so that our numerical results can be more easily reproduced. To derive this system, we started from the equations in physical time, t, rescaled to the dimensionless combinatioñ t ≡ √ 3V 0 2M p t .
(B.1) so thatt = 1 corresponds to the estimated transition time (3.15). 12 This time variable is used to introduce the dimensionless Hubble ratẽ As in the main text, we introduced the dimensionless scalar field ϕ ≡ φ/M p , and we rescale the axion decay constant analogously,f ≡ f /M p . We then rescale the comoving momentum of the gauge modes ask As shown in the previous appendix,k = 1 corresponds to our estimate for the greatest momentum that is excited at the estimated transition timet = 1. Finally, we rescale the gauge field mode functions as Due to the √ 2k factor, the rescaled mode functionÃ has initially magnitude one. Moreover, the last factor rescales away the phase due to the early time solution, so that the mode function is nearly constant in the early time regime (namely, it only varies due to the last term in (3.2), which vanishes at asymptotically early times); this considerably speeds up the numerical evaluation of the gauge modes in the early time / UV regime.
Using the definition (3.18) it is immediate to see that a derivative wrt (rescaled) physical time can be related to a derivative wrt the number of e-folds via d dt =H d dN . We can then write the evolution equations for the scalar field, a second order equation for the scale factor (obtained from a combination of the 00 and of the spatial diagonal part of the background Einstein equations 13 ), and the equation for the gauge field modes, as differential equations in terms of N . Some lengthy but straightforward algebra gives d 2 ϕ dN 2 + 1 − This whole set of equations is written in terms of dimensionless quantities, and it is ready to be integrated numerically (starting from the initial moment N = 0). We note that the two backreaction terms introduce the parameter V 0 /M 4 p . The vector field energy density enters in this system only through the initial condition. This quantity is divergent, being the vacuum energy density of the gauge fields before they undergo any amplification. In principle, one should subtract it via some regularization scheme [59]. However, we only include modes in our numerical evolution for which this term is completely negligible, so, for our practical purposes, this term can simply be disregarded.