# Instabilities of scalar fields around oscillating stars

Taishi Ikeda  Dipartimento di Fisica, “Sapienza” Universitá di Roma, Piazzale
Aldo Moro 5, 00185, Roma, Italy  Vitor Cardoso  CENTRA, Departamento de
Física, Instituto Superior Técnico – IST, Universidade de Lisboa – UL, Avenida
Rovisco Pais 1, 1049-001 Lisboa, Portugal  Niels Bohr International Academy,
Niels Bohr Institute, Blegdamsvej 17, 2100 Copenhagen, Denmark  Miguel Zilhão
CENTRA, Departamento de Física, Instituto Superior Técnico – IST, Universidade
de Lisboa – UL, Avenida Rovisco Pais 1, 1049-001 Lisboa, Portugal

###### Abstract

The behavior of fundamental fields in strong gravity or nontrivial
environments is important for our understanding of nature. This problem has
interesting applications in the context of dark matter, of dark energy physics
or of quantum field theory. The dynamics of fundamental fields has been
studied mainly in static or stationary backgrounds, whereas most of our
Universe is dynamic. In this paper we investigate “blueshift” and parametric
instabilities of scalar fields in dynamical backgrounds, which can be
triggered (for instance) by oscillating stars in scalar-tensor theories of
gravity. We discuss possible implications of our results, which include
constraints on an otherwise hard-to-access parameter space of scalar-tensor
theories.

##  I Introduction

General Relativity (GR) is currently the best description of the gravitational
interaction, and has been successfully tested on different scales Will (2014).
The recent direct detection of gravitational waves (GWs) Abbott  _et al._
(2016a, b, 2020) indicates that even dynamical, strong-field regions are
adequately described by GR (up to the precision probed by current detectors).
In spite of its brilliant status, there are a number of conceptual issues with
GR, ranging from the large-scale description of the cosmos to the fate of
classical singularities in gravitational collapse or the incorporation of
quantum effects, see e.g. Clifton  _et al._ (2012); Berti  _et al._ (2015);
Cardoso and Pani (2019) and references therein. The resolution of some of
these challenges most likely requires that GR be superseded by a more
sophisticated description.

There is currently no single compelling alternative to GR that solves the
above issues without introducing new problems of their own. However, a variety
of modified theories have been proposed, mostly with the view to exploring the
mathematical and physics content of possible contenders to GR. These
frequently include additional degrees of freedom, which might lead to unique
observational signatures. The simplest modification of GR are scalar-tensor
theories, where a new scalar degree of freedom couples to curvature or matter.
Some of these theories arise naturally as possible alternatives, since they
have a well-posed initial value problem and simultaneously evade all known
constraints. Scalars are also a generic prediction of string theory or of
extensions of the standard model of particle physics, and are also natural
candidates for dark energy and dark matter Arvanitaki  _et al._ (2010);
Clifton  _et al._ (2012); Berti  _et al._ (2015); Barack  _et al._ (2019);
Cardoso and Pani (2019).

Table 1: List of some stellar objects which can be prone to the processes described here. We take the radial mode of a neutron star (NS) of mass $1.3M_{\odot}$ and radius $L_{0}=9.7\,{\rm Km}$ described by equation of state A in Ref. Kokkotas and Ruoff (2001). For white dwarfs (WD) we take a degeneracy parameter $1/y_{0}^{2}=0.05$ as quoted in Refs. Schatzman (1961); Singh and Das (1991); Winget and Kepler (2008). The Sun’s description was taken from Refs. Turck-Chieze and Lopes (1993); Turck-Chièze and Lopes (2012). The compactness $\mathcal{C}$ is defined in Eq. (1).  | $\mathcal{C}$ |  $L_{0}$ (km) | $\omega L_{0}/c$ | Radial frequency
---|---|---|---|---
NS | $0.3$ | $10$ | 0.6 | 3 kHz
WD | $10^{-3}$ | $10^{3}$ | 0.004 | 0.2 Hz
Sun | $10^{-6}$ | $7\times 10^{5}$ | 0.029 | 2 mHz

Depending on the coupling to matter, to curvature and on the self-
interactions, such theories and new fundamental fields may give rise to wide
array of new effects, such as the spontaneous “scalarization” of objects
Damour and Esposito-Farese (1993, 1996); Harada (1998); Pani  _et al._ (2011);
Ramazanoğlu and Pretorius (2016); Antoniou  _et al._ (2018); Doneva and
Yazadjiev (2018); Silva  _et al._ (2018), screening mechanisms on
astrophysical scales Khoury and Weltman (2004a, b); Joyce  _et al._ (2015),
etc. Potential astrophysical consequences of these theories were mostly based
on analysis in stationary settings Damour and Esposito-Farese (1993, 1996);
Harada (1998); Pani  _et al._ (2011); Ramazanoğlu and Pretorius (2016);
Cardoso  _et al._ (2020); Khoury and Weltman (2004a, b); Tsujikawa  _et al._
(2009). Time-dependent setups include spontaneous scalarization during the
inspiral and merger of neutron stars Barausse  _et al._ (2013); Palenzuela
_et al._ (2014); Shibata  _et al._ (2014); Khalil  _et al._ (2019), black
holes Silva  _et al._ (2021); East and Ripley (2021), cosmological evolutions
Llinares and Mota (2013, 2013); Hagala  _et al._ (2017) or situations aimed at
understanding well-posedness or other fundamental issues Brito  _et al._
(2014); Silvestri (2011); Nakamura  _et al._ (2021).

Here, we wish to understand possible new phenomena induced by time-periodic
motion, in particular by vibrating stars such as the ones summarized in Table
1. Theories for which a constant scalar is a ground state typically have the
same stationary solutions as GR, making it challenging to tell the two apart.
However, when such stars are disturbed (stochastically, like our Sun, or via
mergers or accretion for a neutron star), they provide a time-periodic
background on which a scalar fluctuation propagates. Indeed, we show that
oscillating stars trigger various instabilities of fundamental fields. Such
instabilities may facilitate constraints on otherwise hard-to-access parameter
space of the theory.

A toy model had been studied in a Minkowski spacetime Wang  _et al._ (2011,
2013). Our results confirm and extend the general picture, but in a nontrivial
way, by dealing also with general relativistic backgrounds and fluctuations
and providing realistic timescales for the mechanism (details on the general
relativistic case are discussed in the appended Supplemental material).

We use geometrized units $c=G=1$ and parameterize stars by their mass $M$,
radius $L_{0}$, and compactness

$\mathcal{C}\equiv\frac{2M}{L_{0}}\,.$ (1)

##  II Setup

Our main purpose is to show that oscillating backgrounds can trigger
instabilities of fundamental fields nonminimally coupled to matter. We will
show this explicitly for the simplest example of a scalar-tensor theory, but
our results are valid for more general setups.

We focus on theories in which a scalar field $\Phi$ is coupled to matter,
described by the Lagrangian density

${\cal
L}=\frac{R[g]}{16\pi}-\frac{1}{2}g^{\mu\nu}\nabla_{\mu}\Phi\nabla_{\nu}\Phi-V(\Phi)+{\cal
L}_{m}\left[\Psi_{m},\tilde{g}\right]\,,$ (2)

where $\tilde{g}_{\mu\nu}=A(\Phi)^{2}g_{\mu\nu}$ is the physical metric and
${\cal L}_{m}$ the Lagrangian describing matter fields $\Psi_{m}$ (these could
be the microscopic, fundamental description leading to notions of density and
pressure for example). The equations of motion for the system are

$\displaystyle R_{\mu\nu}-\frac{1}{2}Rg_{\mu\nu}=8\pi T_{\mu\nu}\,,$ (3a)
$\displaystyle T_{\mu\nu}=\nabla_{\mu}\Phi\nabla_{\nu}\Phi-
g_{\mu\nu}\left(\frac{1}{2}(\nabla\Phi)^{2}+V(\Phi)\right)+T^{(m)}_{\mu\nu}\,,$
$\displaystyle\square\Phi-V^{\prime}(\Phi)=-A(\Phi)^{3}\partial_{\Phi}A(\Phi)\tilde{T}^{(m)}\,.$
(3b)

Here, $\tilde{T}^{(m)}=\tilde{g}^{\mu\nu}\tilde{T}^{(m)}_{\mu\nu}$ while
$T^{(m)}_{\mu\nu}$ and $\tilde{T}^{(m)}_{\mu\nu}$ are the energy-momentum
tensors of the matter fields with respect to metric $g_{\mu\nu}$ and
$\tilde{g}_{\mu\nu}$.

The phenomena we wish to describe are also part of more generic theories: the
fundamental ingredient is a position-dependent effective mass function for the
scalar $\Phi$. Once linearized around a specific background, the right-hand
side of Eq. (3b) is proportional to $\tilde{T}^{(m)}\Phi$. In other words, the
coupling of the scalar to matter gives rise to an effective scalar mass which
depends on the matter content $\tilde{T}^{(m)}\Phi$. This is one key
ingredient of our study. Thus, direct couplings to curvature would also
produce similar effects Silva  _et al._ (2018); Antoniou  _et al._ (2018);
Doneva and Yazadjiev (2018); Ventagli  _et al._ (2020).

For concreteness, we focus exclusively on the following scalar potential and
coupling function,

$V(\Phi)=\frac{\mu_{0}^{2}}{2}\Phi^{2}\,,\qquad
A(\Phi)=e^{\frac{\beta}{2}\Phi^{2}}\,,$ (4)

but our results and methods can be applied to other models. Scalar self-
interactions can play an important role in the nonlinear development of the
instability; here we focus exclusively on the physics at small $\Phi$. Hence
our results describe the early-time dynamics of more general self-interacting
theories. Here, $\mu_{0}$ is the (bare) mass parameter of the scalar field
(note that the mass $m_{\rm s}$ of the field is related to the mass parameter
via $m_{\rm s}=\mu_{0}\,\hbar$ in these units). With this choice of scalar
potential and coupling function, a vanishing scalar field is a solution to the
theory and fluctuations around this solution endow the scalar with an
effective, time and position-dependent mass. Our analysis can also be
generalized to theories with nontrivial, stable scalar profiles, including
popular examples such as spontaneous scalarization Damour and Esposito-Farese
(1993, 1996); Harada (1998); Silva  _et al._ (2018); Antoniou  _et al._
(2018); Doneva and Yazadjiev (2018) and screening mechanisms Khoury and
Weltman (2004a, b); Joyce  _et al._ (2015); Hinterbichler  _et al._ (2011);
Hinterbichler and Khoury (2010).

Let us consider matter fields $\Psi_{m}$ describing a perfect fluid. Focus on
a geometry with vanishing scalar, $\Phi=0$, describing a static star of
(constant, for simplicity) density $\tilde{\rho}_{0}$ and radius $L_{0}$
Shapiro and Teukolsky (1983). The total mass of this solution can be written
as $M=\frac{4\pi}{3}\tilde{\rho}_{0}L_{0}^{3}$.

We ignore the effect of the scalar field on the profile of the fluid. In other
words, we deal only with the Klein-Gordon equation (3b) on this fixed
geometry. Perturbing around a background solution $\Phi_{0}(r)$, one finds
that the stability of such solutions depend dramatically on the value of the
coupling $\beta$ Damour and Esposito-Farese (1993, 1996); Harada (1998); Pani
_et al._ (2011). For concretness, here we take $\beta<0$ with
$\mu^{2}_{0}-\beta\tilde{T}^{(m)}\gg-L_{0}^{-2}$ and such that $\Phi_{0}=0$ is
a stable solution. The assumption of zero background scalar is a conservative
assumption, made to highlight the nontrivial effects of the mechanism we
discuss below, which provide a unique tool to constrain the theory. When the
background scalar is not zero other effects become important. For example, an
oscillating body with a nontrivial background scalar will radiate monopolar
radiation, a topic explored in the past Silvestri (2011); Sotani (2014);
Mendes and Ortiz (2018).

##  III Instabilities of oscillating astrophysical objects

Let us then consider the dynamical behavior of a (nonminimally coupled) scalar
field in a geometry describing a radially oscillating star. For simplicity,
the results shown below assume a Minkowski background and are thus valid only
for Newtonian stars. General Relativity effects change the quantitative but
not the qualitative behavior, and are discussed in the apendded supplemental
material, including an analysis of general-relativistic fluctuations of
compact stars Chandrasekhar (1964); Bardeen  _et al._ (1966); Kokkotas and
Ruoff (2001). Fluctuations around a background value of the scalar field are
described by the Klein-Gordon equation

$\square\Phi=\mu^{2}\Phi\,,$ (5)

with an effective mass

$\mu^{2}=\mu_{0}^{2}-\beta\tilde{T}^{(m)}\,,$

which acts as a position-dependent mass term that depends on the energy
density of the star. If the star oscillates (with a time-dependent radius
$L(t)$), so does the effective mass of the scalar field. For simplicity we
here assume the simple model

$\displaystyle L(t)$ $\displaystyle=L_{0}+\delta L\sin\omega t\,,$ (6)
$\displaystyle\mu^{2}(t,r)$
$\displaystyle=\mu_{0}^{2}+\beta\left(\tilde{\rho}_{0}+\delta\tilde{\rho}(t,r)\sin(\omega
t)\right)\,,$ (7)

where $\delta\tilde{\rho}(t,r)$ is the amplitude of the density oscillation,
and $\tilde{\rho}_{0}$ is the average of the energy density. We assume
$\beta<0$, in such a way that $\mu_{0}^{2}+\beta\tilde{\rho}_{0}\gtrsim 0$ to
avoid tachyonic instabilities. The competition between different mechanisms is
discussed in detail below and in the Supplemental Material, appended at the
end.

Here, we assumed $\delta L\ll L_{0}$ and
$\delta\tilde{\rho}\ll\tilde{\rho}_{0}$. This dependency has two important
aspects:

$\bullet$ Because the effective mass inside the star is smaller than the
exterior bare mass $\mu_{0}$, cf. Eq. (7), long wavelength modes are trapped
inside the star. The oscillating star surface therefore provides reflecting
boundary conditions at the periodically-varying location $L(t)$. For setups
where there is a continuous mass profile, a radially oscillating star
corresponds, nevertheless, to periodically varying field profiles at the
surface, possibly changing on scales smaller than a wavelength.

$\bullet$ The local density oscillations cause a time-varying effective mass
for the scalar, which lends itself to parametric instabilities as we show
below.

In astrophysically realistic stars, the time-dependence of the radius and
density of the star is quite involved. Normal modes of oscillation of a star
will cause a density profile with a sinusoidal-like radial profile as well. We
studied other profiles for the perturbation (in particular
$\delta\tilde{\rho}\propto\sin\omega r$) and fully general-relativistic
settings and find that they lead to similar results as described below and in
the Supplemental material.

We decompose the scalar field into spherical harmonics and evolve each mode,
$\Phi_{lm}(t,r)$, using a fourth-order accurate Runge-Kutta scheme for the
time integration where spatial derivatives are approximated by fourth-order
accurate finite difference stencils. Radiative boundary conditions are imposed
at the boundary of the computational domain, which is not in causal contact
with the star during our numerical evolutions.

We use time-symmetric initial data with a profile parameterized by

$\displaystyle\Phi_{lm}=e^{-\left(\frac{r-r_{0}}{\sigma}\right)^{2}}\,,\qquad\dot{\Phi}_{lm}=0\,,$
(8)

where $\sigma$ and $r_{0}$ denoting the width and initial position of the
scalar field pulse. Since the equation to be solved is linear we set the
initial amplitude of the pulse to unity. We focus our discussion on $l=0$
modes and we fix $r_{0}=0.5L_{0}$ and $\delta L=0.1L_{0}$. We studied
higher-$l$ modes, and found no qualitative difference. The results discussed
below all show fourth-order convergence.

##  IV Results

Consider setups where the effective mass inside [cf. Eq. (7)] the star
satisfies $\mu_{\rm in}^{2}<\mu_{0}^{2}$. Low-energy fluctuations of the
scalar are trapped inside the star and subject to conditions at the surface
which might be prone to blueshift instabilities – the growth of energy in the
field can be understood from a simple particle picture as a cumulative Doppler
effect Dittrich  _et al._ (1994a), which is discussed in greater depth in the
appended supplemental material. Such an analysis indicates the dominant
instability window to be Dittrich  _et al._ (1994a),

$\frac{\pi}{L_{0}+\delta L}<\omega<\frac{\pi}{L_{0}-\delta L}\,.$ (9)

Figure 1: Time evolution of a massive scalar field non-minimally coupled to
geometry, such that it is effectively massless inside the star ($\mu_{\rm
in}=0$) but with bare mass $\mu_{0}L_{0}=100$ and $\delta\tilde{\rho}=0$. The
field “sees” only an oscillating boundary with $\omega L_{0}=3.1$, we froze
density variations inside the star. The field is extracted at $r=0.5L_{0}$ and
the corresponding integrated energy inside the star is shown in the inset.
Here, the initial data is characterized by $\sigma=0.1L_{0}$.

In such cases, we see a transfer of energy from the oscillating star to the
scalar field, increasing the scalar frequency. We can see this behavior in
Figs. 1-2. Upon each reflection at the surface the scalar drifts to higher
frequencies, and after a sufficient amount of reflections it is no longer
confined: the field is finally able to leak to outside the star (and has
frequency $\omega\sim\mu(1+\epsilon),\,\epsilon\ll 1$). This mechanism causes
the energy to grow as $E\sim E_{0}e^{\lambda_{B}t}=E_{0}e^{2n_{\rm
ref}L_{0}\lambda_{B}}$, with $n_{\rm ref}$ the number of reflections and
$\lambda_{B}$ the instability growth rate. Here, $E_{0}\sim\sigma^{-1}$ is the
initial dominant spectral content at low energies, the ones that linger long
enough to be amplified. Thus, the field is able to leak away from the star and
out to spatial infinity when $e^{2n_{\rm
ref}L_{0}\lambda_{B}}\sigma^{-1}>\mu_{0}$. This occurs after a number of
reflections

$\displaystyle n_{\rm
ref}>n_{\ast}\sim\frac{\ln|\mu_{0}\sigma|}{2\lambda_{B}L_{0}}\,.$ (10)

At the critical number of reflections, the total energy confined inside the
star as scalar field is

$\frac{E}{E_{0}}\biggr{\rvert}_{\rm max}=\sigma\mu_{0}\,,$ (11)

this being also the total amount of energy extracted from the star.

Figure 2: Evolution of the integrated energy of a scalar inside an oscillating
star, for different bare mass parameters $\mu_{0}$. The coupling to matter is
such that the field is massless inside the star ($\mu_{\rm in}=0$). The field
“sees” only an oscillating boundary with $\omega L_{0}=3.1$, we froze density
variations inside the star. The initial data is characterized by
$\sigma=0.1L_{0}$.

We have thus far focused only on the effect of oscillating boundaries, while
the interior is non-dynamical. However, real objects will also have an
oscillating density. Thus, the effective mass inside the star will have a
periodic time variability, potentially giving rise to “parametric”
instabilities. Figure 3 summarizes our results for non-vanishing
$\beta\delta\tilde{\rho}$. Two mechanisms now compete, a blueshift and a
parametric mechanism. For small $L_{0}^{2}\delta\tilde{\rho}$, the evolution
of the scalar field is almost identical to that with a vanishing
$\delta\tilde{\rho}L_{0}^{2}$, and the blueshift mechanism dominates. For
large $\delta\tilde{\rho}L_{0}^{2}$, we observe instead that the amplitude of
the scalar field grows while its frequency is barely changing, a clear sign
that we are dealing with a different process. This is in fact a parametric
instability. Contrast this with the blueshift instability cases, where the
scalar field pulse becomes narrower (due to the oscillating boundary) as time
passes, but its amplitude remains constant. Here, the width of the scalar
pulse is roughly constant and the instability is driven by a growth in the
amplitude of the field (due to the oscillating effective density).

Figure 3: Evolution of the energy and scalar field for setups which probe both
an oscillating boundary and a time-varying effective mass inside the star. The
initial data has width $\sigma=0.1L_{0}$, and the star is oscillating with
frequency $\omega L_{0}=3.1$. The (background) effective mass inside the star
is zero, $\mu_{\rm in}^{2}=\mu_{0}^{2}+\beta\tilde{\rho}_{0}\approx 0$, while
the bare mass is $\mu_{0}L_{0}=100$. The time-varying component of the density
is proportional to $\delta\mu^{2}\equiv\beta\delta\tilde{\rho}$. Even though
both the blueshift and parametric instabilities are present, the parametric
mechanism dominates at large $\delta\mu^{2}$.

##  V Application to astrophysical systems

We have shown that pulsating stars or other objects may excite important
instabilities of nonminimally coupled fields. Compact stars, such as neutron
stars, have radial pulsations with frequencies satisfying $\omega L_{0}\sim 1$
for the lowest overtones Kokkotas and Ruoff (2001), and are ideal systems
where such instabilities might be relevant. These change the local density and
distribution of the scalar (which could be one dark matter component) and may
even backreact on the star. A precise description of the evolution of the
instability requires a precise knowledge of stellar oscillations and careful
modeling of the evolution of the scalar in such backgrounds. This is a
challenging program that requires further study.

Notice first that one can apply our results for stars as long as other
dissipation mechanisms are subdominant. Of particular importance are shear
viscosity effects, which in neutron stars have a timescale Cutler and Lindblom
(1987); Barta (2019)

$\tau_{\eta}\approx 100\,\rho^{-5/4}_{17}T_{5}^{2}\left(\frac{L_{0}}{\rm
10km}\right)^{2}\,{\rm s}\,,$ (12)

where $\rho_{17}=\rho/(10^{17}{\rm kg/m^{3}})$, $T_{5}=T/(10^{5}\,{\rm K})$
and $T$ is the neutron star temperature.

The blueshift instability is quenched after

$\displaystyle t_{B}\sim 2L_{0}n_{\ast}$ $\displaystyle=3\times 10^{-2}{\rm
s}\left(\frac{10^{-3}}{\lambda_{B}(\omega)L_{0}}\right)\left(\frac{L_{0}}{10{\rm
km}}\right)$ $\displaystyle\quad{}\times\ln\left|10\frac{\mu}{10^{-10}{\rm
eV}}\frac{\sigma}{10{\rm km}}\right|\,,$ (13)

and on this timescale an energy $E\sim E_{0}\sigma\mu_{0}$ is removed from the
star. This result is not very sensitive to the initial conditions. It is, in
principle, only weakly affected by backreaction, unless $\sigma\mu_{0}$ is an
extremely large number. Our numerical results indicate that
$\lambda_{B}L_{0}\sim 10^{-3}$ is a reasonable estimate for $\delta
L>0.01L_{0}$. The instability window for the blueshift mechanism to work is
tight, however [cf. Eq. (9)]. Only large overtones are affected by it, unless
the star is oscillating nonlinearly with $\delta L\gtrsim 0.1L_{0}$.

Consider now the parametric instability. When
$\beta\sim-\frac{\mu^{2}}{\tilde{\rho}_{0}}$, with $\tilde{\rho}_{0}$ the
temporal average of the density of the star, the relevant dynamics is governed
by the Mathieu equation Abramowitz and Stegun (1972); Bender and Orszag
(1978). This particular case provides a test on our results, and allows to
estimate analytically the timescales involved. For small $\delta\mu^{2}$, the
instability condition amounts to $L_{0}\omega=4\pi/j$ ($j\in\mathbb{Z}$). The
instability rate time scale for $j=1$, for example, is roughly $t_{A}\sim
2\omega/\delta\mu^{2}$, or

$t_{A}\sim\frac{2\omega}{\delta\mu^{2}}\sim 1{\rm s}\left(\frac{10{\rm
km}}{L_{0}}\right)\left(\frac{\omega
L_{0}}{4\pi}\right)\left(\frac{10^{-2}\rho_{17}}{\delta\tilde{\rho}}\right)\left(\frac{-10}{\beta}\right)\,,$

for small $\delta\mu^{2}/\omega^{2}$. These estimates assume that the field is
effectively very light inside the star, which amounts to requiring that

$|\beta|\sim
7\left(\frac{0.3}{\mathcal{C}}\right)\left(\frac{\mu}{10^{-10}{\rm
eV}}\right)^{2}\left(\frac{L_{0}}{10{\rm km}}\right)^{2}\,,$ (14)

but the instabilities discussed here are expected to set in even at small
nonzero effective masses. The star oscillations can be induced by accretion,
tidal effects or even mergers Reisenegger and Goldreich (1994); Chirenti  _et
al._ (2017); Ma  _et al._ (2020). Note that numerical relativity simulations
show that the amplitude of density perturbations during coalesce, for example,
can be significant and of the order of the central density itself Perego  _et
al._ (2019); Bernuzzi  _et al._ (2020). Thus, both mechanisms may act on
timescales short enough to be relevant. In fact, at large couplings $\beta$ –
not yet ruled out for large bare masses – parametric instabilities will be
dominant.

##  VI Final remarks

Nonminimally coupled, massive scalar fields can evade all observational
constraints from the observation of nearly stationary configurations, yet
produce distinct signatures when evolving around oscillating backgrounds. We
have shown that there are at least two possible instability mechanisms: one
blueshifts light scalars inside oscillating stars 111there were important
hints that such an effect could occur. The classical Fermi acceleration
process is one example of energy extraction from a zero-average motion Fermi
(1949). Similar processes exist for fields trapped inside cavities with
periodically-moving boundaries Cooper (1993); Dittrich  _et al._ (1994b);
Moore (1970); Andreata and Dodonov (2000); Dodonov (2001). We showed that
oscillating stars are also prone to such instabilities.; the second mechanism
is of parametric origin, triggered by a periodic oscillation of the star
material. Possible excitation mechanisms for the scalar field could include,
for instance, coherently or stochastically oscillating dark matter.

Both instabilities act on short timescales when compared to viscous timescales
Cutler and Lindblom (1987), and are expected to play a role in neutron star
oscillations. They can backreact on the star – perhaps leading to
gravitational collapse, or (in the blueshift mechanism) simply result in a
leakage of high frequency, high amplitude scalar. If a fraction, or all of
dark matter is made of a scalar component, then these mechanisms can act to
produce overdensities close to neutron stars, providing one more route to
constraining dark matter. Details on observational signatures of these
instabilities require further studies, beyond the scope of this work.

Similar instabilities are expected not only in scalar-tensor theories, but
also in other theories with vectors or spinors Annulli  _et al._ (2019);
Minamitsuji (2020a, b). Scalarized background solutions are not prone to the
type of instabilities discussed here, but star oscillations will lead to
radiation emission Silvestri (2011). We expect that the instabilities
discussed here could describe and affect other systems where a light degree of
freedom is confined to an oscillating background, potentially observable in
condensed matter systems.

Finally, other periodic systems include compact binaries; it can be expected
that similar blueshift mechanisms act on such binaries. Their high degree of
asymmetry makes it more challenging to model, but their immense gravitational
potential energy certainly makes them important candidates.

###### Acknowledgements.

We are grateful to Sebastiano Bernuzzi, Kostas Kokkotas, Shingo Kukita,
Tomohiro Nakamura, and Pantelis Pnigouras for useful correspondence and
advice. V. C. is a Villum Investigator supported by VILLUM FONDEN (grant no.
37766). V.C. acknowledges financial support provided under the European
Union’s H2020 ERC Consolidator Grant “Matter and strong-field gravity: New
frontiers in Einstein’s theory” grant agreement no. MaGRaTh–646597. M.Z.
acknowledges financial support provided by FCT/Portugal through the IF
programme, grant IF/00729/2015. T.I. acknowledges financial support provided
under the European Union’s H2020 ERC, Starting Grant agreement no.
DarkGRA–757480. This project has received funding from the European Union’s
Horizon 2020 research and innovation programme under the Marie Sklodowska-
Curie grant agreement No 690904. We thank FCT for financial support through
Project No. UIDB/00099/2020 and through grants PTDC/MAT-APL/30043/2017 and
PTDC/FIS-AST/7002/2020.

Supplemental Material

We now give more details and expand on the construction outlined in the main
text.

##  VII Setup: nonminimal scalar fields

###  VII.1 The background

We consider matter fields $\Psi_{m}$ describing a perfect fluid. Focus on a
geometry with vanishing scalar, $\Phi=0$, describing a static star of
(constant) density $\tilde{\rho}_{0}$ and radius $L_{0}$. The ADM mass of this
solution can be written as $M=\frac{4\pi}{3}\tilde{\rho}_{0}L_{0}^{3}$ and its
geometry can be expressed as Shapiro and Teukolsky (1983)

$\tilde{g}_{\mu\nu}dx^{\mu}dx^{\nu}=-\tilde{\alpha}^{2}dt^{2}+\tilde{a}^{2}dr^{2}+r^{2}d\Omega_{2}\,,$
(15)

where

$\displaystyle\tilde{\alpha}^{2}=\begin{cases}\left(\frac{3}{2}\left(1-\frac{2M}{L_{0}}\right)^{1/2}-\frac{1}{2}\left(1-\frac{2Mr^{2}}{L_{0}^{3}}\right)^{1/2}\right)^{2}&(r<L_{0})\\\
1-\frac{2M}{r}&(r>L_{0})\end{cases}$

and

$\displaystyle\tilde{a}^{2}=\begin{cases}\left(1-\frac{2Mr^{2}}{L_{0}^{3}}\right)^{-1}&(r<L_{0})\\\
\left(1-\frac{2M}{r}\right)^{-1}&(r>L_{0})\end{cases}\,.$

The profile of the pressure inside the star is

$\displaystyle\tilde{p}_{0}(r)$ $\displaystyle=$
$\displaystyle\tilde{\rho}_{0}\left(\frac{(1-\mathcal{C}\frac{r^{2}}{L_{0}^{2}})^{1/2}-(1-\mathcal{C})^{1/2}}{3(1-\mathcal{C})^{1/2}-(1-\mathcal{C}\frac{r^{2}}{L_{0}^{2}})^{1/2}}\right)\,.$
(16)

We can now linearize the equations of motion around this background. Although
we focus solely on constant-density stars, previous results indicate that the
details of the model are not important Pani  _et al._ (2011). The most
relevant feature is the contribution of the star (trace of) energy tensor
$\tilde{T}=3\tilde{p}-\tilde{\rho}$ to the effective scalar mass, and this
contribution is qualitatively the same for different equations of state.

There are situations when the background scalar is nontrivial, as we explain
below. In such circumstances, it backreacts on the geometry. Although the
formalism to handle backreaction is clear, for simplicity here we ignore the
effect of the scalar field on the profile of the fluid. In other words, we
deal only with the Klein-Gordon equation on a fixed geometry. As a concrete
situation, we consider relativistic stars of constant density as a solution,
parameterized by radius $L_{0}$ and compactness $\mathcal{C}$.

###  VII.2 The dynamical equations

In Sec. X, we deal with full dynamical spacetimes on which the scalar evolves.
For now, we consider simply a time-independent spacetime, since it serves to
convey our main message. Let $\Phi_{0}$ be a static, spherically-symmetric
solution of Eq. (3b), with the metric (15) and scalar potential and coupling
functions given by Eq. (4). $\Phi_{0}$ is then a solution to the equation

$\displaystyle\partial_{r}^{2}\Phi_{0}+\left(\frac{2}{r}+\frac{\partial_{r}\alpha}{\alpha_{0}}-\frac{\partial_{r}a}{a}\right)\partial_{r}\Phi_{0}$
$\displaystyle\qquad+a^{2}\left(\beta\tilde{T}^{(m)}e^{2\beta\Phi_{0}^{2}}\Phi_{0}-\mu_{0}^{2}\Phi_{0}\right)=0\,,$
(17)

where $\alpha$ and $a$ are trivially related to the functions $\tilde{a}$ and
$\tilde{\alpha}$ defined above via

$(\tilde{a},\tilde{\alpha})=(a,\alpha)e^{\frac{\beta}{2}\Phi_{0}^{2}}\,.$

Let us now consider a general fluctuation around this value,

$\displaystyle\Phi(t,r,\theta,\phi)=\Phi_{0}(r)+\epsilon\sum_{l,m}\frac{\psi_{lm}(t,r)}{r}Y_{lm}(\theta,\phi)\,,$
(18)

where $\epsilon$ is a bookkeeping parameter and $Y_{lm}(\theta,\phi)$ are the
usual scalar spherical harmonics. Then, keeping only first order terms in
$\epsilon$, we find from Eq. (3b) that $\psi_{lm}$ is governed by

$-\partial_{t}^{2}\psi-\frac{\alpha^{2}}{a^{2}}\partial_{r}\left(\log{\left(\frac{a}{\alpha}\right)}\right)\partial_{r}\psi+\frac{\alpha^{2}}{a^{2}}\partial_{r}^{2}\psi-
V_{\rm eff}\psi=0\,,$ (19)

where we dropped the $lm$ subscripts for clarity and the effective potential
$V_{\rm eff}$ is given by

$\displaystyle V_{\rm eff}$
$\displaystyle=\frac{\alpha^{2}}{ra^{2}}\left(\frac{\partial_{r}\alpha}{\alpha}-\frac{\partial_{r}a}{a}\right)+\alpha^{2}\frac{l(l+1)}{r^{2}}$
$\displaystyle\quad+\alpha^{2}\left(\mu_{0}^{2}-e^{2\beta\Phi_{0}^{2}}\beta\tilde{T}^{(m)}(1+4\beta\Phi_{0}^{2})\right)\,.$
(20)

The stable background solutions $\Phi_{0}(r)$ depend dramatically on the value
of the coupling $\beta$ Damour and Esposito-Farese (1993, 1996); Harada
(1998); Pani  _et al._ (2011). There is a threshold coupling $\beta_{c}$ that
divides cases where $\Phi_{0}$ is a stable solutions from those where it is
not. For $\beta>\beta_{c}(<0)$, $\Phi_{0}=0$ is a stable solution and the
geometry of the star and all its properties are identical to those in GR. If
$\beta<\beta_{c}$, the trivial scalar profile becomes unstable. In such
conditions, a GR star quickly develops a nontrivial profile $\Phi_{0}$ and is
said, for obvious reasons, to be _scalarized_. In this case, the geometry and
properties of the star change, and it would no longer be described by Eq.
(16), for example. All our techniques below can nonetheless be extended to
scalarized objects, or directly applied if one’s focus is on configurations
where $|\Phi_{0}|\ll 1$ (also known as the decoupling limit).

These scalar-tensor theories are appealing because Damour and Esposito-Farese
(1993, 1996); Harada (1998); Pani  _et al._ (2011); Cardoso  _et al._ (2020)

$\beta_{c}\propto-\left(1+(\mu_{0}R)^{2}\right)\mathcal{C}^{-1}\,.$ (21)

Thus, they predict scalarization in strong field situations, for example in
neutron stars, but not in stars like our Sun. Thus, these theories satisfy
solar-system bounds in a natural way. On the other hand, several observational
predictions regarding neutron stars can be different from GR, and thus
interesting constraints can be imposed. The observation of binary neutron
stars can, in principle, be used to constrain scalarization since scalarized
stars would radiate more energy (thus leading to a more rapid inspiral, as
compared with the standard general relativistic solutions) Damour and
Esposito-Farese (1993, 1996); Harada (1998); Pani  _et al._ (2011); Antoniadis
_et al._ (2013). However, low energy excitations are unable to propagate if
the field is massive. Thus, for fields with

$\displaystyle\mu_{0}\gg 10^{-16}{\rm eV}\,,$ (22)

it is challenging to impose constraints on this mechanism via neutron star
observations Ramazanoğlu and Pretorius (2016).

Figure 4: (Left panel) Effective potential for spherically-symmetric scalar
fluctuations around a relativistic uniform density star with $\mathcal{C}=0.3$
and a trivial background scalar $\Phi_{0}=0$ (with a coupling
$\beta>\beta_{c}$). (Right panel) Same as in the left panel, but now for a
subcritical coupling $\beta<\beta_{c}$, leading to scalarization of the star.
The behavior of the blue and purple lines is hard to see, but their profile is
similar to that of the red and green lines. The corresponding scalar profile
is shown in the inset.

To understand the possible dynamics of the scalar field, it is important to
know the structure of the effective potential $V_{\rm eff}$. This potential is
shown in Fig. 4. The importance of the effective potential is most easily seen
from Eq. (19). For fluctuations of energy $\Omega$,

$\psi\sim e^{-i\Omega t}\Upsilon(r)\,,$ (23)

it takes the form of a Schrödinger-like equation (notice that $\Omega$ has
dimensions of ${\rm mass}^{-1}$ in geometrical units, but is an energy
parameter)

$\frac{\alpha^{2}}{a^{2}}\partial_{r}^{2}\Upsilon-\frac{\alpha^{2}}{a^{2}}\partial_{r}\left(\log{\left(\frac{a}{\alpha}\right)}\right)\partial_{r}\Upsilon+(\Omega^{2}-V_{\rm
eff})\Upsilon=0\,,$ (24)

Thus for small enough energies $\Omega^{2}<V_{\rm eff}$ the fluctuation is
unable to tunnel to infinity if $V_{\rm eff}$ asymptotes to a constant
positive value.

When $\beta<\beta^{c}$, the object is scalarized with a nontrivial profile
$\Phi_{0}\neq 0$ (Fig. 4, right panel). Notice from the shape of the effective
potential that bound states are hard to achieve inside the object in these
cases, since they can easily “tunnel out” (in the language of quantum
mechanics associated with the Schrödinger-like equation (24)).

For non-scalarized solutions, with $\Phi_{0}=0$, Eq. (20) shows that
fluctuations of the scalar are massive with an effective mass

$\mu^{2}=\mu_{0}^{2}-\beta\tilde{T}^{(m)}\,,$ (25)

which effectively acts as a position-dependent mass term. Consider setups
where the effective mass inside the star satisfies $\mu^{2}<\mu_{0}^{2}$ (cf.
Eq. (24) and the left panel of Fig. 4). Then, low-energy fluctuations of the
scalar are trapped inside the star. This aspect plays a critical role in our
analysis and results. Note that for extremely compact objects the trace of the
stress tensor can change sign, and positive couplings may give rise to
spontaneous scalarization as well Mendes (2015). Our results can in principle
be applied to all these different scenarios.

##  VIII Periodic-motion instabilities in one dimension

As we discussed, some classes of scalar-tensor theories are described by a
scalar field with an effective mass that depends on the environment. In
particular, the effective mass $\mu$ of the scalar field inside a star, for
example, can be lighter than its vacuum counterpart $\mu_{0}$. The potential
experienced by the field looks like those depicted in Fig. 4 and bound states
inside the star can appear. Due to this potential barrier, field excitations
with wavelength larger than $\mu^{-1}$ are reflected around the surface of the
object. Thus, the field propagates within a cavity whose size is varying
periodically, and whose density is also varying periodically. We find that
these two aspects give rise two instabilities of different kinds. A varying
cavity size allows for the possibility of energy extraction from the cavity
via “Doppler-like” exchanges, whereas a varying density gives rise to
parametric instabilities.

We describe in details these two types of instability below. Before going into
the more complex case of four-dimensional spacetimes, we consider a simple 1+1
setup in this section.

###  VIII.1 The blueshift instability

Figure 5: A scalar field confined inside a cavity of size $L(t)$. The left
boundary is kept fix and the right boundary (periodically) oscillates in time.

We decouple the two effects mentioned above and consider first a scalar field
of constant mass, confined to a one-dimensional oscillating cavity in flat
space – see Fig. 5. In particular, one cavity boundary is fixed at $x=0$ while
the other is located at $L(t)$, as described by Eq. (6). The scalar field is
then governed by the massive Klein-Gordon equation with Dirichlet boundary
conditions imposed at a time-dependent boundary,

$\displaystyle-\partial_{t}^{2}\Phi+\partial_{x}^{2}\Phi-\mu^{2}\Phi=0\,,$
(26) $\displaystyle\Phi(t,0)=\Phi(t,L(t))=0\,,$

where $\mu$ is the (effective) mass of the scalar field inside the cavity.
This models an oscillating star with a scalar field whose effective mass is
small inside the star and large outside. The system has been investigated
previously in other contexts, and shown to give rise to blue or redshifted
reflections Cooper (1993); Dittrich  _et al._ (1994b). After a large enough
number of reflections, the energy of the field increases without bound for

$\displaystyle\frac{N\pi}{L_{0}+\delta L}<\omega<\frac{N\pi}{L_{0}-\delta
L}\,,$ (27)

where $N$ is an integer Dittrich  _et al._ (1994b). We call this the blueshift
condition and the associated instability the blueshift instability. In the
following we study the mechanism and associated instability in detail.

####  VIII.1.1 Estimate from Doppler effect

Figure 6:  Time evolution of energy of a photon inside a box with an
oscillating boundary, solution of recursion relation (29). The amplitude of
oscillation of the boundary is $\delta L=0.1L_{0}$, the initial position
$x_{\rm init}$ takes different selected values. For oscillation frequencies
$\omega L_{0}=2.8$ (top panel), the energy oscillates, which means that
there’s energy transfer between the photon and the boundary, and they reached
equilibrium. This also means that the blueshift instability condition (27) is
not satisfied. For oscillation frequencies that do satisfy that condition, the
boundary is constantly giving energy to the photon, resulting in an
exponential growth of its energy, as seen in the bottom panel.

The growth of energy in the field can be understood as a cumulative Doppler
effect Dittrich  _et al._ (1994a). To understand this in more depth and from a
different perspective, we consider the following toy model. A massless
particle is inside a perfectly reflecting cavity whose boundary oscillates,
and is emitted at $x=x_{\rm init}$ moving to the right with initial (angular)
frequency $\Omega_{0}$. After one reflection at the (moving) boundary, the
frequency is Doppler-shifted to $\Omega_{1}$,

$\Omega_{1}=\frac{1-v(t_{1})}{1+v(t_{1})}\Omega_{0}\,,$ (28)

where $v(t)=\dot{L}=\delta L\,\omega\cos\omega t$ is the velocity of the
oscillating boundary, and $t_{1}$ the coordinate time when the particle hits
the boundary determined from the equation $t_{1}+x_{\rm init}=L(t_{1})$. Using
the same argument, one can obtain a recursion relation for the frequency
$\Omega_{n_{\rm ref}}$, and therefore for the energy $E_{n_{\rm ref}}$ after
$n_{\rm ref}$ reflections

$\Omega_{n}=\frac{1-v(t_{n_{\rm ref}})}{1+v(t_{n_{\rm ref}})}\Omega_{n_{\rm
ref}-1}\,,$ (29)

where $t_{n_{\rm ref}}$ is determined by

$t_{n_{\rm ref}}-t_{n_{\rm ref}-1}-L(t_{n_{\rm ref}-1})=L(t_{n_{\rm ref}})\,.$
(30)

The recursion relation (29) can be solved numerically, its solution is shown
in Fig. 6 for different choices of parameters fixing $\delta L=0.1L_{0}$. Note
that, for this choice of $\delta L$, condition (27) becomes $2.85<\omega
L_{0}/N<3.49$. We thus choose two representative values for the oscillation
frequency $\omega$, one outside the instability interval and one inside. The
top panel of the figure shows that when the instability criterion is not
satisfied the energy of the particle oscillates in time. In some interactions
with the boundary it gains energy, in some it loses energy, the average gain
is zero. The bottom panel shows that there are conditions that lead to a
constant gain of energy by the particle, and where each reflection is
accompanied by a transfer of energy from the boundary to the bouncing massless
particle. In particular, for certain conditions, the energy grows
exponentially,

$E\propto e^{\lambda_{B}t}\,,$ (31)

and the exponent $\lambda_{B}>0$ is insensitive (or very mildly dependent
only) on the initial position $x_{\rm init}$.

Figure 7: Relation between the instability rate $\lambda$ and the frequency of
the boundary for each $\delta L$.

Figure 7 shows the relation between the instability rate, $\lambda_{B}$, and
the frequency of oscillation of the boundary $\omega$. An approximate relation
for the lowest-frequency instability window in Fig. 7 is

$\lambda_{B}L_{0}\sim a(\omega)\left(\frac{\delta
L}{L_{0}}-\left|\frac{\pi}{\omega L_{0}}-1\right|\right)\,,$ (32)

valid in the interval $\pi/(L_{0}+\delta L)<\omega<\pi/(L_{0}-\delta L)$. The
coefficient $a(\omega)$ has a nontrivial dependence on the boundary
oscillation frequency and is $a\sim 3.5$ for $\omega L_{0}=\pi$, and $a\sim
4.6$ for $\omega L_{0}=0.9\pi$.

####  VIII.1.2 An exact solution

The results above are based on a simple particle picture, compelling but valid
for massless and localized fields only. We now wish to show that this behavior
arises also when solving the wave equation (26). In fact, an exact solution
for certain prescribed motion can be found for massless scalars subjected to
Dirichlet conditions at the boundaries Law (1994). The boundaries are located
at

$x=\\{0,L(t)\\}\,,$ (33)

but – here and only here – we will sacrifice the simple boundary motion (6) in
order to find an exact solution to the problem. In one dimension, the solution
of the wave equation takes the form

$\Phi=e^{-ik\pi{\cal R}(t+x)}-e^{-ik\pi{\cal R}(t-x)}\,,$ (34)

which ensures that the massless KG equation is satisfied. Here, $k$ is an
arbitrary constant and $k\pi{\cal R}$ is the frequency $\Omega$ of the pulse.
To satisfy the BCs one takes

${\cal R}(t+L(t))={\cal R}(t-L(t))+2\,,$ (35)

The solution to this recursion relation solves the problem. For prescribed
motion,

$\displaystyle L(t)$ $\displaystyle=$ $\displaystyle
L_{0}+\frac{L_{0}}{2\pi}\left(\arcsin{\left(\sin\theta\cos\frac{2\pi
t}{L_{0}}\right)}-\theta\right)\,,$ (36) $\displaystyle\theta$
$\displaystyle=$ $\displaystyle\arctan\frac{\epsilon\pi}{L_{0}}\,,$ (37)

where $\epsilon$ describes the amplitude of the motion, the solution can be
found as

${\cal R}(2n_{\rm ref}L_{0}+\zeta)=2n_{\rm
ref}+\frac{1}{2}-\frac{1}{\pi}\arctan\left(\cot\frac{\zeta\pi}{L_{0}}-\frac{2n_{\rm
ref}\epsilon\pi}{L_{0}}\right)\,,$ (38)

where $\zeta$ is a variable in $[-L_{0},L_{0}]$ and $n\geq 1$ is a positive
integer. The branch of $\arctan$ must be selected carefully in order to avoid
a discontinuity. Notice that for small $\epsilon$,

$\frac{L(t)}{L_{0}}=1-\frac{\epsilon}{L_{0}}\sin^{2}\frac{\pi t}{L_{0}}+{\cal
O}(\epsilon/L_{0})^{3}$ (39)

It is easy to show that the system is unstable, precisely because the
frequency content is increasing Law (1994). This construction can be
generalized to a 3+1 setup.

####  VIII.1.3 Time evolutions

The results above are either based on a simple particle picture, or arise
within a very special and prescribed boundary motion. We now wish to show that
this behavior is a generic consequence also when solving the (massive) wave
equation (26). To numerically evolve this equation, we first start by
introducing new coordinates $(T,X)$ defined by

$T(t,x)=t\,,\qquad X(t,x)=\frac{x}{L(t)}\,.$ (40)

In these coordinates, the equation of motion for the scalar field $\phi(T,X)$
take the form

$\displaystyle-\partial_{T}^{2}\Phi+\left(\frac{1}{L^{2}}-\frac{X^{2}\dot{L}^{2}}{L^{2}}\right)\partial_{X}^{2}\Phi+X\left(\frac{\ddot{L}}{L}-\frac{2\dot{L}^{2}}{L^{2}}\right)\partial_{X}\Phi$
$\displaystyle\qquad+2X\frac{\dot{L}}{L}\partial_{T}\partial_{X}\Phi-\mu^{2}\Phi=0\,,$
(41)

Note that we interchangeably use both dots and $\partial_{T}$ to denote time
derivatives with respect to $T$. The boundary conditions (now at _constant_
spatial coordinate $X$) are

$\displaystyle\Phi(T,X=0)=\Phi(T,X=1)=0\,.$ (42)

In the new coordinates, the effect of the moving boundary appears as a shift
vector term. To numerically solve Eq. (41), we use the characteristic
variables $\Phi_{\pm}$ defined by

$\Phi_{\pm}\equiv\frac{1}{2}\left(\dot{\Phi}+\lambda_{\mp}\partial_{X}\Phi\right)\,,\quad\lambda_{\pm}=\pm\frac{1}{L}-X\frac{\dot{L}}{L}\,.$
(43)

The equations of motion are then

$\displaystyle\dot{\Phi}_{\pm}$
$\displaystyle=-\lambda_{\pm}\partial_{X}\Phi_{\pm}-\frac{\mu^{2}}{2}\Phi\,,$
(44a) $\displaystyle\dot{\Phi}$
$\displaystyle=L\left(\lambda_{+}\Phi_{+}-\lambda_{-}\Phi_{-}\right)\,,$ (44b)

with the constraint

$\frac{\partial_{X}\Phi}{L}=-\Phi_{+}+\Phi_{-}\,.$ (45)

In these variables, the boundary conditions become

$\displaystyle\Phi_{+}(T,0)$
$\displaystyle=\frac{\lambda_{-}(T)}{\lambda_{+}(T)}\Phi_{-}(T,0)\,,$ (46a)
$\displaystyle\Phi_{-}(T,1)$
$\displaystyle=\frac{\lambda_{+}(T)}{\lambda_{-}(T)}\Phi_{+}(T,1)\,.$ (46b)

Finally, the energy of the scalar field inside the cavity is given by

$\displaystyle E$
$\displaystyle=\frac{1}{2}\int_{0}^{L(t)}\left((\partial_{t}\Phi)^{2}+(\partial_{x}\Phi)^{2}+\mu^{2}\Phi^{2}\right)dx$
$\displaystyle=\int_{0}^{1}LdX\left\\{(\Phi_{+}^{2}+\Phi_{-}^{2})+\frac{\mu^{2}}{2}\Phi^{2}\right\\}\,.$
(47)

Figure 8: Convergence study of the evolution of the scalar field at fixed
radius for the simulation of Fig. 9. The purple line shows the difference
between results obtained with low ($4h$) and medium ($2h$) resolutions, while
the (dashed) green line shows the difference between results obtained with
medium and high ($h$) resolutions multiplied by 16, the expected factor for
fourth-order convergence.

We evolve Eqs. (44) with boundary conditions (46) using a finite differencing
scheme where spatial derivatives are approximated with 4th-order accurate
upwind stencils and a 4th-order accurate Runge Kutta scheme is used for the
time integration. We have checked that this code convergences with the
expected fourth-order accuracy, as shown in Fig. 8.

In all our simulations we use time-symmetric initial data, where the scalar
field is parameterized by

$\displaystyle\Phi(0,X)=e^{-\left(\frac{r-r_{0}}{\sigma}\right)^{2}}\,,\quad\dot{\Phi}(0,X)=\frac{\dot{L}}{L}X\partial_{X}\Phi\,,$
(48)

where $\sigma$ and $r_{0}$ denote the width, and initial position of the
scalar field pulse. Note that the wave equation is linear, hence the amplitude
can arbitrarily be set to unity. In this work we focus on
$\sigma=0.2L_{0},\,r_{0}=0.5L_{0}$ and $\delta L=0.1L_{0}$. Our results are
summarized in Figs. 9-12. We surveyed other regions in parameter space and
found no qualitative new features with respect to those discussed below.

Figure 9: Evolution of a gaussian pulse of scalar field and its energy
(inset). The field is massless, but confined to a cavity whose right boundary
oscillates at $\omega L_{0}=3.1$ with an amplitude $\delta L=0.1L_{0}$. The
gaussian pulse has width $\sigma=0.2L_{0}$ and is centered at
$r_{0}=0.5L_{0}$. The green dotted line in the inset is the exponential
$\exp{(0.31t/L_{0})}$.

Consider first massless fields, $\mu=0$, inside a cavity with an oscillating
boundary. One typical evolution is shown in Fig. 9. This figure is very clear:
the scalar field evolves in time with a negligible variation in amplitude, but
with a significant evolution of frequency content. Note that each reflection
is visible in the figure. The scalar field becomes “sharper” upon each
reflection, signaling a frequency increase. The inset shows the associated
increase in the total energy inside the cavity. The energy increase is
exponential and the exponent is consistent with the simple estimate from the
blueshift condition (27). Thus, the picture outlined in Sec. VIII.1.1 seems to
describe evolutions of the wave equation. For these parameters, the energy
content inside the cavity grows exponentially, driven by the blueshift
conversion. We also verified that the growth timescale $\lambda$ does not
strongly depend on the initial position $r_{0}$ of the pulse, consistent with
the particle picture (cf. Fig. 6).

Figure 10: Time evolution of the total energy inside the cavity for different
boundary oscillation frequencies $\omega L_{0}$. Notice that the purple line,
which does not satisfy the instability condition, remains bounded.

This interpretation is made stronger with the analysis of different
parameters. Figure 10 shows us that the timescale is sensitive to the boundary
oscillation frequency $\omega L_{0}$ and that at small frequencies the
instability disappears. In fact, our numerical results are very consistent
with the toy model of Sec. VIII.1.1, including the threshold frequency
prediction. For frequencies outside the instability window, our results show a
clear oscillating pattern, again consistent with the previous discussion. To
summarize, there is a clear parallel between our toy model of a massless
particle in Sec. VIII.1.1, and the evolution of the massless wave equation. In
particular, there are regions in parameter space where the energy in the
scalar field inside the cavity grows exponentially with time. For other
parameters, the energy oscillates in time, providing further support to the
interpretation in terms of blue and redshift processes.

Figure 11: Time evolution of the total energy for different values of $\mu
L_{0}$. Observe that the energy remains bounded in the top panel (where the
blueshift condition is not satisfied) whereas in the bottom panel it increases
without bound (where the blueshift condition is satisfied).

Thus far, we studied only massless fields. Figure 11 shows the evolution of
the energy in the cavity for different field mass $\mu$ and different
frequencies $\omega$. Notice that one could expect qualitatively new features
depending on whether the Compton wavelength of the field is larger or smaller
than the cavity size. Surprisingly, our results show clearly that the
instability rate and instability window is only weakly dependent on $\mu$,
even for large $\mu\,L_{0}$.

Figure 12: Time evolution of the scalar field and energy, for a scalar
confined to a cavity of length $L_{0}$. The scalar is massive, with
$\mu^{2}<0$, and the cavity boundary oscillates with $\omega L_{0}=2.8$ (top
panel) and $\omega L_{0}=3.1$ (bottom panel). For small $|L_{0}\mu|$ we
recover the previous results for massless fields; in particular a blueshift
instability at certain frequency ranges (bottom panel) and stability outside
this window. However, at negative large $(L_{0}\mu)^{2}$ a new mechanism sets
in, that of scalarization governed by a tachyonic mode Damour and Esposito-
Farese (1993); Cardoso  _et al._ (2020). (Notice that the tachyonic mode – the
green line – grows exponentially, the scale just doesn’t allow to see any
further high-amplitude oscillation).

Because we are interested also in scalarization, for which the effective mass
parameter $\mu^{2}_{\rm eff}$ can be negative, we studied this case as well,
summarized in Fig. 12. Now, another effect can set in, which affects even
stationary geometries. Namely, scalarization induced by a tachyonic mode
Damour and Esposito-Farese (1993); Cardoso  _et al._ (2020). The tachyonic
instability is described by a rate $L_{0}\lambda_{\rm tachyon}\sim\sqrt{|\mu
L_{0}|^{2}-1}$; thus at large $|L_{0}\mu|$ it can dominate the entire physics.

This is indeed apparent in Fig. 12. The upper panel summarizes the evolution
outside the blueshift instability window. For small $|L_{0}\mu|$ our results
indicate indeed just an oscillatory behavior for the total energy. At large
values of this coupling however, the tachyonic mode sets in. We observer a
similar phenomena inside the blueshift instability window. At large values of
$|L_{0}\mu|$ the tachyon mode takes over the entire process, but for small
couplings we recover all the $\mu=0$ results. Our results indicate that there
is not relevant change in the scalarization threshold induced by the boundary
oscillation.

###  VIII.2 Parametric instabilities

Up to now, we studied the effects of an oscillating boundary alone. However,
whenever an astrophysical object vibrates, it also changes its local density
in a periodic way. To mimic this effect in a one-dimensional setting let us
consider the system

$\displaystyle-\partial_{t}^{2}\Phi+\partial_{x}^{2}\Phi-\mu^{2}\Phi=0\,,$
(49) $\displaystyle\Phi(t,0)=\Phi(t,L_{0})=0\,,$

where $\mu^{2}(t)=\bar{\mu}^{2}+\delta\mu^{2}\sin(\omega t)$. Thus,
$\bar{\mu}^{2}$ corresponds to the average effective mass in the system. The
oscillating piece $\delta\mu^{2}$ mimics the time-varying density. Later on,
when dealing with a realistic setup, a certain condition between
$\delta\mu^{2}$ and $\delta L$ will be necessary to ensure mass conservation.
For now, we proceed with this idealization and without further constraints.

In Fourier space, that is $\Phi(t,x)=\sum_{k}\phi_{k}(t)e^{ikx}$, Eq. (49) can
be written as

$\displaystyle\partial_{\tau}^{2}\phi_{k}+\left(\delta+2\,\epsilon\sin(\tau)\right)\phi_{k}=0\,,$
(50)

where $\phi_{k}$ is the amplitude of the Fourier mode with wavelength $k$, and

$\tau=\omega
t\,,\qquad\delta=\frac{k^{2}+\bar{\mu}^{2}}{\omega^{2}}\,,\qquad\epsilon=\frac{\delta\mu^{2}}{2\,\omega^{2}}\,.$
(51)

Equation (50) is the well-known Mathieu equation Abramowitz and Stegun (1972);
Bender and Orszag (1978), and we refrain from providing further details. For
certain combination of the parameters, Mathieu equation admits unstable
solutions for $\delta\sim j^{2}/4$ with $j\in\mathbb{Z}$ Boskovic  _et al._
(2019); Chen  _et al._ (1996). In particular, for small $\epsilon$ and $j=1$,

$\delta=\frac{1}{4}+\epsilon\,\delta_{1}\,,$ (52)

one finds the instability rate $\lambda_{A}$ (i.e. $\Phi\sim
e^{\lambda_{A}t}$) Boskovic  _et al._ (2019); Chen  _et al._ (1996),

$\lambda_{A}=\omega\epsilon\sqrt{1-\delta_{1}^{2}}\,.$ (53)

This instability differs, fundamentally, from the previous “blueshift” or
Doppler-like mechanism. It is of parametric origin, and it does not lead to no
migration to higher frequencies. Instead, modes in the frequency range given
by Eq. (52) are amplified with the rate (53). Thus, this instability can only
be quenched via nonlinear effects.

##  IX Instabilities of oscillating objects in flat space

As we saw, fluctuations in some background value of the scalar are described
by the Klein-Gordon equation (5) with an effective mass $\mu^{2}=\mu^{2}(t,r)$
that depends on the energy density of the star. If the star oscillates, so
does the effective mass of the scalar field. We model this dependence as in
Eq. (7), where $\delta\tilde{\rho}(t,r)$ is the amplitude of the density
oscillation, and $\tilde{\rho}_{0}$ is the temporal and spatial average of the
energy density. Imposing conservation of the mass $M$ of the star, we find

$\tilde{\rho}_{0}=\frac{3}{4\pi}ML_{0}^{-3}\,,\quad\delta
L(t)=-\int_{0}^{L_{0}}dr\frac{r^{2}}{L_{0}^{2}}\frac{\delta\tilde{\rho}}{\rho_{0}}\,.$
(54)

To summarize, we evolve equation (5) with an effective mass

$\mu^{2}(t,r)=\begin{cases}\mu_{\rm in}^{2}+\delta\mu^{2}\sin(\omega
t)&(r<r_{-})\\\ \mu_{\rm smooth}^{2}(r)&(r_{-}<r<r_{+})\\\
\mu_{0}^{2}&(r>r_{+})\end{cases}$ (55)

where $r_{\pm}=L(t)\pm\frac{1}{2}\Delta{L}$, and $\Delta L$ defines a
thickness of the surface of the star, introduced to make $\mu$ a smooth
function. Here, $\mu_{\rm smooth}^{2}(r)$ is a smooth fifth-order polynomial
connecting the effective mass inside the star, $\mu_{\rm
in}^{2}+\delta\mu^{2}\sin(\omega t)$, to its exterior bare value
$\mu_{0}^{2}$. Note that $\mu_{\rm in}^{2}=\mu_{0}^{2}+\beta\tilde{\rho}_{0}$,
from Eq. (7).

###  IX.1 Numerical setup

Decomposing the scalar field in spherical harmonics, the evolution of the
$(lm)$ mode, $\Phi_{lm}(t,r)$, is governed by

$\left(\partial_{r}^{2}+\frac{2}{r}\partial_{r}-\frac{l(l+1)}{r^{2}}-\mu^{2}\right)\Phi_{lm}-\partial_{t}^{2}\Phi_{lm}=0\,.$
(56)

Regularity of the field at $r=0$ implies that $\Phi_{lm}\propto r^{l}$ for
$l\geq 0$, and we impose radiative boundary conditions at spatial infinity.

As in Sec. VIII.1.3, we numerically solve Eq. (56) using a fourth-order
accurate Runge-Kutta scheme for the time integration where spatial derivatives
are approximated by fourth-order accurate finite difference stencils. We
monitor the energy of the scalar field inside the star, defined as

$\displaystyle E_{lm}$
$\displaystyle=\frac{1}{2}\int_{0}^{L_{0}}\bigg{(}|\dot{\Phi}_{lm}|^{2}+|\partial_{r}\Phi_{lm}|^{2}+\frac{l(l+1)}{r^{2}}|\Phi_{lm}|^{2}$
$\displaystyle\qquad\qquad\quad+\mu^{2}(t,r)|\Phi_{lm}|^{2}\bigg{)}\,r^{2}dr\,,$
(57)

(although we mostly focus on spherically symmetric modes).

We once again use time-symmetric initial data, with a profile parameterized by
(8), where $\sigma$ and $r_{0}$ denoting the width and initial position of the
scalar field pulse. Since the equation to be solved is linear we set the
initial amplitude of the pulse to unity. As we will see below, we will mostly
focus on $l=0$ cases and we fix $r_{0}$ and $\delta L$ to

$r_{0}=0.5L_{0}\,,\qquad\delta L=0.1L_{0}\,.$ (58)

Figure 13: The purple line shows the difference between results obtained with
low ($4h$) and medium ($2h$) resolutions, while the (dashed) green line shows
the difference between results obtained with medium and high ($h$) resolutions
multiplied by 16, the expected factor for fourth-order convergence.

All our results are numerically convergent. An example is shown in Fig. 13,
showing the expected fourth-order convergence.

###  IX.2 Results

####  IX.2.1 Parametric instabilities

Figure 14: Instability rate of oscillating star, when density is allowed to
oscillate (for initial condition with width $\sigma=0.1L_{0}$ but results are
insensitive to this choice). Now, two mechanisms compete: a Doppler like or
blueshift instability and a parametric instability. In this example, with
$\mu_{0}L_{0}=100$, the latter always dominates. Crosses are data points for
$\delta\mu^{2}L_{0}^{2}=16$ and different $\delta L/L_{0}$; blue line is
prediction from Mathieu equation [Eq. (53)] with $k=2\pi/L_{0}$.

Figure 14 shows that the typical growth rate $\lambda_{A}$ of the
amplification instability as a function of $\omega L_{0}$ with
$\delta\mu^{2}L_{0}^{2}=16$ and different $\delta L/L_{0}$. We used the
evolution of the scalar to estimate the rate $\lambda_{A}$ of the instability.
We also numerically evolved Mathieu equation [Eq. (50)] with
$k=\frac{2\pi}{L_{0}}$, $\bar{\mu}=0$ and different $\delta\mu^{2}$, and
obtained the instability timescale, summarized in Fig. 14 (blue line). Note
that $k=\frac{2\pi}{L_{0}}$ is the longest wavelength which can be excited
inside the star and the relevant mode governing the instability. Our results
show that the numerical evolution is in good agreement with the simple
analysis of Sec. VIII.2, leading to a Mathieu equation. For small $\delta L$,
the timescale in the numerical simulation is in good agreement with the
prediction from Eq. (52), and even for $\delta L=0.1L_{0}$ – the difference
between the growth rate and the “Mathieu” prediction is only about 10%. Here,
using the well-known properties Boskovic  _et al._ (2019); Chen  _et al._
(1996) of the Mathieu equation, we can derive the resonance band in small
$\epsilon$, defined in Eq. (52):

$\omega\sim\frac{4\pi}{nL_{0}}\,,$ (59)

where $n$ is an integer. This expression and Fig. 14 indicate that there are
infinite number of instability windows at low frequency. The rate for $n=1$ at
small $\epsilon$ reads

$L_{0}\lambda\sim\frac{\beta L_{0}^{2}\delta\rho}{8\pi}\,,$ (60)

where it is assumed that $L_{0}\lambda\ll 1$. These results were derived
assuming a simple constant profile for $\delta\rho$. We have studied
nontrivial profiles, including sinusoidal shapes (more appropriate for the
description of normal modes) and our results are still compatible with the
above description.

##  X Instability of oscillating stars in General Relativity

We have thus far been considering toy models for instabilities using a flat
space approximation. To study the existence of instabilities in astrophysical
objects we will now consider radial perturbations around a relativistic
constant density star. The corresponding eigenvalue and stability problem was
studied previously in classical works Chandrasekhar (1964, 1964); Bardeen  _et
al._ (1966); Kokkotas and Ruoff (2001). As we will see, the amplification
instability discussed in previous sections also manifests itself in this
setting. This is a highly nontrivial result: the governing equations, as well
and the perturbation density profile are complicated functions, which can
hardly be related to any Mathieu-like equation. Nevertheless, our results show
that such stars are prone to parametric instabilities.

###  X.1 Formulation

Figure 15: The first four normal mode frequencies of a constant-density star
with compactness $\mathcal{C}=0.3$, for different index $\gamma$. We recover
well known results in the literature Chandrasekhar (1964); Bardeen  _et al._
(1966); Kokkotas and Ruoff (2001).

Let us consider perturbations around the metric, density, pressure, and
4-velocity profiles of the relativistic constant density star given by Eq.
(15). Any of these dynamical quantities $X(t,r)=(\alpha,\,a,\,\rho,\,p)$, is
then written as its background, stationary value corresponding to Eq. (15),
plus some small linearized fluctuation,

$X(t,r)=X_{0}(r)+\epsilon\delta X(t,r)\,,$ (61)

where $\epsilon$ is a bookkeeping parameter for the perturbation. We also need
to consider fluid displacements and we write the radial component of the four-
velocity
$u_{r}(t,r)=\epsilon\frac{\tilde{a}_{0}^{2}}{\tilde{\alpha}_{0}}\partial_{t}\xi(t,r)$,
with $\xi(t,r)$ the Lagrangian displacement. From the equation of motion for
fluid one obtains the following equation for $\xi$

$\displaystyle-$
$\displaystyle\frac{\tilde{a}_{0}^{2}}{\tilde{\alpha}_{0}^{2}}\left(\tilde{\rho}_{0}+\tilde{p}_{0}\right)\ddot{\xi}=\frac{4}{r}\tilde{p}_{0}^{\prime}\xi-\frac{1}{\tilde{a}_{0}\tilde{\alpha}_{0}^{2}}\left(\tilde{a}_{0}\tilde{\alpha}_{0}^{3}\frac{\gamma\,\tilde{p}_{0}}{r^{2}}\left(\frac{r^{2}}{\tilde{\alpha}_{0}}\xi\right)^{\prime}\right)^{\prime}$
(62) $\displaystyle+$ $\displaystyle
8\pi\tilde{a}_{0}^{2}\tilde{p}_{0}(\tilde{\rho}_{0}+\tilde{p}_{0})\xi-\frac{1}{\tilde{\rho}_{0}+\tilde{p}_{0}}(\tilde{p}_{0}^{\prime})^{2}\xi\,,$

where primes stand for radial derivatives and $\gamma$ is the adiabatic index
defined as

$\displaystyle\gamma=\frac{\tilde{\rho}+\tilde{p}}{\tilde{p}}\left(\frac{d\tilde{p}}{d\tilde{\rho}}\right)\,.$
(63)

For a constant density star, $\gamma$ is infinite (incompressible fluid);
however, we will assume that it is of (constant) finite value to perform the
perturbation Chandrasekhar (1964); Posada and Chirenti (2019). The density,
pressure, and metric perturbation can be obtained from the Lagrangian
displacement

$\displaystyle\delta\rho$
$\displaystyle=\frac{\tilde{p}_{0}+\tilde{\rho}_{0}}{\gamma}\left(4\pi
r\tilde{a}_{0}^{2}\tilde{p}_{0}(1+\gamma)+\frac{\tilde{\rho}_{0}}{2r\tilde{p}_{0}}(a^{2}-1)-\gamma\xi^{\prime}\right)$
$\displaystyle\quad{}+\frac{\tilde{p}_{0}+\tilde{\rho}_{0}}{2r\gamma}\left(-1-5\gamma+\tilde{a}_{0}^{2}\left(1+\gamma+8\pi
r^{2}\tilde{\rho}_{0}\right)\right)\xi\,,$ (64a) $\displaystyle\delta p$
$\displaystyle=\frac{\tilde{p}_{0}}{\tilde{\rho}_{0}+\tilde{p}_{0}}\gamma\delta\rho\,,$
(64b) $\displaystyle\delta a$ $\displaystyle=-4\pi
r\tilde{a}_{0}(\tilde{p}_{0}+\tilde{\rho}_{0})\xi\,,$ (64c)
$\displaystyle\delta\alpha^{\prime}$
$\displaystyle=-4\pi\tilde{a}_{0}^{4}(\tilde{\rho}_{0}+\tilde{p}_{0})(1+8\pi
r^{2}\tilde{p}_{0})\xi$
$\displaystyle\quad{}+\left(\frac{-1+\tilde{a}_{0}^{2}}{r}+8\pi
r\tilde{a}_{0}^{2}\tilde{p}_{0}-\frac{\tilde{\alpha}_{0}^{\prime}}{\tilde{\alpha}_{0}}\right)\delta\alpha\,.$
(64d)

We focus on the normal mode of the perturbation and assume that
$\xi(t,r)=e^{i\omega t}\zeta_{\omega}(r)=e^{i\omega
t}\frac{\tilde{\alpha}_{0}(r)}{r^{2}}\zeta_{\omega}$. The equation for
$\zeta_{\omega}$ is

$\displaystyle\frac{\zeta_{\omega}^{\prime\prime}}{\zeta_{\omega}}$
$\displaystyle=-\frac{\tilde{a}_{0}^{2}}{\gamma\tilde{\alpha}_{0}^{2}}\left(1+\frac{\tilde{\rho}_{0}}{\tilde{p}_{0}}\right)\omega^{2}+8\pi\frac{\tilde{a}_{0}^{2}}{\gamma}(\tilde{p}_{0}+\tilde{\rho}_{0})+\frac{4\tilde{p}_{0}^{\prime}}{r\gamma\tilde{p}_{0}}$
$\displaystyle\quad{}-\frac{\tilde{p}_{0}^{\prime
2}}{\gamma\tilde{p}_{0}(\tilde{\rho}_{0}+\tilde{p}_{0})}+\left(\frac{2}{r}-3\frac{\tilde{\alpha}_{0}^{\prime}}{\tilde{\alpha}_{0}}-\frac{\tilde{a}_{0}^{\prime}}{\tilde{a}_{0}}-\frac{\tilde{p}_{0}^{\prime}}{\tilde{p}_{0}}-\frac{\gamma^{\prime}}{\gamma}\right)\frac{\zeta_{\omega}^{\prime}}{\zeta_{\omega}}\,.$

The boundary condition for $\xi_{\omega}(r)$ is

$\xi_{\omega}(0)=0\,,\qquad\delta
p(L_{0})\propto\xi_{\omega}^{\prime}(L_{0})=0\,,$ (65)

which corresponds to the following boundary condition for $\zeta_{\omega}(r)$

$\displaystyle\zeta_{\omega}(0)$ $\displaystyle=0\,,$ (66)
$\displaystyle\zeta^{\prime}_{\omega}(L_{0})$
$\displaystyle=-\left(\frac{\tilde{\alpha}_{0}^{\prime}(L_{0})}{\tilde{\alpha}_{0}(L_{0})}-\frac{1}{L_{0}}\right)\zeta_{\omega}(L_{0})\,.$

We solved the eigenvalue problem for $\omega$ with the above boundary
conditions for different adiabatic index $\gamma$ with a shooting method. Our
results give a relation between eigenfrequency $\omega$ and index $\gamma$ for
different overtones, labeled by $n$. They are summarized in Fig. 15, and are
consistent with classical results in the literature Chandrasekhar (1964);
Bardeen  _et al._ (1966); Kokkotas and Ruoff (2001). For the simulation of the
scalar field, we introduce a 5th order polynomial smooth transition function
$W(r;R_{1},R_{2})$ satisfying

$\displaystyle W(r;R_{1},R_{2})$ $\displaystyle=1~{}~{}(r\leq R_{1}),$
$\displaystyle W(r;R_{1},R_{2})$ $\displaystyle=0~{}~{}(r\geq R_{2}),$
$\displaystyle W^{\prime}(R_{1};R_{1},R_{2})$
$\displaystyle=W^{\prime}(R_{2};R_{1},R_{2})=W^{\prime\prime}(R_{1};R_{1},R_{2})$
$\displaystyle=W^{\prime\prime}(R_{2};R_{1},R_{2})=0\,,$

and replace the Lagrangian displacement with the smoothed profile as follows

$\displaystyle\xi_{\omega}(r)$
$\displaystyle\to\xi_{\omega}(r)W(r;R_{1},R_{2})\,.$

With this we construct density, pressure, and metric profile from

$\displaystyle Q(t,r)$ $\displaystyle=Q_{0}(r+\xi(t,r))+\delta
Q(t,r)W(r;R_{1},R_{2})$ $\displaystyle Q_{0}(r)$ $\displaystyle=Q_{\rm
in}(t,r)W(r;R_{1},R_{2})+Q_{\rm out}(1-W(r;R_{1},R_{2}))$

where $Q(t,r)$ denotes the density, pressure, and the metric, and $Q_{0}(r)$
denotes the background profile for each variable.

We are then able to evolve the scalar field $\phi=\frac{\psi}{r}$ through the
evolution equation

$\displaystyle\partial_{t}^{2}\phi=\partial_{t}\ln\left|\frac{\alpha}{a}\right|\partial_{t}\phi+\frac{\alpha^{2}}{a^{2}}\partial_{r}\ln\left|\frac{\alpha}{a}\right|\partial_{r}\phi+\frac{\alpha^{2}}{a^{2}}\left(\frac{2}{r}\partial_{r}\phi+\partial_{r}^{2}\phi\right)-\frac{l(l+1)}{r^{2}}\alpha^{2}\phi+\beta\tilde{T}\alpha^{2}\phi-\mu^{2}_{0}\alpha^{2}\phi\,.$
(67)

For the simulations reported below, we use a momentarily static Gaussian pulse
as initial data parameterized by

$\displaystyle\phi(0,r)=Ae^{-\left(\frac{r-r_{0}}{w}\right)^{2}}\,,\qquad\dot{\phi}(0,r)=0\,,$

where $A,w,r_{0}$ are the initial amplitude, width, and position of the pulse,
respectively. We will always set $r_{0}=0.5L_{0}$.

During the evolution we monitor the scalar field at fixed radius, and the
energy of the scalar field inside star, given by

$\displaystyle E=\frac{1}{2}\int_{0}^{L_{0}}dr\alpha ar^{2}\Big{(}$
$\displaystyle\frac{1}{\alpha^{2}}\dot{\phi}^{2}+\frac{1}{a^{2}}\phi^{\prime
2}+\frac{l(l+1)}{r^{2}}\phi^{2}$
$\displaystyle{}+(\mu_{0}^{2}-\beta\tilde{T})\phi^{2}\Big{)}$ (68)

###  X.2 Results

Figure 16:  Time evolution with $\mu_{0}L_{0}=100$, $\beta=-280000$,
$\delta\rho_{c}\simeq 0.08\tilde{\rho}_{0}$, and different initial width of
the pulse. The initial width $w/L_{0}$ is $0.2$ (upper panel), and $0.4$
(lower panel). The background star is perturbed with the fundamental mode
($n=0$).

We now take an oscillating compact star. The star oscillations could be caused
by accretion or because it was the product of a merger. The amplitude of such
oscillations can in any case be significant Reisenegger and Goldreich (1994);
Chirenti  _et al._ (2017); Ma  _et al._ (2020); Perego  _et al._ (2019);
Bernuzzi  _et al._ (2020). In our simulations we consider a single mode of
oscillation, characterized by the number of nodes $n$. The generalization to
accommodate for different modes is straightforward. The control parameters for
the background spacetime are the compactness of the star $\mathcal{C}$, the
adiabatic index $\gamma$, the central value of the density perturbation
$\delta\rho_{c}$, and $n$. Here, the compactness $\mathcal{C}$ and adiabatic
index $\gamma$ are fixed to $0.3$ and $2.0$, which are representative of
neutron stars Lattimer and Prakash (2001). The most drastic astrophysical
situations exciting large fluctuations of neutron stars concerns the remnant
phase just after a neutron star merger. According to numerical relativity
simulations, the amplitude of the density perturbations can be significant and
of the order of the central density itself Perego  _et al._ (2019); Bernuzzi
_et al._ (2020). Although such large density perturbations lie beyond the
validity range of the linear approximation for the perturbation, we consider
${\rm max}(\delta\rho(t,0))\simeq 0.28\tilde{\rho}_{0},0.08\tilde{\rho}_{0}$,
and use the linear normal mode derived in previous subsection.

Figure 16 shows the time evolution of the scalar field at $r=0.5L_{0}$ with
$\mu_{0}L_{0}=100$, $\beta=-280000$, and $\delta\rho_{c}\simeq
0.08\tilde{\rho}(0)$, and the fundamental density perturbation mode $n=0$.
Note how the scalar has a time-periodic behavior, but is not growing neither
in amplitude nor in frequency. With these parameters we could not observe any
instabilities.

Figure 17:  Time evolution with $\mu_{0}L_{0}=100$, $\beta=-280000$,
$\delta\rho_{c}\simeq 0.28\tilde{\rho}_{0}$, $w/L_{0}=0.2$. We plot radial
modes of density perturbation with $n=0,1,2$. Notice that the density
fluctuation is large, but still smaller than typical values in seen during the
coalescence of two neutron stars Perego  _et al._ (2019); Bernuzzi  _et al._
(2020). Also, for smaller $\beta$, the instability sets in at small density
fluctuations.

However, an instability sets in with either larger values of density
fluctuation or smaller values of $\beta$ (i.e., larger $\beta^{2}$). We note
that for massive fields – our subject here – observational constraints on
$\beta$ are poor or non-existing and any of the values discussed here is
compatible with current observations. For example, Fig. 17 shows the time
evolution of the scalar field at $r=0.5L_{0}$, but now with
$\delta\rho_{c}\simeq 0.28\tilde{\rho}_{0}$. We can compare and contrast the
results with those of Fig. 16. For these larger values of density
fluctuations, a parametric instability sets in. The frequency of the
fundamental density perturbation mode with $\gamma=2.0$ is about $0.36/L_{0}$
and the corresponding period is about $17L_{0}$. This period is clearly
evident in Fig. 17, in the form of a plateau in the evolution of the scalar.
In other words, the growth proceeds upon interaction with the star with a
periodicity dictated by the star eigenfrequency. As is clear from the
evolution equations or the background structure of the star itself, it is far
from trivial that such oscillating stars in General Relativity would trigger
parametric instabilities. Our results confirm it. We note that Fig. 17 refers
to a significant density fluctuation. However, our results are mostly
sensitive to the combination $\beta\tilde{\rho}_{0}$, thus much milder
fluctuations also result in instability, provided $\beta^{2}$ is sufficiently
large.

## References

  * Will (2014) Clifford M. Will, “The Confrontation between General Relativity and Experiment,” Living Rev. Rel. 17, 4 (2014), arXiv:1403.7377 [gr-qc] .
  * Abbott  _et al._ (2016a) B.P. Abbott _et al._ (LIGO Scientific, Virgo), “GW150914: The Advanced LIGO Detectors in the Era of First Discoveries,” Phys. Rev. Lett. 116, 131103 (2016a), arXiv:1602.03838 [gr-qc] .
  * Abbott  _et al._ (2016b) B. P. Abbott _et al._ (LIGO Scientific, Virgo), “Tests of general relativity with GW150914,” Phys. Rev. Lett. 116, 221101 (2016b), [Erratum: Phys. Rev. Lett.121,no.12,129902(2018)], arXiv:1602.03841 [gr-qc] .
  * Abbott  _et al._ (2020) R. Abbott _et al._ (LIGO Scientific, Virgo), “Tests of General Relativity with Binary Black Holes from the second LIGO-Virgo Gravitational-Wave Transient Catalog,” (2020), arXiv:2010.14529 [gr-qc] .
  * Clifton  _et al._ (2012) Timothy Clifton, Pedro G. Ferreira, Antonio Padilla, and Constantinos Skordis, “Modified Gravity and Cosmology,” Phys. Rept. 513, 1–189 (2012), arXiv:1106.2476 [astro-ph.CO] .
  * Berti  _et al._ (2015) Emanuele Berti _et al._ , “Testing General Relativity with Present and Future Astrophysical Observations,” Class. Quant. Grav. 32, 243001 (2015), arXiv:1501.07274 [gr-qc] .
  * Cardoso and Pani (2019) Vitor Cardoso and Paolo Pani, “Testing the nature of dark compact objects: a status report,” Living Rev. Rel. 22, 4 (2019), arXiv:1904.05363 [gr-qc] .
  * 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] .
  * Barack  _et al._ (2019) Leor Barack _et al._ , “Black holes, gravitational waves and fundamental physics: a roadmap,” Class. Quant. Grav. 36, 143001 (2019), arXiv:1806.05195 [gr-qc] .
  * Kokkotas and Ruoff (2001) Kostas. D. Kokkotas and Johannes Ruoff, “Radial oscillations of relativistic stars,” Astron. Astrophys. 366, 565 (2001), arXiv:gr-qc/0011093 .
  * Schatzman (1961) E. Schatzman, “Sur la période de pulsation radiale des naines blanches,” Annales d’Astrophysique 24, 237 (1961).
  * Singh and Das (1991) Harinder P. Singh and M. K. Das, “Radial Oscillations of Finite Temperature White Dwarfs,” Astrophysics and Space Science 186, 117–124 (1991).
  * Winget and Kepler (2008) D. E. Winget and S. O. Kepler, “Pulsating White Dwarf Stars and Precision Asteroseismology,” Ann. Rev. Astron. Astrophys. 46, 157 (2008), arXiv:0806.2573 [astro-ph] .
  * Turck-Chieze and Lopes (1993) Sylvaine Turck-Chieze and I. Lopes, “Toward a unified classical model of the sun: On the sensitivity of neutrinos and helioseismology to the microscopic physics,” Astrophys. J. 408, 347–367 (1993).
  * Turck-Chièze and Lopes (2012) Sylvaine Turck-Chièze and Ilídio Lopes, “Solar-stellar astrophysics and dark matter,” Research in Astronomy and Astrophysics 12, 1107–1138 (2012).
  * Damour and Esposito-Farese (1993) Thibault Damour and Gilles Esposito-Farese, “Nonperturbative strong field effects in tensor - scalar theories of gravitation,” Phys. Rev. Lett. 70, 2220–2223 (1993).
  * Damour and Esposito-Farese (1996) Thibault Damour and Gilles Esposito-Farese, “Tensor - scalar gravity and binary pulsar experiments,” Phys. Rev. D 54, 1474–1491 (1996), arXiv:gr-qc/9602056 .
  * Harada (1998) Tomohiro Harada, “Neutron stars in scalar tensor theories of gravity and catastrophe theory,” Phys. Rev. D 57, 4802–4811 (1998), arXiv:gr-qc/9801049 .
  * Pani  _et al._ (2011) Paolo Pani, Vitor Cardoso, Emanuele Berti, Jocelyn Read, and Marcelo Salgado, “The vacuum revealed: the final state of vacuum instabilities in compact stars,” Phys. Rev. D 83, 081501 (2011), arXiv:1012.1343 [gr-qc] .
  * Ramazanoğlu and Pretorius (2016) Fethi M. Ramazanoğlu and Frans Pretorius, “Spontaneous Scalarization with Massive Fields,” Phys. Rev. D 93, 064005 (2016), arXiv:1601.07475 [gr-qc] .
  * Antoniou  _et al._ (2018) G. Antoniou, A. Bakopoulos, and P. Kanti, “Evasion of No-Hair Theorems and Novel Black-Hole Solutions in Gauss-Bonnet Theories,” Phys. Rev. Lett. 120, 131102 (2018), arXiv:1711.03390 [hep-th] .
  * Doneva and Yazadjiev (2018) Daniela D. Doneva and Stoytcho S. Yazadjiev, “New Gauss-Bonnet Black Holes with Curvature-Induced Scalarization in Extended Scalar-Tensor Theories,” Phys. Rev. Lett. 120, 131103 (2018), arXiv:1711.01187 [gr-qc] .
  * Silva  _et al._ (2018) Hector O. Silva, Jeremy Sakstein, Leonardo Gualtieri, Thomas P. Sotiriou, and Emanuele Berti, “Spontaneous scalarization of black holes and compact stars from a Gauss-Bonnet coupling,” Phys. Rev. Lett. 120, 131104 (2018), arXiv:1711.02080 [gr-qc] .
  * Khoury and Weltman (2004a) Justin Khoury and Amanda Weltman, “Chameleon fields: Awaiting surprises for tests of gravity in space,” Phys. Rev. Lett. 93, 171104 (2004a), arXiv:astro-ph/0309300 .
  * Khoury and Weltman (2004b) Justin Khoury and Amanda Weltman, “Chameleon cosmology,” Phys. Rev. D 69, 044026 (2004b), arXiv:astro-ph/0309411 .
  * Joyce  _et al._ (2015) Austin Joyce, Bhuvnesh Jain, Justin Khoury, and Mark Trodden, “Beyond the Cosmological Standard Model,” Phys. Rept. 568, 1–98 (2015), arXiv:1407.0059 [astro-ph.CO] .
  * Cardoso  _et al._ (2020) Vitor Cardoso, Arianna Foschi, and Miguel Zilhao, “Collective scalarization or tachyonization: when averaging fails,” Phys. Rev. Lett. 124, 221104 (2020), arXiv:2005.12284 [gr-qc] .
  * Tsujikawa  _et al._ (2009) Shinji Tsujikawa, Takashi Tamaki, and Reza Tavakol, “Chameleon scalar fields in relativistic gravitational backgrounds,” JCAP 05, 020 (2009), arXiv:0901.3226 [gr-qc] .
  * Barausse  _et al._ (2013) Enrico Barausse, Carlos Palenzuela, Marcelo Ponce, and Luis Lehner, “Neutron-star mergers in scalar-tensor theories of gravity,” Phys. Rev. D 87, 081506 (2013), arXiv:1212.5053 [gr-qc] .
  * Palenzuela  _et al._ (2014) Carlos Palenzuela, Enrico Barausse, Marcelo Ponce, and Luis Lehner, “Dynamical scalarization of neutron stars in scalar-tensor gravity theories,” Phys. Rev. D 89, 044024 (2014), arXiv:1310.4481 [gr-qc] .
  * Shibata  _et al._ (2014) Masaru Shibata, Keisuke Taniguchi, Hirotada Okawa, and Alessandra Buonanno, “Coalescence of binary neutron stars in a scalar-tensor theory of gravity,” Phys. Rev. D 89, 084005 (2014), arXiv:1310.0627 [gr-qc] .
  * Khalil  _et al._ (2019) Mohammed Khalil, Noah Sennett, Jan Steinhoff, and Alessandra Buonanno, “Theory-agnostic framework for dynamical scalarization of compact binaries,” Phys. Rev. D 100, 124013 (2019), arXiv:1906.08161 [gr-qc] .
  * Silva  _et al._ (2021) Hector O. Silva, Helvi Witek, Matthew Elley, and Nicolás Yunes, “Dynamical Descalarization in Binary Black Hole Mergers,” Phys. Rev. Lett. 127, 031101 (2021), arXiv:2012.10436 [gr-qc] .
  * East and Ripley (2021) William E. East and Justin L. Ripley, “Dynamics of spontaneous black hole scalarization and mergers in Einstein-scalar-Gauss-Bonnet gravity,” (2021), arXiv:2105.08571 [gr-qc] .
  * Llinares and Mota (2013) Claudio Llinares and David Mota, “Releasing scalar fields: cosmological simulations of scalar-tensor theories for gravity beyond the static approximation,” Phys. Rev. Lett. 110, 161101 (2013), arXiv:1302.1774 [astro-ph.CO] .
  * Hagala  _et al._ (2017) R. Hagala, C. Llinares, and D.F. Mota, “Cosmic Tsunamis in Modified Gravity: Disruption of Screening Mechanisms from Scalar Waves,” Phys. Rev. Lett. 118, 101301 (2017), arXiv:1607.02600 [astro-ph.CO] .
  * Brito  _et al._ (2014) Richard Brito, Alexandra Terrana, Matthew Johnson, and Vitor Cardoso, “Nonlinear dynamical stability of infrared modifications of gravity,” Phys. Rev. D 90, 124035 (2014), arXiv:1409.0886 [hep-th] .
  * Silvestri (2011) Alessandra Silvestri, “Scalar radiation from Chameleon-shielded regions,” Phys. Rev. Lett. 106, 251101 (2011), arXiv:1103.4013 [astro-ph.CO] .
  * Nakamura  _et al._ (2021) Tomohiro Nakamura, Taishi Ikeda, Ryo Saito, Norihiro Tanahashi, and Chul-Moon Yoo, “Dynamical Analysis of Screening in Scalar-Tensor Theory,” Phys. Rev. D 103, 024009 (2021), arXiv:2010.14329 [gr-qc] .
  * Wang  _et al._ (2011) Charles H. T. Wang, Paolo M. Bonifacio, Robert Bingham, and J. Tito Mendonca, “Parametric instability induced scalar gravitational waves from a model pulsating neutron star,” Phys. Lett. B 705, 148–151 (2011), arXiv:0905.0706 [gr-qc] .
  * Wang  _et al._ (2013) C. H. T. Wang, A. O. Hodson, A. St. J. Murphy, T. B. Davies, R. Bingham, and J. T. Mendonça, “Dynamical trapping and relaxation of scalar gravitational fields,” Phys. Lett. B 726, 791–794 (2013), arXiv:1301.0269 [gr-qc] .
  * Ventagli  _et al._ (2020) Giulia Ventagli, Antoine Lehébel, and Thomas P. Sotiriou, “Onset of spontaneous scalarization in generalized scalar-tensor theories,” Phys. Rev. D 102, 024050 (2020), arXiv:2006.01153 [gr-qc] .
  * Hinterbichler  _et al._ (2011) Kurt Hinterbichler, Justin Khoury, Aaron Levy, and Andrew Matas, “Symmetron Cosmology,” Phys. Rev. D 84, 103521 (2011), arXiv:1107.2112 [astro-ph.CO] .
  * Hinterbichler and Khoury (2010) Kurt Hinterbichler and Justin Khoury, “Symmetron Fields: Screening Long-Range Forces Through Local Symmetry Restoration,” Phys. Rev. Lett. 104, 231301 (2010), arXiv:1001.4525 [hep-th] .
  * Shapiro and Teukolsky (1983) S. L. Shapiro and S. A. Teukolsky,  _Black holes, white dwarfs, and neutron stars: The physics of compact objects_ (Wiley, New York, 1983).
  * Sotani (2014) Hajime Sotani, “Scalar gravitational waves from relativistic stars in scalar-tensor gravity,” Phys. Rev. D 89, 064031 (2014), arXiv:1402.5699 [astro-ph.HE] .
  * Mendes and Ortiz (2018) Raissa F. P. Mendes and Néstor Ortiz, “New class of quasinormal modes of neutron stars in scalar-tensor gravity,” Phys. Rev. Lett. 120, 201104 (2018), arXiv:1802.07847 [gr-qc] .
  * Chandrasekhar (1964) S. Chandrasekhar, “The Dynamical Instability of Gaseous Masses Approaching the Schwarzschild Limit in General Relativity.” Astrophys. J. 140, 417 (1964).
  * Bardeen  _et al._ (1966) James M. Bardeen, Kip S. Thorne, and David W. Meltzer, “A Catalogue of Methods for Studying the Normal Modes of Radial Pulsation of General-Relativistic Stellar Models,” Astrophys. J. 145, 505 (1966).
  * Dittrich  _et al._ (1994a) J Dittrich, P Duclos, and P Šeba, “Instability in a classical periodically driven string,” Physical Review E 49, 3535 (1994a).
  * Cutler and Lindblom (1987) Curt Cutler and Lee Lindblom, “The Effect of Viscosity on Neutron Star Oscillations,” Astrophys. J. 314, 234 (1987).
  * Barta (2019) Dániel Barta, “Effect of viscosity and thermal conductivity on the radial oscillation and relaxation of relativistic stars,” Classical and Quantum Gravity 36, 215012 (2019), arXiv:1904.00907 [gr-qc] .
  * Abramowitz and Stegun (1972) M. Abramowitz and I. A. Stegun,  _Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables_ (Dover, New York, 1972).
  * Bender and Orszag (1978) C. M. Bender and S. A. Orszag,  _Advanced Mathematical Methods for Scientists and Engineers_ (McGraw-Hill, New York, 1978).
  * Reisenegger and Goldreich (1994) Andreas Reisenegger and Peter Goldreich, “Excitation of Neutron Star Normal Modes during Binary Inspiral,” Astrophys. J. 426, 688 (1994).
  * Chirenti  _et al._ (2017) Cecilia Chirenti, Roman Gold, and M. Coleman Miller, “Gravitational waves from f-modes excited by the inspiral of highly eccentric neutron star binaries,” Astrophys. J. 837, 67 (2017), arXiv:1612.07097 [astro-ph.HE] .
  * Ma  _et al._ (2020) Sizheng Ma, Hang Yu, and Yanbei Chen, “Excitation of f-modes during mergers of spinning binary neutron star,” Phys. Rev. D 101, 123020 (2020), arXiv:2003.02373 [gr-qc] .
  * Perego  _et al._ (2019) Albino Perego, Sebastiano Bernuzzi, and David Radice, “Thermodynamics conditions of matter in neutron star mergers,” Eur. Phys. J. A 55, 124 (2019), arXiv:1903.07898 [gr-qc] .
  * Bernuzzi  _et al._ (2020) Sebastiano Bernuzzi _et al._ , “Accretion-induced prompt black hole formation in asymmetric neutron star mergers, dynamical ejecta and kilonova signals,” Mon. Not. Roy. Astron. Soc. 497, 1488–1507 (2020), arXiv:2003.06015 [astro-ph.HE] .
  * Note (1) There were important hints that such an effect could occur. The classical Fermi acceleration process is one example of energy extraction from a zero-average motion Fermi (1949). Similar processes exist for fields trapped inside cavities with periodically-moving boundaries Cooper (1993); Dittrich  _et al._ (1994b); Moore (1970); Andreata and Dodonov (2000); Dodonov (2001). We showed that oscillating stars are also prone to such instabilities.
  * Annulli  _et al._ (2019) Lorenzo Annulli, Vitor Cardoso, and Leonardo Gualtieri, “Electromagnetism and hidden vector fields in modified gravity theories: spontaneous and induced vectorization,” Phys. Rev. D 99, 044038 (2019), arXiv:1901.02461 [gr-qc] .
  * Minamitsuji (2020a) Masato Minamitsuji, “Stealth spontaneous spinorization of relativistic stars,” Phys. Rev. D 102, 044048 (2020a), arXiv:2008.12758 [gr-qc] .
  * Minamitsuji (2020b) Masato Minamitsuji, “Spontaneous vectorization in the presence of vector field coupling to matter,” Phys. Rev. D 101, 104044 (2020b), arXiv:2003.11885 [gr-qc] .
  * Antoniadis  _et al._ (2013) John Antoniadis _et al._ , “A Massive Pulsar in a Compact Relativistic Binary,” Science 340, 6131 (2013), arXiv:1304.6875 [astro-ph.HE] .
  * Mendes (2015) Raissa F. P. Mendes, “Possibility of setting a new constraint to scalar-tensor theories,” Phys. Rev. D 91, 064024 (2015), arXiv:1412.6789 [gr-qc] .
  * Cooper (1993) J. Cooper, “Long-time behavior and energy growth for electromagnetic waves reflected by a moving boundary,” IEEE Transactions on Antennas and Propagation 41, 1365–1370 (1993).
  * Dittrich  _et al._ (1994b) J. Dittrich, P. Duclos, and P. Šeba, “Instability in a classical periodically driven string,” Phys. Rev. E 49, 3535–3538 (1994b).
  * Law (1994) C.K. Law, “Resonance Response of the Quantum Vacuum to an Oscillating Boundary,” Phys. Rev. Lett. 73, 1931–1934 (1994).
  * Boskovic  _et al._ (2019) Mateja Boskovic, Richard Brito, Vitor Cardoso, Taishi Ikeda, and Helvi Witek, “Axionic instabilities and new black hole solutions,” Phys. Rev. D99, 035006 (2019), arXiv:1811.04945 [gr-qc] .
  * Chen  _et al._ (1996) Lin-Yuan Chen, Nigel Goldenfeld, and Y. Oono, “The Renormalization group and singular perturbations: Multiple scales, boundary layers and reductive perturbation theory,” Phys. Rev. E 54, 376–394 (1996), arXiv:hep-th/9506161 .
  * Chandrasekhar (1964) S. Chandrasekhar, “The Dynamical Instability of Gaseous Masses Approaching the Schwarzschild Limit in General Relativity,” Astrophys. J. 140, 417–433 (1964), [Erratum: Astrophys.J. 140, 1342 (1964)].
  * Posada and Chirenti (2019) Camilo Posada and Cecilia Chirenti, “On the radial stability of ultra compact Schwarzschild stars beyond the Buchdahl limit,” Class. Quant. Grav. 36, 065004 (2019), arXiv:1811.09589 [gr-qc] .
  * Lattimer and Prakash (2001) J. M. Lattimer and M. Prakash, “Neutron star structure and the equation of state,” Astrophys. J. 550, 426 (2001), arXiv:astro-ph/0002232 .
  * Fermi (1949) Enrico Fermi, “On the Origin of the Cosmic Radiation,” Phys. Rev. 75, 1169–1174 (1949).
  * Moore (1970) G. T. Moore, “Quantum theory of the electromagnetic field in a variable‐length one‐dimensional cavity,” Journal of Mathmatical Physics 11, 2679 (1970).
  * Andreata and Dodonov (2000) M A Andreata and V V Dodonov, “Energy density and packet formation in a vibrating cavity,” Journal of Physics A: Mathematical and General 33, 3209–3223 (2000).
  * Dodonov (2001) V. V. Dodonov, “Nonstationary Casimir effect and analytical solutions for quantum fields in cavities with moving boundaries,” Adv. Chem. Phys. 119, 309–394 (2001), arXiv:quant-ph/0106081 .