# Effects of finite sizes of atomic nuclei on shear modulus and torsional
oscillations in neutron stars

Hajime Sotani1,2 , Hajime Togashi3,4 and Masatoshi Takano5,6
1Astrophysical Big Bang Laboratory, RIKEN, Saitama 351-0198, Japan
2Interdisciplinary Theoretical & Mathematical Science Program (iTHEMS), RIKEN,
Saitama 351-0198, Japan
3Department of Physics, Tohoku University, Sendai 980-8578, Japan
4Waseda Research Institute for Science and Engineering, Waseda University,
Tokyo 169-8555, Japan
5Department of Physics, Waseda University, Tokyo 169-8555, Japan
6Institute for Advanced Theoretical and Experimental Physics, Waseda
University, Tokyo 169-8555, Japan<EMAIL_ADDRESS>
###### Abstract

The shear modulus of neutron star matter is one of the important properties
for determining torsional oscillations in neutron stars. We take into account
the effects of finite sizes of spherical nuclei on the shear modulus and
examine the frequencies of crustal torsional oscillations. The shear modulus
decreases owing to the finite-size effect, which in turn decreases the
frequencies of torsional oscillations. In particular, the finite-size effect
becomes more crucial for oscillations with a larger azimuthal quantum number
and for neutron star models with a weaker density dependence of nuclear
symmetry energy. In practice, when one identifies the quasi-periodic
oscillations from a neutron star, where the magnetic effect is negligible,
with crustal torsional oscillations, the finite-size effect can be more
significant at frequencies higher than $\sim 100$ Hz.

###### keywords:

asteroseismology – stars: neutron – stars: oscillations

##  1 Introduction

A neutron star produced through a supernova explosion is one of the most
suitable environments for probing the physics under extreme states. The
density inside neutron stars easily exceeds the standard nuclear density, and
the gravitational and magnetic fields can be much stronger than those observed
in our solar system (Shapiro & Teukolsky, 1983). By observing neutron stars or
their associated phenomena, one may obtain information that cannot be produced
on Earth. For example, the discovery of 2$M_{\odot}$ neutron stars excludes
all the soft equations of state (EOSs), with which the expected maximum mass
does not reach the observed masses (Demorest et al., 2010; Antoniadis et al.,
2013; Cromartie et al., 2020). In addition, the pulsar light curve is affected
by the general-relativistic light-bending effect, which primarily depends on
stellar compactness, i.e., the ratio of the stellar mass to the stellar
radius. Thus, by carefully observing the pulsar light curves, the mass and
radius of neutron stars can be restricted (Riley et al., 2019, 2021; Miller et
al., 2019, 2021). The resultant constraints are still largely uncertain, but
they must be helpful in constraining the EOSs by accumulating the
observational data on the mass and radius for additional objects. The
gravitational wave event arising from the merging of binary neutron stars,
GW170817 (Abbott et al., 2017a, b), also provides information on the
constraint on tidal deformability, which indicates that the $1.4M_{\odot}$
neutron star radius should be less than 13.6 km (Annala et al., 2018).

The oscillation frequencies of neutron stars are another important piece of
information. Since the oscillation spectrum strongly depends on the interior
properties of an object, one may inversely obtain the interior information by
observing the corresponding spectrum. This technique is known as
asteroseismology. The asteroseismology of neutron stars was developed much
earlier in order to interpret various possible neutron star oscillations
detected in the electromagnetic spectrum (e.g., McDermott, Van Horn, & Hansen
(1988)). Even now, the most important results of asteroseismology are those
associated with electromagnetic observations, e.g., quasi-periodic
oscillations (QPOs) from soft-gamma repeaters (SGRs) as we will describe
below. Furthermore, it has been proposed that the neutron star mass, radius,
and EOS can be constrained by observing the gravitational waves from isolated
compact objects (e.g., Andersson & Kokkotas (1996, 1998); Sotani, Tominaga &
Maeda (2001); Sotani, Kohri & Harada (2004); Sotani et al. (2011); Doneva et
al. (2013); Sotani (2021); Sotani & Kumar (2021); Sotani & Dohi (2022)).
Regarding the gravitational asteroseismology of neutron stars, the important
role of $r$-modes in rapidly rotating neutron stars has also been often
discussed (e.g., Andersson & Kokkotas (2001)). Additionally, to understand the
gravitational wave signals appearing in the numerical simulation of a core-
collapse supernova explosion, asteroseismology has recently been adopted
(e.g., Morozova et al. (2018); Sotani et al. (2019); Sotani & Sumiyoshi
(2019); Torres-Forné et al. (2019); Sotani, Takiwaki & Togashi (2021)).
Gravitational waves from isolated neutron stars have never been detected yet,
but once they can be directly observed, one may be able to extract neutron
star properties by asteroseismology. Instead of gravitational waves, however,
several QPO frequencies have been detected in the x-ray afterglow following
giant flares from SGRs, e.g., SGR 1900+14 and SGR 1806-20, whose frequencies
are in the range of a few tenths of Hz up to kHz (Strohmayer & Watts, 2005,
2006). If one assumes that these QPOs come from crustal torsional oscillations
in a neutron star, one can constrain the crust properties (e.g., Steiner &
Watts (2009); Gearheart et al. (2011); Passamonti & Andersson (2012); Sotani
et al. (2012, 2013a, 2013b); Sotani (2014); Sotani, Iida & Oyamatsu (2016,
2017a, 2018, 2019)), even though the emission mechanism of QPOs is not
understood yet. We note that, since the SGRs are magnetars, i.e., highly
magnetized neutron stars, one may have to take into account the magnetic field
effect on torsional oscillations. In fact, the frequencies of crustal
torsional oscillations increase as the field strength increases (Sotani,
Kokkotas & Stergioulas, 2007; Gabler et al., 2013). Then, eventually, the
crustal torsional oscillations disappear and Alfven oscillations become
dominant as the magnetar oscillations. Nevertheless, considering the
observational estimation of the magnetic field strength of SGR 1900+14 and SGR
1806-20 inferred from the spin-down rate (Kouveliotou et al., 1998; Hurley et
al., 1999), the crustal torsional oscillations must still survive.

Crustal torsional oscillations are characterized by the shear modulus in the
crust region. The shear modulus is a property associated with rigidity against
distortion. Therefore, it is mainly determined by the lattice energy in the
crust. In the phase composed of spherical nuclei in the crust, where the
nuclei form a body-centered cubic (bcc) lattice, the shear modulus for zero-
temperature matter has been formulated as a function of the ion number density
and charge of a nucleus (Strohmayer et al., 1991) (see section 2 for details).
However, in deriving the formula of the shear modulus, it is assumed that the
nuclei in the crust are point-like particles. Indeed, this assumption may be
valid for the nuclei in the outer crust, where neutrons are still confined in
the nuclei, whereas the finite-size effect of nuclei must become non-
negligible in the inner crust. Thus, in this paper, we first consider the
shear modulus in the bcc lattice by including the finite-size effect of the
nuclei and then discuss this effect on the frequency of crustal torsional
oscillations in a nonmagnetized neutron star. The so-called pasta phase
composed of nonspherical nuclei also contributes to crustal torsional
oscillations (e.g., Sotani, Iida & Oyamatsu (2018, 2019)), but the thickness
of such a phase is much smaller than that of the phase composed of spherical
nuclei. Therefore, in this paper, we simply consider crustal torsional
oscillations, whose restoring force comes from only the elasticity of the
lattice structure with spherical nuclei. We note that the torsional
oscillations are confined inside the crust region because the shear modulus
vanishes inside the core region, where the matter behaves as a fluid.

##  2 Shear modulus

The elasticity of the neutron star crust is characterized by the shear
modulus. The shear modulus in the bcc lattice composed of spherical nuclei,
which are the main components of the crust region, is approximately described
in the limit of zero temperature as

$\mu_{0}=0.1194\,\frac{n_{i}(Ze)^{2}}{a},$ (1)

where $n_{i}$, $Z$, and $a$ denote the ion number density, charge number of
the ion, and the radius of a Wigner–Seitz cell, respectively, i.e., $4\pi
a^{3}/3=1/n_{i}$ (Strohmayer et al., 1991). We note that this expression has
been derived by averaging the overall direction of infinitesimal displacements
with respect to the shear motion under the assumption that each spherical
nucleus is only a point charge (Ogata & Ichimaru, 1990). In Eq. (1), one can
observe that the shear modulus in the bcc lattice is proportional to the
lattice (or Coulomb) energy. In fact, the shear modulus in two-dimensional
geometry is directly expressed by the Coulomb energy in equilibrium (Pethick &
Potekhin, 1998).

To take into account the size effect of spherical nuclei on the shear modulus,
we assume the situation that the protons distribute uniformly with the density
$n_{\rm p}$ in the region of $r\leq r_{\rm p}$, where $r_{\rm p}$ is the
reference proton radius defined by $4\pi r_{\rm p}^{3}n_{\rm p}/3=Z$. We note
that the neutron radius $r_{\rm n}$ is generally not equivalent to $r_{\rm
p}$, and the difference $r_{\rm n}-r_{\rm p}$ is the neutron skin thickness.
Then, the lattice energy $W_{\rm L}$ of the unit cell within the Wigner–Seitz
approximation for the compressible liquid drop model is written as (Baym,
Bethe & Pethick, 1971)

$W_{\rm L}=-\frac{9Z^{2}e^{2}}{10a}f_{\rm L},$ (2)

where $f_{L}$ is a function of $r_{\rm p}/a$ given as

$f_{\rm L}=1-\frac{1}{3}\left(\frac{r_{\rm p}}{a}\right)^{2}.$ (3)

Here, $W_{\rm L}$ in the limit of $r_{\rm p}=0$, where $f_{L}=1$, becomes the
lattice energy with point-like nuclei, i.e., $f_{L}$ corresponds to a
contribution function of the finite-size effect of spherical nuclei. In Fig.
1, we show the behavior of $f_{\rm L}$ and the volume fraction $w\equiv(r_{\rm
p}/a)^{3}$ as a function of $r_{\rm p}/a$. From this figure, one can observe
that the lattice energy slightly depends on the size of spherical nuclei.
Since the ratio of $r_{\rm p}$ to $a$ (or the volume fraction) generally
increases with density, one can observe the stronger finite-size effect of the
nuclei in the higher-density region. Since the shear modulus is directly
associated with lattice energy as mentioned above, it would be natural to
express the shear modulus including the finite-size effect of spherical
nuclei, $\mu_{f}$, as

$\mu_{f}=\mu_{0}f_{\rm L}.$ (4)

Furthermore, nuclei may deform during the shear motion. If so, one has to take
into account the effect of nuclear deformation on the shear modulus in the
case of the finite-size nuclei. In fact, since the timescale of strong
interactions (associated with the deformation) is much smaller than that of
the Coulomb interaction, nuclei can deform, depending on their position.
However, it is unclear how large the deformation of nuclei is during the shear
motion because the amplitude of shear oscillations is not constrained by the
observations. Thus, in this study, we simply adopt Eq. (4) as the shear
modulus in the case of finite-size nuclei. We note that the phonon
contribution to the shear modulus (Baiko, 2011), the electron screening effect
on the shear modulus (Kobyakov & Pethick, 2013), and the elastic properties of
polycrystals (Kobyakov & Pethick, 2015) have been studied as other effects,
but we simply neglect such effects in this study. We also note that the total
Coulomb energy of a unit cell within the Wigner–Seitz approximation for a
compressible liquid drop model is written as

$E_{\rm Coul}=\frac{3Z^{2}e^{2}}{5r_{\rm p}}+W_{\rm
L}=\frac{3Z^{2}e^{2}}{5r_{\rm p}}f_{\rm C},$ (5)

where $f_{\rm C}$ is a contribution function of the finite-size effect of
spherical nuclei given as

$f_{\rm C}=1-\frac{3}{2}\frac{r_{\rm p}}{a}+\frac{1}{2}\left(\frac{r_{\rm
p}}{a}\right)^{3}.$ (6)

One may think that the shear modulus including the finite-size effect is
constructed with $f_{\rm C}$ instead of $f_{\rm L}$ shown in Eq. (4), where
the dependence of $f_{\rm C}$ on $r_{\rm p}/a$ is stronger than that of
$f_{\rm L}$. In this study, however, we employ $f_{\rm L}$ as shown in Eq.
(4).

Figure 1: The contribution function of the finite-size effect of spherical
nuclei ($f_{\rm L}$) and the volume fraction ($w$) are plotted as functions of
the ratio of the reference proton radius to the radius of a Wigner–Seitz cell,
where the solid and dotted lines correspond to $f_{\rm L}$ and $w$,
respectively.

To examine how the frequencies of torsional oscillations depend on the EOSs,
in this study, we adopt three different EOSs listed in Table 1, i.e., the EOS
based on the relativistic mean-field (Shen EOS) (Shen et al., 1998, 2011), the
EOS with the Skyrme-type interaction (SLy4 EOS) (Douchin & Haensel, 2001), and
the EOS derived by the variational method (Togashi EOS) (Togashi et al.,
2017). The derivation of the reference proton radius $r_{\rm p}$ used in Eq.
(3) for the SLy4 EOS is different from that for the other EOSs. That is, the
SLy4 EOS was constructed within the compressible liquid drop model for non-
uniform matter in the neutron star crust, so $r_{\rm p}$ can be strictly
defined in the model. On the other hand, the Thomas–Fermi approximation was
adopted to construct the Shen and Togashi EOSs for non-uniform matter, where
the diffuseness of the nuclear surface was taken into account. For the Shen
and Togashi EOSs, therefore, $r_{\rm p}$ used in this study is evaluated so
that the root-mean-square proton radius calculated in the Thomas–Fermi
approximation reproduces that in the compressible liquid drop model with
$r_{\rm p}$. In Fig. 2, the ratio of $r_{\rm p}$ to $a$ for different EOSs is
shown as a function of the baryon number density, where the dashed, solid, and
dotted lines denote the results obtained using the Shen, SLy4, and Togashi
EOSs, respectively. The values of $r_{\rm p}$ obtained using the Shen and
Togashi EOSs are reasonable as compared with that obtained using the SLy4 EOS.

We note that the constraint on the $1.4M_{\odot}$ neutron star radius obtained
from GW170817 has excluded the Shen EOS (Annala et al., 2018), but we still
consider this EOS to examine the EOS dependence of the frequencies of
torsional oscillations. In Table 1, we list the EOS parameters, where $K_{0}$
is the incompressibility of the symmetric nuclear matter, $L$ is the density
dependence of nuclear symmetry energy, and $\eta$ is given by
$\eta\equiv(K_{0}L^{2})^{1/3}$ (a suitable combination for expressing low-mass
neutron stars proposed by Sotani et al. (2014)). In addition, we also list the
maximum mass and the transition density from the phase composed of spherical
nuclei to uniform matter, i.e., the density at the bottom of the inner crust.
We note that fission-like instability occurs when the volume fraction $w$
reaches $\sim 1/8$ (Pethick & Ravenhall, 1995), which leads to the phase
transition from the spherical nuclei phase to the so-called pasta phase and
eventually to uniform matter, as the density increases. Even so, since the
thickness of the pasta phase is much smaller than that of the phase composed
of spherical nuclei (Sotani, Iida & Oyamatsu, 2017b), we simply neglect the
existence of the pasta phase in this study. This assumption is expected to be
valid at least for the EOSs adopted in this study for the following reasons:
The crustal torsional oscillations are confined inside the phase of spherical
and cylindrical nuclei, since the slab-like nuclei behave as fluid against
linear perturbation (Pethick & Potekhin, 1998). Furthermore, the frequencies
inside the phase of spherical nuclei are almost the same as those in the phase
of spherical and cylindrical nuclei when one considers the EOS with
$L\,\raisebox{-3.44444pt}{$\stackrel{{\scriptstyle\textstyle>}}{{\sim}}$ }40$
MeV, which was confirmed by Sotani, Iida & Oyamatsu (2018).

Table 1: EOS parameters adopted in this study, the maximum gravitational mass $M_{\rm max}$ for a cold neutron star constructed with the corresponding EOSs, and the density at the bottom edge of the neutron star crust, $n_{sp}$. EOS |  $K_{0}$ (MeV) |  $L$ (MeV) |  $\eta$ (MeV) | $M_{\rm max}/M_{\odot}$ |  $n_{sp}$ (fm-3)
---|---|---|---|---|---
Shen | 281 | 111 | 151 | 2.17 | 0.0602
SLy4 | 230 | 45.9 | 78.5 | 2.05 | 0.0760
Togashi | 245 | 38.7 | 71.6 | 2.21 | 0.0758

In Fig. 2, the ratio of $r_{\rm p}$ to $a$ increases with density, which
becomes more than $\sim 60\%$ at the bottom of the crust independently of the
EOS. One also observes the EOS dependence (or $L$-dependence) in the profile
of $r_{\rm p}/a$, i.e., $r_{\rm p}/a$ for an EOS with larger $L$ tends to be
smaller. This is because the symmetry energy at subnuclear density becomes
smaller for an EOS with a larger $L$, which decreases the charge number of
nuclei. The top panel in Fig. 3 shows the shear modulus with and without the
finite-size effect for the neutron star model constructed with the Togashi
EOS; $\mu_{f}$ becomes smaller than $\mu_{0}$ as the density increases owing
to the finite-size effect. The bottom panel in Fig. 3 shows the relative
deviations of the shear modulus for different EOSs; the deviations amount to
$\sim 15\%$ at around the bottom of the crust. We also note that the EOS
dependence of the relative deviation shown in the bottom panel in Fig. 3 is
very similar to that of $r_{\rm p}/a$, as expected.

Figure 2: The ratio $r_{\rm p}/a$ is shown as a function of the baryon number
density for different EOSs, where the dashed, solid, and dotted lines
correspond to the results obtained using the Shen, SLy4, and Togashi EOSs,
respectively.  Figure 3: Shear moduli with and without the finite-size effect
are respectively shown with the solid and dashed lines as a function of the
baryon number density for the Togashi EOS in the top panel. In the bottom
panel, we show the relative deviation of the shear modulus without the finite-
size effect ($\mu_{0}$) to that with the finite-size effect ($\mu_{f}$) for
different EOSs, where the dashed, solid, and dotted lines correspond to the
results obtained using the Shen, SLy4, and Togashi EOSs, respectively.

##  3 Torsional oscillations

In this study, we simply consider a neutron star without a magnetic field. As
mentioned above, since the torsional oscillations excited in a nonmagnetized
neutron star are confined inside the neutron star crust, we focus only on the
crust region. Spherically symmetric neutron star models are generally
constructed by integrating the Tolman–Oppenheimer–Volkoff (TOV) equations
together with an appropriate EOS from the center up to the surface outward,
for a given central density. On the other hand, to eliminate uncertainties in
the EOS for the higher-density region, which corresponds to the neutron star
core region, in this study, we construct the neutron star crust by integrating
the TOV equations inward from the stellar surface down to the bottom edge of
the crust for selected values of stellar mass and radius, as in Sotani et al.
(2012, 2013a). In such a crust model, we consider the torsional oscillations.

Since the crustal torsional oscillations are of the axial type, the
relativistic Cowling approximation works well, i.e., the metric perturbations
during the oscillations can be omitted. The perturbation equation for the
angular displacement can be derived from the linearized equation of motion
(Schumaker & Thorne, 1983). Since the core region of a neutron star behaves as
a fluid, torsional oscillations are confined only to the crust region. Thus,
the boundary conditions are imposed on the bottom edge of the crust and at the
stellar surface. In practice, we adopt the zero-traction condition at the
bottom edge of the crust and the zero-torque condition at the stellar surface
(Schumaker & Thorne, 1983; Sotani, Kokkotas & Stergioulas, 2007). The concrete
perturbation equation and boundary conditions are shown in Sotani et al.
(2012). Then, the problem to solve becomes an eigenvalue problem.

As shown in Fig. 3, the shear modulus decreases owing to the finite-size
effect. As a result, since the shear velocity given by
$(\mu/(\varepsilon+p))^{1/2}$ with the energy density $\varepsilon$ and
pressure $p$ also decreases, one can expect that the frequencies of torsional
oscillations with the finite-size effect become smaller than those without the
effect. In Fig. 4, we show the relative deviation $\Delta$ given as

$\Delta=\frac{{}_{\ell}f_{n}^{(-)}-{}_{\ell}f_{n}^{(+)}}{{}_{\ell}f_{n}^{(+)}},$
(7)

where ${}_{\ell}f_{n}^{(+)}$ and ${}_{\ell}f_{n}^{(-)}$ respectively
correspond to the frequencies of the $\ell$-th torsional oscillations with and
without the finite-size effect, while $n$ denotes the nodal number in the
eigenfunction. Since the results obtained using the SLy4 EOS are almost the
same as those obtained using the Togashi EOS, hereafter we focus only on the
results obtained using the Shen and Togashi EOSs. From this figure, one can
observe that the frequencies decrease owing to the finite-size effect, as
expected, but the relative deviation $\Delta$ is almost independent of $\ell$.
One also finds that the relative deviation in the fundamental frequencies
$(n=0)$ decreases for the neutron star model with a larger $L$, which is
caused by the $L$-dependence in the shear modulus as shown in Fig. 3, whereas
the relative deviation in overtones does not seem to be very simple because
the frequencies of overtones depend on not only the shear modulus but also the
crust thickness (Hansen & Cioffi, 1980)111We also note that, for strongly
magnetized neutron stars, the frequency modulation due to the fine-size effect
of atomic nuclei as shown in Fig. 4 is much smaller than that due to the
magnetic fields (e.g., Sotani, Kokkotas & Stergioulas (2007); Gabler et al.
(2013)).. However, the fundamental oscillations may be more important, if one
is going to identify the QPO frequencies from a neutron star with crustal
torsional oscillations. In fact, the frequencies observed in magnetar giant
flares, such as SGR 1806-20 or SGR 1900$+$14, range from $\sim 10$ Hz up to
the order of kilohertz, among which the frequencies of about 100 Hz, e.g., 150
Hz in SGR 1806-20 or 155 Hz in SGR 1900+14, are the most suitable ones to be
identified with the fundamental frequencies of the crustal torsional
oscillations for certain values of $\ell$ (see Sotani et al. (2012, 2013a)).
In the following, therefore, we focus on fundamental oscillations.

Figure 4: The relative deviation of frequencies without the finite-size effect
from those with such an effect, defined by Eq. (7), is shown as a function of
$\ell$, where the circles, squares, and diamonds denote the results for the
fundamental oscillations ($n=0$) and the 1st ($n=1$) and 2nd overtones
$(n=2)$, respectively. The left and right panels correspond to the neutron
star model with $1.4M_{\odot}$ and $1.8M_{\odot}$, whereas the top and bottom
panels correspond to the neutron star model constructed from the results
obtained using the Shen and Togashi EOSs, respectively.

Although the relative deviation $\Delta$ is almost independent of $\ell$ as
shown in Fig. 4, the absolute value of the deviation of frequency increases
with $\ell$ because the fundamental frequency itself increases with $\ell$. To
choose the value of $\ell$ above which the finite-size effect is not
negligible, we examine the following two quantities, i.e., $\delta f_{1}$,
which is the deviation between the $\ell$-th and $(\ell-1)$-th fundamental
frequencies without the finite-size effect, and $\delta f_{2}$, which is the
deviation between the $\ell$-th fundamental frequencies without and with the
finite-size effect. Explicitly, they are given by

$\displaystyle\delta f_{1}(\ell)$ $\displaystyle=$
$\displaystyle{}_{\ell}f_{0}^{(-)}-{}_{\ell-1}f_{0}^{(-)},$ (8)
$\displaystyle\delta f_{2}(\ell)$ $\displaystyle=$
$\displaystyle{}_{\ell}f_{0}^{(-)}-{}_{\ell}f_{0}^{(+)}.$ (9)

The resultant values of $\delta f_{1}$ and $\delta f_{2}$ are shown in Fig. 5,
where the filled and open circles respectively correspond to $\delta f_{1}$
and $\delta f_{2}$. As seen in this figure, $\delta f_{2}$ increases with
$\ell$, which means that the finite-size effect becomes more important as
$\ell$ increases. In fact, we show the $\ell$ value above which $\delta
f_{2}/\delta f_{1}$ becomes larger than 0.5 with the vertical dashed line in
Fig. 5 for reference. For the neutron star model constructed using the
Togashi, SLy4, and Shen EOSs, the corresponding $\ell$ values are 12, 13, and
19, respectively. This implies that the finite-size effect becomes more
important for the stellar model with a smaller $L$.

Figure 5: The deviations of the frequencies, $\delta f_{1}$ and $\delta
f_{2}$, are shown as a function of $\ell$, where $\delta f_{1}$ and $\delta
f_{2}$ are given by Eqs. (8) and (9), respectively. The vertical dashed line
denotes the values of $\ell$ at which the ratio of $\delta f_{2}$ to $\delta
f_{1}$, i.e., $\delta f_{2}/\delta f_{1}$, becomes more than 0.5.

##  4 Conclusion

The QPOs observed in the afterglow following magnetar giant flares are
considered to be strongly associated with neutron star oscillations. One of
the possible theoretical explanations of the observations is crustal torsional
oscillations, in which the oscillation frequencies are governed by the shear
modulus. The shear modulus in the bcc lattice composed of spherical nuclei is
well known and formulated under the assumption that the nuclei are point-like
particles. In this study, we focus on the finite-size effect of spherical
nuclei and take it into account the value of the shear modulus. The shear
modulus decreases owing to the finite-size effect, which is more important for
EOSs with a smaller $L$. We also examine the frequencies of crustal torsional
oscillations in a nonmagnetized neutron star to evaluate the finite-size
effect. We find that the finite-size effect becomes more crucial in torsional
oscillations with a larger $\ell$ and in the neutron star models constructed
using the EOSs with a smaller $L$. In practice, if one identifies the observed
QPO frequencies with crustal torsional oscillations, the finite-size effect is
not negligible for QPOs of more than $\sim 100$ Hz. We note that the QPO
frequencies detected in the afterglow of magnetar giant flares can be more
strongly affected by magnetic fields (e.g., Sotani, Kokkotas & Stergioulas
(2007); Gabler et al. (2013)) than by the finite-size effects. Therefore, the
latter effects would be more pronounced if QPOs associated with torsional
oscillations were detected from weakly magnetized neutron stars.

## Acknowledgements

We would like to express our special thanks to H. Shen for providing us with
numerical data on parameters of the Thomas–Fermi calculation. This work was
supported in part by Japan Society for the Promotion of Science (JSPS) KAKENHI
Grant Numbers JP18K13551, JP19KK0354, JP20H04753, JP21H01088, and JP21K13924
and by Pioneering Program of RIKEN for Evolution of Matter in the Universe
(r-EMU).

## Data availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Abbott et al. (2017a) Abbott B. P. et al. (LIGO Scientific Collaboration and Virgo Collaboration), 2017a, Phys. Rev. Lett., 119, 161101
  * Abbott et al. (2017b) Abbott B. P. et al. (LIGO Scientific Collaboration and Virgo Collaboration), 2017b, ApJ, 848, L12
  * Andersson & Kokkotas (1996) Andersson N., Kokkotas K. D., 1996, Phys. Rev. Lett., 677, 4134
  * Andersson & Kokkotas (1998) Andersson N., Kokkotas K. D., 1998, MNRAS, 299, 1059
  * Andersson & Kokkotas (2001) Andersson N., Kokkotas K. D., 2001, Int. J. Mod. Phys. D, 10, 381
  * Annala et al. (2018) Annala E., Gorda T., Kurkela A., Vuorinen A., 2018, Phys. Rev. Lett., 120, 172703
  * Antoniadis et al. (2013) Antoniadis J. et al., 2013, Science, 340, 6131
  * Baiko (2011) Baiko D. A., 2011, MNRAS, 416, 22
  * Baym, Bethe & Pethick (1971) Baym G., Bethe H. A., Pethick C. J., 1971, Nucl. Phys. A, 175, 225
  * Cromartie et al. (2020) Cromartie H. T. et al., 2020, Nat. Astron., 4, 72
  * Demorest et al. (2010) Demorest P., Pennucci T., Ransom S., Roberts M., Hessels J., 2010, Nature, 467, 1081
  * Doneva et al. (2013) Doneva D. D., Gaertig E., Kokkotas K. D., Krüger C., 2013, Phys. Rev. D, 88, 044052
  * Douchin & Haensel (2001) Douchin F., Haensel P., 2001, A&A, 380, 151
  * Gabler et al. (2013) Gabler M., Cerdá-Durán P., Font J. A., Müller E., Stergioulas N., 2013, MNRAS, 430, 1811
  * Gearheart et al. (2011) Gearheart M., Newton W. G., Hooker J., Li B. A., 2011, MNRAS, 418, 2343
  * Hansen & Cioffi (1980) Hansen C., Cioffi D. F., 1980, ApJ, 238, 740
  * Hurley et al. (1999) Hurley K. et al., 1999, ApJ, 510, L111
  * Kobyakov & Pethick (2013) Kobyakov D., Pethick C. J., 2013, Phys. Rev. C, 87, 055803
  * Kobyakov & Pethick (2015) Kobyakov D., Pethick C. J., 2015, MNRAS, 449, L110
  * Kouveliotou et al. (1998) Kouveliotou C. et al., 1998, Nature, 393, 235
  * McDermott, Van Horn, & Hansen (1988) McDermott P. N., Van Horn H. M., Hansen C. J., 1988, ApJ, 325 725
  * Miller et al. (2019) Miller M. C. et al., 2019, ApJ, 887, L24
  * Miller et al. (2021) Miller M. C. et al., 2021, ApJ, 918, L28
  * Morozova et al. (2018) Morozova V., Radice D., Burrows A., Vartanyan D., 2018, ApJ, 861, 10
  * Ogata & Ichimaru (1990) Ogata S., Ichimaru S., 1990, Phys. Rev. A, 42, 4867
  * Passamonti & Andersson (2012) Passamonti A., Andersson N., 2012, MNRAS, 419, 638
  * Pethick & Potekhin (1998) Pethick C. J., Potekhin A. Y., 1998, Phys. Lett. B, 427, 7
  * Pethick & Ravenhall (1995) Pethick C. J., Ravenhall D. G., 1995, Annu. Rev. Nucl. Part. Sci., 45, 429
  * Riley et al. (2019) Riley T. E. et al., 2019, ApJ, 887, L21
  * Riley et al. (2021) Riley T. E. et al., 2021, ApJ, 918, L27
  * Schumaker & Thorne (1983) Schumaker B. L., Thorne K. S., 1983, MNRAS, 203, 457
  * Shapiro & Teukolsky (1983) Shapiro S. L., Teukolsky S. A., in Black Holes, White Dwarfs, and Neutron Stars (A Wiley-Interscience Publication, 1983)
  * Shen et al. (1998) Shen H., Toki H., Oyamatsu K., Sumiyoshi K., 1998, Nucl. Phys. A, 637, 435
  * Shen et al. (2011) Shen H., Toki H., Oyamatsu K., Sumiyoshi K., 2011, ApJS, 197, 20
  * Sotani (2014) Sotani H., 2014, Phys. Lett. B, 730, 166
  * Sotani (2021) Sotani H., 2021, Phys. Rev. D, 103, 123015
  * Sotani & Kumar (2021) Sotani H., Kumar B., 2021, Phys. Rev. D, 104, 123002
  * Sotani & Dohi (2022) Sotani H., Dohi A., 2022, Phys. Rev. D, 105, 023007
  * Sotani, Iida & Oyamatsu (2016) Sotani H., Iida K., Oyamatsu K., 2016, New Astron., 43, 80
  * Sotani, Iida & Oyamatsu (2017a) Sotani H., Iida K., Oyamatsu K., 2017a, MNRAS, 464, 3101
  * Sotani, Iida & Oyamatsu (2017b) Sotani H., Iida K., Oyamatsu K., 2017b, MNRAS, 470, 4397
  * Sotani, Iida & Oyamatsu (2018) Sotani H., Iida K., Oyamatsu K., 2018, MNRAS, 479, 4735
  * Sotani, Iida & Oyamatsu (2019) Sotani H., Iida K., Oyamatsu K., 2019, MNRAS, 489, 3022
  * Sotani et al. (2014) Sotani H., Iida K., Oyamatsu K., Ohnishi A., 2014, Prog. Theor. Exp. Phys. 2014, 051E01
  * Sotani, Kohri & Harada (2004) Sotani H., Kohri K., Harada T., 2004, Phys. Rev. D, 69, 084008
  * Sotani, Kokkotas & Stergioulas (2007) Sotani H., Kokkotas K. D., Stergioulas N., 2007, MNRAS, 375, 261
  * Sotani et al. (2019) Sotani H., Kuroda T., Takiwaki T., Kotake K., 2019, Phys. Rev. D, 99, 123024
  * Sotani et al. (2012) Sotani H., Nakazato K., Iida K., Oyamatsu K., 2012, Phys. Rev. Lett., 108, 201101
  * Sotani et al. (2013a) Sotani H., Nakazato K., Iida K., Oyamatsu K., 2013a, MNRAS, 428, L21
  * Sotani et al. (2013b) Sotani H., Nakazato K., Iida K., Oyamatsu K., 2013b, MNRAS, 434, 2060
  * Sotani & Sumiyoshi (2019) Sotani H., Sumiyoshi K., 2019, Phys. Rev. D, 100, 083008
  * Sotani, Takiwaki & Togashi (2021) Sotani H., Takiwaki T., Togashi H., 2021, Phys. Rev. D, 104, 123009
  * Sotani, Tominaga & Maeda (2001) Sotani H., Tominaga K., Maeda K. I., 2001, Phys. Rev. D, 65, 024010
  * Sotani et al. (2011) Sotani H., Yasutake N., Maruyama T., Tatsumi T., 2011, Phys. Rev. D, 83, 024014
  * Steiner & Watts (2009) Steiner A. W., Watts A. L., 2009, Phys. Rev. Lett., 103, 181101
  * Strohmayer et al. (1991) Strohmayer T., van Horn H. M., Ogata S., Iyetomi H., Ichimaru S., 1991, ApJ, 375, 679
  * Strohmayer & Watts (2005) Strohmayer T. E., Watts A. L., 2005, ApJ, 632, L111
  * Strohmayer & Watts (2006) Strohmayer T. E., Watts A. L., 2006, ApJ, 653, 593
  * Togashi et al. (2017) Togashi H., Nakazato K., Takehara Y., Yamamuro S., Suzuki H., Takano M., 2017, Nucl. Phys. A, 961, 78
  * Torres-Forné et al. (2019) Torres-Forné A., Cerdá-Durán P., Passamonti A., Obergaulinger M., Font J. A., 2019, MNRAS, 482, 3967


# Scalar perturbations of Galileon cosmologies in the mechanical approach
within the late Universe

Jan Novák<EMAIL_ADDRESS>Physics Department, Technical University
of Liberec, Studentská 1402/2, 461 17, Liberec, Czech republic

###### Abstract

We investigate the Universe at the late stage of its evolution and inside the
cell of uniformity 150-370 Mpc. We consider the Universe to be filled at these
scales with dust like matter, a minimally coupled Galileon field and
radiation. We use the mechanical approach. Therefore, the peculiar velocities
of the inhomogeneities as well as the fluctuations of the other perfect fluids
can be considered nonrelativistic. Such fluids are said to be coupled because
they are concentrated around the inhomogeneities. We investigate the
conditions under which the physical Galileon field, i.e. compatible with
results of the latest gravitational wave experiments GW150914, GW151226,
GW170104, GW170814, GW170817 and GW170608, can become coupled. We know that at
the background level coupled scalar fields behave as a two-component perfect
fluids: one which mimics a network of frustrated cosmic string and an another
one which corresponds to an effective cosmological constant. We found a
correction for the Galileon field, which behaves like a matter component.

Galileon cosmology, mechanical approach, scalar perturbations, late-time
acceleration

###### pacs:

Valid PACS appear here

##  I Introduction

The $\Lambda$CDM model is consistent with observational data but the energy
scale of dark energy is too low Weinberg ; Starobinsky . Therefore, this
cosmological constant is not compatible with the cosmological constant
originated from the vacuum energy in a quantum field theory. Already many
models have been proposed to explain the present accelerated expansion of the
Universe. If the history of particle physics is any guide, then one can assume
that the dark energy is due to a new field. Within this setup, the most
popular ones are the models based on scalar fields, Amendola .

Dark energy models based on minimaly coupled scalar field are named
quintessence Lopez . In this case equation of state (EoS) parameter
$w\in(-1,0)$. It is as well possible to have phantom scalar fields where
$w<-1$. Finally quintom models are based on scalar fields, where there is a
crossing of the value $w=-1$. In addition, a dynamical EoS parameter, $w$, can
help to solve the coincidence problem.

Another possibility to describe the late-time acceleration of the Universe is
to modify the law of gravity from general relativity (GR) at large distances
Faraoni ; Capozziello . There have been two main approaches. The first one is
to introduce a Lagrangian for gravity constructed from the Ricci, Riemann and
metric tensors. Examples are $f(R)$ theories or Gauss-Bonnet gravities. The
second possibility is to build up the laws of gravity from higher dimensional
models that realize the cosmic acceleration through gravitational leakage to
extra dimensions. The DGP (Dvali-Gabadadze-Poratti) model belongs to this
class111 Modified gravity models as a mean to explain the late-time
acceleration of the universe need to be constructed to recover GR behavior in
the regions of high density for the consistency with local gravity
experiments. There are two possible ways to do it: the chameleon mechanism and
Vainshtein mechanism, Ref. Khoury ; Babichev ;  Chimento .

Mostly inspired by the DGP model, the authors of Ref. Rattazzi  derived the
five Lagrangians that lead to field equations invariant under the Galileon
symmetry $\partial_{\mu}\phi\rightarrow\partial_{\mu}\phi+b_{\mu}$ (the
vectorial parameter $b_{\mu}$ corresponds to a constant shift) in a Minkowski
spacetime:

$\displaystyle L_{1}=M^{3}\phi,\ L_{2}=(\nabla\phi)^{2},\
L_{3}=(\square\phi)(\nabla\phi)^{2}/M^{3},$ $\displaystyle
L_{4}=(\nabla\phi)^{2}\big{[}2(\nabla\phi)^{2}-2\phi_{;\mu\nu}\phi^{;\mu\nu}-R(\nabla\phi)^{2}/2\big{]}/M^{6},$
$\displaystyle
L_{5}=(\nabla\phi)^{2}\big{[}(\nabla\phi)^{3}-3\square\phi\phi_{;\mu\nu}\phi^{;\mu\nu}+2\phi_{;\mu}^{\
\nu}\phi_{;\nu}^{\ \rho}\phi_{;\rho}^{\
\mu}-6\phi_{;\mu}\phi^{;\mu\nu}\phi^{;\rho}G_{\nu\rho}\big{]}/M^{9},M^{3}=M_{Pl}H_{0}^{2}$

The scalar field $\phi$ that respects the Galileon symmetry is named the
Galileon field Felice . Each of the five terms give origin to a second order
differential equation of motion, which keep the theory free from unstable
spin-2 ghost degrees of freedom. If we carry the analysis in a curved
spacetime, we need to change the Lagrangians to their covariant form222Some of
the terms $L_{i}$ are already disfavoured by the results of the experiment
GW170817, Miguel , but the term which we will consider in our work is fully
viable..

It is of great importance to suggest a mechanism which can verify the
viability of Galileon models. The theory of cosmological perturbations is a
powerful tool to investigate cosmological models Mukhanov . In fact, we will
study the Universe at the late stages of its evolution and deep inside the
cell of uniformity. In this setup, there are discrete inhomogeneities in this
cell, because galaxies, group and clusters of galaxies were already formed and
can be considered as discrete sources for the gravitational potential. It was
shown in previous works that in this case the mechanical approach Eingorn ;
Zhuk ; Alexander  is a powerful tool to study the scalar perturbations. It
enables us to get the gravitational potential and to describe the motions of
galaxies.

The hydrodynamical approach is a good tool to investigate the growth of
structure of the early universe. It works well in the linear approximation.
However, it became inapplicable at the strongly nonlinear regime. It starts
for $z$ of few dozens. However, on much bigger scales matter becomes on
average homogeneous and isotropic with matter in the form of the perfect
fluid. It is important to define theoretically at which scales we should
define the transition from the highly inhomogeneous mechanical distribution to
the smooth hydrodynamical one.

The Universe is filled with inhomogeneously distrubuted discrete structures at
the scale 150-370 MPc. The mechanical approach enables us to obtain the
gravitational potential for an arbitrary number of randomly distributed
inhomogeneities. We can investigate the relative motion of galaxies and the
formation of the Hubble flow with the expression for the gravitational
potential. The mechanical approach works well for the $\Lambda$CDM model,
where the peculiar velocities of the inhomogeneities could be considered as
negligibly small when we compare them with the speed of light. Additionally,
we consider scales deep inside the cell of uniformity. Consequently, we can
drop the peculiar velocities at first order of approximation. Such models were
generalized also to the case of cosmologies with different perfect fluids,
which can play the role of dark energy and dark matter Burgazli ; Brilenkov ;
Akarsu ; Morais ; Perfect .

The fluctuations of these additional perfect fluids form also their own
inhomogeneities. It is supposed in the mechanical approach that the velocities
of displacement of such inhomogeneities are of the order of the peculiar
velocities of the inhomogeneities corresponding to dust like matter. These
types of inhomogeneities are coupled to each other in the sense mentioned in
Kumar . This means that for the considered models, we investigate the possible
existence of such coupled fluids. They can as well play the role of dark
matter as shown in Perfect .

In the present paper, we consider a cosmological model with a Galileon field
minimally coupled to gravity. The Universe is also filled with dust like
matter and radiation. We study the theory of scalar perturbations, Ref.
Albarran , for such models and obtain a condition under which the
inhomogeneities of the dust like matter and the inhomogeneities of the scalar
field can be coupled to each other. We show that this condition imposes rather
strong restrictions on the scalar field itself. The coupled scalar field
behaves at the background level as a three component perfect fluid: a
cosmological constant, a term which mimics a network of frustrated cosmic
strings and a further component which behaves as matter. Though the
gravitational wave experiments constrain the Galileon cosmologies, here we
provide a further tool based on the mechanical approach to choose the viable
theories.

The work is structered followingly: we present the action, which we consider
in this article, in the first section; We study the mechanical approach in the
next part III. And we present the results in the final part IV.

##  II Summary of the Galileon models

Given that non-gravitational coupling between the Galileon field and gravity
are physically allowed, we will simply consider the following self-interaction
term for the Galileon field:

$S_{I}=\alpha\int_{M}\sqrt{|g|}\square\phi\
\partial_{\mu}\phi\partial^{\mu}\phi\ d^{4}x,$

On the previous expression, $\alpha$ is a small parameter, which measure the
deviation from the model of a minimally coupled scalar field and has units of
volume. We intend to demonstrate the typical behavior for Galileon models on
this Lagrangian within the mechanical approach. First of all, we will obtain
the tensor of energy momentum for this Lagrangian as:

$\displaystyle T_{\mu\nu}=\frac{2}{\sqrt{-g}}\frac{\delta S_{I}}{\delta
g^{\mu\nu}},$ (1)

where $\frac{\delta S}{\delta g^{\mu\nu}}$ is the variational derivative of
the action. We then get Quiros ,

$\displaystyle T_{\mu\nu}=\alpha\ 2\partial_{\mu}\phi\ \partial_{\nu}\phi\
\square\phi-2\alpha\ \phi_{;(\mu}\ (\partial_{\rho}\phi\
\partial^{\rho}\phi)_{;\nu)}+\ \alpha\ g_{\mu\nu}\ \phi_{;\lambda}\
(\partial_{\rho}\phi\ \partial^{\rho}\phi)^{;\lambda}.$ (2)

We would like to use the following perturbed metric in our computations
according to Kumar :

$\displaystyle ds^{2}=a^{2}\big{[}(1+2\Phi)\ d\eta^{2}-(1-2\Psi)\
\gamma_{ij}dx^{i}dx^{j}\big{]}$ (3)

For simplicity we assume spatially flat solutions
$\gamma_{ij}=\gamma^{ij}={\hbox{1\kern-1.49994pt\vrule
width=0.8pt,depth=-0.5pt}}_{3\times 3}$, therefore
$\gamma=\gamma_{ij}\gamma^{ij}=3$.

##  III Mechanical approach

Now we can write the tensor of energy-momentum for whole action, when we
include also the minimally coupled scalar field:

$\displaystyle
S=\int_{M}\sqrt{|g|}\big{(}\frac{1}{2}\partial_{\rho}\phi\partial^{\rho}\phi-V(\phi)+\alpha\square\phi\partial_{\mu}\phi\partial^{\mu}\phi\big{)}\
d^{4}x.$

The characteristic term in Galileon cosmologies is the term $\square\phi$, so
we need to compute $\square\phi$ up to first order on the perturbative setup.
We use the notation $\phi=\phi_{c}+\varphi$, where $\phi_{c}$ is the
background term and $\varphi$ the perturbed quantity (all derivatives will be
with respect to the conformal time):

$\displaystyle\square\phi=\frac{\phi_{c}^{\prime\prime}}{a^{2}}-\frac{a^{\prime}\phi_{c}^{\prime}}{a^{3}}+\frac{a^{\prime}\phi_{c}^{\prime}}{a^{3}}\gamma+\big{[}-\frac{2}{a^{2}}\phi_{c}^{\prime\prime}\Phi+\frac{\varphi^{\prime\prime}}{a^{2}}-\frac{\Delta_{ij}\varphi}{a^{2}}-\frac{\gamma}{a^{2}}\phi_{c}^{\prime}(\Psi^{\prime}+2\Phi\frac{a^{\prime}}{a})-$
$\displaystyle-\frac{a^{\prime}}{a}\frac{\gamma}{a^{2}}\varphi^{\prime}-\frac{\varphi^{\prime}a^{\prime}}{a^{3}}+\frac{4\Phi\phi_{c}^{\prime}a^{\prime}}{a^{3}}-\frac{\phi_{c}^{\prime}}{a^{3}}(\Phi^{\prime}a+2a^{\prime}\Phi)-\square_{ij}\varphi\frac{1}{a^{2}}\big{]},$
(4)

where we have defined $\Delta_{ij}\varphi\equiv\varphi_{,ij}\gamma^{ij}$ and
$-\varphi_{,k}\Gamma^{k}_{ij}\gamma^{ij}=\square_{ij}\varphi-\Delta_{ij}\varphi$.
So, we can finally compute $T^{0}_{0}$ for $S_{I}$:

$\displaystyle\bar{T}^{0}_{0}+\delta T^{0}_{\ 0}=2\alpha\
\frac{1}{a^{2}}(\phi_{c}^{\prime})^{2}\big{(}\frac{\phi_{c}^{\prime\prime}}{a^{2}}+\frac{\gamma}{a^{2}}\phi_{c}^{\prime}\frac{a^{\prime}}{a}-\frac{a^{\prime}\phi_{c}^{\prime}}{a^{3}}\big{)}-2\alpha\frac{1}{a^{4}}(\phi_{c}^{\prime})^{2}.(\phi_{c}^{\prime\prime}-2\phi_{c}^{\prime}H+\phi_{c}^{\prime}H)+$
$\displaystyle+2\alpha\frac{(\phi_{c}^{\prime})^{2}}{a^{2}}\big{[}-\frac{2}{a^{2}}\phi_{c}^{\prime\prime}\Phi+\frac{\varphi^{\prime\prime}}{a^{2}}-\frac{\Delta_{ij}\varphi}{a^{2}}+\frac{\gamma}{a^{2}}(-\phi_{c}^{\prime}\Psi^{\prime}-2\Phi\frac{a^{\prime}}{a}\phi_{c}^{\prime}-\frac{\varphi^{\prime}a^{\prime}}{a})-\frac{\varphi^{\prime}a^{\prime}}{a^{3}}+\frac{4\Phi\phi_{c}^{\prime}a^{\prime}}{a^{3}}-\frac{\phi_{c}^{\prime}}{a^{3}}(\Phi^{\prime}a+2a^{\prime}\Phi)+$
$\displaystyle+\frac{1}{a^{2}}(\square_{ij}-\Delta_{ij})\varphi\big{]}-4\alpha\frac{\Phi}{a^{2}}(\phi_{c}^{\prime})^{2}(\frac{\phi_{c}^{\prime\prime}}{a^{2}}-\gamma\frac{\phi_{c}^{\prime}a^{\prime}}{a^{3}}-\frac{a^{\prime}\phi_{c}^{\prime}}{a^{3}})+4\alpha\varphi^{\prime}\frac{\phi_{c}^{\prime}}{a^{2}}(\frac{\phi_{c}^{\prime\prime}}{a^{2}}-\gamma\frac{\phi_{c}^{\prime}a^{\prime}}{a^{3}}-\frac{a^{\prime}\phi_{c}^{\prime}}{a^{3}})+4\alpha\frac{\Phi(\phi_{c}^{\prime})^{2}\phi_{c}^{\prime\prime}}{a^{4}}-$
$\displaystyle-2\alpha\frac{\varphi^{\prime}\phi_{c}^{\prime}\phi_{c}^{\prime\prime}}{a^{4}}-2\alpha\frac{\phi_{c}^{\prime}\phi_{c}^{\prime\prime}\varphi^{\prime}}{a^{4}}-2\alpha\frac{(\phi_{c}^{\prime})^{2}\varphi^{\prime\prime}}{a^{4}}+4\alpha\frac{(\phi_{c}^{\prime})^{2}\phi_{c}^{\prime\prime}\Phi}{a^{4}}-4\alpha\frac{a^{\prime}(\phi_{c}^{\prime})^{3}\Phi}{a^{5}}+2\alpha\frac{a^{\prime}\varphi^{\prime}(\phi_{c}^{\prime})^{2}}{a^{5}}+2\alpha\frac{(\phi_{c}^{\prime})^{3}}{a^{2}}(\frac{\Phi^{\prime}}{a^{2}}-2\frac{\Phi
a^{\prime}}{a^{3}})=$
$\displaystyle=2\alpha\gamma\frac{(\phi_{c}^{\prime})^{2}}{a^{4}}(\phi_{c}^{\prime}\Phi^{\prime}+4\Phi\frac{a^{\prime}}{a}\phi_{c}^{\prime}-3\varphi^{\prime}\frac{a^{\prime}}{a})-2\alpha\frac{(\phi_{c}^{\prime})^{2}}{a^{4}}\square_{ij}\varphi$
(5)

And we must do the variation of the tensor of energy momentum:

$\displaystyle\delta{T}^{0}_{\
0}=2\alpha\frac{(\phi_{c}^{\prime})^{2}}{a^{2}}\big{[}-\frac{2}{a^{2}}\phi_{c}^{\prime\prime}\Phi+\frac{\varphi^{\prime\prime}}{a^{2}}-\frac{\Delta_{ij}\varphi}{a^{2}}+\frac{\gamma}{a^{2}}(-\phi_{c}^{\prime}\Psi^{\prime}-2\Phi\frac{a^{\prime}}{a}\phi_{c}^{\prime}-\frac{\varphi^{\prime}a^{\prime}}{a})-\frac{\varphi^{\prime}a^{\prime}}{a^{3}}+\frac{4\Phi\phi_{c}^{\prime}a^{\prime}}{a^{3}}-\frac{\phi_{c}^{\prime}}{a^{3}}(\Phi^{\prime}a+2a^{\prime}\Phi)+$
$\displaystyle+\frac{1}{a^{2}}(\square_{ij}-\Delta_{ij})\varphi\big{]}-4\alpha\frac{\Phi}{a^{2}}(\phi_{c}^{\prime})^{2}(\frac{\phi_{c}^{\prime\prime}}{a^{2}}-\gamma\frac{\phi_{c}^{\prime}a^{\prime}}{a^{3}}-\frac{a^{\prime}\phi_{c}^{\prime}}{a^{3}})+4\alpha\varphi^{\prime}\frac{\phi_{c}^{\prime}}{a^{2}}(\frac{\phi_{c}^{\prime\prime}}{a^{2}}-\gamma\frac{\phi_{c}^{\prime}a^{\prime}}{a^{3}}-\frac{a^{\prime}\phi_{c}^{\prime}}{a^{3}})+4\alpha\frac{\Phi(\phi_{c}^{\prime})^{2}\phi_{c}^{\prime\prime}}{a^{4}}-$
$\displaystyle-2\alpha\frac{\varphi^{\prime}\phi_{c}^{\prime}\phi_{c}^{\prime\prime}}{a^{4}}-2\alpha\frac{(\phi_{c}^{\prime})^{2}\varphi^{\prime\prime}}{a^{4}}-2\alpha\frac{\phi_{c}^{\prime}\phi_{c}^{\prime\prime}\varphi^{\prime}}{a^{4}}+4\alpha\frac{(\phi_{c}^{\prime})^{2}\phi_{c}^{\prime\prime}\Phi}{a^{4}}-4\alpha\frac{a^{\prime}(\phi_{c}^{\prime})^{3}\Phi}{a^{5}}+2\alpha\frac{a^{\prime}\varphi^{\prime}(\phi_{c}^{\prime})^{2}}{a^{5}}+2\alpha\frac{(\phi_{c}^{\prime})^{3}}{a^{2}}(\frac{\Phi^{\prime}}{a^{2}}-2\frac{\Phi
a^{\prime}}{a^{3}})=$
$\displaystyle=2\alpha\gamma\frac{(\phi_{c}^{\prime})^{2}}{a^{4}}(\phi_{c}^{\prime}\Phi^{\prime}+4\Phi\frac{a^{\prime}}{a}\phi_{c}^{\prime}-3\varphi^{\prime}\frac{a^{\prime}}{a})-2\alpha\frac{(\phi_{c}^{\prime})^{2}}{a^{4}}\square_{ij}\varphi$
$\displaystyle\delta T^{i}_{j}=\alpha
g^{i}_{j}\frac{\phi_{c}^{\prime}}{a^{6}}\big{[}4\phi_{c}^{\prime\prime}\varphi^{\prime}a^{2}+2\varphi^{\prime\prime}\phi_{c}^{\prime}a^{2}-6\phi_{c}^{\prime}\varphi^{\prime}aa^{\prime}-8\Phi\phi_{c}^{\prime}a^{2}\phi_{c}^{\prime\prime}+8\Phi(\phi_{c}^{\prime})^{2}a^{\prime}a-2(\phi_{c}^{\prime})^{2}a^{2}\Phi^{\prime}\big{]}$
$\displaystyle\delta T^{0}_{\ i}=2\alpha\
g^{00}\phi^{\prime}_{c}\phi_{,i}\square\phi-\alpha\
g^{00}\phi^{\prime}_{c}(\phi_{c,\rho}\phi_{,\sigma}g^{\sigma\rho})_{,i}-\alpha\varphi_{,i}g^{00}({\phi_{,\rho}}\phi_{,\sigma}g^{\sigma\rho})_{,0}=$
$\displaystyle=\frac{2}{a^{2}}\alpha\phi^{\prime}_{c}\varphi_{,i}(\frac{\phi_{c}^{\prime\prime}}{a^{2}}-\frac{\gamma}{a^{2}}\phi_{c}^{\prime}\frac{a^{\prime}}{a}-\frac{a^{\prime}\phi_{c}^{\prime}}{a^{3}})-\alpha
g^{00}(\phi_{c}^{\prime}+\varphi^{\prime})(\phi_{,\rho
i}\phi_{,\sigma}g^{\sigma\rho}+\phi_{,\rho}\phi_{,\sigma
i}g^{\sigma\rho}+\phi_{,\rho}\phi_{,\sigma}g^{\sigma\rho}_{\ \ ,i})-$
$\displaystyle-\alpha g^{00}\varphi_{,i}(\phi_{,\rho
0}\phi_{,\sigma}g^{\sigma\rho}+\phi_{,\rho}\phi_{,\sigma
0}g^{\sigma\rho}+\phi_{,\rho}\phi_{,\sigma}g^{\sigma\rho}_{\ \ ,0})=$
$\displaystyle=\frac{2\alpha}{a^{2}}\phi_{c}^{\prime}\varphi_{,i}(\frac{\phi_{c}^{\prime\prime}}{a^{2}}-\frac{\gamma}{a^{2}}\phi_{c}^{\prime}\frac{a^{\prime}}{a}-\frac{a^{\prime}\phi_{c}^{\prime}}{a^{3}})-2\alpha\frac{(\phi_{c}^{\prime})^{2}}{a^{4}}\varphi_{,0i}+2\alpha\frac{\Phi_{,i}}{a^{4}}(\phi_{c}^{\prime})^{3}+2\alpha\frac{\varphi_{,i}}{a^{2}}\big{[}-\phi_{c,00}\frac{\phi_{c,0}}{a^{2}}+(\phi_{c,0})^{2}\frac{a^{\prime}}{a^{3}}\big{]},$
(6)

where $g^{\sigma\rho}_{\ \ ,0}$ is a derivative of contravariant metric with
respect to the conformal time.

So, the pertubed part of the Einstein equations read Kumar :

$\displaystyle\Delta\Phi-3H(\Phi^{\prime}+H\Phi)+3K\Phi=\frac{\kappa}{2}a^{2}(\delta\epsilon_{dust}+\delta\epsilon_{rad})+\frac{\kappa}{2}\big{[}-(\phi_{c}^{\prime})^{2}\Phi+\phi_{c}^{\prime}\varphi^{\prime}+a^{2}\frac{dV}{d\phi}(\phi_{c})\varphi+$
$\displaystyle+2\alpha\gamma\frac{(\phi_{c}^{\prime})^{2}}{a^{2}}(\phi_{c}^{\prime}\Phi^{\prime}-3\varphi^{\prime}\frac{a^{\prime}}{a})-2\alpha\frac{(\phi_{c}^{\prime})^{2}}{a^{2}}\square_{ij}\varphi\big{]},$
(7)

$\displaystyle\partial_{i}\Phi^{\prime}+H\partial_{i}\Phi=\frac{\kappa}{2}\big{[}\phi_{c}^{\prime}\partial_{i}\varphi+2\alpha\phi_{c}^{\prime}\partial_{i}\varphi(\frac{\phi_{c}^{\prime\prime}}{a^{2}}-\frac{\gamma}{a^{2}}\phi_{c}^{\prime}\frac{a^{\prime}}{a}-\frac{a^{\prime}}{a}\phi_{c}^{\prime})-2\alpha\frac{(\phi_{c}^{\prime})^{2}}{a^{2}}\varphi_{,0i}+\frac{2\alpha\Phi_{,i}}{a^{2}}(\phi_{c}^{\prime})^{3}+2\alpha\frac{\varphi_{,i}}{a^{2}}[-\phi_{c,00}\phi_{c,0}+(\phi_{c,0})^{2}\frac{a^{\prime}}{a}\big{]},$
(8)

$\displaystyle\frac{2}{a^{2}}\big{\\{}\Phi^{\prime\prime}+3H\Phi^{\prime}+\Phi(2\frac{a^{\prime\prime}}{a}-H^{2}-K)\big{\\}}=\kappa\big{\\{}\delta
p_{rad}-\frac{(\phi_{c}^{\prime})^{2}\Phi}{a^{2}}+\frac{\phi_{c}^{\prime}\varphi^{\prime}}{a^{2}}-\frac{dV}{d\phi}\phi_{c}\varphi-\frac{\alpha\phi_{c}^{\prime}}{a^{6}}[4\phi_{c}^{\prime\prime}\varphi^{\prime}a^{2}+2\varphi^{\prime\prime}\phi_{c}^{\prime}a^{2}-6\phi_{c}^{\prime}\varphi^{\prime}aa^{\prime}-$
$\displaystyle-8\Phi\phi_{c}^{\prime}a^{2}\phi_{c}^{\prime\prime}+8\Phi(\phi_{c}^{\prime})^{2}a^{\prime}a-2(\phi_{c}^{\prime})^{2}a^{2}\Phi^{\prime}]\big{\\}}.$
(9)

We obtain the field equation for the field $\phi$ by carrying the variation of
(II) and the result reads:

$\displaystyle-2\alpha(\square\phi)^{2}+2\alpha\nabla^{\mu}\nabla^{\nu}\phi\nabla_{\mu}\nabla_{\nu}\phi+2\alpha\nabla^{\mu}\phi\nabla^{\nu}\phi
R_{\mu\nu}-\square\phi-\frac{dV}{d\phi}=0$ (10)

After plugging in the definition of the energy-momentum tensor, we have that

$\displaystyle\bar{\epsilon}_{\varphi}=-8\alpha\frac{(\phi_{c}^{\prime})^{3}a^{\prime}}{a^{5}}+\frac{2\alpha\phi_{c}^{\prime\prime}}{a^{4}}(\phi_{c}^{\prime})^{2}+\frac{1}{2a^{2}}(\phi_{c}^{\prime})^{2}+V(\phi_{c}),$
(11)

$\displaystyle\bar{p}_{\varphi}=-6\alpha\frac{\phi_{c}^{\prime\prime}\phi_{c}^{\prime
2}}{a^{4}}+6\alpha\frac{a^{\prime}(\phi_{c}^{\prime})^{3}}{a^{5}}+\frac{1}{2a^{2}}(\phi_{c}^{\prime})^{2}-V(\phi_{c}),$
(12)

where $\bar{\epsilon}_{\varphi}$ and $\bar{p}_{\varphi}$ are the energy
density of the Galileon field at the background level.

Consequently, Friedmann and Raychaudhuri equations read:

$H^{2}=\frac{\kappa
a^{2}}{3}\big{[}\bar{\epsilon}_{DUST}+\bar{\epsilon}_{RAD}+\frac{1}{2a^{2}}(\phi_{c}^{\prime})^{2}+V(\phi_{c})-8\frac{\alpha
a^{{}^{\prime}}(\phi_{c}^{\prime})^{3}}{a^{5}}+2\alpha(\phi_{c}^{\prime})^{2}\frac{\phi_{c}^{\prime\prime}}{a^{4}}\big{]}-K$
(13)

$H^{\prime}=-\frac{\kappa
a^{2}}{6}\big{[}\bar{\epsilon}_{DUST}+2\bar{\epsilon}_{RAD}-16\alpha\frac{(\phi_{c}^{\prime})^{2}\phi_{c}^{\prime\prime}}{a^{4}}+10\alpha\frac{a^{{}^{\prime}}}{a^{5}}(\phi_{c}^{\prime})^{3}+2\frac{(\phi_{c}^{\prime})^{2}}{a^{2}}-2V(\phi_{c})\big{]}$
(14)

We continue now in two steps: we compute the derivative of
$\Phi^{\prime}+H\Phi$ and we plug the second Einstein equation (8) to the
third Einstein equation (9):

$\Phi^{\prime}+H\Phi=\frac{\kappa}{2}\big{\\{}\phi_{c}^{\prime}\varphi+2\alpha\phi_{c}^{\prime}\varphi(-\frac{\gamma\phi_{c}^{\prime}a^{\prime}}{a^{3}})-2\frac{\alpha(\phi_{c}^{\prime})^{2}a^{\prime}\varphi}{a^{3}}-2\alpha\frac{(\phi_{c}^{\prime})^{2}}{a^{2}}\varphi^{\prime}+2\alpha\frac{\Phi}{a^{2}}(\phi_{c}^{\prime})^{3}+2\alpha\frac{\varphi}{a^{2}}[(\phi_{c}^{\prime})^{2}\frac{a^{\prime}}{a}]\big{\\}}$
(15)

$\displaystyle\Phi^{\prime\prime}+H^{\prime}\Phi+H\Phi^{\prime}=\frac{\kappa}{2a^{4}}\\{\phi_{c}^{\prime\prime}a^{4}\varphi+\phi_{c}^{\prime}a^{4}\varphi^{\prime}-12\alpha\phi_{c}^{\prime\prime}\phi_{c}^{\prime}aa^{\prime}\varphi-2\alpha(\phi_{c}^{\prime})^{2}\varphi^{\prime}a^{\prime}a-6\alpha(\phi_{c}^{\prime})^{2}a^{\prime\prime}a\varphi+$
$\displaystyle+18\alpha(\phi_{c}^{\prime})^{2}(a^{\prime})^{2}\varphi-2\alpha(\phi_{c}^{\prime})^{2}\varphi^{\prime\prime}a^{2}-4\alpha\phi_{c}^{\prime}\phi_{c}^{\prime\prime}\varphi^{\prime}a^{2}+2\alpha\Phi^{\prime}(\phi_{c}^{\prime})^{3}a^{2}+6\alpha\Phi(\phi_{c}^{\prime})^{2}\phi_{c}^{\prime\prime}a^{2}-4\alpha\Phi(\phi_{c}^{\prime})^{3}a^{\prime}a\\}$
(16)

So after plugging to the third Einstein equation (9):

$\displaystyle\frac{2}{a^{2}}\big{\\{}\frac{\Phi
a^{2}}{3}(\frac{-3}{2}\bar{\epsilon}_{DUST}-2\bar{\epsilon}_{RAD}+6\alpha\frac{(\phi_{c}^{\prime})^{2}\phi_{c}^{\prime\prime}}{a^{4}}+3\alpha\frac{a^{\prime}(\phi_{c}^{\prime})^{3}}{a^{5}}-\frac{3}{2}\frac{(\phi_{c}^{\prime})^{2}}{a^{2}})+$
$\displaystyle+H[\phi_{c}^{\prime}\varphi-6\alpha\phi_{c}^{\prime}\varphi\frac{\phi_{c}^{\prime}a^{\prime}}{a^{3}}-2\alpha\frac{(\phi_{c}^{\prime})^{2}\varphi^{\prime}}{a^{2}}+2\alpha\Phi\frac{(\phi_{c}^{\prime})^{3}}{a^{2}}]+\frac{1}{2}\big{[}\frac{1}{a^{4}}(\phi_{c}^{\prime\prime}a^{4}\varphi+\phi_{c}^{\prime}\varphi^{\prime}a^{4}-12\alpha(\phi_{c}^{\prime\prime})a\phi_{c}^{\prime}a^{\prime}\varphi-$
$\displaystyle-2\alpha(\phi_{c}^{\prime})^{2}\varphi^{\prime}a^{\prime}a-4\alpha\phi_{c}^{\prime}\phi_{c}^{\prime\prime}a^{2}\varphi^{\prime}-6\alpha(\phi_{c}^{\prime})^{2}a^{\prime\prime}a\varphi+18\alpha(\phi_{c}^{\prime})^{2}(a^{\prime})^{2}\varphi-2\alpha(\phi_{c}^{\prime})^{2}\varphi^{\prime\prime}a^{2}+2\alpha(\phi_{c}^{\prime})^{3}a^{2}+$
$\displaystyle\big{[}-H\Phi+\frac{\kappa}{2}\big{(}\phi_{c}^{\prime}\varphi-6\alpha\phi_{c}^{\prime}\varphi(\frac{\phi_{c}^{\prime}a^{\prime}}{a^{3}})-2\alpha\frac{(\phi_{c}^{\prime})^{2}}{a^{2}}\varphi^{\prime}+2\alpha\Phi\frac{(\phi_{c}^{\prime})^{3}}{a^{2}})+6\alpha\Phi(\phi_{c}^{\prime})^{2}\phi_{c}^{\prime\prime}a^{2}-4\alpha\Phi(\phi_{c}^{\prime})^{3}a^{\prime}a\big{]}\big{\\}}+(\phi_{c}^{\prime})^{2}\frac{\Phi}{a^{2}}+\frac{dV}{d\phi}\varphi-\frac{\phi_{c}^{\prime}\varphi^{\prime}}{a^{2}}+$
$\displaystyle+\alpha\frac{\phi_{c}^{\prime}}{a^{6}}(4\phi_{c}^{\prime\prime}\varphi^{\prime}a^{2}+2\varphi^{\prime\prime}\phi_{c}^{\prime}a^{2}-6\phi_{c}^{\prime}\varphi^{\prime}aa^{\prime}-8\Phi\phi_{c}^{\prime}\phi_{c}^{\prime\prime}a^{2}+8\Phi(\phi_{c}^{\prime})^{2}a^{\prime}a-$
$\displaystyle-2(\phi_{c}^{\prime})^{2}a^{2}\big{\\{}-H\Phi+\frac{\kappa}{2}\big{[}\phi_{c}^{\prime}\varphi+2\alpha\phi_{c}^{\prime}\varphi\big{(}-3\frac{\phi_{c}^{\prime}a^{\prime}}{a^{3}}-2\frac{\alpha(\phi_{c}^{\prime})^{2}}{a^{2}}\varphi^{\prime}+2\alpha\Phi\frac{(\phi_{c}^{\prime})^{3}}{a^{2}}\big{)}\big{]}\big{\\}}=\frac{\delta\epsilon_{rad}}{3}$
(17)

We further simplify this equation as:

$\displaystyle\frac{\delta\epsilon_{rad}}{3}=-\frac{1}{3a^{6}}[(-6a^{\prime}\phi_{c}^{\prime}a^{3}-3\frac{dV}{d\phi}a^{6}-3\phi_{c}^{\prime\prime}a^{4}+12\alpha\phi_{c}^{\prime\prime}\phi_{c}^{\prime}aa^{\prime}+18\alpha(\phi_{c}^{\prime})^{2}a^{\prime\prime}a-18\alpha(\phi_{c}^{\prime})^{2}(a^{\prime})^{2}\varphi-30\alpha\Phi(\phi_{c}^{\prime})^{3}a^{\prime}a-$
$\displaystyle-6\alpha\Phi(\phi_{c}^{\prime})^{2}\phi_{c}^{\prime\prime}a^{2}+(36\alpha(\phi_{c}^{\prime})^{2}a^{\prime}a)\varphi^{\prime}+3\Phi
a^{6}\bar{\epsilon}_{dust}+4\Phi a^{6}\bar{\epsilon}_{rad}]$ (18)

So, after inserting the equation of motion, we get the following equation:

$\displaystyle\frac{\delta\epsilon_{rad}}{3}=\frac{-1}{3a^{6}}[(36\alpha
H\phi_{c}^{\prime}\phi_{c}^{\prime\prime}+18\alpha
H^{2}(\phi_{c}^{\prime})^{2}a^{2}+18\alpha
H^{\prime}(\phi_{c}^{\prime})^{2}a^{2}++12\alpha\phi_{c}^{\prime\prime}a\phi_{c}^{\prime}a^{\prime}+18\alpha(\phi_{c}^{\prime})^{2}a^{\prime\prime}a-18\alpha(\phi_{c}^{\prime})^{2}(a^{\prime})^{2})\varphi-$
$\displaystyle-30\alpha\Phi(\phi_{c}^{\prime})^{3}a^{\prime}a-6\alpha\Phi(\phi_{c}^{\prime})^{2}\phi_{c}^{\prime\prime}a^{2}+36\alpha(\phi_{c}^{\prime})^{2}a^{\prime}a\varphi^{\prime}+3\Phi
a^{6}\bar{\epsilon}_{dust}+4\Phi a^{6}\bar{\epsilon}_{rad}]$

We can plug this expression into the first Einstein equation (7):

$\displaystyle\frac{2}{\kappa}\Delta\Phi+\frac{6K}{\kappa}\Phi-a^{2}(\frac{\delta\rho
c^{2}}{a^{3}}+3\bar{\rho}\frac{\Phi}{a^{3}})-36\alpha
H\frac{(\phi_{c}^{\prime})^{3}\Phi}{a^{2}}-\frac{6\alpha}{a^{2}}\Phi(\phi_{c}^{\prime})^{2}\phi_{c}^{\prime\prime}+6\frac{\kappa\alpha^{2}(\phi_{c}^{\prime})^{6}\Phi}{a^{4}}+(\phi_{c}^{\prime})^{2}\Phi+3\Phi
a^{2}\bar{\epsilon}_{dust}+4\Phi a^{2}\bar{\epsilon}_{rad}-$
$\displaystyle-\big{[}3H\phi_{c}^{\prime}+72\frac{\alpha}{a^{2}}H\phi_{c}^{\prime}\phi_{c}^{\prime\prime}+36\alpha(\phi_{c}^{\prime})^{2}\frac{a^{\prime\prime}}{a^{3}}+a^{2}\frac{dV}{d\phi}\big{]}\varphi+\big{[}\frac{3\alpha(\phi_{c}^{\prime})^{4}}{a^{2}}-18\frac{\alpha^{2}H(\phi_{c}^{\prime})^{5}\kappa\varphi}{a^{4}}\big{]}+(\frac{60\alpha}{a^{2}}H(\phi_{c}^{\prime})^{2}-\phi_{c}^{\prime})\varphi^{\prime}-$
$\displaystyle-\frac{6\alpha^{2}(\phi_{c}^{\prime})^{5}}{a^{4}}\kappa\varphi^{\prime}+2\alpha\frac{(\phi_{c}^{\prime})^{2}\square_{ij}\varphi}{a^{2}}=0$
(19)

##  IV Results

Similarly to Kumar  we suppose that $\Omega=\Omega(r)$ (where $\Phi=\Omega/r$
and $\Omega$ is a function of $a$ and spatial coordinates), which means that
$\phi_{c}^{\prime}=const.$:

$\phi_{c}(\eta)=\beta\eta+\omega,$

where $\omega$ and $\beta$ are constants. Then we get from EoM (IV) that

$6\alpha\beta^{2}a^{\prime\prime}+2a^{2}a^{\prime}\beta+\frac{dV}{d\phi}a^{5}=0.$

Now

$\displaystyle V(\eta)=\frac{\beta^{2}}{a^{2}}+V_{\infty}+\alpha f(a)$ (20)

and the goal is to obtain the dependence $f(\eta)$ in this previous relation,
when we know the dependence $V(a)=\frac{\beta^{2}}{a^{2}}+V_{\infty}$ for the
pure scalar field, Kumar . So we take the equation

$\displaystyle
6\alpha\beta^{3}\frac{a^{\prime\prime}}{a^{3}}+2\frac{a^{\prime}}{a}\beta^{2}+a^{2}V^{\prime}=0$
(21)

and we take the derivate of the previous equation (20) with respect to $\eta$
and obtain:

$\displaystyle
V^{\prime}a^{2}=-2\beta^{2}\frac{a^{\prime}}{a}+\alpha\frac{df}{da}a^{\prime}a^{2}$
(22)

Now we combine these two equations:

$\displaystyle
6\alpha\beta^{3}\frac{a^{\prime\prime}}{a}+\alpha\frac{df}{da}a^{\prime}a^{4}=0$
(23)

We use that $\frac{a^{\prime\prime}}{a}=H^{\prime}+H^{2}$ and substitute $H$
and $H^{\prime}$ from (13) and (14). We finally obtain that:

$\displaystyle
2\alpha\beta^{3}\kappa[\frac{\bar{\rho}c^{2}}{2a^{3}}+3\frac{\beta^{2}}{2a^{2}}+2V_{\infty}+2\alpha
f(a)-13\alpha\frac{a^{\prime}\beta^{3}}{a^{5}}]+\alpha\frac{df}{da}a^{\prime}a^{2}=0$
(24)

We should solve this equation for $f(a)$ and we use the Taylor expansion. We
must find the solution of the equation (13) for $H$:

$\displaystyle 0=H^{2}+\frac{8\kappa\alpha\beta^{3}}{3a^{2}}H-\frac{\kappa
a^{2}}{3}[\bar{\epsilon}_{dust}+\bar{\epsilon}_{rad}+\frac{3\beta^{2}}{2a^{2}}+\alpha
f(a)+V_{\infty}]$ (25)

This is a quadratic equation. We ignore the negative solution like unphysical
and we obtain for the positive one the following expression:

$\displaystyle
H(\alpha)=-\frac{4\alpha\beta^{3}\kappa}{3a^{2}}+\sqrt{\frac{16}{9}(\frac{\alpha\beta^{3}\kappa}{a^{2}})^{2}+\kappa\frac{a^{2}}{3}(\bar{\epsilon}_{DUST}+\bar{\epsilon}_{RAD}+\frac{3\beta^{2}}{2a^{2}}+\alpha
f(a)+V_{\infty})}$ (26)

We expand this expression in $\alpha$ and we get:

$\displaystyle H(\alpha)|_{\alpha=0}=\sqrt{\frac{\kappa
a^{2}}{3}(\bar{\epsilon}_{DUST}+\bar{\epsilon}_{RAD}+\frac{3\beta^{2}}{a^{2}}+V_{\infty})}+[\frac{-4\kappa\beta^{3}}{3a^{2}}+\frac{\kappa
a^{2}f(a)}{6\sqrt{\frac{\kappa
a^{2}}{3}(\bar{\epsilon}_{DUST}+\bar{\epsilon}_{RAD}+\frac{3\beta^{2}}{2a^{2}}+V_{\infty}}}]+O(\alpha^{2})$
(27)

We obtain for the equation (24) in the first order approximation

$\displaystyle
0=2\alpha\beta^{3}\kappa[\frac{\bar{\rho}c^{2}}{2a^{3}}+\frac{3\beta^{2}}{2a^{2}}+2V_{\infty}]+\alpha\frac{df(a)}{da}\sqrt{\frac{\kappa
a^{2}}{3}(\bar{\epsilon}_{dust}+\bar{\epsilon}_{rad}+\frac{3\beta^{2}}{2a^{2}}+V_{\infty})a^{3}}$
(28)

We see after some simplifications that

$\frac{df(a)}{da}=\frac{-2\beta^{3}\kappa[\bar{\rho}c^{2}{2}+\frac{3\beta^{2}}{2}a+2V_{\infty}a^{3}]}{a^{5}\sqrt{H_{0}^{2}\Omega_{r}+\frac{\kappa\bar{\rho}c^{2}}{3}a+\frac{\kappa\beta^{2}}{2}a^{2}+\frac{V_{\infty}\kappa
a^{4}}{3}}}$ (29)

So $f(a)$ behaves like matter:

$\displaystyle f(a)\sim\frac{1}{a^{3}}$ (30)

##  V Conclusion

We were studying Galileon cosmologies in this work. We used the mechanical
approach, because we worked deep inside the cell of uniformity 150-370 Mpc and
in the late universe. We filled the Universe with minimally coupled Galileon
field and also dust like matter in the form of discrete distributed galaxies
and groups of galaxies. We further included radiation.

All types of inhomogeneities have non-relativistic velocities. Different types
of inhomogeneities do not not run away massively during the universe evolution
in this case. Fluctuations of the energy density of such perfect fluids are
usually concentrated around the inhomogeneities of dust like matter. Therefore
we call those perfect fluids coupled Perfect . They can even screen the
gravitational potential of galaxies Alvina  or they can play a role of dark
matter flattening the rotation curves of dwarf galaxies Brilenkov . In the
present work, we have investigated the possibility for Galileon scalar field
to be coupled with Galaxies in the late universe. For such Galileon scalar
fields to exist, we have shown that they have to meet certain conditions. At
the background level, such scalar field behaves as a three component perfect
fluid: a network, which mimics a behavior of frustrated cosmic strings with
the EoS parameter $w=-1/3$, a cosmological constant and some matter component,
(30). This is the main result of our work.

##  VI Acknowledgement

I would like to thank Mariam Bouhmadi-López for discussions on the topic and
for the possibility to stay at the university UPV/EHU in Bilbao in 2017 and
2018.

## References

  * (1) S. Weinberg, The cosmological constant problem, Rev.Mod.Phys. 61, 1 , 1989, arXiv: 0005265 [astro-ph]
  * (2) V. Sahni, A. Starobinsky, The case for a positive cosmological Lambda-term, Int.J.Mod.Phys. D9: 373-444, 2000, arXiv: 9904398 [astro-ph]
  * (3) L. Amendola, S. Tsujikawa, Dark energy, Cambridge University Press, 2010
  * (4) M. Bouhmadi-López, K. Sravan Kumar, J. Marto, J. Morais and A. Zhuk, K-essence model from the mechanical approach point of view: coupled scalar field and the late cosmic acceleration, JCAP 07, 050, 2016, arXiv: 1605.03212v2 [gr-qc]
  * (5) S. Capozziello, V. Faraoni, Beyond Einstein Gravity: A Survey of Gravitational Theories for Cosmology and Astrophysics, Fundam.Theor.Phys.170, 2010
  * (6) S. Capozziello, M. De Laurentis, Extended Theories of Gravity, Physics Reports 509 (4-5): 167-321, 2011, arXiv: 1108.6266 [gr-qc]
  * (7) M. Bouhmadi-López, L. P. Chimento, Phys. Rev. D 82, 103506, arXiv:1007.4141 [astro-ph.CO]
  * (8) A. Nicolis, R. Rattazzi, E. Trincherini, The galileon as a local modification of gravity, Phys.Rev. D 79, 064036, 2009, arXiv:0811.2197 [hep-th]
  * (9) A. De Felice, S. Tsujikawa, Generalized Galileon cosmology, Phys. Rev. D84: 124029, 2011, arXiv: 1008.4236 [hep-th]
  * (10) J. Khoury, A. Weltman, Chameleon cosmology, Phys.Rev. D 69, 044026, 2004, arXiv: 0309411 [astro-ph]
  * (11) E. Babichev, C. Deffayet, An introduction to the Vainshtein mechanism, Class. Quantum Grav. 30, 184001, 2013, arXiv:1304.7240 [gr-qc]
  * (12) J. M. Ezquiaga, M. Zumalacárregui, Dark energy after GW170817: dead ends and the road ahead, Phys. Rev. Lett. 119, 251304, 2017, arXiv:1710.05901v2 [astro-ph.CO]
  * (13) V. F. Mukhanov, H. A. Feldman, R. H. Brandenberger, Theory of cosmological perturbations, Physics Reports 215, 203, 1992
  * (14) M. Eingorn, First-order cosmological perturbations engendered by point-like masses, Astroph. J. 825, no.2, 84, arXiv: 1509.03835 [gr-qc]
  * (15) M. Eingorn, A. Zhuk, Remarks on mechanical approach to observable Universe, JCAP 05, 024, 2014, arXiv: 1309.4924 [astro-ph.CO]
  * (16) M. Eingorn, A. Zhuk, Hubble flows and gravitational potentials in observable Universe, JCAP 09, 026, 2012, arXiv: 1205.2384 [astro-ph.CO]
  * (17) A. Burgazli, M. Eingorn, A. Zhuk, Rigorous theoretical constraint on constant negative EoS parameter w and its effect for the late Universe, Eur. Phys. J. C 75, 118, 2015
  * (18) M. Brilenkov, M. Eingorn, L. Jenkovszky, A. Zhuk, Scalar perturbations in cosmological models quark nuggets, Eur. Phys. J. C 74, 3011, 2014, arXiv: 1310.4540 [astro-ph.CO]
  * (19) Ö. Akarsu, M. Bouhmadi-López, M. Brilenkov, R. Brilenkov, M. Eingorn, A. Zhuk, Are dark energy models with variable EoS parameter w compatible with the late inhomogeneous Universe?, JCAP 07, 038, 2015, arXiv: 1502.04693 [gr-qc]
  * (20) M. Bouhmadi-López, M. Brilenkov, R. Brilenkov, J. Morais, A. Zhuk, Scalar perturbations in the late Universe: viability of the Chaplyagin gas model, JCAP 12, 037, 2015 [gr-qc]
  * (21) A. Zhuk, Perfect fluids coupled to inhomogeneities in the late Universe, Gravitation and Cosmology, 22, 159, 2016, arXiv: 1601.01939 [gr-qc]
  * (22) A. Burgazli, A. Zhuk, J. Morais, M. Bouhmadi-López, K. Sravan Kumar: Coupled scalar fields in the late Universe, The mechanical approach and the late cosmic acceleration, JCAP 09, 045, 2016, arXiv: 1512:03819 [gr-qc]
  * (23) I. Albarran, M. Bouhmadi-López, J. Morais, Cosmological perturbations in an effective and genuinely phantom dark energy Universe, Phys. Dark Univ. 16, 94-108, 2017, arXiv: 1611.00392 [astro-ph.CO]
  * (24) I. Quiros, R. García-Salcedo, T. Gonzalez, F. Antonio Horta-Rangel, J. Saavedra, Brans-Dicke Galileon and the variational principle, arXiv:1605.00326v2 [gr-qc]
  * (25) M. Eingorn, A. Kudinova, A. Zhuk, Dynamics of astrophysical objects against the cosmological background, JCAP 04, 010, 2013, arXiv:1211.4045 [astro-ph.CO]
  * (26) M. Eingorn, J. Novák, A. Zhuk, arXiv: 1401.5410, f(R) gravity: scalar perturbations in the late Universe, The European physical journal C 74, 3005, 2014, arXiv:1401.5410 [astro-ph.CO]
  * (27) A. Burgazli, M. Eingorn and A. Zhuk, Rigorous theoretical constraint on constant negative EoS parameter $w$ and its effect for the late Universe, Eur. Phys. J. C 75 (2015) 118, arXiv:1301.0418 [astro-ph.CO]
  * (28) M. Brilenkov, M. Eingorn, L. Jenkovszky and A. Zhuk, Scalar perturbations in cosmological models with quark nuggets, Eur. Phys. J. C 74 (2014) 3011 arXiv:1310.4540 [astro-ph.CO]


#  exoplanet: Gradient-based probabilistic inference for exoplanet data &
other astronomical time series

Daniel Foreman-Mackey  Rodrigo Luger  Eric Agol  Thomas Barclay  Luke G. Bouma
Timothy D. Brandt  Ian Czekala  Trevor J. David  Jiayin Dong  Emily A. Gilbert
Tyler A. Gordon  Christina Hedges  Daniel R. Hey  Brett M. Morris  Adrian M.
Price-Whelan  Arjun B. Savel

††margin:  DOI:
https://doi.org/10.21105/joss.03285https://doi.org10.21105/joss.03285 Software
• https://github.com/openjournals/joss-
reviews/issues/3285https://doi.orgReview • https://github.com/exoplanet-
dev/exoplanethttps://doi.orgRepository •
https://doi.org/10.5281/zenodo.5006965https://doi.orgArchive Editor:
https://arfon.orghttps://doi.orgArfon Smith
Reviewers: • https://github.com/grburgesshttps://doi.org@grburgess •
https://github.com/benjaminpopehttps://doi.org@benjaminpope Submitted: 04 May
2021
Published: 22 June 2021 License
Authors of papers retain copyright and release the work under a Creative
Commons Attribution 4.0 International License
(http://creativecommons.org/licenses/by/4.0/https://doi.orgCC BY 4.0).

## Summary

exoplanet is a toolkit for probabilistic modeling of astronomical time series
data, with a focus on observations of exoplanets, using PyMC3 (Salvatier et
al., 2016). PyMC3 is a flexible and high-performance model-building language
and inference engine that scales well to problems with a large number of
parameters. exoplanet extends PyMC3’s modeling language to support many of the
custom functions and probability distributions required when fitting exoplanet
datasets or other astronomical time series.

While it has been used for other applications, such as the study of stellar
variability (e.g., Gillen et al., 2020; Medina et al., 2020), the primary
purpose of exoplanet is the characterization of exoplanets (e.g., Gilbert et
al., 2020; Plavchan et al., 2020) or multiple-star systems (e.g., Czekala et
al., 2021) using time-series photometry, astrometry, and/or radial velocity.
In particular, the typical use case would be to use one or more of these
datasets to place constraints on the physical and orbital parameters of the
system, such as planet mass or orbital period, while simultaneously taking
into account the effects of stellar variability.

## Statement of need

Time-domain astronomy is a priority of the observational astronomical
community, with huge survey datasets currently available and more forthcoming.
Within this research domain, there is significant investment into the
discovery and characterization of exoplanets, planets orbiting stars other
than our Sun. These datasets are large (on the scale of hundreds of thousands
of observations per star from space-based observatories such as _Kepler_ and
_TESS_), and the research questions are becoming more ambitious (in terms of
both the computational cost of the physical models and the flexibility of
these models). The packages in the _exoplanet_ ecosystem are designed to
enable rigorous probabilistic inference with these large datasets and high-
dimensional models by providing a high-performance and well-tested
infrastructure for integrating these models with modern modeling frameworks
such as PyMC3. Since its initial release at the end of 2018, exoplanet has
been widely used, with 64 citations of the Zenodo record (Foreman-Mackey et
al., 2021) so far.

## The _exoplanet_ software ecosystem

Besides the primary exoplanet package, the _exoplanet_ ecosystem of projects
includes several other libraries. This paper describes, and is the primary
reference for, this full suite of packages. The following provides a short
description of each library within this ecosystem and discusses how they are
related.

  * •

exoplanet111https://github.com/exoplanet-dev/exoplanet is the primary library,
and it includes implementations of many special functions required for
exoplanet data analysis. These include the spherical geometry for computing
orbits, some exoplanet-specific distributions for eccentricity (Kipping,
2013a; Van Eylen et al., 2019) and limb darkening (Kipping, 2013b), and
exposure-time integrated limb-darkened transit light curves.

  * •

exoplanet-core222https://github.com/exoplanet-dev/exoplanet-core provides
efficient, well-tested, and differentiable implementations of all of the
exoplanet-specific operations that must be compiled for performance. These
include an efficient solver for Kepler’s equation (based on the algorithm
proposed by Raposo-Pulido & Peláez, 2017) and limb darkened transit light
curves (Agol et al., 2020). Besides the implementation for PyMC3, exoplanet-
core includes implementations in numpy (Harris et al., 2020) and jax (Bradbury
et al., 2018).

  * •

celerite2333https://celerite2.readthedocs.io, is an updated implementation of
the _celerite_ algorithm444https://celerite.readthedocs.io (Foreman-Mackey,
2018; Foreman-Mackey et al., 2017) for scalable Gaussian Process regression
for time series data. Like exoplanet-core, celerite2 includes support for
numpy, jax, and PyMC3, as well as some recent generalizations of the
_celerite_ algorithm (Gordon et al., 2020).

  * •

pymc3-ext555https://github.com/exoplanet-dev/pymc3-ext, includes a set of
helper functions to make PyMC3 more amenable to the typical astronomical data
analysis workflow. For example, it provides a tuning schedule for PyMC3’s
sampler (based on the method used by the Stan project and described by
Carpenter et al., 2017) that provides better performance on models with
correlated parameters.

  * •

rebound-pymc3666https://github.com/exoplanet-dev/rebound-pymc3 provides an
interface between _REBOUND_ (Rein & Liu, 2012), _REBOUNDx_ (Tamayo et al.,
2020), and PyMC3 to enable inference with full N-body orbit integration.

## Documentation & case studies

The main documentation page for the _exoplanet_ libraries lives at
https://docs.exoplanet.codeshttps://doi.orgdocs.exoplanet.codes where it is
hosted on https://readthedocs.orghttps://doi.orgReadTheDocs. The tutorials
included with the documentation are automatically executed on every push or
pull request to the GitHub repository, with the goal of ensuring that the
tutorials are always compatible with the current version of the code. The
celerite2 project has its own documentation page at
https://celerite2.readthedocs.iohttps://doi.orgcelerite2.readthedocs.io, with
tutorials that are similarly automatically executed.

Alongside these documentation pages, there is a parallel “Case Studies”
website at
https://gallery.exoplanet.codeshttps://doi.orggallery.exoplanet.codes that
includes more detailed example use cases for exoplanet and the other libraries
described here. Like the tutorials on the documentation page, these case
studies are automatically executed using GitHub Actions, but at lower cadence
(once a week and when a new release of the exoplanet library is made) since
the runtime is much longer. Figure 1 shows the results of two example case
studies demonstrating some of the potential use cases of the exoplanet
software ecosystem.

Figure 1: Some examples of datasets fit using exoplanet. The full analyses
behind these examples are available on the “Case Studies” page as Jupyter
notebooks. (left) A fit to the light curves of a transiting exoplanet observed
by two different space-based photometric surveys: Kepler and TESS. (right) The
phase-folded radial velocity time series for an exoplanet observed from
different observatories with different instruments, fit simultaneously using
exoplanet.

## Similar tools

There is a rich ecosystem of tooling available for inference with models such
as the ones supported by exoplanet. Each of these tools has its own set of
strengths and limitations and we will not make a detailed comparison here, but
it is worth listing some of these tools and situating exoplanet in this
context.

Some of the most popular tools in this space include (and note that this is
far from a comprehensive list!) EXOFAST (J. Eastman et al., 2013; J. D.
Eastman et al., 2019), radvel (Fulton et al., 2018), juliet (Espinoza et al.,
2019), exostriker (Trifonov, 2019), PYANETI (Barragán et al., 2019),
allesfitter (Günther & Daylan, 2021), and orbitize (Blunt et al., 2020).
Similar tools also exist for modeling observations of eclipsing binary
systems, including JKTEBOP (Southworth et al., 2004), eb (Irwin et al., 2011),
and PHOEBE (Conroy et al., 2020). These packages all focus on providing a
high-level interface for designing models and then executing a fit. In
contrast, exoplanet is designed to be lower level and more conceptually
similar to tools like batman (Kreidberg, 2015), PyTransit (Hannu Parviainen,
2015), ldtk (H. Parviainen & Aigrain, 2015), ellc (Maxted, 2016), starry
(Luger et al., 2019), or Limbdark.jl (Agol et al., 2020), which provide the
building blocks for evaluating the models required for inference with
exoplanet datasets. In fact, several of the higher-level packages listed above
include these lower-level libraries as dependencies, and our hope is that
exoplanet could provide the backend for future high-level libraries.

As emphasized in the title of this paper, the main selling point of exoplanet
when compared to other tools in this space is that it supports differentiation
of all components of the model and is designed to integrate seamlessly with
the aesara (Willard et al., 2021) automatic differentiation framework used by
PyMC3. It is worth noting that aesara was previously known as Theano (Theano
Development Team, 2016), so these names are sometimes used interchangeably in
the PyMC3 or exoplanet documentation777More information about this distinction
is available at https://docs.exoplanet.codes/en/stable/user/theano/. This
allows the use of modern inference algorithms such as No U-Turn Sampling
(Hoffman & Gelman, 2014) or Automatic Differentiation Variational Inference
(Kucukelbir et al., 2017). These algorithms can have some computational and
conceptual advantages over inference methods that do not use gradients,
especially for high-dimensional models. The computation of gradients is also
useful for model optimization; this is necessary when, say, searching for new
exoplanets, mapping out degeneracies or multiple modes of a posterior, or
estimating uncertainties from a Hessian. Care has been taken to provide
gradients which are numerically stable, and more accurate and faster to
evaluate than finite-difference gradients.

## Acknowledgements

We would like to thank the Astronomical Data Group at Flatiron for listening
to every iteration of this project and for providing great feedback every step
of the way.

This research was partially conducted during the _Exostar19_ program at the
_Kavli Institute for Theoretical Physics_ at UC Santa Barbara, which was
supported in part by the National Science Foundation under Grant No. NSF
PHY-1748958.

Besides the software cited above, exoplanet is also built on top of ArviZ
(Kumar et al., 2019) and AstroPy (Astropy Collaboration et al., 2018, 2013).

## References

Agol, E., Luger, R., & Foreman-Mackey, D. (2020). Analytic Planetary Transit
Light Curves and Derivatives for Stars with Polynomial Limb Darkening. _The
Astronomical Journal_ , _159_(3), 123.
https://doi.org/10.3847/1538-3881/ab4fee

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., Günther, H. M.,
Lim, P. L., Crawford, S. M., Conseil, S., Shupe, D. L., Craig, M. W.,
Dencheva, N., Ginsburg, A., VanderPlas, J. T., Bradley, L. D., Pérez-Suárez,
D., de Val-Borro, M., Aldcroft, T. L., Cruz, K. L., Robitaille, T. P.,
Tollerud, E. J., … Astropy Contributors. (2018). The Astropy Project: Building
an Open-science Project and Status of the v2.0 Core Package. _The Astronomical
Journal_ , _156_ , 123. https://doi.org/10.3847/1538-3881/aabc4f

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., Greenfield, P.,
Droettboom, M., Bray, E., Aldcroft, T., Davis, M., Ginsburg, A., Price-Whelan,
A. M., Kerzendorf, W. E., Conley, A., Crighton, N., Barbary, K., Muna, D.,
Ferguson, H., Grollier, F., Parikh, M. M., Nair, P. H., … Streicher, O.
(2013). Astropy: A community Python package for astronomy. _Astronomy &
Astrophysics_, _558_ , A33. https://doi.org/10.1051/0004-6361/201322068

Barragán, O., Gandolfi, D., & Antoniciello, G. (2019). PYANETI: a fast and
powerful software suite for multiplanet radial velocity and transit fitting.
_Monthly Notices of the Royal Astronomical Society_ , _482_ , 1017–1030.
https://doi.org/10.1093/mnras/sty2472

Blunt, S., Wang, J. J., Angelo, I., Ngo, H., Cody, D., De Rosa, R. J., Graham,
J. R., Hirsch, L., Nagpal, V., Nielsen, E. L., Pearce, L., Rice, M., & Tejada,
R. (2020). orbitize!: A Comprehensive Orbit-fitting Software Package for the
High-contrast Imaging Community. _The Astronomical Journal_ , _159_(3), 89.
https://doi.org/10.3847/1538-3881/ab6663

Bradbury, J., Frostig, R., Hawkins, P., Johnson, M. J., Leary, C., Maclaurin,
D., Necula, G., Paszke, A., VanderPlas, J., Wanderman-Milne, S., & Zhang, Q.
(2018). _JAX: Composable transformations of Python+NumPy programs_ (Version
0.2.5) [Computer software]. http://github.com/google/jax

Carpenter, B., Gelman, A., Hoffman, M. D., Lee, D., Goodrich, B., Betancourt,
M., Brubaker, M., Guo, J., Li, P., & Riddell, A. (2017). Stan: A Probabilistic
Programming Language. _Journal of Statistical Software_ , _76_(1), 1–32.
https://doi.org/10.18637/jss.v076.i01

Conroy, K. E., Kochoska, A., Hey, D., Pablo, H., Hambleton, K. M., Jones, D.,
Giammarco, J., Abdul-Masih, M., & Prša, A. (2020). Physics of Eclipsing
Binaries. V. General Framework for Solving the Inverse Problem. _The
Astrophysical Journal Supplement Series_ , _250_(2), 34.
https://doi.org/10.3847/1538-4365/abb4e2

Czekala, I., Ribas, Á., Cuello, N., Chiang, E., Macías, E., Duchêne, G.,
Andrews, S. M., & Espaillat, C. C. (2021). A Coplanar Circumbinary
Protoplanetary Disk in the TWA 3 Triple M Dwarf System. _The Astrophysical
Journal_ , _912_(1), 6. https://doi.org/10.3847/1538-4357/abebe3

Eastman, J. D., Rodriguez, J. E., Agol, E., Stassun, K. G., Beatty, T. G.,
Vanderburg, A., Gaudi, B. S., Collins, K. A., & Luger, R. (2019). EXOFASTv2: A
public, generalized, publication-quality exoplanet modeling code. _arXiv
e-Prints_ , arXiv:1907.09480. http://arxiv.org/abs/1907.09480

Eastman, J., Gaudi, B. S., & Agol, E. (2013). EXOFAST: A Fast Exoplanetary
Fitting Suite in IDL. _Publications of the Astronomical Society of the
Pacific_ , _125_(923), 83. https://doi.org/10.1086/669497

Espinoza, N., Kossakowski, D., & Brahm, R. (2019). juliet: a versatile
modelling tool for transiting and non-transiting exoplanetary systems.
_Monthly Notices of the Royal Astronomical Society_ , _490_(2), 2262–2283.
https://doi.org/10.1093/mnras/stz2688

Foreman-Mackey, D. (2018). Scalable Backpropagation for Gaussian Processes
using Celerite. _Research Notes of the American Astronomical Society_ ,
_2_(1), 31. https://doi.org/10.3847/2515-5172/aaaf6c

Foreman-Mackey, D., Agol, E., Ambikasaran, S., & Angus, R. (2017). Fast and
Scalable Gaussian Process Modeling with Applications to Astronomical Time
Series. _The Astronomical Journal_ , _154_ , 220\.
https://doi.org/10.3847/1538-3881/aa9332

Foreman-Mackey, D., Savel, A., Luger, R., Czekala, I., Agol, E., Price-Whelan,
A., Hedges, C., Gilbert, E., Barclay, T., Bouma, L., & Brandt, T. D. (2021).
_exoplanet-dev/exoplanet_ (Version 0.4.5) [Computer software]. Zenodo.
https://doi.org/10.5281/zenodo.1998447

Fulton, B. J., Petigura, E. A., Blunt, S., & Sinukoff, E. (2018). RadVel: The
Radial Velocity Modeling Toolkit. _Publications of the Astronomical Society of
the Pacific_ , _130_(986), 044504. https://doi.org/10.1088/1538-3873/aaaaa8

Gilbert, E. A., Barclay, T., Schlieder, J. E., Quintana, E. V., Hord, B. J.,
Kostov, V. B., Lopez, E. D., Rowe, J. F., Hoffman, K., Walkowicz, L. M.,
Silverstein, M. L., Rodriguez, J. E., Vanderburg, A., Suissa, G., Airapetian,
V. S., Clement, M. S., Raymond, S. N., Mann, A. W., Kruse, E., … Winters, J.
G. (2020). The First Habitable-zone Earth-sized Planet from TESS. I.
Validation of the TOI-700 System. _The Astronomical Journal_ , _160_(3), 116.
https://doi.org/10.3847/1538-3881/aba4b2

Gillen, E., Briegal, J. T., Hodgkin, S. T., Foreman-Mackey, D., Van Leeuwen,
F., Jackman, J. A. G., McCormac, J., West, R. G., Queloz, D., Bayliss, D.,
Goad, M. R., Watson, C. A., Wheatley, P. J., Belardi, C., Burleigh, M. R.,
Casewell, S. L., Jenkins, J. S., Raynard, L., Smith, A. M. S., … Vines, J. I.
(2020). NGTS clusters survey - I. Rotation in the young benchmark open cluster
Blanco 1. _Monthly Notices of the Royal Astronomical Society_ , _492_(1),
1008–1024. https://doi.org/10.1093/mnras/stz3251

Gordon, T. A., Agol, E., & Foreman-Mackey, D. (2020). A Fast, Two-dimensional
Gaussian Process Method Based on Celerite: Applications to Transiting
Exoplanet Discovery and Characterization. _The Astronomical Journal_ ,
_160_(5), 240. https://doi.org/10.3847/1538-3881/abbc16

Günther, M. N., & Daylan, T. (2021). Allesfitter: Flexible Star and Exoplanet
Inference from Photometry and Radial Velocity. _The Astrophysical Journal
Supplement Series_ , _254_(1), 13. https://doi.org/10.3847/1538-4365/abe70e

Harris, C. R., Millman, K. J., Walt, S. J. van der, Gommers, R., Virtanen, P.,
Cournapeau, D., Wieser, E., Taylor, J., Berg, S., Smith, N. J., Kern, R.,
Picus, M., Hoyer, S., Kerkwijk, M. H. van, Brett, M., Haldane, A., Río, J. F.
del, Wiebe, M., Peterson, P., … Oliphant, T. E. (2020). Array programming with
NumPy. _Nature_ , _585_(7825), 357–362.
https://doi.org/10.1038/s41586-020-2649-2

Hoffman, M. D., & Gelman, A. (2014). The No-U-Turn sampler: adaptively setting
path lengths in Hamiltonian Monte Carlo. _Journal of Machine Learning
Research_ , _15_(1), 1593–1623. https://jmlr.org/papers/v15/hoffman14a.html

Irwin, J. M., Quinn, S. N., Berta, Z. K., Latham, D. W., Torres, G., Burke, C.
J., Charbonneau, D., Dittmann, J., Esquerdo, G. A., Stefanik, R. P., Oksanen,
A., Buchhave, L. A., Nutzman, P., Berlind, P., Calkins, M. L., & Falco, E. E.
(2011). LSPM J1112+7626: Detection of a 41 Day M-dwarf Eclipsing Binary from
the MEarth Transit Survey. _The Astrophysical Journal_ , _742_(2), 123.
https://doi.org/10.1088/0004-637X/742/2/123

Kipping, D. M. (2013a). Parametrizing the exoplanet eccentricity distribution
with the beta distribution. _Monthly Notices of the Royal Astronomical
Society_ , _434_ , L51–L55. https://doi.org/10.1093/mnrasl/slt075

Kipping, D. M. (2013b). Efficient, uninformative sampling of limb darkening
coefficients for two-parameter laws. _Monthly Notices of the Royal
Astronomical Society_ , _435_ , 2152–2160.
https://doi.org/10.1093/mnras/stt1435

Kreidberg, L. (2015). batman: BAsic Transit Model cAlculatioN in Python.
_Publications of the Astronomical Society of the Pacific_ , _127_(957), 1161.
https://doi.org/10.1086/683602

Kucukelbir, A., Tran, D., Ranganath, R., Gelman, A., & Blei, D. M. (2017).
Automatic Differentiation Variational Inference. _Journal of Machine Learning
Research_ , _18_(14), 1–45. http://jmlr.org/papers/v18/16-107.html

Kumar, R., Carroll, C., Hartikainen, A., & Martin, O. A. (2019). ArviZ a
unified library for exploratory analysis of Bayesian models in Python. _The
Journal of Open Source Software_. https://doi.org/10.21105/joss.01143

Luger, R., Agol, E., Foreman-Mackey, D., Fleming, D. P., Lustig-Yaeger, J., &
Deitrick, R. (2019). starry: Analytic Occultation Light Curves. _The
Astronomical Journal_ , _157_ , 64. https://doi.org/10.3847/1538-3881/aae8e5

Maxted, P. F. L. (2016). ellc: A fast, flexible light curve model for detached
eclipsing binary stars and transiting exoplanets. _Astronomy & Astrophysics_,
_591_ , A111. https://doi.org/10.1051/0004-6361/201628579

Medina, A. A., Winters, J. G., Irwin, J. M., & Charbonneau, D. (2020). Flare
Rates, Rotation Periods, and Spectroscopic Activity Indicators of a Volume-
complete Sample of Mid- to Late-M Dwarfs within 15 pc. _The Astrophysical
Journal_ , _905_(2), 107. https://doi.org/10.3847/1538-4357/abc686

Parviainen, Hannu. (2015). PYTRANSIT: fast and easy exoplanet transit
modelling in PYTHON. _Monthly Notices of the Royal Astronomical Society_ ,
_450_(3), 3233–3238. https://doi.org/10.1093/mnras/stv894

Parviainen, H., & Aigrain, S. (2015). LDTK: Limb Darkening Toolkit. _Monthly
Notices of the Royal Astronomical Society_ , _453_(4), 3821–3826.
https://doi.org/10.1093/mnras/stv1857

Plavchan, P., Barclay, T., Gagné, J., Gao, P., Cale, B., Matzko, W., Dragomir,
D., Quinn, S., Feliz, D., Stassun, K., Crossfield, I. J. M., Berardo, D. A.,
Latham, D. W., Tieu, B., Anglada-Escudé, G., Ricker, G., Vanderspek, R.,
Seager, S., Winn, J. N., … Zilberman, P. (2020). A planet within the debris
disk around the pre-main-sequence star AU Microscopii. _Nature_ , _582_(7813),
497–500. https://doi.org/10.1038/s41586-020-2400-z

Raposo-Pulido, V., & Peláez, J. (2017). An efficient code to solve the Kepler
equation. Elliptic case. _Monthly Notices of the Royal Astronomical Society_ ,
_467_(2), 1702–1713. https://doi.org/10.1093/mnras/stx138

Rein, H., & Liu, S.-F. (2012). REBOUND: an open-source multi-purpose N-body
code for collisional dynamics. _Astronomy & Astrophysics_, _537_ , A128.
https://doi.org/10.1051/0004-6361/201118085

Salvatier, J., Wiecki, T. V., & Fonnesbeck, C. (2016). Probabilistic
programming in python using PyMC3. _PeerJ Computer Science_ , _2_ , e55.
https://doi.org/10.7717/peerj-cs.55

Southworth, J., Maxted, P. F. L., & Smalley, B. (2004). Eclipsing binaries in
open clusters - II. V453 Cyg in NGC 6871. _Monthly Notices of the Royal
Astronomical Society_ , _351_(4), 1277–1289.
https://doi.org/10.1111/j.1365-2966.2004.07871.x

Tamayo, D., Rein, H., Shi, P., & Hernandez, D. M. (2020). REBOUNDx: a library
for adding conservative and dissipative forces to otherwise symplectic N-body
integrations. _Monthly Notices of the Royal Astronomical Society_ , _491_(2),
2885–2901. https://doi.org/10.1093/mnras/stz2870

Theano Development Team. (2016). Theano: A Python framework for fast
computation of mathematical expressions. _arXiv e-Prints_ , _abs/1605.02688_.
http://arxiv.org/abs/1605.02688

Trifonov, T. (2019). _The Exo-Striker: Transit and radial velocity interactive
fitting tool for orbital analysis and N-body simulations_ (p. ascl:1906.004).

Van Eylen, V., Albrecht, S., Huang, X., MacDonald, M. G., Dawson, R. I., Cai,
M. X., Foreman-Mackey, D., Lundkvist, M. S., Silva Aguirre, V., Snellen, I., &
Winn, J. N. (2019). The Orbital Eccentricity of Small Planet Systems. _The
Astronomical Journal_ , _157_(2), 61. https://doi.org/10.3847/1538-3881/aaf22f

Willard, B. T., Osthege, M., Ho, G., Vieira, R., Wiecki, T., Foreman-Mackey,
D., Chaudhari, K., Legrand, N., Kumar, R., Lao, J., Abril-Pla, O., Fonnesbeck,
C., Goldman, R. P., & Gorelli, M. (2021). _pymc-devs/aesara_ (Version 2.0.7)
[Computer software]. Zenodo. https://doi.org/10.5281/zenodo.4695331