Chameleons in the Early Universe: Kicks, Rebounds, and Particle Production

Chameleon gravity is a scalar-tensor theory that includes a non-minimal coupling between the scalar field and the matter fields and yet mimics general relativity in the Solar System. The scalar degree of freedom is hidden in high-density environments because the effective mass of the chameleon scalar depends on the trace of the stress-energy tensor. In the early Universe, when the trace of the matter stress-energy tensor is nearly zero, the chameleon is very light, and Hubble friction prevents it from reaching the minimum of its effective potential. Whenever a particle species becomes non-relativistic, however, the trace of the stress-energy tensor is temporarily nonzero, and the chameleon begins to roll. We show that these"kicks"to the chameleon field have catastrophic consequences for chameleon gravity. The velocity imparted to the chameleon by the kick is sufficiently large that the chameleon's mass changes rapidly as it slides past its potential minimum. This nonadiabatic evolution shatters the chameleon field by generating extremely high-energy perturbations through quantum particle production. If the chameleon's coupling to matter is slightly stronger than gravitational, the excited modes have trans-Planckian momenta. The production of modes with momenta exceeding 1e7 GeV can only be avoided for small couplings and finely tuned initial conditions. These quantum effects also significantly alter the background evolution of the chameleon field, and we develop new analytic and numerical techniques to treat quantum particle production in the regime of strong dissipation. This analysis demonstrates that chameleon gravity cannot be treated as a classical field theory at the time of Big Bang Nucleosynthesis and casts doubt on chameleon gravity's viability as an alternative to general relativity.


I. INTRODUCTION
Light scalar fields are of great interest in cosmology because they arise in many explanations for the current acceleration of the expansion of the Universe [1][2][3][4][5]. It is challenging for these models to evade the stringent experimental limits on fifth forces within the Solar System and the laboratory, however, because no theory has been constructed that both explains current cosmological observations and forbids interactions between the scalar field and Standard-Model particles. Axionic quintessence models come the closest, because they possess a shift symmetry that forbids a direct coupling to the stress-energy tensor of Standard-Model particles, but they still interact with photons [6]. In all other cases, we must reconcile ourselves to a coupling between matter and the scalar field; unless such a coupling is forbidden, we must include it in our theory as it will be generated by quantum effects. Problematically, the presence of a new light scalar field coupled to matter usually implies the existence of a new long-range fifth force, and no new forces have been seen in either laboratory experiments or Solar-System tests of general relativity. The precision of these experiments constrains the strength of any new force to be many or-ders of magnitude weaker than gravity [7]. In a simple Yukawa model, this constraint forces the energy scale that controls the strength of the coupling between the scalar field and matter to be many orders of magnitude above the Planck scale. Such a large energy scale is almost impossible to justify in any reasonable effective field theory.
In 2003 Khoury and Weltman proposed chameleon gravity, which contains a scalar degree of freedom whose potential can provide the vacuum energy required for cosmic acceleration [8,9]. The chameleon scalar field is a light field that interacts with the matter fields of the Standard Model, but it possesses a dynamical mechanism to hide these interactions in dense environments. The chameleon's potential function contains non-linear terms that, when combined with the chameleon's coupling to matter fields, make the chameleon's effective mass dependent on its environment. The chameleon is heavy in dense environments, which suppresses its ability to mediate a fifth force. Such dynamical mechanisms for suppressing fifth forces are known as screening mechanisms; the chameleon mechanism is one of only three known screening mechanisms capable of making scalar-tensor gravitational theories compatible with experimental constraints on fifth forces [10]. Furthermore, the chameleon mechanism is essential for f (R) theories of modified gravity, which generate cosmic acceleration by making the gravitational Lagrangian a non-trivial function of the Ricci scalar [11]. Such a theory can be rewritten as a metric theory with the standard Einstein-Hilbert action and an additional scalar field that couples to matter in the same way as the chameleon [12]. In order for an f (R) theory to successfully pass observational tests, this scalar field must have a potential function that allows it to employ the chameleon mechanism [13][14][15][16].
In the original chameleon theory, and in f (R) gravity, the coupling between the chameleon field and the matter fields was assumed to have gravitational strength. 1 Later it was found that much stronger couplings are also allowed [17,18]; when the coupling to matter is stronger, the screening mechanism is also stronger, and the scalar field can still be hidden from fifth-force experiments. We can search for strongly coupled chameleons in high-precision low-energy photon experiments [19][20][21][22][23][24], with ultra-cold neutrons [25][26][27], in precision atomic measurements [28], in Casimir force experiments [29,30], with dark-matter direct-detection experiments [31], and in particle colliders [32,33]. It has also been suggested to look for strongly coupled chameleons produced in the Sun [34,35] and to seek chameleon signatures in observations of stars and galaxies [34][35][36][37][38][39], the cosmic microwave background [40][41][42], and the 21 cm power spectrum [43]. These experiments exploit the fact that strongly coupled chameleons interact strongly with matter particles and photons in vacuum, so if an experiment or an astrophysical observation is targeted at a diffuse environment, it has the potential to see a chameleon signal. A number of these experiments have been purposely designed to look for chameleons [23,30,31,35], while other results come from exploiting measurements made for other purposes.
Chameleons with gravitational-strength couplings in vacuum are harder to detect directly and are best sought by searching for deviations from general relativity. Constraints on these theories come from laboratory searches for modifications of gravity [44,45] and from astrophysical observations, including constraining the effects of the chameleon on the formation of structure and the current matter power spectrum [46][47][48][49], on weak lensing [50], and on the evolution of stars [51]. We will summarize the best constraints on chameleon theories in Section II B; for our purposes, the essential constraint is that the chameleon potential must have a steep section in which small changes in the chameleon field (∆φ ≃ 0.01 eV) lead to significant changes in the chameleon potential and its derivatives.
The chameleon potential was designed to provide the chameleon screening mechanism and does not originate from fundamental physics. An approach to constructing a chameleon model from a KKLT compactification of string theory was discussed in Ref. [52] and extended in Refs. [53,54]; previous attempts to embed the chameleon into UV-complete theories were unsuccessful [55,56]. In the absence of a UV-complete theory, chameleon gravity is usually treated as an effective field theory that should only be trusted at relatively low energies. Quantum corrections to this theory have largely been ignored, even though one-loop corrections to the chameleon potential can be significant in laboratory environments [57], and oscillations of the curvature scalar in f (R) gravity can lead to particle production [58,59].
In a recent letter [60], we exposed an additional quantum instability in chameleon gravity: the chameleon's behavior just prior to the time of Big Bang Nucleosynthesis (BBN) triggers catastrophic quantum effects that transfer most of the chameleon's energy to perturbations with momenta greater than 10 7 GeV. Increasing the strength of the chameleon's coupling to matter increases the energies of the generated perturbations, and chameleons with matter couplings that are moderately stronger than gravitational interactions experience trans-Planckian excitations. In this work, we provide a more detailed treatment of this phenomenon, including the derivations that were omitted from Ref. [60]. We also extend our analysis to power-law chameleon potentials and find that they generate even more energetic perturbations than the exponential potentials studied in our earlier work.
The chameleon's behavior in the early Universe was first investigated in Ref. [61]. Of particular concern is how much the chameleon scalar field evolves between the time of BBN and the present day. Large variations in the chameleon's value can be interpreted as large variations in particle masses, and yet we know that particle masses at the time of BBN do not significantly differ from the masses that we measure today [e.g. 62,63]. Ref. [61] found that the chameleon is driven toward its current value prior to BBN, thus ensuring that the nucleon masses are sufficiently close to their observed values that BBN is unaffected, regardless of the chameleon's initial value (but also see Ref. [64]). Although the chameleon is usually light while the Universe is radiation dominated, the field is able to overcome Hubble friction and approach its present-day value because it becomes momentarily heavy whenever the Universe's temperature equals the mass of a particle species in equilibrium with the radiation bath. We will discuss how these mass thresholds dramatically perturb the dynamics of the chameleon field in Section II; in summary, they kick the chameleon scalar field closer to the minimum of its effective potential, thus enabling it to approach the value it holds today.
We will show that these kicks are generally too effective; the chameleon reaches the minimum of its effective potential with a large velocity (φ ≫ MeV 2 ) and climbs up the steep part of its potential. Since the chameleon potential changes significantly when the chameleon value changes by 0.01 eV, these large velocities lead to rapid changes in the chameleon's effective mass, which generate perturbations via quantum particle production. These perturbations have sufficiently high energies that they push chameleon gravity outside its low-energy regime of validity, and quantum correc-tions dominate the chameleon's potential. Therefore, the chameleon's evolution during BBN cannot be understood using only a low-energy effective field theory, which casts doubt on chameleon gravity's viability. Previous studies of the chameleon's evolution prior to and during BBN [61,64] treated the chameleon purely classically and consequently missed these important quantum effects.
We begin by reviewing chameleon gravity in Section II; we focus on the shape of the chameleon potential and the chameleon's dynamics in a radiation-dominated Universe. In Section III, we present a novel solution to the chameleon's equations of motion in the presence of the aforementioned mass-threshold kicks. We apply this solution in Section IV, where we consider how the chameleon responds to the kicks generated by Standard-Model particles, and we calculate the chameleon's velocity when it reaches the minimum of its effective potential. In Section V, we show that these velocities lead to non-adiabatic changes in the chameleon's effective mass, and we investigate the resulting particle production both numerically and analytically. Finally, we discuss the implications of our results and conclude in Section VI. Appendices A and B provide further details about the kicks from Standard-Model particles and the chameleon field's evolution at high temperatures, and we review the fundamental theory of quantum particle production in an expanding Universe in Appendix C.

II. CHAMELEON GRAVITY
In chameleon gravity, the spacetime metricg µν that appears in the matter Lagrangian is a conformal rescaling of the metric g * µν that solves Einstein's equations: where β is a dimensionless coupling constant, φ is the chameleon field, and M Pl = (8πG) −1/2 . The action for this theory can be written as where g * is the determinant of the metric g * µν , R * is its Ricci scalar, V (φ) is the chameleon potential, and S m is the action for the matter fields. If we define T * µν ≡ (−2/ √ −g * )δS m /δg µν * , then varying this action with respect to g * µν yields the Einstein field equation with a stress-energy tensor equal to the sum of T * µν and the stress-energy tensor for the chameleon field. Consequently, g * µν and T * µν are respectively called the Einsteinframe metric and stress-energy tensor. Meanwhile, the matter fields couple to the Jordan-frame metricg µν . We can also define a Jordan-frame stress-energy tensor: T µν ≡ (−2/ √ −g)δS m /δg µν . The conformal relationship betweeng µν and g * µν implies that T µ * ν = e 4βφ/M PlT µ ν .
Therefore, if the matter fields are perfect fluids with density ρ and pressure p, ρ * = e 4βφ/M Plρ , and the equation of state parameter w ≡ p/ρ is the same in both frames.
A. The chameleon potential The choice of the potential V (φ) is crucial to the success of the chameleon mechanism, as it is the nonlinearities of the potential that allow the chameleon field to hide from fifth-force experiments. Here we give a brief review of the properties required of a chameleon potential. Due to the coupling between the chameleon field and matter, the chameleon explores a wide region of its potential, and its evolution is driven by the ambient energy distribution. First, it is necessary for there to be a region of the potential that is close to being flat; if we want the chameleon to be cosmologically relevant today, then the field must have a very small mass in cosmological environments. Second, the potential must have a steep section to provide a barrier that limits the chameleon's excursion from its cosmological value in the interior of the objects used in fifth-force searches. Such objects include the Earth, the Moon, the Sun, and laboratory test masses.
To study such situations it is easiest to assume that the source object is spherically symmetric, static, and composed of non-relativistic matter. Then the equation of motion governing the behavior of the chameleon in the Einstein frame is Therefore, the chameleon's evolution is governed by an effective potential This effective potential has a minimum at φ min , where V ′ (φ min ) = −βρ * /M Pl , which implies that φ min depends on the ambient matter density ρ * . At this minimum, the effective mass of the chameleon field is For the chameleon mechanism to operate, m 2 must be a positive and monotonically increasing function of the density ρ * for all relevant densities. Since φ min increases as the density decreases, V ′′′ (φ) must be negative over the range of accessible φ values. It is not sufficient for the effective mass of the chameleon to simply increase as the density increases; screening the fifth force mediated by the chameleon requires the mass to increase sharply. As discussed above, V (φ) must be chosen so that, when the density approaches the cosmological background density, the chameleon's Compton wavelength approaches cosmological distance scales. To show that the mass must be much greater inside a test mass, we consider the condition for the scalar potential well generated by a source to be shallower than that source's gravitational well. This condition, known as the thin-shell condition, guarantees that the fifth force mediated by the scalar field will be much weaker than the gravitational force [8,9]. If the source object has mass M c and radius R, then the thin-shell condition is where φ c and φ ∞ are the positions of the minimum of the effective potential inside and far outside the source object, respectively. Since V ′′ (φ) must be positive over the relevant range of φ values, To get a rough estimate of a bound on the mass of the chameleon inside the source object [m 2 c = V ′′ (φ c )], we Taylor expand the right-hand side of this inequality. On rearranging, and assuming ρ ∞ ≪ ρ c , we find Although this inequality only approximates the thin-shell condition, it illuminates the essential component of the chameleon screening mechanism: inside a source, the chameleon is too massive to carry a force beyond the source's boundary. Since 1/R ≫ H for all sources, there must be a region of the chameleon potential that is much steeper than the region probed by cosmological densities. This feature of the chameleon potential will be crucial in the discussion that follows.

B. Current constraints on chameleon models
Chameleon gravity is constrained by laboratory experiments, gravitational tests in the Solar System, and astrophysical observations. The best current bounds on the coupling parameter β come from laboratory experiments that study diffuse systems in a vacuum. The chameleon behaves as a very light scalar field within a laboratory vacuum, and none of its effects are screened. Consequently, precision measurements of very diffuse systems can detect signatures of the chameleon. The best current constraint comes from measurements of atomic energy levels in hydrogen [28], which would be perturbed by the existence of a new chameleon force. These measurements set an upper bound on the coupling parameter: β 10 14 . For reference, gravitational-strength coupling corresponds to β ∼ O(1), and f (R) gravity theories have β = 1/ √ 6.
Constraints on the energy scale that controls the chameleon potential, which we will denote M , are more model dependent. The best constraints come from laboratory searches for fifth forces and from Casimir experiments. The Eöt-Wash experiment currently provides the best constraints on weakly coupled chameleon theories (β ∼ < 20) [44,45]. Over a wide range of the parameter space, however, these fifth-force experiments are not sensitive to the chameleon. In these experiments, thin plates are often used to shield electromagnetic forces, but unfortunately, these plates often shield the chameleon forces too. Experimental searches for Casimir effects do not use such shielding; since they look for forces between parallel plates held very closely together, they can be very constraining for chameleon theories. The constraints on M typically depend on the value of β and the precise form of the chameleon potential [18]. It is possible to make some general statements, however: for power-law potentials of the form V (φ) = M 4+n φ −n and β 0.1, searches for Casimir effects and fifth forces constrain M ≤ 0.01 eV. For 0.1 β 10 5 the constraints on M are typically much stronger in specific models; for full details of the constraints on specific choices of the chameleon potential, we refer the reader to Refs. [18,44,45].
In our analysis, we will consider both the power-law potential and the exponential potential considered in previous studies of the chameleon's cosmological evolution [61]: with n > 0. If M ∼ 0.001 eV, this potential provides the vacuum energy required to drive cosmic acceleration at late times. In low-density environments, φ ≫ M , and this potential is effectively a power-law potential plus a constant. Therefore, it is subject to the same laboratory constraints as power-law potentials: M < 0.01 eV. This potential meets all the requirements discussed in the previous subsection: it is nearly flat when φ ≫ M ; it is steep when φ ∼ < M ; V ′′ (φ) is always positive; and V ′′′ (φ) is always negative. Constraints on this potential were analyzed in detail in Refs. [44,45], where it was found that if M is chosen to be the dark energy scale, constraints from the Eöt-Wash experiment demand fairly large values for n and β: roughly n ∼ > 10 for β ∼ 1, and n 4 for β ∼ 10. Finally, there is an astrophysical constraint on the present-day cosmological mass of the chameleon (m ∞ in the notation of the previous section) [65,66]. Since the chameleon force must be screened inside galaxies, galaxies must satisfy the thin-shell condition given by Eq. (6). Satisfying this constraint requires m ∞ /H 0 10 3 , which corresponds to 1/m ∞ 10 Mpc. For both the power-law potential [V (φ) = M 4+n φ −n ] and the potential given by Eq. (9), this constraint implies that M ∼ < 0.07β 2/3 MeV for n = 2, M ∼ < 2β 3/4 GeV for n = 4, and M ∼ < 1000β 4/5 GeV for n = 6.

C. A cosmological chameleon
If the Universe is homogeneous and isotropic, then the Einstein-frame and Jordan-frame metrics are FRW metrics with conformally related scale factors (ã = e βφ/M Pl a * ) and proper times (dt = e βφ/M Pl dt * ). In the Jordan frame, the matter fields do not interact with the scalar fields, so the matter stress-energy is conserved:∇ µT µ ν = 0 andρ ∝ã −3(1+w) . It follows that : the energy density in radiation is proportional to a −4 * in the Einstein frame, but the Einstein-frame energy density in matter is not proportional to a −3 * . While T µ * ν is not conserved in the Einstein frame, the sum of T µ * ν and the stress-energy tensor for the scalar field φ is conserved.
Varying Eq. (2) with respect to φ gives the chameleon equation of motion where a dot represents differentiation with respect to Einstein proper time t * and H * ≡ȧ * /a * . In the second line, we have evaluated the trace of the stress-energy tensor T µ * µ : ρ * R is the Einstein-frame energy density of the cosmic radiation bath, f m is the Einstein-frame density of nonrelativistic matter divided by ρ * R , and Σ ≡ (ρ * R − 3p * R )/ρ * R , where p * R is the Einstein-frame pressure of the radiation bath. Since both f m and Σ are ratios of elements of the stress-energy tensor, the conformal relationship betweenT µ ν and T µ * ν implies that these ratios are the same in the Einstein frame and the Jordan frame, so f m and Σ may be evaluated using Jordan-frame energy densities and pressures. As in Eq. (4), we can use Eq.(11) to define an effective potential for the chameleon, which is minimized when φ = φ min . If M ∼ 0.001 eV, φ min < M in the pre-BBN Universe.
If the radiation bath only consisted of photons, then Σ would be zero. In the early Universe, however, several massive particles were in thermal equilibrium with the photons, and we include the energy densities of these particles in ρ * R . When the temperature of the radiation is much larger than the mass of the particle, these particles are relativistic, and their contribution to Σ is zero. As the radiation cools, the particles' pressure decreases faster than their energy density and their contribution to Σ increases. When the temperature is much less than the mass of the particle, the particles are Boltzmann suppressed and their contribution to Σ decreases again. Therefore, each species of massive particles makes a contribution to Σ that peaks when the temperature of the radiation bath is nearly equal to the mass of the particle. In Appendix A, we evaluate (ρ i − 3p i )/ρ R for a particle with mass m i and g i degrees of freedom that is in thermal equilibrium with a radiation bath at temper-ature T J [67][68][69]: J ] −1 is the number of relativistic degrees of freedom. In the denominator of the integrand, the + sign applies to fermions, and the − sign applies to bosons. While m i ≃ T J and Σ = 0, the chameleon experiences a force that drives it to smaller φ values; Σ effectively "kicks" the chameleon.
To numerically solve the chameleon's equation of motion, we will need to specify how the Jordan-frame temperature depends on a * and φ. It is useful to consider the entropy density of the radiation bath, s R = (ρ R +p R )/T J , and to define g * S ≡ s R [(2π 2 /45)T 3 J ] −1 . Entropy conservation in the Jordan frame implies that g * S (T J )ã 3 T 3 J is constant, which gives us an (implicit) expression for T J in terms of Einstein-frame variables, including the values of a * and φ at some fixed time [φ(a * ,i ) ≡ φ i ], and the Jordan-frame temperature T J,i at that same time: (13) We evaluate g * S (T J ) for the Standard Model particle spectrum (see Appendix A for details), and then we numerically invert the function f (T J ) = T J × [g * S (T J )] 1/3 to obtain T J . The initial temperature T J,i is chosen so that Σ + f m ≪ 1 for T J > T J,i . We expect that the chameleon is at rest prior to the onset of the kicks, because any velocity it may have obtained during reheating would be damped by Hubble friction [61]. In this case, the chameleon will remain at rest while Σ + f m ≪ 1, so its subsequent evolution does not depend on the specific value of T J,i .
We assume a rather generic initial condition for φ: M ≪ φ i ∼ < M Pl . For φ ∼ < M , the chameleon potential is very steep, and if φ is less than φ min in the early Universe, the chameleon will quickly roll to larger field values, where it will eventually stick due to Hubble friction. This evolution was demonstrated explicitly in Ref. [61]; if φ i is less than φ min prior to BBN, then the driving term from the chameleon potential dominates over the frictional Hubble term, and the field rolls until it stops at a value where Ω (i) φ < 1 is the initial fraction of the Universe's energy density in the chameleon field. For the purposes of our analysis, φ i = φ stop in this scenario, and thus we expect φ i ∼ < M Pl if φ ≃ φ min during inflation, when φ min ≪ M . It is also conceivable that φ had an initial value that far exceeded M . We still restrict our analysis to φ i < M Pl because we expect that matter loop corrections will significantly renormalize the bare chameleon potential for larger field values. We note, however, that our analysis is applicable to larger values of φ i if one is willing to consider φ i > M Pl in spite of these concerns. In fact, we will frequently set φ i > M Pl in Section IV in order to illustrate how far the chameleon rolls during the kicks.
A change of variables facilitates our analysis of the chameleon's evolution. We define a new time variable: p ≡ ln(a * /a * ,i ), and we use a prime to denote differentiation with respect to p. We also define a dimensionless scalar field ϕ ≡ φ/M Pl . With these definitions, and using the Friedmann equations in the Einstein frame, the equation of motion for φ [Eq. 11] becomes which we can use to eliminate H 2 * from Eq. (15). We can further simplify the chameleon's equation of motion by noting that Σ ∼ < 0.1 and f m ∼ < 10 −6 prior to BBN (see Section IV). Therefore, we can approximate 1 + f m ≃ 1 + Σ ≃ 1, which leaves Finally, to close the system of equations, we note that Eq. (13) implies that which explicitly shows that, if there are no changes in the number of relativistic degrees of freedom, ρ * R ∝ a −4 * as expected.
In the next two sections, we use this system of equations to examine the evolution of the chameleon field during the radiation-dominated era. In the absence of massive particles, so that both Σ and f m are zero, analytically solving Eq. (17) for M ≪ φ i ∼ < M Pl reveals that any initial field velocity possessed by the chameleon will quickly damp away to zero, and the field will freeze at a value φ > φ min that is fixed by its initial conditions [61]. We will see that the kicks dislodge the chameleon and send it rolling toward the minimum of its effective potential.

III. THE SURFING SOLUTION
In this section, we analyze the chameleon's response to the kicking function Σ analytically, and we expose a new solution to its equation of motion. First, we simplify Eq. (17) by noting that V (φ) ≪ ρ * R while φ ∼ > M in the early Universe. Second, we assume that φ ≫ φ min , as is the case for the range of initial conditions that we consider (M ≪ φ i M Pl ), so that the dV /dϕ driving term in Eq. (17) is negligible compared to the driving term from Σ. Finally, we neglect the background density of non-relativistic matter f m . These approximations simplify the system and allow us to analytically examine the dynamics of the kicks. However, it is not necessary to make such assumptions in order to study the evolution of the chameleon kicks numerically, and the numerical results of the next section do not require such assumptions. Working with the simplified system, the chameleon equation of motion (17) is approximately In deriving this equation, we have assumed that ϕ ′ < √ 6, which is equivalent to assuming that the chameleon's kinetic energy is less than the critical density of the Universe.
Recall that the Jordan-frame temperature T J is given by Eq. (13) as a function of the Einstein-frame scale factor, the scalar field value, and the number of relativistic degrees of freedom: The dependence of the Jordan-frame temperature on the chameleon field is important because it allows the existence of a novel solution to the chameleon equation of motion, which we call the surfing solution: where p s is the time at which the surfing behavior begins, and the value of the field at this time is where λ is a constant. Inserting this ansatz for ϕ into Eq. (20) yields which implies that the Jordan-frame temperature is constant while the chameleon follows the surfing solution.
We call this constant value of T J the surfing temperature (T s ). Since the surfing solution has ϕ ′′ = 0 and Therefore, for a given kick function Σ, the surfing temperature is determined by β, and then λ is set by Eq. (23). On the surfing solution, the value of the scalar field decreases with time p, so eventually the field value approaches φ min , and the bare scalar potential is no longer negligible. At this point our approximations break down, and the surfing solution ceases to exist. It is important to stress, however, that once the chameleon reaches the surfing solution, it will remain on that solution until the scalar field gets close to the minimum of its effective potential. The name "surfing solution" was chosen because the chameleon field surfs the wave of the kick function all the way to the minimum of its effective potential.
While the chameleon is surfing, the Jordan-frame scale factorã = e βϕ a * remains constant, and the Jordan-frame Universe is static. The Einstein-frame scale factor does continue to increase, but this expansion has no observable effects on particles or Jordan-frame energy densities. The surfing solution effectively pauses the evolution of the Universe from the time the Jordan-frame temperature reaches the surfing temperature to the time when φ ≃ φ min . Despite the interruption in the Jordan-frame expansion, we stress that time in both the Jordan and Einstein frames continues to evolve forward during the surfing phase.
The surfing solution is relevant for the cosmological evolution of the chameleon field only if it is an attractor in the space of solutions. Otherwise, the field only surfs the kicks if an unlikely fine-tuning of initial conditions occurs. It can be seen that the surfing solution is reached regardless of the initial value of ϕ by noting that, while V (ϕ) can be neglected, the equations of motion are invariant under the transformations for constant C. Therefore, all initial values for ϕ are equivalent up to a time translation; changing ϕ i or, equivalently, changing T J,i changes the value of p s and λ but does not change the existence of the surfing solution or the surfing temperature. To see that the surfing solution is an attractor as the initial field velocity is also varied, it is easiest to solve the field equations numerically and plot the phase portraits, as shown in Fig. 1. We take numerical values of the parameters that approximate the kick coming from the electron and positron (discussed in more detail in Appendix A): g i = 4 and g * S = g * = 10.75. With these parameters, the maximum value of the kick function is Max[Σ(T )] ≈ 0.062, so surfing solutions exist for all β 2.3. In Fig. 1, we show phase portraits both for β = 2, for which a surfing solution does not exist, and β = 3, which has a surfing solution. We consider initial field velocities with |ϕ ′ i | ∼ < √ 3, as is required to prevent the scalar field's kinetic energy from dominating the Universe (Ωφ ∼ < 0.5). When generating these figures, we have omitted V (φ), so there is nothing preventing ϕ from going negative.
The top panel of Fig. 1, with β = 2, shows that the kick moves the chameleon field to smaller values if ϕ ′ i ∼ < 1, and in all cases, the field's velocity goes to zero after the kick passes. In the bottom panel, we see that all chameleons with ϕ ′ i ∼ < 1 (corresponding to Ωφ ∼ < 1/6) end up on the surfing solution, with ϕ ′ = −1/β. These plots show that, when the surfing solution exists, it is an attractor in the space of solutions for all except the largest positive values of ϕ ′ i . Although we only display phase portraits for one particular choice of parameter values, this behavior persists for surfing solutions over the whole range of parameter space for Σ.

IV. KICKS FROM THE STANDARD MODEL
We now consider how the chameleon responds to the kicks generated by the Standard-Model particle spectrum. As described in Appendix A, we evaluate Σ by summing contributions from the particles in the Standard Model, including a Higgs particle with a mass of 125 GeV. Figure 2 shows the resulting Σ(T J ). We see that the individual kicks from different particles are not distinct events; instead, the Standard-Model particles produce four "combo-kicks." Each combo-kick has a larger amplitude than the previous kicks because each particle's contribution to Σ is suppressed by a factor of 1/g * (T J ) (see Eq. 12), and the number of relativistic degrees of freedom decreases as the Universe cools. The discontinuity between the second and third combo-kicks arises from the QCD phase transition, which we assume happens instantaneously at a temperature of 170 MeV. The longest pause between kicks occurs prior to the last kick, when Σ reaches a minimum value of 0.00026 at a temperature of 7.4 MeV. Even at this temperature, Σ is much larger than the matter fraction f m . Moreover, even if we assume that the current matter content of the Universe, including dark matter, is decoupled and nonrelativistic at all temperatures, f m ≪ Σ for all temperatures greater than 50 keV. Therefore, f m does not affect the chameleon's evolution during the kicks, and we do not consider it further.
Our calculation of Σ(T ) at temperatures above 100 MeV is an incomplete treatment that provides a minimal value for the kick function. First, it underestimates Σ(T ) during the QCD phase transition. Lattice QCD calculations indicate that contributions from other hadrons and interactions between fields cause Σ to increase sharply during the QCD phase transition, reaching values between 0.2 and 0.4 [70][71][72]. Due to the discrepancies between different lattice QCD calculations of Σ, we choose to neglect these additional contributions. We note, however, that this additional peak in the kick function would lower the minimal value of β required for the surfing solution, and it would enhance the impact velocity of the chameleon field for models that reach the minimum of their effective potential at temperatures below 400 MeV. Second, we do not include contributions from particles beyond the Standard Model; if nothing else, the dark matter particle should contribute to Σ. Third, we neglect the electroweak phase transition and use the particle spectrum given in Appendix A at high temperatures. Ref. [72] showed that this approximation differs only slightly from a one-loop treatment of electroweak thermodynamics [73] for temperatures less than 100 GeV, and the approximation is accurate within an order of magnitude for temperatures between 100 GeV and 300 GeV. At these high temperatures, the potential contributions from beyond-Standard Model particles dwarf our calculation of Σ(T ), making our neglect of the electroweak phase transition irrelevant. Finally, we do not include the QCD trace anomaly [74,75]. In the perturbative regime of QCD (i.e. energies above 100 GeV), the QCD trace anomaly implies that (ρ * R − 3p * R )/ρ * R ≃ 0.001 even if all components of the plasma are relativistic. At temperatures less than the electroweak phase transition, this contribution to Σ is much smaller than the contributions from the Standard-Model kicks, and it does not significantly affect the chameleon's evolution if β ∼ < 5. For more strongly coupled chameleons, we show in Appendix B that the trace anomaly significantly increases the chameleon's velocity toward the minimum of its effective potential.
In Section III, we showed that chameleons with β > 1/(3Σ max ), where Σ max is the maximum value of Σ during the kick, will "surf" the kick and approach the minimum of the effective potential with a velocity dφ/d ln a * = −M Pl /β. The amplitude of the first combo- kick (due to the top quark and the W, Z, and Higgs bosons) implies that all chameleons with β > 3.05 can surf this kick, but numerically solving Eq. (17) reveals that the chameleon reaches the surfing solution during the first kick only if β ≥ 3.07. If β < 3.07, then the first combo-kick will push the chameleon toward the potential minimum, but as the Jordan-frame temperature cools, Σ will decrease, and the chameleon will eventually come to a halt at a new position. When we consider nonsurfing chameleons later in this section, we will compute how far the chameleon moves during a kick that it cannot surf. For now, let us assume that the value of the chameleon prior to the first combo-kick was sufficiently large that φ ≫ φ min after the passage of all the kicks that the chameleon cannot surf. In that case, the chameleon will surf subsequent kicks if β > 1/(3Σ max ) for these kicks. Chameleons with β > 1.82 can surf the final kick, which occurs when the electrons and positrons become non-relativistic. Since this kick has the largest amplitude of the four combo-kicks we consider, chameleons with β < 1.82 cannot surf. As previously mentioned, however, the peak in Σ due to the QCD phase transition would extend the surfing solution to smaller values of β.
The existence of the surfing solution guarantees that all chameleons with β > 1.82 will be kicked to φ ∼ < φ min , regardless of their initial field values. If β ≥ 3.07, then the chameleon surfs the first combo-kick, and it approaches the minimum of its effective potential with a constant velocity. If 1.82 < β < 3.07, then the chameleon's velocity at φ ≃ φ min depends on its value prior to the kicks (φ i ). If φ i is greater than the displacement caused by the previous kicks, then the chameleon will reach the FIG. 3: The value of |φ| and Ts when φ ≪ M Pl for chameleons that follow the surfing solution described in Section III. This scenario is unavoidable if β ≥ 3.07, and it also applies to smaller values of β if the initial value of φ is sufficiently large. The computation of these velocities does not include the contribution to Σ from the QCD trace anomaly and also ignores the electroweak phase transition. minimum of its effective potential by surfing the first kick for which β > 1/(3Σ max ), and it too will approach φ min with a constant velocity. For all surfing chameleons, the chameleon's velocity depends only on where T s is the temperature in the Jordan frame during the surfing phase: Σ(T s ) = 1/(3β 2 ). In deriving this equation, we assumed that exp[4βφ/M Pl ] ≃ 1 so that we could equate ρ * R to the radiation density in the Jordan frame. Therefore, this equation is only applicable when φ min ∼ < φ ≪ M Pl . Figure 3 showsφ and T s for the surfing solution given the Standard-Model particle spectrum and no contribution from the QCD trace anomaly. We see thatφ(β) is a smooth function for β ≥ 3.07; these chameleons surf the first kick. Chameleons with 1.82 < β < 3.07 values must wait for Σ(T J ) to equal 1/(3β 2 ) near the peaks of subsequent kicks. Since there are no surfing solutions in the gaps between the kicks, T s (β) andφ(β) are discontinuous for β < 3.07. If β > 4.5, then the surfing temperature T s exceeds 150 GeV, and the chameleon is sensitive to the details of the electroweak phase transition and the QCD trace anomaly. Since Σ decreases as the temperature increases for T J > 100 GeV, adding the trace anomaly will increase the value of T s that satisfies Σ(T s ) = 1/(3β 2 ) for a given value of β. From Eq. (27), we see thatφ ∝ T 2 s , so the trace anomaly increases the value of |φ| during the surfing phase. The trace anomaly can have a more profound impact if its contribution implies that Σ > 1/(3β 2 ) at all temperatures. This scenario is discussed in Appendix B, where we show that a constant high-temperature plateau in Σ with Σ > 1/(3β 2 ) implies that the temperature in the Jordan frame increases as the chameleon rolls toward the minimum of the effective potential, leading to a drastic increase in |φ| compared to the values shown in Fig. 3. Furthermore, the inclusion of kicks from massive particles beyond the Standard Model will also increase Σ at high temperatures, and consequently, T s and |φ|. Therefore, the |φ| values shown in Fig. 3 should generally be considered lower bounds.
As previously mentioned, chameleons with 1.82 < β < 3.07 can only surf the second, third, or fourth kicks if the earlier kicks leave φ ≫ φ min . The numerical solution to the chameleon equation of motion for β = 1.83 confirms that the chameleon surfs the last kick; it rolls 3.5M Pl toward φ min during the first three kicks, and then its velocity reaches a value of dφ/d ln a * = −M Pl /β = −0.54M Pl near the peak of the last kick. The chameleon maintains that velocity until the surfing solution is no longer valid (φ ∼ < φ min ). If β = 1.99, just below the threshold for surfing the third kick, then the chameleon rolls 13.5M Pl toward φ min before it begins to surf the last kick. Chameleons with 2.00 < β < 2.73 will surf the third kick, and they roll between 1.5M Pl (for β = 2.0) and 25.5M Pl (for β = 2.73) during the first two kicks. Finally, chameleons with 2.74 ≤ β ≤ 3.06 will surf the second kick. The displacement of the chameleon due to the first kick depends on the electroweak phase transition; the higher the temperature at which the top quark becomes massive, the larger the displacement. For all 2.74 ≤ β ≤ 3.06, however, the chameleon will roll more than M Pl during the first kick, even if Σ = 0 at temperatures greater than 150 GeV. In summary, chameleons with 1.83 ≤ β < 3.07 can only surf if the initial value of the chameleon significantly exceeds M Pl . For smaller values of φ i , the chameleon will reach φ ≃ φ min before it can surf. Figure 4 shows the evolution of the chameleon field for several values of β. If β < 1.82, the chameleon experiences four rolling episodes, corresponding to the four combo-kicks produced by the Standard-Model particles, and then it stops rolling when T J ∼ < 10 −5 GeV and Σ ≃ 0. In Fig. 4, we purposefully chose large values for ϕ i so that ϕ ≫ ϕ min during all four kicks. In this case, the total displacement in the chameleon field produced by Σ(T J ) does not depend on ϕ i . We refer to this displacement as ∆ϕ, and it is a function of β alone.
Earlier treatments of the chameleon's response to the kicks [61] obtained an analytic estimate of ∆ϕ by neglecting the Jordan-frame temperature's dependence on ϕ. If we assume that e β(ϕi−ϕ) ≃ 1 in Eq. (13), then Eq. (11) may be integrated twice to obtain ϕ(p): (28) where the subscript "1" indicates that this is the firstorder solution, derived assuming that β(ϕ i − ϕ) ≃ 0. In deriving this expression, we assumed that ϕ ≫ ϕ min during all four kicks so that we can neglect V ′ (φ) in Eq. (11). If we also neglect changes in the number of relativistic degrees of freedom and take T J = T J,i e −p , we can analytically evaluate ∆ϕ. First consider a single kick, with Σ(T J ) given by Eq. (12) for one species with mass m i and g i degrees of freedom. If the chameleon is initially at rest, where it is best to evaluate g * at T J = m i /2.45 for fermions and T J = m i /2.30 for bosons because τ ± reaches its maximum at these temperatures. Remarkably, this integral can be evaluated analytically: where the top number applies to fermions, and the bottom number applies to bosons. Since ϕ 1 (p) is linearly dependent on Σ(T J ), we can obtain the total field displacement by summing over all the contributions from the Standard Model. We find that ∆ϕ = −1.45β. Meanwhile, numerically evaluating Eq. (28) gives ∆ϕ 1 = −1.57β; the difference arises because the evaluation of ∆ϕ 1 fully accounts for the changes in the number of relativistic degrees of freedom when evaluating Σ(T J ) and . Since the cooling of the Jordan-frame radiation is slowed by the energy injected by annihilating particles (T J (p) > T J,i e −p after the first kick), each kick lasts a little longer (in terms of the Einstein clock p) than it does if changes in g * are neglected. The extra duration of the kicks in the Einstein frame leads to a slightly larger displacement of the chameleon field.
A similar effect implies that ∆ϕ 1 will not accurately describe the chameleon's displacement if β∆ϕ ∼ > 1. In this case, exp[β(ϕ i − ϕ)] is significantly greater than one during the kicks. Consequently, T J (p) > T J,i e −p after the first kick [see Eq. (13)], and Σ is nonzero for a larger range of p values. This dilation of Σ is illustrated in Fig. 5, which shows Σ(p) using both T J [φ(p), p] and T J (p) = T J (φ = φ i , p) for β = 1.5. The longer duration of the kicks in the Einstein frame increases ∆ϕ; for example, if β = 1.5, ∆ϕ 1 = −2.355, but Fig. 4 shows that ∆ϕ = −3.5. A more accurate analytical estimate for ∆ϕ can be obtained by iterating the solutions to Eq. (28): (32) Figure 6 shows ∆ϕ 1 , ∆ϕ 2 , ∆ϕ 3 , and ∆ϕ 4 for different values of β and compares them to the full numerical solution for ∆ϕ. We see that the analytical approximations always underestimate ∆ϕ, but the first-order approximation ∆ϕ 1 is accurate to within 10% for β ∼ < 0.9, and ∆ϕ 4 is accurate to within 10% for β ∼ < 1.7. As β approaches 1.83, |∆ϕ| increases rapidly and higher-order analytical approximations are required. The surfing solution, which corresponds to ∆ϕ = −∞, is the extension of this pattern.
We have shown that all chameleons with β ≥ 1.83 or ϕ i ≤ ∆ϕ will reach ϕ ∼ < ϕ min during the kicks. There is an additional constraint, however: the chameleon must satisfy ϕ < 0.1/β prior to the last kick to ensure that Einstein-frame particle masses do not vary by more than 10% between now and BBN [61]. The last kick is the most powerful kick; it alone displaces the chameleon by at least |∆ϕ 1 | = 0.56β. Therefore, if β > 0.43, then requiring that ϕ < 0.1/β at the onset of BBN implies that the last kick will take the chameleon to ϕ < ϕ min . For smaller values of β, the chameleon can avoid reaching φ ∼ < φ min while satisfying the BBN constraint only if  corresponds to f (R) gravity, then this condition implies 0.23 ∼ < (ϕ i − xβ) ∼ < 0.245, which is a very finely tuned initial condition! We conclude that the chameleon can only avoid being kicked to ϕ ∼ < ϕ min if β < 0.43 and ϕ i is within the limited range given by Eq. (33).
Having established that the kicks almost always take the chameleon to φ min , we now consider the chameleon's velocity when it reaches the minimum of its effective potential: for ρ * R ≫ V (φ). Since φ min ≪ M Pl , ρ * R at impact is nearly equal to the Jordan-frame radiation density: ρ ≡ (π 2 /30)g * (T J )T 4 J . During the kicks, ϕ ′ never exceeds 0.55, soφ ∼ T 2 J ϕ ′ , where T J is evaluated when φ = φ min . Figure 7 shows how the chameleon's velocity when φ ≪ M Pl depends on β and the chameleon's initial value φ i . In most cases, chameleons with larger initial values will reach φ min with smaller velocities because the Jordan frame will have longer to cool prior to impact. Surfing chameleons are an exception to this rule because T J ∼ > T s when they reach φ min , regardless of the initial field value. For a given value of ϕ i , increasing β always increases the impact velocity; for nonsurfing chameleons, increasing β generally increases both |ϕ ′ | and T J at impact, while for surfing chameleons, increasing β increases the surfing temperature T s (see Fig. 3), which more than compensates for the reduction in |ϕ ′ |. The key result of this section is that |φ| ≫ M 2 when φ ≃ φ min : at impact, T J ∼ > 0.5 MeV and ϕ ′ is nearly always > 0.02, so |φ| > 5 × 10 −9 GeV 2 in all but a few finely tuned cases. Moreover, |φ| is usually much larger; for instance, a surfing chameleon with β ≥ 3.07 has |φ| > 4000 GeV 2 when it reaches the minimum of its effective potential.

V. PARTICLE PRODUCTION
In the previous section, we derived the velocity imparted to the chameleon by the kicks, and we showed that |φ| ≫ M 2 when the chameleon reaches the minimum of its effective potential (φ ≃ φ min ). While φ ∼ > φ min , we can neglect the chameleon's bare potential V (φ), but V (φ) becomes important when the chameleon rolls to smaller values. The chameleon will climb up its bare potential until it exhausts its kinetic energy, and then it will roll back toward the minimum of its effective potential. During this rebound, V ′′ (φ) changes rapidly, and the chameleon condensate cannot adjust its mass adiabatically. Instead, we will show that the rapid changes in mass excite high-energy perturbations in the chameleon field.

A. A first look at the rebound
Since the rebound occurs when φ ∼ < φ min , we must first determine the value of φ min during the kicks. From Eq. (11), we see that the value of φ min in the early Universe is determined by the value of ρ * R Σ. Specifically, we can make the approximation ρ * R ≃ρ because φ min is always much smaller than M Pl . After the electroweak phase transition,ρΣ monotonically decreases as the Universe cools, so φ min moves to larger values as the kicks progress. At the peak of the final kick, T J = 0.16 MeV andρΣ = 32 g cm −3 , which exceeds the mean density of the Earth. Since the Earth must satisfy the thin-shell condition discussed in Section II A, the chameleon's rebound after it is kicked past φ min will sample the steep part of the chameleon's potential. For the exponential potential given by Eq. (9) and inverse-power-law poten- 16 MeV. For example, if n = 2 and β = 2, then φ min increases from 0.14M to 0.62M between the peaks of the first and the last kicks if the potential is exponential, and it increases from (5.3 × 10 −9 )M to 0.26M if the potential is an inverse power law. Since the chameleon potential diverges as φ approaches zero, the change in φ during the rebound must be less than M . Later in our analysis, we will see that quantum particle production starts while φ min < φ ∼ < M , so we will broaden our definition of the rebound to include all times while φ ∼ < M . We can estimate the duration of the rebound as ∆t ∼ M/φ, whereφ is the velocity imparted by the kicks.
In the previous section we showed thatφ ∼ T 2 J ; it follows that the duration of the rebound is much smaller than the Hubble time: H * ∆t ∼ M/M Pl . Consequently, the expansion of the Universe will not affect the evolution of the chameleon during the rebound, and ρ * R Σ will not change significantly during the rebound. During the rebound, the matter coupling induces a change in the chameleon's velocity ∆φ ∼ −(β/M Pl )ρ * R Σ∆t. Since ρ * R ∼ T 4 J and ∆t ∼ M/T 2 J , the fractional change in the chameleon's velocity is minuscule: ∆φ/φ ∼ M/M Pl . A more detailed analysis of the chameleon's equation of motion while φ > φ min for both surfing and nonsurfing chameleons confirms these estimates; while the chameleon moves a distance ∆φ ∼ M , neither Hubble friction nor the chameleon's coupling to matter changes its velocity significantly (∆φ/φ ∼ M/M Pl in all cases). Therefore, we will neglect both Hubble friction and the chameleon's coupling to matter while analyzing the rebound. We assume that the chameleon starts at φ ∼ M with a velocitẏ φ M , which is determined by the chameleon's value prior to the kicks and the value of β, as shown in Fig. 7.
In Appendix C, we show that a plane-wave perturbation in the chameleon field with a comoving wavenumber k has an effective mass where τ is conformal time,φ is the spatially averaged value of the chameleon field, and we have dropped the subscript * on the scale factor a because we will work exclusively in the Einstein frame throughout this section. We will also no longer use the variable p = ln(a * /a i ), and primes will denote differentiation with respect to the function's argument. We show in Appendix C that perturbations in the chameleon field are excited when ω ′ k (τ )/ω 2 k ∼ > 1. While the Universe is radiation dominated, a ′′ (τ ) = 0, and We recall that V ′ eff (φ) ≡ V ′ (φ) + (β/M Pl )ρ * R Σ(T J ), and the contribution from the bare potential dominates when φ ≪ φ min , by definition. The matter coupling contributes to higher derivatives of the effective potential through the dependence of T J on φ: from Eq. (13), we see that dT J /dφ ≃ −(β/M Pl )T J . Since differentiation of the bare potential introduces a factor of M −1 while differentiation of the matter coupling term introduces a factor of M −1 Pl , the relative importance of the matter coupling to V ′′ eff (φ) is suppressed by a factor of M/M Pl compared to its relative contribution to V ′ eff (φ). Therefore, the bare potential dominates V ′′ eff (φ) even if φ ∼ > φ min ; for our fiducial exponential potential with n = 2 and M = 10 −3 eV, the bare potential dominates V ′′ eff (φ) for φ ∼ < M and T J ∼ < 1 TeV. The matter coupling's contribution to V ′′′ eff (φ) is suppressed by an additional factor of M/M Pl , so our fiducial bare potential dominates V ′′′ eff (φ) for φ ∼ < 10 6 M and T J ∼ < 1 TeV.
Since the bare potential dominates both terms in Eq. (37), we can neglect the matter coupling when evaluating the relative importance of these terms. For φ ∼ M , HV ′′ (φ) term is negligible. It follows that We will see that the physical wavenumbers (k phys = k/a) of the perturbations that are excited when the chameleon rebounds are much larger than V ′′ eff (φ) while φ ∼ > M , which implies that V ′′ eff (φ) only significantly contributes to the denominator when throughout the rebound, the chameleon's coupling to matter has no impact on the adiabaticity condition and cannot affect particle production. Furthermore, since we have already shown that the matter coupling has a negligible effect onφ while the chameleon rolls ∆φ ∼ < M , we can neglect the chameleon's coupling to matter entirely while analyzing particle production during the rebound. Figure 8 illustrates how the chameleon field rebounds off its bare potential, with V (φ) given by Eq. (9) with n = 2 and M = 10 −3 eV. The initial condition iṡ φ M = −0.1 GeV 2 when φ = 2M , and we have neglected both Hubble friction and the chameleon's coupling to matter when solving Eq. (11). As expected, the chameleon rolls toward zero until V (φ) =φ 2 M /2 and then it quickly turns around and rolls back with the same speed. During this rebound, the effective mass of the field [m 2 eff = V ′′ (φ)] changes dramatically over a time ∆t ≃ 0.005M/|φ M |, momentarily reaching values greater than 10 14 GeV. Figure 8 also shows how the adiabatic ratio evolves during the rebound for several values of the physical wavenumber k phys ; we see that ω ′ k (τ )/ω 2 k ∼ > 1 for k phys ∼ < 400|φ M |/M . Therefore, we expect particle production for modes with k phys ∼ < (∆t) −1 , where ∆t is the timescale over which V ′′ (φ) changes significantly.
When we formulate an analytic model for the rebound in Section V B, we will derive an expression for ∆t. For now, we simply note that ∆t < M/|φ M |, so we expect modes with k phys ∼ > |φ M |/M to be excited during the rebound. From Eq. (C12), we see that the energy density in perturbations per logarithmic interval in k is where n k is the mode occupation number defined in Eq. (C15). In Appendix C, we show how ω ′ k (τ )/ω 2 k ∼ > 1 implies that n k ∼ 1. Therefore, we expect that the rebound will excite modes with k phys ∼ > |φ M |/M with ρ k ∼ (φ M /M ) 4 . If we compare this energy density to the initial energy in the chameleon field, ρ i ≃φ 2 M /2, we find that ρ k /ρ i ∼φ 2 M /M 4 . In the previous section, we found that the kicks impart a velocity to the chameleon that greatly exceeds M 2 , so this naive calculation yields ρ k /ρ i ≫ 1. Therefore, the modes with ω ′ k (τ )/ω 2 k ∼ > 1 are too energetic to have n k ∼ 1 given the energy available to the chameleon field. If these modes are excited during the rebound, this heuristic treatment implies that they must have n k ≪ 1, and even then, we expect them to absorb a significant fraction of the chameleon's energy. Consequently, we must include the backreaction of the perturbations when analyzing the chameleon's evolution. In the next section, we consider this backreaction in detail, and we show that it significantly alters the chameleon's trajectory during the rebound.

B. Analytical model for the rebound
To analyze the excitation of perturbations during the chameleon field's rebound off its bare potential, we first write the chameleon field as whereφ(τ ) is the spatial average of the field. We insert this expression into the chameleon's equation of motion and Taylor expand V ′ (φ) aroundφ to obtain As discussed in the previous section, the chameleon's coupling to matter is irrelevant during the rebound. Therefore, we neglect the matter coupling in this equation, but we note that it could be reinserted by replacing V (φ) with V eff (φ). To find the equation of motion forφ, we take the spatial average of this equation. Since δφ = δφ 3 = 0, we are left witḧ The δφ 2 term in this equation represents the first-order backreaction of the perturbations on the spatially averaged field. We will show that the inclusion of this term is sufficient to ensure that energy is conserved during the rebound. Since our primary aim is to understand how the transfer of energy to the perturbations affects the evolution of the spatially averaged field, we neglect the higher-order terms in Eq. (42). We note, however, that these terms represent higher-order corrections to the chameleon's evolution and may not be negligible, especially if the first-order backreaction term becomes large compared to V ′ (φ).
To evaluate δφ 2 , we expand δφ(τ, x) in terms of creation and annihilation operators as shown in Eq. (C5), and then we take the vacuum expectation value of δφ 2 . As discussed in Appendix C, we regularize the resulting expression by subtracting terms associated with the vacuum state [76], which leaves In Appendix C, we show how the mode functions φ k (τ ) may be expressed in terms of Bogoliubov coefficients α k and β k . Inserting Eq. (C9) into the expression for δφ 2 gives We can use this expression to show how the firstorder backreaction term ensures conservation of energy. The regularized energy density in perturbations ρ fluct is given by Eq. (C14). If we neglect the expansion of the Universe and take a to be constant, then Since n k = |β k | 2 , we can use Eq. (C10b) to evaluateṅ k : In the limit that a is constant, Eq. (36) implies thaṫ ω k = a 2 V ′′′ (φ)φ/(2ω k ). Inserting both of these expressions into Eq. (45) and comparing with Eq. (44) yields The energy density in the spatially averaged field is ρ =φ 2 /2 + V (φ), which implies that where the last line follows from Eq. (42) with H = 0. We see that the sum ρ fluct +ρ is constant on time scales that are short compared to H −1 . Therefore, when energy is transferred to perturbations, the first-order backreaction term in Eq. (42) ensures that an equal amount of energy is extracted from the spatially averaged field.
To probe the backreaction of the perturbations on the spatially averaged field further, we return to Eq. (44), which expresses δφ 2 in terms of the Bogoliubov coefficients α k and β k . In the previous section, we showed that the perturbation modes that we expect to be excited during the rebound (k ∼ |φ M |/M ) must have n k ≪ 1. Since n k = |β k | 2 , the normalization condition for the Bogoliubov coefficients (|α k | 2 − |β k | 2 = 1) demands that |α k | ≃ 1 when n k ≪ 1. In this regime of perturbative particle production, we can obtain an approximate solution for β k by taking α k = 1 in Eq. (C10b) [77]: where we have chosen τ = 0 to correspond to some time before particle production begins. Inserting this expression into Eq. (44) and taking α k = 1 yields Given that |β k | ≪ 1, we expect the |β k | 2 term to be much smaller than the α k β * k term that generates the second term in the integrand of Eq. (50). One may be concerned, however, that this second term involves the integral of an oscillating function and may therefore be suppressed relative to the |β k | 2 term. However, our approximate solution for β k (τ ) implies that so |β k | 2 also contains a cosine integral. Therefore, we may safely assume that the |β k | 2 term makes a negligible contribution to δφ 2 .
Since we are only interested in time scales that are much shorter than the Hubble time, we can further simplify this expression for δφ 2 by taking the scale factor a to be constant and using Eq. (37) to evaluate ω ′ k . With these simplifications, we obtain where we have defined ω 2 k,phys ≡ (ω k /a) 2 = k 2 phys +V ′′ (φ). For the rest of this section, we will neglect the expansion of the Universe and deal only with physical wavenumbers evaluated at the time of the chameleon's rebound. To make the expressions less cluttered, we will omit the "phys" subscript on both k and ω k in subsequent equations. For consistency, we will also drop the Hubble friction term from Eq. (42); recall from SectionV A that this friction term has a negligible impact on the evolution of φ during the rebound.
We now have a new equation for the evolution of the spatially averagedφ field that includes the first-order backreaction:φ where D(t) ≡ (1/2)V ′′′ [φ(t)] δφ 2 with δφ 2 given by Eq. (52) is the "dissipation" term that expresses how energy is transferred from the spatially averaged field to the perturbations. The dissipation term D(t) is non-Markovian; it depends on the entire history of the chameleon's evolution up to t and therefore has "memory." The non-Markovian nature of the dissipation term was also highlighted in an earlier analysis of particle production in scalar field theory [78], which used the "in-in" formalism to calculate dissipation from particle production in φ 4 theory. If ω k is assumed to be constant, then our D(t) is identical to their dissipation term (see Eq. (27) in Ref. [78]), which demonstrates how our perturbative Bogoliubov technique can simplify particle production calculations. Ref. [78] showed that there is no Markovian limit for the dispersion term in φ 4 theory. However, the fact that V ′′′ (φ) sharply increases as φ decreases for chameleon potentials implies that, prior to the rebound, the non-Markovian integral in Eq. (52) is dominated by values of t ′ that are just slightly smaller than t, and we will show that this feature allows us to derive a local approximation for D(t).
The dominance of the t ′ ∼ < t portion of the time integral in Eq. (52) also ensures that the cosine integral in this expression is positive, which allows us to extract the qualitative behavior of D(t). Asφ rolls to smaller values prior to the rebound,φ is nearly constant (as shown in Fig. 8), and D(t) is proportional toφ with a positive coefficient. At this stage, D(t) acts like a drag term, and it slows the chameleon down. However, unlike a drag term, D(t) is actually an integral overφ, which implies that it does not vanish whenφ = 0. Instead, D(t) will remain negative during and shortly after the rebound. Consequently, D(t) acts like a new potential term during the rebound: we will show that D ≃ V ′ D (φ) for some V D (φ).
To derive this new "dissipative potential," we must further simplify our expression for D(t). First, we assume that the modes that are excited have k 2 ≫ V ′′ (φ) throughout the rebound, 2 which allows us to set ω k (t) ≃ k in Eq. (52). Second, we impose an infrared (IR) cutoff on the cosine integral in Eq. (52) and consider only k ≥ k IR . This IR cut-off makes the separation ofφ and δφ explicit; modes with k < k IR are absorbed intō φ, while modes with k ≥ k IR are considered perturbations. Clearly, we must choose k IR to be smaller than the wavenumbers of the modes that we expect to be excited during the rebound. However, we do not want to make k IR arbitrarily small because we only solve a linearized equation for δφ [Eq. (C4)], while the equation of motion forφ is nonlinear. Therefore, decreasing k IR implies that we are neglecting more nonlinear effects. When we calculate particle production numerically in the next section, we will see that taking different values of k IR can be used to determine the importance of nonlinear field interactions. For now though, we simply assume that k IR ∼ < (∆t) −1 , where ∆t ∼ < M/|φ M | is the duration of particle production during the rebound.
With these simplifications, we find that where Ci(x) ≡ − ∞ x (cos y)/y dy. We now take advantage of the fact that V ′′′ (φ) sharply increases as φ decreases, which implies that the integral over t ′ will be dominated by a limited range of values that are just slightly smaller than t. Moreover, Ci(x) is divergent for small x, which further enhances the contribution to the integral from small values of (t − t ′ ). Numerical evaluation of Eq. (54) using theφ(t) solution found numerically in the next section confirms that restricting (t − t ′ ) ∼ < 0.1M/|φ M | does not significantly change D(t) near the rebound. Therefore, we can assume that k IR (t − t ′ ) < 1 and use the approximation where γ E ≃ 0.577 is Euler's constant, to obtain where t min = t − 0.1M/|φ M |. We then integrate Eq. (56) by parts, and we make the approximation again taking advantage of the fact that contributions from t ′ ≃ t dominate the integral, to obtain Since V ′′ (φ) increases sharply as φ decreases, V ′′ φ (t) ≫ V ′′ φ (t min ) as the chameleon climbs its bare potential and turns around. Also, during the rebound, ln [2k IR (t − t min )] changes only slightly, so we may approximate it as constant. We then find that for some constant κ. Comparing the numerical evaluation of D(t) as given by Eq. (54) to the approximation given by Eq. (59) indicates that κ increases from ∼0.02 to ∼ 0.05 during the rebound for both exponential and power-law potentials. Therefore, we expect that a value of κ in this range will accurately approximate D(t) during the rebound; we will see in the next section that this is indeed the case. The approximate expression for D(t) during the rebound given by Eq. (59) shows that the backreaction of the perturbations on the evolution ofφ effectively adds a new term to the chameleon potential: with 0.02 < κ < 0.05. For both exponential and powerlaw potentials, V ′ D (φ) > V ′ (φ) while φ < M , so the dissipative potential will dominate the field's evolution during the rebound. This dominance of the first-order backreaction term is concerning, for it indicates that higher-order contributions to the backreaction are probably not negligible and signals a breakdown of perturbation theory. Our primary aim in this Section, however, is to understand how the extraction of energy from the spatially averaged field affects particle production during the rebound. Since the first-order backreaction captures this energy transfer, we can use our analysis of the first-order backreaction to gain insight into how the evolution of φ is affected by particle production.
Since the dissipative potential is dominant when φ < M , it determines the minimum value ofφ during the rebound, which we denoteφ ta : V D (φ ta ) ≡φ 2 M /2. It follows from Eq. (60) that the maximum value of the chameleon's effective mass during the rebound is much smaller than the wave numbers of the excited modes (k ∼ > |φ M |/M ): V ′′ (φ ta ) = |φ M |/ √ κ ≪φ 2 M /M 2 given that |φ M | ≫ M 2 . Therefore, the backreaction prevents m eff from reaching the extremely large values seen in Fig. 8, and our approximation that ω k ≃ k in Eq. (52) is justified. Furthermore, the field will turn around before the adiabatic ratio ω ′ k (τ )/ω 2 k exceeds unity, which will keep n k ≪ 1 for the excited modes.
Next, we consider how the dissipative potential affects the duration of the rebound. Thus far, we have used ∆t ∼ < M/|φ M | to estimate the duration of the rebound. However, Fig. 8 shows that this definition of ∆t overestimates the duration of the change in m eff , which is the timescale that determines which perturbation modes will be excited. For the parameters shown in Fig. 8, M/|φ M | = 10 −11 GeV −1 , but m eff changes significantly and |ω ′ k |/ω 2 k exceeds unity for a much shorter time: ∼ 5 × 10 −14 GeV −1 . Therefore, we need to refine our calculation of ∆t. We also need to incorporate our new understanding of the evolution of the spatially averaged field. Whileφ ∼φ ta during the rebound, the equation of motion forφ, including the first-order backreaction, is approximatelÿ where we have defined m 2 D ≡ V ′′ D (φ ta ). If we set t = 0 when the field turns around (φ =φ ta ), the solution to Eq. (61) is We can now calculate how long it takes for m 2 eff to change significantly: where the last line follows from Eq. (62). (We are interested in the change in m 2 eff , as opposed to m 2 D , because m 2 eff governs the behavior of the perturbations.) To account for the change in m 2 eff asφ approachesφ ta and as φ rolls away fromφ ta , we multiply Eq. (63) by 2 when evaluating ∆t. Our estimated value of the wave number of the most energetic excited mode (k ex ) is then whereφ ta satisfies V D (φ ta ) ≡φ 2 M /2. If the chameleon has an exponential bare potential, then evaluating Eq. (64) gives This potential is more general than those we have considered previously, because it allows for two mass scales.
To avoid fine-tuning, we will assume that M s /M v is of order unity and we will continue to assume that M v ∼ 10 −3 eV to maintain a connection to dark energy. Numerically solving V D (φ ta ) ≡φ 2 M /2 for 2 ≤ n ≤ 10 and 10 −6 < |φ M |/GeV 2 < 10 6 reveals that (φ ta /M s ) n ∼ < 0.05 in all cases, so we may neglect the O φn ta /M n s term in Eq. (66). We then use V D (φ ta ) ≡φ 2 M /2 to obtain Since V D (φ ta ) ≡φ 2 M /2 with an exponential potential is a transcendental equation, we cannot obtain an exact algebraic expression forφ ta /M s . For a limited range of |φ M | values, however, we can approximateφ ta as where c n is an constant of order unity that depends on the range of |φ M | values under consideration. Inserting this expression into Eq. (67) yields where b n = (c n /2 1/n ) n+1 is also an order-unity constant. 3 For 10 −6 < |φ M |/GeV 2 < 10 6 , b n ranges from 0.25 for n = 2 to 0.38 for n = 10, and Eq. (69) differs from Eq. (67) by less than 7%. Equation (69) shows that k ex is rather insensitive to both κ and n; varying κ within the range 0.01 ≤ κ ≤ 0.06 has nearly no impact on k ex , and increasing n from 2 to 10 changes k ex by less than 25% for 10 −6 < |φ M |/GeV 2 < 10 6 .
If the chameleon has a power-law bare potential, then evaluating Eq. (64) gives For this potential, we can algebraically solve V D (φ ta ) ≡φ 2 M /2 forφ ta to obtain 3 Equation (69) differs from the equation for kex given in our earlier work [60] by a factor of 1/ √ 2 because we initially used the time required for |V ′′ (φ) − V ′′ (φta)|/V ′′ (φta) = 0.5 as the starting point of our derivation of kex. We later realized that |V ′′ (φ) − V ′′ (φta)|/V ′′ (φta) = 1 gave a better fit to the numerical results obtained Section V C, so we modified our expression for kex.
Fortunately, Eq. (72) indicates that k ex is still relatively insensitive to order-unity changes in κ. For powerlaw potentials, however, k ex depends very strongly on n, with larger n values giving smaller values for k ex . Furthermore, comparing Eq. (72) to Eq. (69) reveals that exponential potentials have smaller k ex values than power-law potentials. In both cases, our earlier estimate, k ex ≃ |φ M |/M , significantly underestimates k ex by missing factors related to (M s /φ ta ), which is much greater than unity. Steeper potentials generally have smaller k ex values becauseφ turns around at a larger value.
To summarize the key results of this section, we evaluated the first-order backreaction of the perturbations on the spatially averaged fieldφ and found that the dynamics ofφ during the rebound are governed by a new "dissipative" potential given by Eq. (60). This new potential forces the chameleon field to turn around much earlier than the solution without backreaction predicted, which prevents ω ′ k /ω 2 k from exceeding unity. Consequently, the occupation numbers of the excited modes remain very small. We then used the dissipative potential to estimate the duration of the rebound ∆t. Assuming that the rebound will excite perturbation modes with k ∼ < (∆t) −1 , we derived predictions for the wave number of the most energetic excited mode (k ex ) for both exponential and power-law potentials. We found that k ex ≫ |φ M |/M in both cases, which indicates that the backreaction does not prevent the transfer of energy to extremely energetic modes. Therefore, we now have two reasons to expect that chameleon gravity will suffer a computational breakdown during the rebound; the rebound will excite modes that lie far beyond the expected limits of effective field theory, and if the theory can be trusted during the rebound, the backreaction of these modes will also significantly alter the evolution of the spatially averaged field.

C. Numerical computation of particle production
We test our analytical analysis of the rebound by numerically solving the linearized perturbation equations and the spatially averaged equation with first-order backreaction. In our numerical analysis, we neglect the expansion of the Universe and take the scale factor to be constant. We define "physical" creation and annihilation operatorsâ k,phys = a 3/2â k that obey the commutation relation and express δφ(t, r = a x) in terms of these operators: +â † k,phys φ * k,phys (t)e −ik phys ·r .
Comparing Eq. (74) to Eq. (C5) reveals that φ k,phys = √ aφ k . If we take a to be constant, then Eq. (C7) implies φ k,phys + ω 2 k,phys φ k,phys = 0, where ω 2 k,phys ≡ k 2 phys + V ′′ [φ(t)]. We solve this equation for φ k,phys for several logarithmically spaced k phys values with k IR < k phys < k max . As in the previous section, k IR separates the long-wavelength perturbations that are included inφ from the shorter-wavelength perturbations that compose δφ, and it is chosen to be smaller than the wavenumbers of the modes we expect to be excited during the rebound (k IR ≪ k ex ). We also do not expect modes with k phys ≫ k ex to be excited, so we choose a value of k max that is much larger than k ex . The number of k values we sample depends on the ratio k max /k IR and is chosen so that the interval between log k values is ∼0.05.
To evaluate ω 2 k,phys in Eq. (75), we have to solve Eq. (53) forφ(t): We evaluate δφ 2 at each time step by converting Eq. (44) into physical variables and then using the φ k,phys solutions to compute the integral (77) With the φ k,phys solutions, we can also evaluate the occupation number n k = |β k | 2 : which corresponds to Eq. (C15) if a is constant. We also use Eq. (39) to evaluate ρ k for each φ k,phys solution. When we numerically solve Eqs. (75) and (76), we must choose initial conditions forφ(t) and all the φ k,phys (t) functions. We initially setφ = 2M witḣ φ =φ M , whereφ M is chosen from the range of velocities shown in Figure 7. Since most of the particle production occurs nearφ ta ∼ < 0.1M , the final spectrum of perturbations is insensitive to the initial value ofφ, provided that it is greater than 0.5M . We neglect the tiny amount of particle production that occurs beforeφ = 2M and initially set n k = 0. If n k = 0 at some initial time t i , Eqs. (C9) and (C10) imply that We do not solve directly for φ k,phys (t) in our numerical analysis because it is impossible to accurately evaluate Eq. (78) when n k ≪ 1. Instead, we express φ k,phys (t) as FIG. 9: The evolution of the spatially averaged chameleon field, including the backreaction from particle production, as it rebounds off its bare potential for the same parameters as in Figure 8. The top panel shows that the field value turns around at a much larger value than predicted by the classical solution, and the middle panel shows that the effective mass is confined to much smaller values. The bottom panel shows that the earlier turn-around ensures that the adiabatic ratio is always much less than unity.
Equation (75) with an initial condition given by Eq. (79) requires thatθ k (t) = W k (t). Equation (75) also provides an evolution equation for W k (t), and demanding that n k = 0 at the initial time requires W k (t i ) = ω k,phys (t i ) andẆ k (t i ) = 0. We then define a new function ̟(t) ≡ W k (t) − ω k,phys (t) that describes the deviation of φ k,phys (t) from Eq. (79), and we numerically solve Eq. (75) for the evolution of ̟(t). When the φ k,phys terms in Eq. (78) are expressed as functions of ̟, the result is 1/2+f (̟,̟). Therefore, we can evaluate n k directly from ̟(t) without having to numerically subtract the vacuum contribution (the 1/2 term in Eq. 78), thus avoiding the numerical errors introduced by subtracting two numbers with values that far exceed their difference. Figure 9 shows the numerically computed evolution of the spatially averaged chameleon field after starting at φ = 2M with a velocityφ M = −0.1 GeV 2 . The potential V (φ) is given by Eq. (9) with n = 2 and M = 10 −3 eV. Comparing this figure to the classical solution shown in Fig. 8 for the same V (φ) and initial conditions reveals how profoundly the evolution ofφ is affected by the transfer of energy to perturbations. As predicted in the previous section, the field turns around when V D (φ) =φ 2 M /2, which gives a much larger value forφ ta than the classical V (φ ta ) =φ 2 M /2. Consequently,φ turns around before its effective mass exceeds a GeV and before the adiabatic ratio for k ≃ (∆t) −1 exceeds unity, in stark contrast to the classical evolution depicted in Fig. 8. Since the adiabatic ratio is always very small, we expect n k ≪ 1 as well. Figure 10 shows the evolution of n k for the same wavenumbers; we see that n k increases dramatically as the field rebounds, and then it maintains a constant value as the field rolls out to larger values. We also see that n k ≪ 1; as expected, the excited modes are so energetic that n k ≪ 1 is required to conserve energy.
Even though n k ≪ 1, the fluctuations still contain a significant fraction of the chameleon's energy. Figure 9 illustrates that the rebound is not elastic;φ rolls out with a smaller velocity because some energy has been transferred to the fluctuations. Figure 11 shows the evolution of ρ fluct , as defined by Eq. (C14), forφ M = −100 GeV 2 . We see that all of the chameleon's energy is transferred to fluctuations at the rebound (at t = 0), but then some of that energy is returned to the spatially averaged field; all values ofφ M share this basic behavior. The postrebound transfer of energy from the fluctuations toφ is a manifestation of the dissipative potential V D (φ); since the backreaction of the perturbations onφ acts as a new potential immediately after the rebound, it can accelerateφ as it rolls out to larger values. The amount of energy returned toφ after the rebound depends on k IR ; Figure 11 shows that smaller values of k IR lead to less final energy in fluctuations. Since, k IR determines which modes are treated linearly, this dependence on k IR indicates that the final value of ρ fluct depends on nonlinear interactions that are not included in our analysis. Therefore, we cannot determine how much energy is transferred to perturbations during the rebound. This limitation is disappointing, but not surprising; as we discussed in the previous section, the dominance of the dissipative potential V D (φ) over V (φ) during the rebound foretold that our linear analysis with only a first-order backreaction would be insufficient to determine the chameleon's final state. However, the fact that increasing k IR (and thus including more nonlinear effects) increases the final value of ρ fluct indicates that nonlinear interactions are unlikely to prevent the transfer of energy to fluctuations.
Although the total energy transferred to fluctuations depends on k IR , the energy spectrum of the fluctuations is more robust. Figure 12 shows the post-rebound fluctuation energy density per logarithmic interval in k, as defined in Eq. (39), for bothφ M = −100 GeV 2 anḋ φ M = −20000 GeV 2 . For both cases, the spectra are shown for several values of k IR . As expected, the rebound generates a spectrum of fluctuations that is rather sharply peaked at a specific wavelength, with larger |φ M | values exciting more energetic fluctuations. Figure 12 shows that the amplitude of the fluctuation spectrum depends on the value of k IR , but the shape of the spectrum does not. We conclude that the basic characteristics of the fluctuation spectrum, particularly which wave numbers receive the most energy, does not depend on nonlinear effects. As discussed in Section V B, the duration of the rebound determines which fluctuation modes are excited, and we see no evidence that nonlinear effects change the rebound's basic timescale. On the contrary, Fig. 9 illustrates that even the first-order backreaction does not significantly alter the duration of the rebound. Furthermore, the analytic calculation of k ex in Section V B successfully predicts the peak in the fluctuation spectrum; Eq. (67) with κ = 0.03 gives k ex = 2.4 × 10 16 GeV forφ M = −100 GeV 2 and k ex = 5.6 × 10 18 GeV forφ M = −20000 GeV 2 . The turn-around value ofφ is also independent of k IR ; for all the spectra shown in Fig. 12, V D (φ ta ) ≃φ 2 M /2. We conclude that nonlinear effects only become important after the rebound, when the generated fluctuations begin to interact. Since both φ ta and k ex are determined by the dynamics ofφ prior to the rebound, they are insensitive to k IR .
Having established that the fluctuation spectra are insensitive to nonlinear effects, we begin a more extensive comparison of the analytic model developed in Section V B to the numerical results. Figure 13 shows how the fluctuation spectra depend on bothφ M and V (φ). The top panel shows spectra for an exponential potential with n = 2, while the bottom panel shows spectra for both exponential and power-law potentials with various values of n. In all cases, the fluctuation spectrum is sharply peaked at a wave number we call k peak , and k peak increases asφ M increases. We also see that k peak does not vary much with n for exponential potentials, but k peak is very sensitive to n for power-law potentials, with smaller n values producing larger k peak values. Finally, we see that power-law potentials produce larger k peak values than exponential potentials. All of these findings are consistent with the predictions of Section V B, and we see that the rebound does indeed excite fluctuations with extremely high energies.
We quantitatively test our analytic model for the rebound in Fig. 14, which shows the turn-around value of φ for exponential potentials with n = 2, 4 and 10 as a function ofφ M . The solid lines show V D (φ ta ) =φ 2 M /2; the values of κ were chosen to match the numerical results, and all lie within the range predicted in Section V B (0.02 ∼ < κ ∼ < 0.05). The dotted curve directly beneath each solid line shows V (φ ta ) =φ 2 M /2. We see that particle production does forceφ to turn around before . In all cases with exponential potentials, kIR = 10|φM |/M , which implies that 0.03kex ∼ < kIR ∼ < 0.06kex, and for the power-law spectra kIR ≃ 0.05kex. V (φ ta ) =φ 2 M /2, and the value ofφ ta can be predicted by the dissipative potential (up to the choice of κ). Figure  15 shows k peak for the same potentials and compares it to k ex given by Eq. (67); in all cases k ex is within 10% of k peak . This level of agreement between k ex and k peak is rather shocking; the derivation of k ex was based on the rough estimate that the rebound would excite modes with k ∼ (∆t) −1 , where ∆t was the duration of the rebound, and we expected it to differ from k peak by a factor of order unity. Instead, we see that Eq. (67) and Eq. (69) accurately predict k peak , so we can use these expressions to predict the fluctuation spectra generated by other exponential potentials. We see from Eq. (67) that the value of n has a mixed impact on k ex ; k ex is proportional to n, but increasing n also slightly increasesφ ta , which decreases k ex . As a result, n does not significantly affect k peak , as seen in Fig. 15. For n ≤ 10, changing n changes k ex by less than 25% for 10 −6 GeV 2 < |φ M | < 10 6 GeV 2 . Figure 16 shows that the analytic model developed in Section V B is equally successful for power-law potentials. We see that V D (φ ta ) =φ 2 M /2 with κ = 0.025 successfully predicts the turn-around value ofφ, and k ex given by Eq. (71) matches k peak to within 16%. Power-law potentials are more numerically challenging than exponential potentials because of the small values ofφ ta . Consequently, we can only numerically study potentials with n ≥ 6, and we are restricted to a smaller range ofφ M values. However, the agreement between k ex and k peak shown in Fig. 16 indicates that Eq. (71) provides an accurate estimate of the wave numbers of perturbations excited by rebounds off more general power-law potentials. In particular, potentials with smaller values of n give significantly larger values for k ex ; for example, if n = 2, then k ex exceeds the reduced Planck mass for |φ M | ∼ > 2 GeV 2 . In summary, numerically solving Eqs. (75) and (76) for both power-law and exponential potentials confirms the predictions of Section V B; the evolution ofφ during the rebound is governed by the dissipative potential V D (φ), and the rebound excites extremely energetic fluctuations, with shallower potentials generating higher-energy fluctuations. When combined with the values ofφ obtained in Section IV, both of these facts indicate that chameleon gravity experiences a catastrophic breakdown of calculability just prior to BBN.

VI. SUMMARY AND DISCUSSION
Chameleon gravity runs into trouble in the early Universe because it attempts to unite two radically different energy scales: the MeV-GeV scale of the Standard Model and the meV scale of dark energy. The same coupling to the trace of the stress-energy tensor that enables the chameleon to evade astronomical and labora-tory constraints on fifth forces also makes the chameleon susceptible to excitation whenever the trace of the stressenergy tensor changes, as occurs when particles become non-relativistic in the early Universe. Meanwhile, the chameleon's ability to evade detection is also contingent on the presence of a different energy scale in the chameleon's potential: V (φ/M ) where M ∼ < 0.01 eV. The possibility of using the chameleon's potential to drive the current epoch of cosmic acceleration leads us to consider even smaller values for M : M ≃ 0.001 eV. We have shown that the combination of these disparate energy scales leads to a breakdown of calculability just prior to BBN. The chameleon's coupling to matter accelerates the field to GeV-scale velocities, causing its effective mass to change rapidly and leading to the production of extremely energetic fluctuations that violate the limits of Effective Field Theory. Moreover, the production of these fluctuations significantly alters the chameleon's evolution, leaving its state during BBN unknown.
The chameleon's difficulties begin when the temperature of the radiation bath in the early Universe falls below the mass of a particle that is in thermal equilibrium. At that time, the trace of the stress-energy tensor momentarily increases because the pressure of the massive particles decreases faster than their density. During this transition, the trace of the stress-energy tensor is sufficiently large that the chameleon's coupling to it overcomes Hubble friction and forces the chameleon field to roll down the slope of its effective potential. Earlier treatments of chameleon cosmology used these kicks to the chameleon field to move the chameleon field from its expected initial value (M ≪ φ i ∼ < M Pl ) to the minimum of its effective potential (φ min ≃ M ) prior to BBN [61,64].
We have shown that the chameleon field does not just reach the minimum of its effective potential; the field rolls past it with a very large velocity. Unlike earlier work, our analysis includes the effect the chameleon's evolution has on the expansion of the Universe in the Jordan frame. As the chameleon rolls, it slows the Jordan-frame expansion, which extends the duration of the kicks and enhances their impact. If the chameleon coupling to matter is slightly stronger than gravitational (β ∼ > 1.8), this effect halts the expansion of the Jordan frame entirely until the chameleon reaches the minimum of its effective potential. We call this novel solution to the chameleon's equation of motion the "surfing solution," and it guarantees that all chameleon fields with β ∼ > 2 will reach the minima of their effective potentials with velocities |φ| ∼ > 10 −3 GeV 2 , regardless of their initial value. In general, nearly all chameleons with sub-Planckian initial values reach their potential minima with |φ| ∼ > MeV 2 ; only weakly coupled chameleons (β ∼ < 0.4) with finely tuned initial conditions avoid this fate. Moreover, our calculation of |φ| uses a minimal prescription for the kicks that only includes contributions from Standard-Model particles. The inclusion of the QCD trace anomaly [74,75], interactions during the QCD phase transition [72], or ad-  67) for an exponential potential (Eq. 9) with M = 10 −3 eV and n = 2. Increasing n does not significantly alter the fluctuation spectrum. In the white region, φi is sufficiently large that the field does not reach the minimum of its effective potential. In the region marked "BBN excluded," φ > 0.1M Pl /β when the temperature is 1 MeV, which spoils the success of BBN.
ditional particles would give the chameleon field an even larger velocity when it reaches the minimum of its effective potential.
When the chameleon rolls past the minimum of its effective potential with |φ| ≫ M 2 , it climbs up the steep portion of its bare potential, where φ ∼ < M . This steep region is required to give the chameleon a large mass in dense environments, thus screening the fifth force it generates. After the kicks, however, the steepness of V (φ) becomes a severe liability. Since V ′′ (φ) changes sharply for small displacements (∆φ < M ), the chameleon's velocity after the kicks causes its effective mass to change rapidly. These non-adiabatic changes in mass trigger the excitation of fluctuations with wave numbers k ≃ (∆t) −1 , where ∆t < M/|φ| is the timescale over which the mass changes. Since |φ| ≫ M 2 , these fluctuations are extremely energetic, with k ≫ 10 7 GeV in most cases. The classical evolution of the chameleon field as it rolls up and rebounds off the steep part of its potential predicts that these fluctuations should be generated with occupation numbers of order unity, but there is insufficient energy in the chameleon field to generate such high-energy particles. Therefore, the backreaction of the fluctuations must significantly alter the evolution of the chameleon field.
To account for the transfer of energy to the fluctuations, we added a first-order backreaction term to the chameleon's equation of motion that ensures that energy is conserved during the rebound. An analysis of this backreaction term revealed that it effectively introduces a new term to the chameleon's potential that dominates over V (φ) during the rebound. This "dissipative potential" does not significantly change the timescale of the rebound, but it halts the chameleon's climb up its potential while the occupation numbers of the excited modes are still very small. Nevertheless, numerical calculations confirmed that a large fraction of the chameleon's energy is transferred to fluctuations with k ≃ (∆t) −1 during the rebound, which significantly alters the trajectory of the spatially averaged field.
Furthermore, these numerical calculations revealed that a simple derivation of ∆t using the dissipative potential accurately predicts the fluctuation wave numbers that receive the most energy for a given value ofφ after the kicks. These wave numbers are shown in Figure  17 for an exponential chameleon potential. Power-law potentials generate fluctuations with even higher wave numbers. Thus we see that nearly all chameleons generate extremely high-energy perturbations during the rebound, with strongly coupled chameleons generating trans-Planckian fluctuations. Since we cannot trust chameleon gravity up to such high energies, we cannot predict what impact these fluctuations may have on BBN.
Our calculation of the excitation of perturbations during the rebound is not a complete treatment. First, we only included the first-order backreaction of the perturbations on the evolution of the spatially averaged field. This backreaction term is a one-loop correction to the chameleon's equation of motion, and our omission of higher-order backreaction terms implies that we are neglecting higher-order quantum corrections. During the rebound, the first-order backreaction term dominates over the bare potential, which indicates that the higherorder backreaction terms are not insignificant. Since the "dissipative potential" generated by the first-order backreaction is responsible for the transfer of energy from the perturbations to the spatially averaged field after the rebound, higher-order backreaction terms would probably affect the final distribution of energy. Nevertheless, the first-order backreaction term alone provides insight into how the transfer of energy to the perturbations affects the rebound. Since the inclusion of the backreaction term does not change the timescale of the rebound, it does not affect which modes are excited. The backreaction strongly affects the evolution of the spatially averaged field, however, and it suppresses the occupation numbers of the excited modes.
Second, we linearized the equation of motion for the perturbation modes, which neglects interactions between modes with different wavelengths. We explored the importance of these nonlinear interactions by varying the IR cutoff of our perturbations, which shifts the distinction between the perturbations and the spatially averaged field. Since the equation of motion for the spatially averaged field is not linearized, the IR cutoff determines which nonlinear effects are neglected. Changing the IR cutoff strongly affects the final state of the chameleon field because it varies how much energy is transferred to the perturbations. Therefore, a full nonlinear treatment of the perturbation equations and their backreaction is required to determine how the chameleon field evolves after the rebound. However, changing the IR cutoff does not alter the position of the peak in the perturbation spectrum. The timescale of the rebound determines which modes are excited, and this timescale is not affected by nonlinear interactions. Therefore, the utility of a fully nonlinear treatment of the perturbations would be limited by our ignorance of physics at energy scales much greater than 1000 TeV.
Another source of uncertainty in our calculations is the chameleon potential V (φ). The large velocity imparted to the chameleon field by the kicks implies that the chameleon reaches very small values during the rebound; even when its journey is curtailed by the transfer of energy to perturbations, the field reaches values that are smaller than the values obtained in the densest astrophysical objects. There are reasons to distrust the chameleon potential at these small values of φ: the one-loop Coleman-Weinberg correction to V (φ) is much larger than V (φ) itself, and the dimensionless quartic coupling d 4 V /dφ 4 exceeds unity. These facts alone support our conclusion that chameleon gravity suffers a computational breakdown after the kicks, independently of the severe violations of adiabaticity that trigger particle production.
These flaws are commonplace in chameleon gravity, however; even for moderate densities (ρ ≃ 10 g cm −3 ), d 4 V /dφ 4 > 1 for nearly all values of n and β, and the oneloop quantum corrections dominate the potential for several chameleon models [57]. Most analyses of chameleon gravity ignore these difficulties, implicitly assuming that V (φ) is protected from quantum effects in some way. We have extended chameleon gravity the same privilege by assuming that V (φ) is valid for all values of φ during the rebound, and we have demonstrated that a computational breakdown still occurs prior to BBN. Yet, it may be possible to avoid this calamity by changing V (φ) for small φ values. It would be interesting to search for wellbehaved potentials that enable the chameleon screening mechanism while avoiding particle production after the kicks. Any such potential would have to be more complicated than the exponential and power-law potentials that we considered, and we leave this investigation for future work.
Another interesting avenue for further investigation is to determine if any other modified gravity theories suffer from a similar computational breakdown in the early Universe. Coupled dark energy theories with quintessence fields that couple to dark matter but not baryons [79,80] may also be susceptible; if the dark matter particle is a thermal relic, then the quintessence field will be kicked when the dark matter particle becomes nonrelativistic. Since these coupled dark energy theories do not need the chameleon mechanism to evade laboratory and Solar System constraints on fifth forces, the scalar's potential function is less constrained than the potential functions we considered. Nevertheless, the fact that the quintessence field is supposed to drive cosmic acceleration implies that its potential function involves at least one mass scale that is much smaller than the mass of the dark matter particle. Furthermore, fifth forces exclusive to the dark sector can be constrained by tidal streams [81,82] and anisotropies in the cosmic microwave background [83], which require β ∼ < 0.07 in the absence of a screening mechanism. Therefore, these theories can only include gravitational-strength couplings if they employ the chameleon mechanism, which would require a chameleon scalar potential.
In general, chameleon gravity provides a cautionary tale about the dangers of uniting dark energy and highenergy physics in a single theory; the extreme hierarchy of energy scales can produce surprising and uncontrollable effects. The key ingredients that make chameleon gravity vulnerable are 1) a field that couples to the stress-energy tensor and is initially displaced from the minimum of its potential, and 2) a potential function for which V ′′ (φ) changes significantly when the scalar field rolls by a small amount (∆φ ≪ GeV). Given these two features, we expect that the field will be accelerated to a GeV-scale velocity in the early Universe, and then that large velocity will induce rapid changes in the field's mass, leading to the excitation of high-energy particles. An easy way to prevent a violation of adiabaticity may be to impose a shift symmetry; if V (φ) is insensitive to the value of φ, it does not matter how largeφ becomes. Therefore, the perils faced by chameleon gravity in the early Universe provide additional motivation for including a shift symmetry in any scalar-tensor theory that attempts explain the current epoch of cosmic acceleration. thermal equilibrium are where g is the number of degrees of freedom for the particle species; m is the particle's mass; T is the temperature of the radiation bath; and the + sign in the denominator applies to fermions, while the − sign applies to bosons. We define ρ R to be the sum of the the energy densities for all particles that are in thermal equilibrium with the radiation bath in the early Universe (including neutrinos), and we define g * ≡ ρ R [(π 2 /30)T 4 ] −1 . To compute the kick function Σ ≡ (ρ R − 3p)/ρ R , we first evaluate ρ − 3p for each particle species: where we have introduced u = E/T as the integration variable in the last line. Dividing this expression by ρ R = g * (π 2 /30)T 4 yields the contribution to Σ from a single particle species, as given by Eq. (12): (A5) In Section IV, we evaluate Σ(T ) for the Standard Model by summing the contributions from the particle species listed in Table I. Prior calculations of the kick function assumed that g * was constant during each kick and computed g * for each particle species by summing the contributions of all particles with m ≤ m i : where T i is the temperature of the particle species. At temperatures greater than 1 MeV, T i = T for all particles included in ρ R , but at lower temperatures, the neutrinos decouple from the photon bath, and T ν = T . This computation of g * overestimates g * during the kick because it treats the particle responsible for the kick as if it were relativistic throughout the kick. To avoid underestimating Σ in this way, we compute g * (T ) by numerically evaluating ρ(T ) for all the particles in Table I and adding these energy densities to the energy densities of the relativistic particles. Prior to the QCD phase transition, the relativistic bosons are gluons (g = 16) and photons (g = 2), and the relativistic fermions are light quarks (g = 36 for u, d, s), muons, electrons, and neutrinos (g = 6). Therefore, g * (T ) smoothly decreases from 106.75 to 61.75 prior to the QCD phase transition. We treat the QCD phase Contributions to Σ fermions bosons particle g m (GeV) particle g m (GeV) before QCD phase transition top  12  172  Higgs 1  125  bottom 12  4.2  Z  3  91  charm 12 1.3 W ± 6 80 tau 4 1.8 after QCD phase transition muon 4 0.106 π 0 1 0.140 electron 4 5.11 × 10 −4 π ± 2 0.135 transition as an instantaneous event that occurs at a temperature of 170 MeV. Below this temperature, the quarks and the gluons are bound into hadrons, and the only particles that contribute to g * are pions, muons, electrons, neutrinos and photons. Consequently, g * discontinuously changes from 61.75 to 17.25 when T = 170 MeV, which generates the discontinuity in Σ(T ) seen in Figure 2. After the QCD phase transition, g * (T ) decreases smoothly from 17.25 to 10.75 as the temperatures decreases from 170 MeV to 10 MeV.
The calculation of Σ(T ) for T ∼ < 1 MeV is complicated by the decoupling of the neutrinos from the radiation bath. After they decouple, the neutrinos' temperature is proportional to 1/a, which causes it to differ from the temperature of the radiation bath during the final kick. To evaluate the neutrino temperature T ν at temperatures below 1 MeV, we employ the conservation of entropy to obtain T ν = g * S (T ) 10.75 where g * S is the total entropy density of all particles in thermal equilibrium with the radiation bath divided by (2π 2 /45)T 3 . We evaluate g * S (T ) by computing the entropy density, s = (ρ + P )/T , from Eqs. (A1) and (A2) in the same way that we used Eq. (A1) to evaluate g * (T ). In addition to providing the neutrino temperature, this computation of g * S (T ) is used to numerically solve Eq. (13) for T J (a * ). Once we know how T ν /T changes as the temperature cools below 1 MeV, we can use Eq. (A6) to evaluate the neutrinos' contribution to g * (T ) during the final kick, during which g * (T ) decreases from 10.75 to 3.36. This evolution of g * significantly enhances the amplitude of the final kick compared to earlier calculations that assumed a fixed value of g * = 10.75.
The energy density of the fluctuations is defined as Like the zero-point energy of an infinite tower of harmonic oscillators, this integral is divergent. We can expose this divergence by substituting Eq. (C9) for φ k and using Eqs. (C10) and (C11), which gives Thus we see that |β k | 2 = 0 corresponds to the lowest possible energy density, and the divergent portion of the integral corresponds to the infinite energy of this ground state. Therefore, the |β k | 2 = 0 state is interpreted as the vacuum state, and we regularize ρ fluct by subtracting the energy density of the vacuum state [84], which gives where we have introduced the occupation number n k : We see from Eq. (C13) that n k (τ ) = |β k (τ )| 2 .
If we start from a vacuum state, then Eq. (C10b) implies that we will remain in a vacuum state while the adiabaticity condition, is satisfied. Conversely, particle production occurs when the effective mass of the perturbations varies sufficiently rapidly that ω ′ k (τ )/ω 2 k ≃ 1. In that case, Eq. (C10b) states that the vacuum state will evolve to a state with n k = |β k | 2 = 0.