# A novel integrated Sachs-Wolfe effect from Early Dark Energy

Tristan L. Smith  Department of Physics and Astronomy, Swarthmore College,
Swarthmore, PA 19081, USA  Center for Cosmology and Particle Physics,
Department of Physics, New York University, New York, NY 10003, USA  John T.
Giblin Jr  Department of Physics, Kenyon College, 201 N College Rd, Gambier,
OH 43022, USA  Department of Physics/CERCA/Institute for the Science of
Origins, Case Western Reserve University, Cleveland, OH 44106-7079 – USA
Mustafa A. Amin  Department of Physics and Astronomy, Rice University,
Houston, TX 77005, USA  Mary Gerhardinger  Department of Physics, Kenyon
College, 201 N College Rd, Gambier, OH 43022, USA  Center for Particle
Cosmology, Department of Physics and Astronomy, University of Pennsylvania,
Philadelphia, Pennsylvania 19104, USA  Ericka Florio  Department of Physics,
Kenyon College, 201 N College Rd, Gambier, OH 43022, USA  Department of
Mathematics and Theoretical Physics, University of Cambridge, Wilberforce Rd,
Cambridge UK, CB3 0WA  Matthew Cerep  Department of Physics and Astronomy,
Swarthmore College, Swarthmore, PA 19081, USA  Department of Physics and
Astronomy, West Virginia University, P.O. Box 6315, Morgantown, WV 26506, USA
Center for Gravitational Waves and Cosmology, West Virginia
University,Chestnut Ridge Research Building, Morgantown, WV 26505, USA  Shar
Daniels  Department of Physics and Astronomy, Swarthmore College, Swarthmore,
PA 19081, USA  Department of Physics and Astronomy, University of Delaware,
Newark, Delaware 19716, USA

###### Abstract

We study the nonlinear effects of minimally coupled, massless, cosmological
scalar fields on the cosmic microwave background (CMB). These fields can
exhibit post-recombination parametric resonance and subsequent nonlinear
evolution leading to novel contributions to the gravitational potential. We
compute the resulting contributions to the CMB temperature anisotropies
through the time-variation of the gravitational potential (i.e., the
integrated Sachs-Wolfe (ISW) effect). We find that fields that constitute 5%
of the total energy density and become dynamical at $z_{c}\simeq 10^{4}$ can
produce marginally observable ISW signals at multipoles $\ell\simeq 2000$.
Fields that become dynamical at earlier times and/or have initial
displacements at a flatter part of their potential, produce ISW contributions
that are significantly larger and at higher multipoles. We calculate these
dynamics and the resulting evolution of gravitational perturbations using
analytic estimates alongside detailed nonlinear lattice simulations, which
couple scalar fields and cosmological fluids to a perturbed metric. Finally,
we discuss the possibility of detecting these features with future high-
resolution CMB observations.

##  I Introduction

Scalar fields are ubiquitous in cosmology, from the inflaton Guth (1981) that
phenomenologically gives a mechanism for the flatness and structure of our
observed universe, to quintessence Caldwell  _et al._ (1998) which provides a
dynamical explanation for the current epoch of accelerated expansion. In
between these two epochs, there have been a plethora of proposed cosmological
scalar fields, with a range of masses, motivated by many beyond-the-standard
model scenarios from inflationary model building Baumann and McAllister
(2015), to string theory Svrcek and Witten (2006); Arvanitaki  _et al._
(2010); Kane  _et al._ (2015); Stott  _et al._ (2017), to attempts to explain
the seemingly fine-tuned nature of the current epoch of accelerated expansion
Griest (2002); Linder and Smith (2011); Kamionkowski  _et al._ (2014). Indeed,
adding scalar fields is often the first line of attack to confront
inconsistencies in our cosmological scenarios. After that, however, it is
imperative that these scenarios be forced to predict how their existence will
impact our cosmological observations.

In this paper we consider the consequences of a single cosmological scalar
field, which is minimally coupled, with a canonical kinetic term, and a
potential that goes as $V=\lambda\varphi^{4}$ about its minimum. With these
restrictions, the scalar field will be held fixed at some initial displacement
by Hubble friction. When the Hubble parameter drops below some critical value,
the field will become dynamical and eventually oscillate about the minimum of
the potential. Such a field is often referred to as ‘Early Dark Energy’ (EDE).

There is a rich literature exploring the linear and nonlinear dynamics in
cosmological scalar fields both in the very early and late universe. In the
context of the post-recombination universe, a linearized analysis of the
resonant growth of scalar field fluctuations, and its impact gravitational
potentials was explored earlier in Amin  _et al._ (2012). In that case,
however, the scalar field was not always subdominant in energy density. A
linearized analysis of resonant growth of perturbations in certain oscillating
dark energy models was investigated in Johnson and Kamionkowski (2008). In
terms of numerical simulations, the resonant growth of perturbations, and
nonlinear dynamics in scalar fields for $\lambda_{n}\varphi^{n}$ potentials
has been studied in detail in the early universe context Lozanov and Amin
(2017, 2018); Khlebnikov and Tkachev (1996).111In general, there is a large
literature on how parametric resonance and subsequent nonlinear dynamics may
also play a central role in thermalizing the energy stored in the inflaton at
the end of inflation, initiating the epoch of radiation domination. See Amin
_et al._ (2014); Lozanov (2019) for recent reviews.

Indirect constraints on cosmological scalar fields which are dynamically
relevant between inflation and today can be placed using observations of the
cosmic microwave background (CMB) and clustering of matter (through galaxy
clustering and/or weak lensing). Given that the field is minimally coupled,
its effects on observables come from how it modifies the expansion history and
how its perturbations gravitate. Changes to the homogeneous expansion leads to
changes to Hubble friction, modifying the growth of dark matter perturbations.
Given that the field oscillates about an approximately quartic minimum, non-
resonant modes have significant pressure support which prevents them from
growing. This, in turn, limits their impact on cosmological observables. For
fields that become dynamical during radiation domination non-resonant modes,
at the linear level, act like a perfect, $w=-1/3$, fluid, with a maximum
fractional energy density of
$\rho_{\varphi}/(\rho_{\varphi}+\rho_{\gamma}+\rho_{\nu})\lesssim 5\%$ Poulin
_et al._ (2019), where $\rho_{\gamma}$ is the energy density of the CMB
photons and $\rho_{\nu}$ is the energy density of three standard ultra-
relativistic neutrino species.

The self-interaction from the $\lambda\varphi^{4}$ potential allows the
fluctuations in such fields to grow significantly through parametric
resonance, whereby linear perturbations grow exponentially by efficiently
extracting energy from the oscillating homogeneous background field. Ref.
Smith  _et al._ (2020) noted that such a process occurs in fields which may
resolve the Hubble tension. The underlying reason is that for such potentials,
the resonant wave numbers do not flow out of the resonance band as the
universe expands.222As discussed in Ref. Smith  _et al._ (2020), this implies
that the typical potential used to resolve the Hubble tension, which has
$V\propto\varphi^{6}$ around its minimum, does not experience significant
effects through parametric resonance. The growth rate of field fluctuations
eventually exceeds the Hubble rate (even if the field is subdominant in energy
density), leading to nonlinear, spatially inhomogeneous evolution of the
scalar field.

In the post-recombination $\Lambda$CDM universe, and during matter domination,
the gravitational potential $\Psi$ is determined primarily by the dark matter
perturbations, with $\dot{\Psi}\simeq 0$. The significant resonant growth of
scalar field perturbations, even if its background energy density is
subdominant, can lead to a significant contribution to $\dot{\Psi}$. This
evolution can leave an imprint on the CMB via the Integrated Sachs-Wolfe
effect (ISW). It is this effect that we explore in detail in this paper.333For
fields which become dynamical well before matter-radiation equality, they may
make a non-negligible contribution to gravitational potential at the resonant
scale. The fact that the scalar field energy density is subdominant limits its
impact.

Through both detailed simulations and analytic scaling equations, we show that
this novel ISW contribution is larger the earlier the field becomes nonlinear.
The time at which nonlinear evolution starts, in turn, is earlier if the field
either becomes dynamical earlier and/or starts at a location where the
potential is flatter. As a result of this we find that scalar fields which may
address the Hubble tension become nonlinear too late to produce an observable
ISW. For fields that become dynamical earlier/start on a flatter part of their
potential, the resulting ISW can dominate over the primary CMB power spectrum,
with a peak at very small angular scales which cannot be probed by current CMB
measurements, but which may be accessible to proposed CMB telescopes such as
CMB-HD Sehgal  _et al._ (2019); Aiola  _et al._ (2022).

The paper is organized as follows. In Section II, we describe our fiducial
model and present some details of our linear and nonlinear analyses,
highlighting its important features. We also describe the metric perturbations
and fluid descriptions used in our nonlinear analysis to model a realistic
universe. In Section III we describe our numerical methods. Section IV, we
present the results of our linear analysis, showing the evolution of the
homogeneous mode of EDE, as well as a comparison between our linear and
nonlinear results which demonstrate the presence of parametric resonance and
validate the code used for our nonlinear analysis. We then outline our
calculations of the ISW from the nonlinear simulation. In Section V we present
the resulting ISW contributions for several different scalar field models. In
Section VI, we summarize our results and discuss their significance with
regards to falsifying the EDE model. Note that we work in units where
$c=\hbar=k_{B}=1$ and where the reduced Planck mass is $m_{\rm pl}\equiv(8\pi
G)^{-1/2}$.

##  II The Dynamics of A Cosmological Scalar Field

The action for a minimally coupled scalar field, $\varphi$, with a canonical
kinetic term, is given by

$\displaystyle S=\int d^{4}x\sqrt{-g}\left[\frac{m_{\rm
pl}^{2}}{2}R-\frac{1}{2}\partial_{\mu}\varphi\partial^{\mu}\varphi-V(\varphi)+\mathcal{L}_{{\Lambda\textrm{CDM}}}\right],$
(1)

where $\mathcal{L}_{{\Lambda\textrm{CDM}}}$ includes contributions from the
cosmological constant, cold-dark matter and the rest of the Standard Model. We
work in the Conformal Newtonian gauge; a perturbed, conformal time foliation
of FLRW without vector and tensor perturbations,

$\displaystyle ds^{2}=a^{2}[-(1+2\Phi)d\tau^{2}+(1-2\Psi)d\mathbf{x}\cdot
d\mathbf{x}],$ (2)

where $a$ is the scale factor. Throughout our work here, we assume that there
is no significant gravitational anisotropic stress, so that the two Bardeen
potentials are the same, $\Psi=\Phi$. The gravitational potential $\Psi$ is
evaluated using Einstein’s equations linearized in $\Psi$ (but not the field
$\varphi$).444The consistency of these assumptions can be tested using full
nonlinear gravitational methods Giblin  _et al._ (2019); however, we do not
expect any significant corrections to our results when we relax these
assumptions Giblin and Tishue (2019). With this metric, the equation of motion
for this field is

$\displaystyle\ddot{\varphi}-\left(1+4\Psi\right)\nabla^{2}\varphi+2\left(\mathcal{H}-2\dot{\Psi}\right)\dot{\varphi}$
(3) $\displaystyle+a^{2}\left(1+2\Psi\right)\frac{dV}{d\varphi}$
$\displaystyle=0$

where an over-dot is a partial derivative with respect to conformal time,
$\mathcal{H}\equiv(da/d\tau)/a$, and $\mathcal{H}$ includes contributions from
both the spatially averaged scalar field, as well as the
${\Lambda\textrm{CDM}}$ components. That is,

$\displaystyle\mathcal{H}^{2}=\frac{a^{2}}{3m_{\rm
pl}^{2}}\left[\bar{\rho}_{\varphi}+\bar{\rho}_{{\Lambda\textrm{CDM}}}\right],$
(4)

where

$\bar{\rho}_{\Lambda{\rm CDM}}=3H_{0}^{2}m_{\rm
pl}^{2}\left(1-\Omega_{m}+\frac{\Omega_{m}}{a^{3}}+\frac{\Omega_{r}}{a^{4}}\right),$
(5)

$\Omega_{m}$ and $\Omega_{r}$ is the total matter (cold dark matter plus
baryons) and total radiation (photons plus neutrinos) density parameters and a
bar indicates the quantity is spatially averaged. The scalar field stress
energy tensor is given by

$(T_{\varphi})^{\mu}_{\
\nu}=\nabla^{\mu}\varphi\nabla_{\nu}\varphi-\frac{1}{2}\left[\nabla^{\alpha}\varphi\nabla_{\alpha}\varphi+2V(\varphi)\right]\delta^{\mu}_{\
\nu}.$ (6)

For concreteness, we take the potential associated with the scalar field to be

$\displaystyle
V(\varphi)=m^{2}f^{2}\left(1-\cos\frac{\varphi}{f}\right)^{\\!2}\,,$ (7)

which, in the limit where $\varphi/f\ll 1$, resembles a massless scalar
$V(\varphi)\simeq(m^{2}/f^{2})\varphi^{4}/4$. Our choice of the precise form
of $V(\varphi)$ is for convenience, and to make contact with earlier work
Smith  _et al._ (2020). The scaling equations derived in this Section lead us
to conclude that any potential which has a quartic minimum, and flattens to a
constant at large field values, will have a qualitatively similar
phenomenology.

Note that the effective mass near the minimum of the potential $m^{2}_{\rm
eff}\equiv d^{2}V/d\varphi^{2}\simeq 3m^{2}(\varphi/f)^{2}$ which is different
from the mass parameter $m^{2}$. In the following discussion equations using
the ‘$\sim$’ symbol are missing factors of order unity whereas those with
‘$\simeq$’ retain those factors but are still approximate.

###  II.1 Background evolution

It will prove useful to develop a set of approximate equations which summarize
the background evolution of the scalar field. Taking the homogeneous limit of
Eq. (3) we have

$\ddot{\bar{\varphi}}+2\mathcal{H}\dot{\bar{\varphi}}+a^{2}V_{,\varphi}=0,$
(8)

where we have introduced the notation $V_{,\varphi}\equiv dV/d\varphi$.
Roughly, when the conformal Hubble parameter is large enough, the friction
term dominates and we have

$\dot{\bar{\varphi}}\sim-\frac{a^{2}}{\mathcal{H}}V_{,\varphi}.$ (9)

During this ‘slow-roll’ evolution the equation of state parameter associated
with the field is given by

$w_{\varphi}\simeq-1+\frac{\dot{\bar{\varphi}}^{2}}{2a^{2}V}\sim-1+\frac{a^{2}}{\mathcal{H}^{2}(a)}\frac{V_{,\varphi}^{2}}{V}.$
(10)

Since $w_{\varphi}>-1$, once the field becomes dynamical the critical scale
factor, $a_{c}$, is roughly given by

$\frac{a_{c}^{2}}{\mathcal{H}^{2}(a_{c})}\frac{V_{,\varphi}^{2}}{V}\sim 1.$
(11)

If we write $V(\bar{\varphi})=m^{2}f^{2}E(\theta\equiv\bar{\varphi}/f)$, then
this condition can be written terms of the mass parameter:

$m\sim\frac{\mathcal{H}(a_{c})}{a_{c}}\frac{\sqrt{E}}{|E_{,\theta}|}\bigg{|}_{\theta=\theta_{i}},$
(12)

where we evaluate the potential-dependent term at the initial field
displacement, $\bar{\varphi}_{i}=\theta_{i}f$.

The fractional contribution the field makes to the total energy density when
it becomes dynamical can now be written

$\frac{\bar{\rho}_{\varphi}(a_{c})}{\bar{\rho}_{\rm
tot}(a_{c})}\simeq\frac{a_{c}^{2}V(\bar{\varphi}_{i})}{m_{\rm
pl}^{2}\mathcal{H}^{2}(a_{c})}\sim\frac{f^{2}}{m_{\rm
pl}^{2}}\left(\frac{E}{E_{,\theta}}\right)^{2}\bigg{|}_{\theta=\theta_{i}}.$
(13)

After the field becomes dynamical it quickly starts to oscillate with a cycle-
averaged equation of state parameter $\langle w_{\varphi}\rangle=1/3$ with an
envelope that decreases as Poulin  _et al._ (2019)

$\varphi_{\rm env}(a)=\varphi_{c}\frac{a_{c}}{a}.$ (14)

Note that empirically we find that $\varphi_{c}\simeq 0.7\varphi_{i}$.

###  II.2 Linear perturbations

If the field were a spectator during inflation then it would generically have
both adiabatic and isocurvature initial conditions Hlozek  _et al._ (2018);
Smith  _et al._ (2020). Here we will ignore the isocurvature
perturbations.555As discussed in Ref. Smith  _et al._ (2020), the amplitude of
the isocurvature perturbations scales with the tensor to scalar ratio, $r$.
Given current upper limits on $r$ from measurements of the $B$-mode
polarization, the EDE isocurvature contribution is, in general, too small to
impact current measurements. The superhorizon adiabatic initial conditions
during slow-roll can be found analytically, giving Smith  _et al._ (2020)

$\frac{\delta\varphi(a\leq
a_{c})}{f}\sim\frac{m^{2}a^{2}}{\mathcal{H}^{2}(a)}E_{,\theta}\bigg{|}_{\theta=\theta_{i}}\zeta_{\rm
ad}(\vec{k}).$ (15)

Outside of resonant phenomena (discussed below), linear perturbations can be
characterized by the cycle-averaged non-adiabatic sound speed Poulin  _et al._
(2019)

$\langle
c_{s}^{2}\rangle=\frac{2a^{2}\varpi^{2}(a)+k^{2}}{6a^{2}\varpi^{2}(a)+k^{2}},$
(16)

where $\varpi(a)$ is the angular frequency of the background field. The
detailed evolution of $\varpi$ is unimportant, since at all times these
perturbations will have significant pressure support, leading to
$\delta\rho_{\varphi}/\bar{\rho}_{\varphi}\simeq{\rm constant}$ on subhorizon
scales.

###  II.3 Parametric resonance

The fluctuations of the scalar field are unstable and undergo exponential
growth in a narrow band of wavenumbers, as discussed in detail in Ref. Smith
_et al._ (2020) for the EDE context. We summarize the main results here.
Ignoring the effects of gravity, the Fourier modes of field perturbations at
linear order satisfy

$\displaystyle\delta\ddot{\varphi}_{k}+2\mathcal{H}\delta\dot{\varphi}_{k}+\left[k^{2}+a^{2}V^{\prime\prime}(\bar{\varphi})\right]{\delta\varphi}_{k}=0\,.$
(17)

The $a^{2}V^{\prime\prime}(\bar{\varphi})\propto a^{2}\bar{\varphi}^{2}$ term
provides a time-dependent, approximately periodic, oscillatory contribution to
the effective frequency.

Soon after the field becomes dynamical, the amplitude of the oscillations in
the background field are damped so that $\varphi/f\ll 1$ and the potential can
be approximated as $V\propto\varphi^{4}$. When the potential is well-
approximated by a power-law the perturbations evolve as

$\displaystyle\delta\varphi_{k}(a)\simeq\delta\varphi_{k}(a_{c})\frac{a_{c}}{a}\exp\left[\int_{a_{c}}^{a}\frac{\Re[\mu_{k}]}{\mathcal{H}}d\ln
a\right].$ (18)

The perturbations grow exponentially fast around the resonant wavenumber
$k_{\rm res}$, with a width $\Delta k_{\rm res}$, and a Floquet exponent
$\mu_{k_{\rm res}}$ given by Smith  _et al._ (2020)

$\displaystyle k_{\rm res}\simeq
1.27m\left(\frac{\varphi_{c}}{f}\right)a_{c}\,,$ (19)
$\displaystyle\gamma_{\rm res}\equiv\frac{\Delta k_{\rm res}}{k_{\rm
res}}\simeq 0.17\,,$ $\displaystyle\Re[\mu_{k_{\rm res}}]\simeq
0.036m\left(\frac{\varphi_{c}}{f}\right)a_{c}\,.$

The above discussion allows us to derive an approximate equation for the
scale-factor when the resonant wavenumber becomes nonlinear, $a_{\rm nl}$. In
the following we will assume that $a_{\rm nl}\ll 1$, ensuring that the
universe is filled with just matter and radiation. Using the approximate,
resonantly growing solution (18), we can estimate when the perturbations
become nonlinear:

$\displaystyle\frac{\gamma_{\rm res}k_{\rm
res}^{3}}{2\pi^{2}}P_{{\frac{\delta\varphi}{\varphi_{{\scaleto{\rm
env}{1.8pt}}}}}}(k_{\rm res},a_{\rm nl})\equiv\gamma_{\rm
res}\Delta^{2}_{{\frac{\delta\varphi}{\varphi_{{\scaleto{\rm
env}{1.8pt}}}}}}(k_{\rm res},a_{\rm nl})\simeq 1,$ (20)

where
$\langle\delta\varphi_{\bm{k}}(a)\delta\varphi_{\bm{k}^{\prime}}(a)\rangle=(2\pi)^{3}P_{\delta\varphi/\varphi_{\rm
env}}(k,a)\delta({\bm{k}}-\bm{k^{\prime}})$. Solving for $a_{\rm nl}$ we have

$\displaystyle\frac{a_{\rm nl}}{a_{c}}$ $\displaystyle\simeq$ $\displaystyle
1-\frac{14\ln[\gamma_{\rm
res}\Delta^{2}_{{\frac{\delta\varphi}{\varphi_{{\scaleto{\rm
env}{1.8pt}}}}}}(k_{\rm res},a_{c})]}{ma_{c}/\mathcal{H}(a_{c})\theta_{c}}$
(21) $\displaystyle+$ $\displaystyle\frac{a_{c}}{a_{c}+a_{\rm
eq}}\left(\frac{14\ln[\gamma_{\rm
res}\Delta^{2}_{{\frac{\delta\varphi}{\varphi_{{\scaleto{\rm
env}{1.8pt}}}}}}(k_{\rm
res},a_{c})]}{2ma_{c}/\mathcal{H}(a_{c})\theta_{c}}\right)^{2},$

where $a_{\rm eq}\equiv\Omega_{r}/\Omega_{m}$ and is the scale factor at which
the matter and radiation energy densities are equal. Since fluctuations in
$\varphi$ are still linear at $a_{c}$, both terms are positive and give
$a_{\rm nl}>a_{c}$.

Since $k_{\rm res}\sim ma_{c}\sim\mathcal{H}(a_{c})$ we can approximately
compute $\Delta_{\delta\varphi/\varphi_{\rm env}}(k_{\rm res},a_{c})$ using
the superhorizon solution given in Eq. (15) which allows us to write

$\displaystyle\Delta_{{\frac{\delta\varphi}{\varphi_{{\scaleto{\rm
env}{1.8pt}}}}}}(k_{\rm res},a_{c})$ $\displaystyle\simeq$ $\displaystyle
A_{s}^{1/2}\frac{m^{2}a_{c}^{2}}{\mathcal{H}^{2}(a_{c})}\frac{E_{,\theta}(\theta_{i})}{\theta_{c}},$
(22) $\displaystyle\sim$ $\displaystyle\frac{A_{s}^{1/2}}{p_{i}},$

where $p_{i}\equiv d\ln E/d\ln\theta|_{\theta=\theta_{i}}$ is the effective
power-law index of the potential at the initial field displacement and
$\Delta^{2}_{\zeta}(k)\simeq A_{s}$ assuming scale invariant initial
conditions.

###  II.4 The ISW contribution

The ISW can be calculated using a line of sight integral Seljak and
Zaldarriaga (1996)

$\displaystyle T_{\rm ISW}(\hat{n})$ $\displaystyle=$ $\displaystyle
2\int_{\tau_{i}}^{\tau_{0}}d\tau\frac{\partial\Psi(\tau,\vec{x})}{\partial\tau}\bigg{|}_{\vec{x}=(\tau_{0}-\tau)\hat{n}},$
(23)

where $\tau_{i}$ is some initial conformal time and $\tau_{0}$ is the
conformal time today. The angular power spectrum due to the ISW effect is

$\displaystyle C_{l}^{\rm ISW}$ $\displaystyle=$ $\displaystyle
8\pi\int_{\tau_{i}}^{\tau_{0}}d\tau\int_{\tau_{i}}^{\tau_{0}}d\tau^{\prime}\int
d\ln k\Delta^{2}_{\dot{\Psi}}(k;\tau,\tau^{\prime})$ $\displaystyle\times$
$\displaystyle j_{l}(k[\tau_{0}-\tau])j_{l}(k[\tau_{0}-\tau^{\prime}]),$

where $j_{l}(x)$ is a spherical Bessel function. This allows us to see that
since $d\tau=d\ln a/\mathcal{H}$, the ISW contribution scales as $C_{l}^{\rm
ISW}\propto\Delta^{2}_{\dot{\Psi}/\mathcal{H}}(\ell\simeq k[\tau_{0}-\tau_{\rm
nl}])$.

The scalar field contribution to the rate of change of the Newtonian potential
on subhorizon scales ($k\gg\mathcal{H}$) can be approximated by

$\nabla^{2}\dot{\Psi}_{\varphi}\,\simeq\,\frac{a^{2}}{2m_{\rm
pl}^{2}}\nabla^{j}T^{0}_{\ j}.$ (25)

where we have taken the subhorizon limit since $k_{\rm res}$ is well within
the horizon at $a_{\rm nl}$. In Fourier space, and linearizing around
$\bar{\varphi}$, we have

$\frac{\dot{\Psi}_{\varphi}}{\mathcal{H}}\simeq\frac{1}{\mathcal{H}m_{\rm
pl}^{2}}\dot{\bar{\varphi}}\delta\varphi,$ (26)

which has a power spectrum

$\displaystyle\Delta^{2}_{\dot{\Psi}_{\varphi}/\mathcal{H}}(k,a)$
$\displaystyle\simeq$ $\displaystyle\frac{\dot{\bar{\varphi}}^{2}}{m_{\rm
pl}^{4}\mathcal{H}^{2}}\varphi_{\rm
env}^{2}\Delta^{2}_{{\frac{\delta\varphi}{\varphi_{{\scaleto{\rm
env}{1.8pt}}}}}},$ $\displaystyle\simeq$
$\displaystyle\frac{\bar{\rho}_{\varphi}(a)}{\bar{\rho}_{\rm
tot}(a)}\theta_{c}^{2}\frac{f^{2}}{m_{\rm
pl}^{2}}\left(\frac{a_{c}}{a}\right)^{2}\Delta^{2}_{{\frac{\delta\varphi}{\varphi_{{\scaleto{\rm
env}{1.8pt}}}}}}(k,a).$

The maximum contribution to the ISW effect from the scalar field will be at
$k=k_{\rm res}$ and $a=a_{\rm nl}$, where
$\Delta^{2}_{{\frac{\delta\varphi}{\varphi_{{\scaleto{\rm
env}{1.8pt}}}}}}(k_{\rm res},a_{\rm res})\simeq 1$, giving

$\Delta^{2}_{\dot{\Psi}_{\varphi}/\mathcal{H}}(k_{\rm res},a_{\rm
nl})\simeq\frac{\bar{\rho}_{\varphi}(a_{\rm nl})}{\bar{\rho}_{\rm tot}(a_{\rm
nl})}\left(\frac{f}{m_{\rm pl}}\frac{a_{c}}{a_{\rm nl}}\theta_{c}\right)^{2}.$
(28)

Assuming that $a_{\rm nl}>a_{\rm eq}$,666In order to produce a measurable ISW
contribution we must have $a_{\rm nl}>a_{\rm rec}>a_{\rm eq}$ so this
assumption is required to produce an observable signal. we have
$\bar{\rho}_{\varphi}/\bar{\rho}_{\rm tot}\propto 1/a$, and using Eq. (13) we
can write this as

$\displaystyle\Delta^{2}_{\dot{\Psi}_{\varphi}/\mathcal{H}}(k_{\rm res},a_{\rm
nl})$ $\displaystyle\sim$
$\displaystyle\left(\frac{\bar{\rho}_{\varphi}(a_{c})}{\bar{\rho}_{\rm
tot}(a_{c})}\right)^{2}p_{i}^{2}$ $\displaystyle\times$
$\displaystyle\frac{{\rm max}(a_{c},a_{\rm eq})}{a_{\rm
nl}}\left(\frac{a_{c}}{a_{\rm nl}}\right)^{2}.$

Figure 1: The scale factor when the field becomes nonlinear, $a_{\rm nl}$ (top
panel), the factor that determines the amplitude of the ISW contribution [see
Eq. (28)] (middle panel), and the resonant wavenumber, all as a function of
the critical redshift, $z_{c}$, (bottom panel). The stars in the middle panel
correspond the four models that we show in Fig. 6. These locations explicitly
show that decreasing $\theta_{i}$ decreases the ISW peak; however, a larger
value of $z_{c}$ increases the strength of the signal. Note that when
computing this figure we use Eqns. (21), (28, and solve for the exact
background field dynamics to determine $f/m_{\rm pl}$ $m$, and $\theta_{c}$.

We are now in a position to anticipate what aspects of the cosmological scalar
field will control its contribution to the ISW effect. First note that the
analytic approximations, Eqns. (28) and (II.4), indicate that the relevant
dynamics depend on the shape of the potential around the initial field
displacement (i.e., Eqns. (12), (13), (22)), and on the fact that the
potential goes as $V\propto\varphi^{4}$ around its minimum. In other words,
the exact shape of the potential between the initial field displacement and
the minimum has a subdominant effect on the relevant field dynamics. Fixing
$a_{c}$ and $\bar{\rho}_{\varphi}(a_{c})/\bar{\rho}_{\rm tot}(a_{c})$, Eq.
(II.4) tells us that the amplitude of the ISW contribution has a strong
dependence the ratio $a_{c}/a_{\rm nl}$, with a smaller $a_{\rm nl}\rightarrow
a_{c}$ leading to a larger ISW contribution.

As we now argue, the overall ISW contribution increases as the power-law slope
of the potential at the initial field displacement decreases. First note that
Eq. (II.4) tells us that the ISW contribution scales as $\propto
p_{i}^{2}/a_{\rm nl}^{3}$. Next note that
$[m\theta_{c}a_{c}/\mathcal{H}(a_{c})]^{-1}\sim
p_{i}\sqrt{E(\theta_{c})/\theta_{c}^{4}}$ from Eq. (12), so that Eq. (21)
tells us $a_{\rm nl}\propto p_{i}$ or $a_{\rm nl}\propto p_{i}^{2}$, if the
first or second term dominates, respectively. Therefore, the ISW contribution
will roughly scale as $\propto 1/p_{i}^{n}$, with $1\lesssim n\lesssim 4$.

Fig. 1 shows how the different quantities that determine the amplitude of the
ISW contribution depends on both the critical redshift, $z_{c}$, and the
initial field displacement, $\theta_{i}$. This figure shows that for a fixed
value of $\bar{\rho}_{\varphi}(a_{c})/\bar{\rho}_{\rm tot}(a_{c})=0.05$, the
amplitude of the ISW is set by the quantity plotted in the middle panel which
shows that we can achieve a similar amplitude by increasing $z_{c}$ and/or
increasing $\theta_{i}$. Note that the pairs of model parameters which have a
similar amplitude in the middle panel of Fig. 1 (indicated by the stars) also
have similar ISW contributions, as shown in the right panel of Fig. 6.

##  III Linear and Nonlinear Numerical Methods

In all cases we split the ${\Lambda\textrm{CDM}}$ content into separate
contributions from matter, radiation and dark energy. Each of these
contributions will contribute to the background evolution of the scale factor,
via Eq. (5). Most simulations that study the parametric decay of scalar
condensates contain only scalar fields and self-consistently calculate the
evolution of the homogeneous spacetime. Sometimes these simulations will
either impose an expansion history, as was an option in Ref. Felder and
Tkachev (2008), or add a homogeneously diluting component that contributes to
the self-consistent evolution, e.g. Ref. Weiner  _et al._ (2020). In this work
we will additionally track the inhomogeneities of matter and radiation using a
fluid treatment. Since we will keep perturbation in the these fluids to linear
order, we will consider these perturbations in momentum space following the
method derived in Ref. Ma and Bertschinger (1995). For each fluid, $i$, with
(constant) equation of state, $w_{i}\equiv p/\rho$, the fluid perturbation is

$\displaystyle\delta_{i}\equiv\delta\rho_{i}/\bar{\rho_{i}},$ (30)

which we keep throughout the simulation as a momentum-space quantity. The
variable $\theta$ is defined by

$\left(\bar{\rho}+\bar{P}\right)\theta\equiv ik^{j}\delta T^{0}_{j},$ (31)

where the contributions to $\theta$ from each species, $\theta_{i}$, are just
the divergences of the fluid velocities,
$\theta_{i}=ik^{j}\partial_{j}{v_{i}}$, and

$\sum_{i}\left(\bar{\rho}_{i}+\bar{P}_{i}\right)\theta_{i}\equiv\left(\bar{\rho}+\bar{P}\right)\theta.$
(32)

We then have a set of evolution equations for the fluid variables, where $i$
can be either matter or radiation,

$\displaystyle\dot{\delta}_{i}$ $\displaystyle=-(1+w)(\theta_{i}-3\dot{\Psi})$
(33) $\displaystyle\dot{\theta}_{i}$
$\displaystyle=-\frac{\dot{a}}{a}(1-3w)\theta_{i}+\frac{w}{1+w}k^{2}\delta_{i}+k^{2}\Psi.$
(34)

Of course, both the field and fluid equations require us to know the metric
perturbations. These can be found from the linearized Einstein Equations,

$\displaystyle{\Psi}$
$\displaystyle=-\frac{1}{k^{2}}\left[\tilde{S}_{1}-\frac{1}{k^{2}}\frac{\dot{a}}{a}\left(\tilde{S_{2}}-\frac{3}{2m_{\rm
pl}^{2}}a^{2}\sum_{\rm species}(\bar{\rho}+\bar{P})\theta\right)\right]$ (35)
$\displaystyle\dot{{\Psi}}$
$\displaystyle=-\frac{1}{3k^{2}}\left[-\frac{1}{k^{2}}\left(\tilde{S_{2}}-\frac{3}{2m_{\rm
pl}^{2}}a^{2}\sum_{\rm
species}(\bar{\rho}+\bar{P})\theta\right)+\frac{\dot{a}}{a}{\Psi}\right]$ (36)

where the two quantities $\tilde{S}_{1}$ and $\tilde{S}_{2}$ are the Fourier
Transforms of

$\displaystyle S_{1}=\frac{1}{2m_{\rm
pl}^{2}}\delta\rho_{\varphi}=\frac{1}{2m_{\rm
pl}^{2}}\left(\frac{1}{2}\dot{\varphi}^{2}+\frac{1}{2}\left(\nabla\varphi\right)^{2}+V(\varphi)-\bar{\rho}_{\varphi}\right)$
(37)

and

$\displaystyle S_{2}=\frac{3}{2m_{\rm
pl}^{2}}\partial_{i}(\partial_{0}\varphi\partial_{i}\varphi)=\frac{3}{2m_{\rm
pl}^{2}}\left[\nabla\varphi\cdot\nabla\dot{\varphi}+\dot{\varphi}\nabla^{2}\varphi\right],$
(38)

which are the contributions to the Poisson equations from the field.

In both the linear and nonlinear analyses we solve the Friedmann constraint,
Eq. (4), alongside the fluid equations of motion Eq. (33) and Eq. (34) and the
two Poisson equations for $\Psi$, Eq. (35), and $\dot{\Psi}$, Eq. (36). In our
linear analysis, we additionally separate out the field average from its
fluctuations, $\varphi=\bar{\varphi}+\delta\varphi$, and solve for the
dynamics of the homogeneous mode, $\bar{\varphi}$, via Eq. (8) separately from
the perturbations, $\delta\varphi$, see equation Eq. (17). We complete this
analysis entirely in momentum space.

In our nonlinear analysis we utilized a modified version of GABE (Grid And
Bubble Evolver) Child and Giblin (2012) that simulates the scalar field in
configuration space according to Eq. (3) while still evolving the fluid
variables in momentum space. The sources, (37) and (38), are calculated in
configuration space, then Fourier transformed to allow us to invert Eqns. (35)
and (36) to find the gravitational perturbations. These perturbations, $\Psi$
and $\dot{\Psi}$ are stored both in momentum space (to be used in the
evolution equations for the fluids) and inverse Fourier transformed into
configuration space (to be used in the evolution equation for the field). This
procedure allows us to solve for all linearized quantities in momentum space,
while allowing us to treat the field evolution nonlinearly and fully resolve
the dynamics of the configuration-space metric perturbations.

In order to generate initial conditions for our lattice simulations, we
numerically solve the set of coupled differential equations in a simplified
Einstein-Boltzmann hierarchy and approximate recombination as instantaneous.
For this analysis, we include a tightly coupled baryon-photon fluid, CDM,
neutrinos, and the scalar field. The neutrinos are treated as a perfect fluid
(i.e. their anisotropic stress vanishes, $\sigma_{\nu}=0$). We have also
evolved the system with free-streaming neutrinos with a Boltzmann hierarchy
that is truncated at the third moment (while still using the approximation
$\Phi=\Psi$), using the proscription outlined in Ref. Ma and Bertschinger
(1995), and found no difference in the resulting field dynamics.

In the lattice simulations we treat the CDM/baryons as a single matter fluid,
and the photons/neutrinos as a single radiative constituent. We chose the
following parameters for all of our simulations: $A_{s}=2\times 10^{-9}$,
$h=0.67$, $\Omega_{M}=\Omega_{c}+\Omega_{b}=0.314$,
$\Omega_{R}=\Omega_{\gamma}+\Omega_{\nu}=9.16\times 10^{-5}$. We evolve the
system from an initial scale factor, $a_{i}$, is set to be small enough so
that the field dynamics are linear, to $a_{f}=1/30$, using a time step of
$\Delta t=L/N/100$, where $L$ is the length of an edge of the simulation box
and $N$ is the number of pixels on one side of the box. Our fiducial lattice
size is $N^{3}=256^{3}$ and we choose $L$ to ensure that the resonant
wavenumber, $k_{\rm res}=2\pi/\lambda_{\rm res}$, is well inside of the box,
$L=20\lambda_{\rm res}$. We have confirmed that using smaller boxes do not
alter our results.

##  IV Resonance and nonlinear evolution of the Scalar Field

###  IV.1 The homogeneous mode

In this Subsection we will focus on two models–where we set
$\bar{\rho}_{\varphi}(a_{c})/\bar{\rho}_{\rm tot}(a_{c})=0.05$,
$z_{c}=10^{4.5}$, and choose either $\theta_{i}=1.5$ or $\theta_{i}=3$.

We start by comparing the homogeneous evolution of the scalar field in our
lattice simulation with our linear calculations. The field $\varphi$ starts
out as roughly homogeneous, and starts oscillating when $a=a_{c}$, where
$a_{c}$ is determined by
$a_{c}^{2}V^{\prime\prime}(\varphi(a_{c}))\sim\mathcal{H}^{2}(a_{c})$. As the
field enters a period of coherent oscillations, we expect that the
contribution to $\bar{\rho}_{\varphi}$ ceases to look dark energy-like
(approximately constant) and starts to look radiation-like – since the minimum
is massless. For our linear analysis, the homogeneous mode will continue to
oscillate about its minimum, decaying only due to Hubble friction. However, in
the nonlinear analysis, we expect the homogeneous mode to show signs of
earlier decay when the field exits the linearized regime.

Fig. 2 shows a comparison of the homogeneous evolution of the scalar field
between the linear and nonliear simulations. We can see that, in each case,
there exists a $z_{\rm nl}$ at which the homogenous mode in the fully
nonlinear simulation starts to decay away from the linear solution. This is
due to the transfer of energy from the homogeneous mode to the $k>0$ modes,
indicating the presence of resonance from the nonlinear self-coupling.

Figure 2: The background evolution of the scalar field field for
$z_{c}=10^{4.5}$ from our nonlinear simulations (blue) compared to a
background-only calculation (red) with $\theta_{i}=3$ (top panel) and
$\theta_{i}=1.5$ (bottom panel). The decrease in amplitude in the nonlinear
simulation comes from the coupling between the homogeneous mode and the
perturbations and indicates that the field dynamics are nonlinear.

This also gives us an opportunity to test the accuracy of the analytic
equations presented in Sec. II.4. The top panel of Fig. 1 shows that for
$z_{c}=10^{4.5}$, $\theta_{i}=3.0$ has $z_{\rm nl}\simeq 500$ and for
$\theta_{i}=1.5$ has $z_{\rm nl}\simeq 120$. Comparing this to the values of
$z_{\rm nl}$ denoted in Fig. 2 ($z_{\rm nl}=300$ and $z_{\rm nl}=70$,
respectively) shows that the analytic formulae provide a good approximation to
the redshift of nonlinearity within a factor of order unity.

Figure 3: A comparison between the dimensionless power spectra for different
quantities at a later time in the simulation (in this case $z_{c}=10^{4.5}$
and $\theta_{i}=1.5$). For each panel, the blue dashed line shows the power
spectrum from a linear analysis at $z=500$, while the red (solid) curves come
from the $z=500$ slice from a run initialized at $z=1000$. The resonant wave
number, $k_{\rm res}\simeq 0.26\ {\rm Mpc}^{-1}$, is indicated by the dotted
line in the top two panels.

###  IV.2 The inhomogeneous dynamics

We begin by validating our simulations by showing consistency between our
linear and nonlinear analyses at a time when the nonlinear simulations have
significantly evolved, but where we still expect to see excellent agreement
with the linear analysis. Fig. 3 shows how well the nonlinear simulations
reproduce the mode-by-mode behavior of the radiation and matter fluids, the
field and the gravitational perturbations when compared to the linear
simulations. It is clear that both the field and fluid dynamics are being
solved correctly at the linear level. Small differences between the two
spectra at large scales are due to sample variance and binning in the
nonlinear simulations. Note that the Newtonian potential, $\Psi$, is mainly
sourced by the matter, and its time derivative is mainly sourced by the
(oscillating) radiation content.

In order to further make a connection between the full simulation and the
analytic expressions in Sec. II.4, we plot the evolution of the resonant
wavenumber in Fig. 4. As expected, once the field becomes nonlinear the
perturbations remain relatively constant. The resulting contribution to the
ISW is shown in the bottom panel of Fig. 4. There we can see that this
contribution peaks at $z_{\rm nl}$ and, as expected from the middle panel of
Fig. 1, the peak with $\theta_{i}=3$ is higher than $\theta_{i}=1.5$. We can
make this comparison more quantitative by noting that Eq. (II.4) gives
$\Delta^{2}_{\dot{\Psi}/\mathcal{H}}|_{k=k_{\rm res};\theta_{i}=3}=1.5\times
10^{-8}$ and $\Delta^{2}_{\dot{\Psi}/\mathcal{H}}|_{k=k_{\rm
res};\theta_{i}=1.5}=8\times 10^{-10}$ which is within a factor of a few of
the results from the simulation in the bottom panel of Fig. 4, which give
$3.5\times 10^{-9}$ and $2\times 10^{-10}$, respectively.

Figure 4: The evolution of the scalar field dimensionless power spectrum
normalized by the field’s envelope at the resonant wavenumber (top panel). We
can see that once the field becomes nonlinear the field’s dimensionless power
spectrum is relatively constant. The dimensionless power spectrum of
$\mathcal{H}^{-1}\partial\Psi/\partial\tau$ at the resonant wavenumber. As
discussed in Sec. II.4, this quantity measures the contribution that the
scalar field makes to the late-time ISW effect (bottom panel). In both panels
the simulation with $\theta_{i}=3$ is shown in blue and $\theta_{i}=1.5$ in
red.

###  IV.3 Calculating the ISW contribution

The contribution to the ISW is calculated by computing Eq. (23) along several
lines of sight through the simulation box; a cartoon of the procedure is given
in Fig. 5. Slices are taken along the line of sight through the box, from
$x=0$ to $x=L$. We use periodic boundary conditions, and so after $x=L$ our
next slice is taken at $x=0$, and we integrate over the box again. We compute
this integral from some initial redshift $z_{i}\gg z_{\rm nl}$ up to a time
that is late enough, $z_{f}$, so that the ISW contribution at the final
redshift slice is negligible (as shown in Fig. 4 the scalar field’s
contribution to $\dot{\Psi}$ decreases rapidly after $z_{\rm nl}$). We note
that repeatedly traversing the same box will introduce some spurious
correlations. However, given that the light-crossing time is larger than the
dynamical time-scales for $k_{\rm res}\gg 1/L$, we expect they will be small
for modes well within the simulation box.

Figure 5: A cartoon describing the line-of-sight integration used to compute
the scalar field’s ISW contribution, given in Eq. (23).

Figure 6: The ISW contribution from a scalar field that addresses the Hubble
tension with $\theta_{i}=2.18$ (red) and $\theta_{i}=3.14$ (pink) compared to
the primary CMB (black) (left panel). The ISW contribution from four different
cosmological scalar field models (right panel). Unlike the left panel, the
chosen parameters are not constrained to resolve the Hubble tension. The
dashed orange curve in the right panel shows the residual foreground
contribution required to achieve the science goals of CMB-HD.

Once we have evaluated (23), we end up with a partial-sky map of the
contribution to the ISW from $\dot{\Psi}$ from our nonlinear simulations. We
convert the spatial location of each ray at $z_{f}$ to an angular position on
our sky and compute the resulting 2D Fourier transform of the angular map to
determine the ISW contribution to the CMB temperature power spectrum. We have
confirmed that the power spectra we compute is insensitive to moving the final
redshift slice closer to today.

##  V The scalar field ISW contribution

We are now in a position to compute the ISW contribution for this cosmological
scalar field. In order to convert the temperature power spectra to units of
$\mu{\rm K}^{2}$, we take $T_{\rm CMB}=2.7255\ {\rm K}$ Fixsen (2009).

###  V.1 The contribution from a scalar field that address the Hubble tension

It is of interest to determine whether a scalar field which may address the
Hubble tension may include nonlinear dynamics which will produce a significant
contribution to the late-time ISW effect. The best-fit parameters for such a
field is given in the Appendix of Ref. Smith  _et al._ (2020):
$\log_{10}z_{c}=3.52$, $\rho_{\varphi}(a_{c})/\rho_{\rm tot}(a_{c})=0.09$,
$\theta_{i}=2.18$ which corresponds to $m=6\times 10^{-28}\ {\rm eV}$ and
$f/m_{\rm pl}=0.2$. These parameters lead to $k_{\rm res}=0.043\ {\rm
Mpc}^{-1}$ and $z_{\rm nl}\simeq 3.3$, which imply that the ISW contribution
will peak at $l_{\rm peak}\simeq k_{\rm res}(\tau_{0}-\tau_{\rm nl})\simeq
300$.

As shown in Eq. (II.4), the overall contribution to the ISW effect is
primarily determined by the ratio $a_{c}/a_{\rm nl}$, and in the case of a
field that resolves the Hubble tension this factor is $\sim 10^{-3}$, which
already gives us an indication that the ISW from this model will be very
small. Indeed, as shown in the left panel of Fig. 6, we can see that the
resulting power spectrum is at most $\sim 10^{-2}\ (\mu{\rm K})^{2}$, well
below the cosmic variance limit. Also note, as we anticipated, the
contribution shows a peak at $l\sim 100$.

The main reason that this model makes such a small ISW contribution is the
wide gap between $a_{c}$ and $a_{\rm nl}$. Since $a_{c}\sim a_{\rm eq}$ in
this case, the ISW contribution scales as $(a_{c}/a_{\rm nl})^{3}$ leading to
a suppression of order $10^{-9}$. Fig. 1 shows that, with $z_{c}$ fixed in
order to resolve the Hubble tension, the only way to move these two redshifts
closer is by increasing $\theta_{i}$ (and thereby move towards a flatter part
of the potential). The pink curve in the left panel of Fig. 6 shows how the
ISW contribution grows in amplitude and moves to a smaller scale as
$\theta_{i}\rightarrow\pi$. In this limit we also have
$V^{\prime\prime}(\varphi)<0$, raising the possibility of a tachyonic
instability. However, as shown in Ref. Smith  _et al._ (2020), this
instability is only present for isocurvature perturbations, and even in that
case, leads to a relatively small enhancement of the perturbations. Note that,
for this specific form of the scalar field potential, in order to resolve the
Hubble tension we must have $\theta_{i}\lesssim 2.5$ at the 95% confidence
level Smith  _et al._ (2020). However, other forms for the scalar field
potential (such as $\alpha$-attractors Braglia  _et al._ (2020)) may both
resolve the Hubble tension and allow the field to start at a flatter part of
its potential.

###  V.2 The contribution from a general scalar field

If we allow the scalar field parameters to vary we can produce a measurable
ISW effect. As discussed in the previous Section, the main way to boost the
ISW contribution is to move $z_{\rm nl}$ closer to $z_{c}$. The top panel of
Fig. 1 shows that the most effective way to do this is to increase $z_{c}$.
The bottom pane of Fig. 1 shows that by increasing $z_{c}$ we will also cause
the peak of the ISW contribution to shift to smaller scales, i.e., larger
multipoles.

Fig. 6 (right panel) shows the ISW generated by four possible cosmological
scalar fields. The parameters were chosen so as to produce contributions with
similar amplitudes by exchanging a lower value of $z_{c}$ for a larger value
of $\theta_{i}$ (note that these four models are marked in the middle panel of
Fig. 1 by the four stars). All four models have
$\bar{\rho}_{\varphi}(a_{c})/\bar{\rho}_{\rm tot}(a_{c})=0.05$ (saturating the
95% upper limit from the linear effects of these scalar fields on current
cosmological data).

The power spectra show two clear peaks, which correspond to the ‘fundamental’
and to the ‘first harmonic’ of the resonant wavenumber. At smaller scales the
nonlinear scalar field ISW effect produces a scale-invariant tail. The rise of
the power spectra at the lowest multipoles is, at least in part, due to edge
effects in the simulation, and thus the behavior at these multipoles should
not be taken as a physical result of our analysis.

Although current CMB measurements do not have the angular resolution and
sensitivity to detect these features, future observatories will. The dotted
line in Fig. 6 shows the error bars due to sample variance and instrumental
noise associated with CMB-HD Sehgal  _et al._ (2019); Aiola  _et al._ (2022),
a proposed high-resolution CMB ground-based observatory which would image half
of the sky. Note that the noise curve does not include contributions from
residual uncertainty from the subtraction of foreground contributions. CMB-HD
will need to remove foregrounds to a level that ensures the instrument noise
is comparable to or higher than the residual foreground contribution for
$\ell>7,500$ in order to achieve its CMB lensing science goals Han and Sehgal
(2022), corresponding to a foreground level of about $\simeq 0.3\ \mu{\rm
K}^{2}$.

##  VI Conclusions

In this paper we have calculated the ISW effect from subdominant massless
scalar fields, which can potentially play a role in addressing the Hubble
tension, in the late universe. We present both analytic estimates as well as
high-resolution lattice simulations. We not only included the full nonlinear
evolution of the scalar field, but also included the evolution of (linearized)
radiation and matter fluid perturbations, as well as metric perturbations
sourced by the inhomogeneous scalar field and the fluids. From these nonlinear
simulations, we have evaluated the ISW contribution to photons traveling
through the simulation, validated our results against CLASS in the appropriate
limit, and compared these results, for a set of test cases, to the sensitivity
of CMB-HD, a proposed CMB observatory which is designed to make resolution
measurements.

As anticipated by our scaling equations, the full nonlinear simulations show
that the amplitude of the ISW contribution is mainly determined by the ratio
$a_{c}/a_{\rm nl}$. This is shown through a dramatic increase in the strength
of the ISW effect when choosing parameters that ensure nonlinearity sets in at
an earlier time. The impact can also be enhanced by assuming larger initial
field displacements, which lead to stronger, earlier resonance.

While there is wide theoretical motivation for studying extra scalar degrees
of freedom, a recent and popular invocation of such a model has been to ease
the emerging Hubble tension Poulin  _et al._ (2023) via an EDE field. From
what we show here, early dark energy with a potential of the form
$V=m^{2}f^{2}(1-\cos\varphi/f)^{2}$, that can resolve the Hubble tension will
produce a contribution to the ISW that is unlikely to be constrained or
validated in upcoming CMB experiments. We make this choice as a concrete
proof-of-concept, although there has been some indication that such a
potential can be realized from high energy physics McDonough and Scalisi
(2022); Cicoli  _et al._ (2023). However, models which use a different form of
the potential– such as $\alpha$-attractors Kallosh and Linde (2013); Kallosh
_et al._ (2013); Galante  _et al._ (2015); Braglia  _et al._ (2020)– may
resolve the Hubble tension just as well and at an initial field value where
the potential is flatter producing a larger ISW effect. The flattening of the
potentials in such models can also bring in additional non-perturbative
dynamical effects Lozanov and Amin (2018) at early times. We leave an
exploration of different scalar field potentials to future work.

Although we have established that for some model parameters the EDE ISW is
larger than the noise in an instrument like CMB-HD (see Fig. 6), it remains to
be seen whether this signal can be distinguished from expected foreground
contamination at these small angular scales. We note that the specific
harmonic structure of the EDE ISW contribution, relative to the smooth
foreground power spectra, may aide in its detectability.

In addition to the effect on the evolution of the gravitational potentials,
the nonlinear fragmentation of the scalar field also source gravitational
waves at frequencies of order $k_{\rm res}/a_{\rm nl}$. The post recombination
production of gravitational waves at such small frequencies might lead to
additional $B$-mode polarization at reionization (e.g., Ref. Geller  _et al._
(2021)). It may also be possible to see the ISW contribution using other
probes (i.e., 21cm observations Raccanelli  _et al._ (2016)), and would be
worth considering in the future.

Cosmological scalar fields provide a rich phenomonology which touches all
aspects of cosmology, from inflation to the current epoch of accelerated
expansion, from linear to resonant nonlinear dynamics. Here we have shown that
scalar fields may make a novel contribution to the ISW effect, imparting
characteristic features in the very small-scale CMB temperature anisotropies.
This adds further motivation to build a CMB observatory targeting the very
small-scale anisotropies (such as CMB-HD) which will give us access to new
ways to understand the origin and evolution of the universe.

###### Acknowledgements.

We thank Kaloian Lozanov for significant initial involvement, including
setting up preliminary scalar-field lattice simulations and calculations of
the corresponding ISW contribution. We thank Neelima Sehgal, Dongwong Han, and
Amanda Macinnis, for discussions about CMB-HD. We thank Yacine Ali-Haïmoud,
Daniel Grin and Zachary Weiner for prompt, detailed and insightful feedback on
the manuscript. The numerical work presented here utilized the Strelka
Computing Cluster, which is run by Swarthmore College, and the Vera Rubin
Cluster at Kenyon College, which is supported by the National Science
Foundation (NSF) and the Kenyon College Department of Physics. T.L.S. is
supported by the NSF, Grant No. AST-2009377 and the Research Corporation.
J.T.G., M.G., and E.F. are supported by the NSF, Grant No. PHY-2013718. M.G.
is supported by the U.S. Department of Energy, Office of Science, Office of
Advanced Scientific Computing Research, Department of Energy Computational
Science Graduate Fellowship under Award Number DE-SC0023112. M. A. is
supported by a NASA ATP-Theory Grant 80NSSC20K0518. This report was prepared
as an account of work sponsored by an agency of the United States Government.
Neither the United States Government nor any agency thereof, nor any of their
employees, makes any warranty, express or implied, or assumes any legal
liability or responsibility for the accuracy, completeness, or usefulness of
any information, apparatus, product, or process disclosed, or represents that
its use would not infringe privately owned rights. Reference herein to any
specific commercial product, process, or service by trade name, trademark,
manufacturer, or otherwise does not necessarily constitute or imply its
endorsement, recommendation, or favoring by the United States Government or
any agency thereof. The views and opinions of authors expressed herein do not
necessarily state or reflect those of the United States Government or any
agency thereof.

## References

  * Guth (1981) Alan H. Guth, “The Inflationary Universe: A Possible Solution to the Horizon and Flatness Problems,” Phys. Rev. D 23, 347–356 (1981).
  * Caldwell  _et al._ (1998) R. R. Caldwell, Rahul Dave, and Paul J. Steinhardt, “Cosmological imprint of an energy component with general equation of state,” Phys. Rev. Lett. 80, 1582–1585 (1998), arXiv:astro-ph/9708069 .
  * Baumann and McAllister (2015) Daniel Baumann and Liam McAllister,  _Inflation and String Theory_, Cambridge Monographs on Mathematical Physics (Cambridge University Press, 2015) arXiv:1404.2601 [hep-th] .
  * Svrcek and Witten (2006) Peter Svrcek and Edward Witten, “Axions In String Theory,” JHEP 06, 051 (2006), arXiv:hep-th/0605206 .
  * Arvanitaki  _et al._ (2010) Asimina Arvanitaki, Savas Dimopoulos, Sergei Dubovsky, Nemanja Kaloper, and John March-Russell, “String Axiverse,” Phys. Rev. D 81, 123530 (2010), arXiv:0905.4720 [hep-th] .
  * Kane  _et al._ (2015) Gordon Kane, Kuver Sinha, and Scott Watson, “Cosmological Moduli and the Post-Inflationary Universe: A Critical Review,” Int. J. Mod. Phys. D 24, 1530022 (2015), arXiv:1502.07746 [hep-th] .
  * Stott  _et al._ (2017) Matthew J. Stott, David J. E. Marsh, Chakrit Pongkitivanichkul, Layne C. Price, and Bobby S. Acharya, “Spectrum of the axion dark sector,” Phys. Rev. D 96, 083510 (2017), arXiv:1706.03236 [astro-ph.CO] .
  * Griest (2002) Kim Griest, “Toward a possible solution to the cosmic coincidence problem,” Phys. Rev. D66, 123501 (2002), arXiv:astro-ph/0202052 [astro-ph] .
  * Linder and Smith (2011) Eric V. Linder and Tristan L. Smith, “Dark Before Light: Testing the Cosmic Expansion History through the Cosmic Microwave Background,” JCAP 04, 001 (2011), arXiv:1009.3500 [astro-ph.CO] .
  * Kamionkowski  _et al._ (2014) Marc Kamionkowski, Josef Pradler, and Devin G. E. Walker, “Dark energy from the string axiverse,” Phys. Rev. Lett. 113, 251302 (2014), arXiv:1409.0549 [hep-ph] .
  * Amin  _et al._ (2012) Mustafa A. Amin, Phillip Zukin, and Edmund Bertschinger, “Scale-Dependent Growth from a Transition in Dark Energy Dynamics,” Phys. Rev. D 85, 103510 (2012), arXiv:1108.1793 [astro-ph.CO] .
  * Johnson and Kamionkowski (2008) Matthew C. Johnson and Marc Kamionkowski, “Dynamical and Gravitational Instability of Oscillating-Field Dark Energy and Dark Matter,” Phys. Rev. D 78, 063010 (2008), arXiv:0805.1748 [astro-ph] .
  * Lozanov and Amin (2017) Kaloian D. Lozanov and Mustafa A. Amin, “Equation of State and Duration to Radiation Domination after Inflation,” Phys. Rev. Lett. 119, 061301 (2017), arXiv:1608.01213 [astro-ph.CO] .
  * Lozanov and Amin (2018) Kaloian D. Lozanov and Mustafa A. Amin, “Self-resonance after inflation: oscillons, transients and radiation domination,” Phys. Rev. D 97, 023533 (2018), arXiv:1710.06851 [astro-ph.CO] .
  * Khlebnikov and Tkachev (1996) S. Yu. Khlebnikov and I. I. Tkachev, “Classical decay of inflaton,” Phys. Rev. Lett. 77, 219–222 (1996), arXiv:hep-ph/9603378 .
  * Amin  _et al._ (2014) Mustafa A. Amin, Mark P. Hertzberg, David I. Kaiser, and Johanna Karouby, “Nonperturbative Dynamics Of Reheating After Inflation: A Review,” Int. J. Mod. Phys. D 24, 1530003 (2014), arXiv:1410.3808 [hep-ph] .
  * Lozanov (2019) Kaloian D. Lozanov, “Lectures on Reheating after Inflation,” (2019), arXiv:1907.04402 [astro-ph.CO] .
  * Poulin  _et al._ (2019) Vivian Poulin, Tristan L. Smith, Tanvi Karwal, and Marc Kamionkowski, “Early Dark Energy Can Resolve The Hubble Tension,” Phys. Rev. Lett. 122, 221301 (2019), arXiv:1811.04083 [astro-ph.CO] .
  * Smith  _et al._ (2020) Tristan L. Smith, Vivian Poulin, and Mustafa A. Amin, “Oscillating scalar fields and the Hubble tension: a resolution with novel signatures,” Phys. Rev. D 101, 063523 (2020), arXiv:1908.06995 [astro-ph.CO] .
  * Sehgal  _et al._ (2019) Neelima Sehgal _et al._ , “CMB-HD: An Ultra-Deep, High-Resolution Millimeter-Wave Survey Over Half the Sky,” (2019), arXiv:1906.10134 [astro-ph.CO] .
  * Aiola  _et al._ (2022) Simone Aiola _et al._ (CMB-HD), “Snowmass2021 CMB-HD White Paper,” (2022), arXiv:2203.05728 [astro-ph.CO] .
  * Giblin  _et al._ (2019) John T. Giblin, James B. Mertens, Glenn D. Starkman, and Chi Tian, “Limited accuracy of linearized gravity,” Phys. Rev. D 99, 023527 (2019), arXiv:1810.05203 [astro-ph.CO] .
  * Giblin and Tishue (2019) John T. Giblin and Avery J. Tishue, “Preheating in Full General Relativity,” Phys. Rev. D 100, 063543 (2019), arXiv:1907.10601 [gr-qc] .
  * Hlozek  _et al._ (2018) Renée Hlozek, David J. E. Marsh, and Daniel Grin, “Using the Full Power of the Cosmic Microwave Background to Probe Axion Dark Matter,” Mon. Not. Roy. Astron. Soc. 476, 3063–3085 (2018), arXiv:1708.05681 [astro-ph.CO] .
  * Seljak and Zaldarriaga (1996) Uros Seljak and Matias Zaldarriaga, “A Line of sight integration approach to cosmic microwave background anisotropies,” Astrophys. J. 469, 437–444 (1996), arXiv:astro-ph/9603033 .
  * Felder and Tkachev (2008) Gary N. Felder and Igor Tkachev, “LATTICEEASY: A Program for lattice simulations of scalar fields in an expanding universe,” Comput. Phys. Commun. 178, 929–932 (2008), arXiv:hep-ph/0011159 .
  * Weiner  _et al._ (2020) Zachary J. Weiner, Peter Adshead, and John T. Giblin, “Constraining early dark energy with gravitational waves before recombination,” (2020), arXiv:2008.01732 [astro-ph.CO] .
  * Ma and Bertschinger (1995) Chung-Pei Ma and Edmund Bertschinger, “Cosmological perturbation theory in the synchronous and conformal Newtonian gauges,” Astrophys. J. 455, 7–25 (1995), arXiv:astro-ph/9506072 .
  * Child and Giblin (2012) Hillary L. Child and John T. Giblin, Jr., “Gravitational Radiation from First-Order Phase Transitions,” JCAP 10, 001 (2012), arXiv:1207.6408 [astro-ph.CO] .
  * Fixsen (2009) D. J. Fixsen, “The Temperature of the Cosmic Microwave Background,” Astrophys. J. 707, 916–920 (2009), arXiv:0911.1955 [astro-ph.CO] .
  * Braglia  _et al._ (2020) Matteo Braglia, William T. Emond, Fabio Finelli, A. Emir Gumrukcuoglu, and Kazuya Koyama, “Unified framework for early dark energy from $\alpha$-attractors,” Phys. Rev. D 102, 083513 (2020), arXiv:2005.14053 [astro-ph.CO] .
  * Han and Sehgal (2022) Dongwon Han and Neelima Sehgal, “Mitigating foreground bias to the CMB lensing power spectrum for a CMB-HD survey,” Phys. Rev. D 105, 083516 (2022), arXiv:2112.02109 [astro-ph.CO] .
  * Poulin  _et al._ (2023) Vivian Poulin, Tristan L. Smith, and Tanvi Karwal, “The Ups and Downs of Early Dark Energy solutions to the Hubble tension: a review of models, hints and constraints circa 2023,” (2023), arXiv:2302.09032 [astro-ph.CO] .
  * McDonough and Scalisi (2022) Evan McDonough and Marco Scalisi, “Towards Early Dark Energy in String Theory,” (2022), arXiv:2209.00011 [hep-th] .
  * Cicoli  _et al._ (2023) Michele Cicoli, Matteo Licheri, Ratul Mahanta, Evan McDonough, Francisco G. Pedro, and Marco Scalisi, “Early Dark Energy in Type IIB String Theory,” (2023), arXiv:2303.03414 [hep-th] .
  * Kallosh and Linde (2013) Renata Kallosh and Andrei Linde, “Universality Class in Conformal Inflation,” JCAP 07, 002 (2013), arXiv:1306.5220 [hep-th] .
  * Kallosh  _et al._ (2013) Renata Kallosh, Andrei Linde, and Diederik Roest, “Superconformal Inflationary $\alpha$-Attractors,” JHEP 11, 198 (2013), arXiv:1311.0472 [hep-th] .
  * Galante  _et al._ (2015) Mario Galante, Renata Kallosh, Andrei Linde, and Diederik Roest, “Unity of Cosmological Inflation Attractors,” Phys. Rev. Lett. 114, 141302 (2015), arXiv:1412.3797 [hep-th] .
  * Geller  _et al._ (2021) Michael Geller, Sida Lu, and Yuhsin Tsai, “B modes from postinflationary gravitational waves sourced by axionic instabilities at cosmic reionization,” Phys. Rev. D 104, 083517 (2021), arXiv:2104.08284 [hep-ph] .
  * Raccanelli  _et al._ (2016) Alvise Raccanelli, Ely Kovetz, Liang Dai, and Marc Kamionkowski, “Detecting the integrated Sachs-Wolfe effect with high-redshift 21-cm surveys,” Phys. Rev. D 93, 083512 (2016), arXiv:1502.03107 [astro-ph.CO] .