11institutetext: Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France
11email<EMAIL_ADDRESS>
# Spiral-wave-driven accretion in quiescent dwarf novæ

M. Van den Bossche 11 G. Lesur 11 G. Dubus 11

(Received XXX; accepted April 26, 2023; in original form )

###### Abstract

Context. In dwarf novæ (DNe) and low-mass X-ray binaries, the tidal potential
excites spiral waves in the accretion disc. Spiral-wave-driven accretion may
be important in quiescent discs, where the angular momentum transport
mechanism has yet to be identified. Previous studies were limited to
unrealistically high temperatures for numerical studies or to specific regimes
for analytical studies.

Aims. We perform the first numerical simulation of spiral-wave-driven
accretion in the cold temperature regime appropriate to quiescent discs, which
have Mach numbers $\gtrsim 100$.

Methods. We used the new GPU-accelerated finite-volume code Idefix to produce
global hydrodynamics 2D simulations of the accretion discs of DN systems with
a sufficiently fine spatial resolution to capture the short scale-height of
cold, quiescent discs with Mach numbers ranging from 80 to 370.

Results. Running the simulations on timescales of tens of binary orbits shows
transient angular momentum transport that decays as the disc relaxes from its
initial conditions. We find the angular momentum parameter $\alpha$ drops to
values of $\ll 10^{-2}$, too weak to drive accretion in quiescence.

###### Key Words.:

CV discs – accretion – compact binary – hydrodynamics

##  1 Introduction

Dwarf novæ (DNe), which are binary systems composed of a white dwarf accreting
from a low-mass companion that fills its Roche lobe, display outbursts due to
changes in the flow of matter through the accretion disc surrounding the white
dwarf (Smak, 1971; Osaki, 1974). Like low-mass X-ray binaries composed of a
black hole or neutron star accreting from a low-mass companion, DNe cycle
between an outburst state and a quiescent state, albeit on a much shorter
recurrence timescale of the order of a month.

The disc instability model (DIM; e.g. Lasota 2001; Hameury 2020 for reviews)
proposes a mechanism to explain the luminosity variation of these systems,
which is related to a gas-ionisation and opacity hysteresis cycle. During
quiescence, the disc gas is mostly neutral and marginally optically thick
($\tau\approx 1$). During this phase, matter from the companion accumulates
until a surface density threshold is reached in the disc. At this point, the
strong dependence of the opacity on temperature as a result of hydrogen
ionisation leads to a thermal runaway that heats the disc. Eventually, the
disc ends up hot, fully ionised, and optically thick: the outburst state. In
this state, the disc has a higher accretion rate onto the white dwarf than the
mass fed at the outer disc from the secondary. The disc density therefore
decreases until its density reaches a lower threshold where hydrogen starts to
recombine. The typical recurrence time for these systems is around 40 days in
total, with about a month for quiescence and a week for the outburst phase
(Osaki, 1996).

The transport of matter is driven by the redistribution of angular momentum in
the disc, historically described as a turbulent viscosity $\nu=\alpha c_{\rm
s}H,$ where $c_{\rm s}$ is the local sound speed and $H$ the local scale
height Shakura & Sunyaev (1973). The parameter $\alpha$ depends on the physics
that governs the transport of angular momentum. Assuming a radiatively
efficient thin disc, it is possible to derive the value of $\alpha$ required
to reproduce the outburst cycles. Mineshige & Osaki (1983), Meyer & Meyer-
Hofmeister (1984), Smak (1984), and Martin et al. (2019), amongst others,
derived an $\alpha$ parameter value of 0.1 to 0.3 during outburst and 0.01
during quiescence, based on observations.

The magneto-rotational instability (MRI, Velikhov 1959, Chandrasekhar 1961,
Balbus & Hawley 1991) is now widely accepted as the source of turbulent
angular momentum transport during the hot, outburst phase (Latter & Papaloizou
2012, Hirose et al. 2014). However, the MRI fails to describe the angular
momentum transport process during quiescence. During this colder phase, the
plasma is too weakly ionised, and so MRI-driven turbulence is expected to be
largely reduced (Gammie & Menou, 1998). This conjecture has been confirmed
through local shearing box models (Scepi et al., 2018), which show that the
molecular resistivity of such cold plasma suppresses MRI-driven angular
momentum transport. Consequently, we do not yet understand the angular
momentum transport mechanism that acts during the quiescent phase.

An alternative route to angular momentum transport in MRI-stable discs is
spiral shock waves. Spiral waves are known to be excited by the tidal
potential of the binary, and patterns observed in the orbital light curves or
Doppler tomograms of cataclysmic variable (CV) accretion discs are consistent
with their presence (e.g. Pala et al., 2019; Ruiz-Carmona et al., 2020).
Studies using linear (Savonije & Papaloizou, 1983; Savonije et al., 1994; Xu &
Goodman, 2018) and non-linear self-similar analysis (Spruit, 1987; Hennebelle
et al., 2016) find that the angular momentum transport driven by spiral shocks
is strongly correlated with temperature (typically with $\alpha\propto
T^{3/2}$). This arises because the temperature controls the disc Mach number
(see Sect. 2), which determines the opening angle of spiral shocks and, in
turn, the angular momentum flux. Therefore, in low-temperature quiescent discs
with Mach numbers of the order of several hundred we expect spiral shock to
drive very weak angular momentum transport.

However, the spiral scenario was recently revived by global simulations of
close binary systems (Ju et al., 2016; Pjanka & Stone, 2020) that emphasise
the important role of spiral density waves in the dynamics of these objects.
Due to the high computational cost of global models, these simulations focused
on hot discs, with Mach numbers of the order of a few tens. Simulating high
Mach numbers is challenging, as this requires high resolution to accurately
capture the shallow angle between the tightly wound spiral and the Keplerian
flow (Matsuda et al., 1990; Rozyczka & Spruit, 1993; Blondin, 2000; Makita et
al., 2000; Kley et al., 2008). Simulations are also typically limited in the
radial range and time span over which they follow the evolution of the spiral
waves. It has therefore not yet been feasible to determine from direct global
simulations whether or not spiral shock waves are a viable angular momentum
transport mechanism in quiescence, when the Mach number can reach $900$ at the
inner edge of the disc (from cold branch temperature estimates of the DIM).

In this paper, we propose to tackle this question using a new GPU-accelerated
finite-volume code, Idefix, which allows us to explore these quiescent
regimes. We present the first 2D global hydrodynamics model with a
sufficiently fine resolution to resolve the spiral waves in a cold quiescent
disc. The code is also fast enough for us to explore timescales well beyond
transients due to relaxation from the initial state. In Sect. 2, we review the
analytical formulation and present the code used to solve the equations. In
Sect. 3, we show two-dimensional simulations of isothermal discs with
realistically high Mach numbers. We discuss the consequences of our findings
in Sect. 4 and present conclusions in Sect. 5.

##  2 Methods

In this paper, we use the Idefix111Stable version and documentation can be
found at https://github.com/idefix-code/idefix (Lesur et al., 2023) finite-
volume astrophysics code to solve the hydrodynamic Euler equations in polar
coordinates $(R,\phi)$. More details on the code are given in Sect. 2.3.

###  2.1 Analytical formulation

We solve the Euler equations for hydrodynamics:

$\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\bm{v})=0,$ (1)
$\frac{\partial}{\partial
t}(\rho\bm{v})+\nabla\cdot(\rho\bm{v}\bm{v}^{T})=-\nabla p-\rho\nabla\Psi,$
(2)

for a fluid of density $\rho$, velocity $\bm{v}$, with pressure $p$ and in a
gravitational potential $\Psi$. The above system is closed with the ideal gas
equation of state:

$p=c_{\text{s}}^{2}\rho.$ (3)

The runs presented here are isothermal: the temperature and sound speed
$c_{\text{s}}$ are constant over time and space. Here, since we only present
2D simulations, the volume density $\rho$ and the pressure $p$ are vertically
integrated. We write these vertically integrated quantities $\Sigma$ and $P$,
respectively.

We work in a rotating frame centred on the white dwarf in which the secondary
star is fixed. For our binary systems, the tidal potential $\Psi$ in the
rotating frame is given by

$\Psi(R,\phi)=-\frac{GM_{\text{WD}}}{R}-\frac{GM_{\text{s}}}{\sqrt{R^{2}+a^{2}-2aR\cos(\phi)}}\\\
+\frac{GM_{\text{s}}}{a^{2}}R\cos(\phi),$ (4)

where $M_{\text{WD}}$ and $M_{\text{s}}$ are the masses of the white dwarf and
secondary star, respectively. $G$ is the gravitational constant, and $a$ is
the binary separation. The third term of this potential accounts for the non-
inertial reference frame. The inertial forces are added by the solver. In
practice, the centrifugal force is added as a radial source term
$\Sigma{v_{\phi}^{2}}/{R}$ with $v_{\phi}$ being the total azimuthal velocity
(i.e. taking both the rotation of the reference frame and a possible advection
velocity). The Coriolis force is included as a modified inter-cell momentum
flux to guarantee angular momentum conservation at machine precision (Mignone
et al., 2012).

###  2.2 Grid and units

For all simulations in this paper, we use a logarithmic grid in the radial
direction, such that $\Delta R\propto R$ and a uniform grid in the azimuthal
direction. The former allows a finer resolution of the inner edge of the disc,
where the spiral shock structures are smaller (Spruit, 1987). We set the inner
boundary of the disc at $r_{0}=0.01\,a$, which corresponds to the expected
radius of the white dwarf. The outer radius of the integration domain is
chosen to be the radial position of the Lagrange point $L_{1}$ of the binary
system. With these definitions, the outer radius changes with the mass ratio
$q=M_{\text{s}}/M_{\text{WD}}$. However, we choose the resolution such that
going to a further radius only appends points to the grid, with no
modification of the inner region (this concerns the runs $\mathrm{Ma}=250$
from $q=0.3$ to $q=0.1$ and $q=0.7$).

In this paper, we choose the time unit as the binary period
$T_{0}=\frac{2\pi}{\omega}=2\pi\sqrt{\frac{a^{3}}{G(M_{\text{s}}+M_{\text{WD}})}}$.
The length unit is the binary separation $a$. Following Ju et al. (2016),
these units are scaled on the dwarf nova system SS Cygni, where $T_{0}$ is
approximately 6.6 hours and $a$ is 1.37$\times 10^{11}$ cm (Bitner et al.,
2007).

We define the Mach number at the inner edge as

$\mathrm{Ma}=\frac{\Omega_{K}(r_{0})r_{0}}{c_{\text{s}}}\approx
364\left(\frac{M_{\rm
WD}}{\text{1}\,M_{\odot}}\right)^{1/2}\left(\frac{10^{9}\rm\,cm}{r_{0}}\right)^{1/2}\left(\frac{10^{4}\rm\,K}{T}\right)^{1/2},$
(5)

where $\Omega_{K}(R)=\sqrt{{GM_{\text{WD}}}/{R^{3}}}$ is the Keplerian angular
frequency around the white dwarf. We note that this definition depends on the
chosen inner boundary of the simulation; for example, in an isothermal setup
$\mathrm{Ma}\propto\sqrt{r_{0}}$. This must be taken into account when
comparing with previous work. In particular, Ju et al. (2016) set their inner
boundary further out at $r_{0}=0.02\,a$.

###  2.3 Algorithm

The simulations presented in this paper all use the Idefix code, which is a
finite-volume conservative Godunov grid-based code with a structure similar to
Pluto (Mignone et al., 2007). The code is written in C++17 and uses the Kokkos
portability library (Edwards et al., 2014) for many-core shared memory
parallelisation and MPI for distributed memory parallelisation, allowing for
high performance on most available architectures.

We use the HLLC Riemann problem solver (Harten et al., 1983; Toro et al.,
1994). We also use the Fargo algorithm (Masset, 2000) with a Keplerian
advection velocity to speed up the integration, as implemented in Mignone et
al. (2012). We verified that we are able to reproduce the work of Ju et al.
(2016) with this code (see details in Appendix B).

Table 1: Runs presented in this paper. $\begin{array}[]{ccc}\hline\cr\vskip
3.0pt plus 1.0pt minus 1.0pt\cr\mathrm{Ma}(r_{0})&\text{Resolution
}(N_{R}\times N_{\phi})&q\\\ \vskip 3.0pt plus 1.0pt minus
1.0pt\cr\hline\cr\vskip 3.0pt plus 1.0pt minus 1.0pt\cr 80&1081\times
1024&0.3\\\ 140&1081\times 1024&0.3\\\ 250&1081\times 1024&0.3\\\
370&1081\times 1024&0.3\\\ \hline\cr 250&1119\times 1024&0.1\\\ 250&1042\times
1024&0.7\\\ 250\textrm{WKZ}&1081\times 1024&0.3\\\ 550\textrm{HR}&4096\times
4096&0.3\\\ \vskip 3.0pt plus 1.0pt minus 1.0pt\cr\hline\cr\end{array}$

###  2.4 Boundary conditions and density floor

In the azimuthal direction, we use periodic boundary conditions while we use
‘Keplerian outflow’ boundary conditions in the radial direction. These last
boundary conditions are identical to standard ‘outflow’ boundaries, with the
exception of the azimuthal velocity, which is set to its Keplerian value in
the ghost zones. We additionally use a density floor to limit the density
contrast in the simulation. If the density of one cell drops below the
threshold $\Sigma_{\mathrm{floor}}=10^{-6}$, then this value is used as
density instead. The density floor is implemented in a total momentum
conserving fashion, so as to prevent the injection of angular momentum by the
procedure.

###  2.5 Initial conditions

We take the initial profile of the disc to have a uniform surface density
$\Sigma=1$. The initial velocity profile is Keplerian with
$v_{\phi}=R\Omega_{\text{K}}=\sqrt{GM_{\text{WD}}/R}$ and no initial radial
velocity. This state is not an equilibrium state, because the initial velocity
only takes into account the gravitational influence of the central white
dwarf, not the companion.

We pretruncate the disc at the outer edge to shorten the time needed for the
outer disc to settle down. The initial truncation radius is chosen to be
slightly larger than the expected $r_{\text{max}}$ of Paczynski (1977). The
truncation is performed using the following mask applied to the initial
density profile:

$\mu(R)=\frac{1}{2}\left(1-\tanh\left(\frac{R-r_{T}}{\delta_{T}}\right)\right),$
(6)

where $r_{T}$ is the chosen truncation radius and $\delta_{T}=0.01$ is the
chosen truncation width. We note that this truncation is not perfect as it is
axisymmetric, unlike the predicted last stable disc orbit (Paczynski, 1977).
At larger radii, we enforce the previously described density floor. We
typically wait for the system to relax over a few binary orbital periods
before measuring any quantity.

###  2.6 Averaging methods

In this work, the dynamical timescale is $\simeq$ 1/10,000${}^{\text{th}}$ of
the binary period, while we are following the evolution of the system for
dozens of binary orbits. To capture the evolution in a computationally
efficient way, we used an ‘on-the-fly’ time and azimuthal averaging scheme: we
sample the simulation every 1/10,000${}^{\text{th}}$ of a binary orbit and
average over 100 of these samples, producing an output each
1/100${}^{\text{th}}$ of a binary orbit. This is sufficient to capture the
rapid inner-edge dynamics.

##  3 Results

In this section, we present globally isothermal models with realistically low
temperatures (high Mach numbers).

###  3.1 Resolution

We chose the resolution such that the expected spiral wave is well resolved
everywhere in the disc. Linear theory predicts that the wavelength at the
inner edge of the disc is (Savonije et al., 1994)

$\displaystyle\frac{\lambda}{R}$
$\displaystyle\simeq\frac{2\pi}{\sqrt{m^{2}-1}}\sqrt{\frac{R}{r_{0}}}\mathrm{Ma}^{-1}(r_{0}),$
(7)

where $m$ is the azimuthal wave number. For the $m=2$ mode that is expected to
be dominant Savonije et al. (1994), a resolution of $N_{R}=1000$ gives a
minimum of 2.5 radial points per wavelength for our maximum Mach number
$\mathrm{Ma}=370$ (Fig. 6). We consider this sufficient to resolve the entire
spiral in the disc. We use about as many points in the azimuthal direction as
in the radial direction to keep a cell aspect ratio of
$\frac{Rd\phi}{dR}\simeq 1.6$, that is, of order unity. Therefore, the
resolution of our isothermal runs is 1081 radial points and 1024 azimuthal
points (Table 1). The run with $\mathrm{Ma}=250$ and $q=0.1$ has a larger
radial extension and consequently larger $N_{R}$ to keep the grid spacing
unchanged in the inner disc. The opposite is true for the run with $q=0.7$.

(a) $\mathrm{Ma}=80$

(b) $\mathrm{Ma}=140$

(c) $\mathrm{Ma}=250$

(d) $\mathrm{Ma}=370$

Figure 1: Density maps after ten binary orbits for isothermal runs with
different Mach numbers.

(a) $\mathrm{Ma}=80$

(b) $\mathrm{Ma}=140$

(c) $\mathrm{Ma}=250$

(d) $\mathrm{Ma}=370$

Figure 2: Azimuthal spectrum of the relative density fluctuation distribution
for modes $m=1$ (red solid line) and $m=2$ (green dashed line) after ten
binary orbits.

###  3.2 Spiral waves

As expected from previous works (Sect. 1), we see that spiral waves develop in
the accretion disc due to the tidal potential (Fig. 1). To get a better grasp
of the excited modes, we show the Fourier decomposition of the density in Fig.
2. The plotted values correspond to the amplitude of the Fourier expansion of
the relative density fluctuation $|\tilde{\Sigma}_{m}(R)|$ smoothed over 50
radial points (moving average), defined from

$\tilde{\Sigma}(R,\phi)=\frac{\Sigma(R,\phi)-\langle\Sigma\rangle_{\phi}(R)}{\langle\Sigma\rangle_{\phi}(R)}=\sum_{m=0}^{\infty}\tilde{\Sigma}_{m}(R)e^{2i\pi
m\phi},$ (8)

where
$\langle\Sigma\rangle_{\phi}(R)=\frac{1}{2\pi}\int_{0}^{2\pi}\Sigma(R,\phi)d\phi$
is the azimuthal average of the surface density. Figure 2 shows strong $m=1$
and $m=2$ spiral modes are present in all runs. The $m=2$ mode predicted by
linear theory dominates only at low Mach numbers and in the outer disc. The
simulations show that the $m=2$ mode is initially excited at the outer
boundary, propagates through the disc, and is reflected when it reaches the
inner edge. The strong $m=1$ mode becomes visible in density maps after this
reflection together with the formation of an eccentric inner cavity.

To investigate the impact of the chosen inner boundary conditions on the
appearance of the $m=1$ mode, we performed an additional run at
$\mathrm{Ma}=250$ and $q=0.3$ implementing a wave-killing zone close to the
inner edge extending from $r_{0}$ to $r_{\text{WKZ}}$. In this zone, we relax
all physical quantities to the initial Keplerian state of the disc, with
constant surface density and no radial velocity. The associated relaxation
timescale is set to one-tenth of the orbital time $\Omega_{K}^{-1}(r_{0})$ at
the inner boundary. For a wave-killing zone to be effective, one needs to
choose its spatial extension to be a few wavelengths of the mode one wishes to
neutralise (see figure 6). In our case, we first take
$r_{\text{WKZ}}=1.2r_{0}$, so that the $m=2$ spiral mode is cancelled at the
inner edge. However, this had no significant impact either on the dynamic or
on the eccentric cavity at the inner edge. Using a more extended wave-killing
zone with $r_{\text{WKZ}}=2r_{0}$ yields the same results. We conclude that
there is more at play than a mere reflection exciting the $m=1$ mode, most
likely because we do not fall within the linear theory regime (see Sect. 4).

###  3.3 Measured $\alpha$

We measure the effective angular momentum transport parameter $\alpha$
(Shakura & Sunyaev, 1973) in our simulations as

$\alpha(R)=\frac{\langle F_{R}(v_{\phi})\rangle_{t,\phi}}{\langle\rho
c_{\text{s}}^{2}\rangle_{t,\phi}},$ (9)

where $F_{R}(v_{\phi})$ is the inter-cell angular momentum flux in the radial
direction $R,$ which can be written $F_{R}(v_{\phi})=\rho
v^{\prime}_{r}v^{\prime}_{\phi}$, where
$v^{\prime}_{\phi}=v_{\phi}-v_{\text{K}}$ is the departure from the Keplerian
velocity. The averages are computed over the whole azimuthal domain
$\phi\in\left[0,2\pi\right]$, and over a few binary orbits. The exact time
span over which the average is carried out will always be specified in the
following. We use the inter-cell angular momentum fluxes computed by the
Riemann solver in the integration loop rather than the reconstructed cell-
centred values. This allows us to better capture the small-scale dynamics. We
compare our method to measure $\alpha$ with that of Ju et al. (2016) in
Appendix C.

Figure 3: Measured angular momentum transport parameter. Top panel: Measured
angular momentum transport parameter $\alpha$ azimuthally averaged and
smoothed over three binary orbits for the run $\mathrm{Ma}=250$ and $q=0.3$.
Here we plot the absolute value of the transport parameter, and note that its
sign changes in the spiral waves (oscillations at $R\gtrsim 0.1$) and is
mostly positive everywhere else. Bottom panel: Azimuthally and radially (over
$0.02\leq R\leq 0.1$) averaged $\alpha$, smoothed over three binary orbits for
three runs.

Figure 3 shows the time evolution of $\alpha$ in the disc. In all simulations,
we first see a transient regime at early times ($t\lesssim 20$) with a high
$\alpha$ that slowly decays. At later times ($t\gtrsim 40$), the transport
parameter in the disc drops to $\alpha\ll 10^{-2}$ in all simulations,
regardless of the mass ratio $q$ or the Mach number. The relaxation timescale
appears to be dependent on the Mach number. For example, the simulation with
$\mathrm{Ma}=370$ shows a slower relaxation than the $\mathrm{Ma}=250$ run
(bottom panel of Fig. 3). The state reached after 80 orbits may still not be a
steady state and $\alpha$ is still decreasing. We did not go beyond this
stage, as $\alpha$ is already much lower than required by the DIM.

We caution that regions of very high $\alpha$ at the inner and outer edges are
unphysical regions where the density floor is triggered because of the growth
of an eccentric inner cavity (Sect. 3.2) or the tidal truncation of the outer
disc. We also note that matter accumulates at the outer radius of the disc at
lower temperatures. This is a result of relaxation of the disc truncation from
its initial state, which assumes a Keplerian velocity profile and an
analytical guess for the outer radius. This residual matter is not efficiently
redistributed by the spiral waves given the weaker angular momentum transport
at lower temperatures.

Figure 4: Measured angular momentum transport parameter $\alpha$ at late times
(between $t=80$ and $t=89.9$) for different Mach numbers with $q=0.3$,
averaged over time and azimuth.

Figure 4 shows the radial distribution of $\alpha$ at late times for
simulations with different Mach numbers, keeping $q=0.3$. The plotted values
are both azimuthally- and time-averaged. The strong oscillations around
$R=0.2-0.3$ reflect the effect of the tidal torque on the spiral waves. Here,
unlike in Ju et al. (2016), the torque term of the angular momentum
conservation equation does not exactly balance the Reynolds stress term. As a
consequence, the measured accretion rate and the measured $\alpha$ are not
exactly zero, but rather oscillate, which creates regions of slightly negative
accretion rate. The efficiency of transport decreases with disc temperature,
and in all cases we get $\alpha\ll 10^{-2}$. We do not find evidence for a
radial steepening of $\alpha$ with higher Mach numbers (e.g. Blondin, 2000).

Figure 5: Measured angular momentum transport parameter $\alpha$ at late times
(between $t=80$ and $t=89.9$) for different mass ratios with
$\mathrm{Ma}=250$, averaged over time and azimuth.

Figure 5 shows how $\alpha$ changes with varying mass ratio $q$, keeping
$\mathrm{Ma}=250$. Again, we obtain $\alpha\ll 10^{-2}$ at late times.
However, unlike the runs at $q=0.3$ and $q=0.7$, the transport appears to
stabilise at a higher value of $\alpha\approx 2\times 10^{-3}$ for $q=0.1$.
This is consistent with a possible coupling with the 3:1 resonance, which is
inside the disc for $q<0.3$ (Lubow, 1991).

###  3.4 Convergence study

In order to make sure that the resolution of our runs is sufficiently fine, we
performed a higher-resolution run of our coolest setup. For this simulation,
we increase the resolution by a factor of four, with a grid of size
$4096\times 4096$. This resolution is much finer than what is required to
resolve the $m=2$ spiral wave from linear theory. However, the theoretical
estimate may not hold given that we are far from the linear regime (see Sect.
4). We observe the same spiral wave structure in the high-resolution run as in
the low-resolution run, albeit with a slower decrease in $\alpha$. This is
likely a consequence of the decreased numerical diffusivity at higher
resolution. The relaxation timescales that we observe in our runs (Fig. 3) are
therefore partly controlled by the numerical diffusivity. We reach $\alpha\ll
10^{-2}$ after 100 binary orbits, a value consistent with the lower-resolution
runs, confirming that these low values of $\alpha$ are not an artefact of
resolution.

##  4 Discussion

###  4.1 One-armed spirals

Our simulations show that the $m=1$ modes are excited concomitantly with $m=2$
spiral waves. Previous works indicated that $m=1$ modes could be excited in a
close binary system provided (1) that the mass ratio is low enough ($q<0.3$)
for the 3:1 resonance radius to lie inside the disc and (2) that the disc is
viscous, so as to couple the $m=2$ and $m=3$ modes to excite the $m=1$ modes
(Heemskerk, 1994; Stehle, 1999; Kornet & Rozyczka, 2000; Kley et al., 2008).
Here, we satisfy neither of these criteria: we do not include viscosity and
our results hold for both $q=0.3$ and $q=0.7$. We do not observe a strong
global eccentricity growth of the disc as Kley et al. (2008) do: the
eccentricity mode is limited to the inner disc region, keeping the overall
shape of the outer disc circular (see Appendix D). The emergence of this inner
eccentric cavity may be of relevance to observations that suggest the presence
of a truncated disc (e.g. Balman & Revnivtsev 2012).

Our results also hold when including a wave-killing zone at the inner boundary
to dampen the $m=2$ reflection. This suggests a less efficient, non-linear
coupling is at play in our simulations that may have been missed in lower-
resolution simulations with lower Mach numbers. Following the analytical
development from Savonije et al. (1994), we can show that linear theory is
expected to fail for the realistic temperature regimes that we explore here.
From the general solution derived in this latter paper, the amplitude of the
density perturbation of the solution $\Sigma^{\prime}$ compared to the
unperturbed density profile $\Sigma_{0}$ can be quantified as

$\sigma(r)=\frac{\Sigma^{\prime}}{\Sigma_{0}}=\frac{3}{4}\frac{1}{1+q}\left(\frac{\omega}{\Omega_{K}}\right)^{2}\text{Ma}^{2},$
(10)

where $\omega$ is the binary angular frequency. For a quiescent disc
temperature of 3,000 K and with typical masses for the white dwarf and its
companion star, we obtain $\sigma(r_{0})\simeq 1$ at the inner edge. As
$\sigma(r)\propto r^{2}$, non-linear effects increase at larger radii. The
solution derived by Savonije et al. (1994) is based on a perturbative WKB
solution where $\sigma\ll 1$ and with the assumption that $r\ll a$. The
emergence of a one-armed spiral wave could not have been predicted by linear
theory, as the regime probed here strongly deviates from this approximation.

###  4.2 Angular momentum transport

The initial high value of the $\alpha$ parameter mentioned in Sect. 3.3 is
likely a consequence of our initial condition, which assumes circular
Keplerian orbits for the gas. This is clearly not an equilibrium when a
secondary star is present. Relaxation from these conditions takes more than 20
orbits. This is longer than for example the time over which Ju et al. (2016)
ran their simulations, meaning that their values of $\alpha$ may still be
impacted by their initial conditions (which are similar to ours). However, the
relaxation time is only a few binary orbits at the lower Mach numbers of their
simulations.

On long timescales, our simulations show that $\alpha$ is always much smaller
than the typical $\alpha\approx 0.01$ in quiescence required by the comparison
of DIM models to observed light curves. The values we reach are of the order
of $\alpha\approx 10^{-4}$ and are still decaying when we end the simulation.
The one exception is for $q=0.1$ where transport stabilises to a larger
$\alpha\approx 10^{-3}$, probably because the 3:1 resonance allows waves
inside the disc to be excited. This can only apply to systems with $q<0.3$,
such as the WZ Sge subclass (Kato, 2015), while DNe show mass ratios in the
range of $0.1\lesssim q\lesssim 0.8$ (Otulakowska-Hypka et al., 2016).

In our simulations with the highest Mach number, the relaxation time is
comparable to the total integration time of 100 binary orbits, which is on the
order of the recurrence time between two normal outbursts in DNe systems
(Otulakowska-Hypka et al., 2016). Spiral waves therefore might still play a
role in transporting angular momentum in quiescence if the relaxation to low
$\alpha$ is slow enough. There are two caveats that make this scenario
unlikely. First, our convergence run indicates that the relaxation time is
likely to depend on the numerical resolution, and so deducing a reliable decay
timescale from the simulation is not straightforward. Second, this would
depend critically on how the disc relaxes from MRI-driven transport in
outburst. Our initial conditions assume an isothermal disc with a flat surface
density profile. A constant-temperature disc is a reasonable approximation for
a DIM quiescent disc. However, the surface density profile is expected to
increase with radius ($\Sigma\propto R$) when the disc enters quiescence due
to the dependence on radius of the critical densities associated with the
temperature hysteresis (Lasota, 2001). This may introduce differences with our
flat $\Sigma$ profile and, more generally, induce a complex interaction in the
time evolution of the $\alpha$ and $\Sigma$ density profiles.

###  4.3 RWI/KHI spiral excitation

We also verified that the spiral-driven accretion that we observe is indeed a
consequence of the presence of a secondary star. Lesur et al. (2015) showed
that an infalling envelope around the disc could trigger a mixture of a Rossby
Wave instability (Lovelace et al., 1999) and centrifugal instability at the
disc outer edge, creating spiral perturbations that propagate inwards. The
outer edge of our disc fulfils the instability criteria of these two
instabilities, in part because of the density floor that induced an inflow
from the ‘void’ surrounding the disc, but also because of the sharp truncation
caused by the secondary.

To check the impact of these edge instabilities, we simulated the same discs
without the companion, keeping the disc truncated density profile. We find
that in this case, these instabilities only drive transport $\alpha\ll
10^{-4}$ in the disc bulk and are therefore of negligible importance in the
present work (see Fig. 12).

However, it is possible that the tidal potential of the secondary star selects
and excites a specific mode of the edge instabilities. This could be a
possible mechanism to explain the unexpected apparition of $m=1$ spiral modes.

###  4.4 Expected instabilities in 3D

Several instabilities expected in 3D are absent from our 2D approach, among
which the vertical shear instability (VSI, Nelson et al., 2013) and the
convective over-stability (COV, Klahr & Hubbard, 2014). The existence of these
two instabilities is mostly controlled by the dimensionless cooling timescale
$\Omega_{K}\tau_{\mathrm{cool}}$ where $\tau_{\mathrm{cool}}$ measures the
typical thermal relaxation timescale resulting from the heating/cooling
equilibrium. In the case of CVs, we expect $\Omega_{K}\tau_{\mathrm{cool}}\ll
1$, a regime in which we expect the VSI to be present but not baroclinic
instabilities (Lesur et al., 2022).

The question is then how much angular momentum transport is to be expected
from the VSI. Because the VSI is driven by the vertical shear of the flow, it
depends strongly on the disc thickness. Using high-resolution 3D
hydrodynamical simulations, Manger et al. (2020) proposes a scaling for the
VSI $\alpha\lesssim(H/R)^{2.6}$, which implies $\alpha<10^{-5}$ in a quiescent
DNe disc subject to VSI turbulence. This means that VSI turbulence, if
present, is probably negligible as an angular momentum transport process in
these systems.

###  4.5 Accretion stream impact and thermodynamics

The contribution of the matter flux coming from the companion star (called
accretion stream) is expected to be the strong heating at the ‘hot spot’ where
it meets the disc. As a consequence, we did not include it in this isothermal
study. Modelling the thermodynamics of the disc would allow us to know the
heating efficiency of the spiral shocks presented in this work. Previous works
such as that of Ju et al. (2016) showed that without any cooling, and using
$\gamma=5/3,$ the disc heats up quickly. This suggests that some parts of the
spiral-heated disc might still be unstable to MRI even during quiescence.

##  5 Conclusion

We present the first 2D hydrodynamics simulations of a DN disc in realistic
quiescence temperature regimes. Our main results can be summarised as follows.

  1. 1.

The dominant spiral mode is a one-armed $m=1$ mode once the disc has settled.
Linear theory, which predicts an $m=2$ mode, fails at low temperatures (high
Mach number), as expected.

  2. 2.

Spiral-wave-driven accretion in realistic temperature regimes is much lower
than the values of $\alpha$ required by the DIM. Running the simulations
clearly shows that, regardless of the initial transport, there is decays as
the disc relaxes from its initial conditions.

We conclude that hydrodynamic spiral waves appear unable to explain accretion
in quiescent dwarf novæ and low mass X-ray binaries (whose outer discs sample
similar densities and temperatures).

###### Acknowledgements.

The authors would like to thank the anonymous referee for their comments that
helped them improve this work. GL acknowledges support from the European
Research Council (ERC) under the European Union Horizon 2020 research and
innovation program (Grant agreement No. 815559 (MHDiscs)). The simulations and
data reduction used in this article were performed using the GRICAD
infrastructure (https://gricad.univ-grenoble-alpes.fr), which is supported by
the Grenoble research communities, as well as the Jean-Zay supercomputer at
IDRIS http://www.idris.fr/jean-zay/. This work was granted access to the HPC
resources of IDRIS under the allocation 2021-A0120402231 made by GENCI. Data
reduction for this paper was performed by making extensive use of the SciPy
(Virtanen et al., 2020), NumPy (Harris et al., 2020), pandas (Reback et al.,
2022) and Matplotlib (Hunter, 2007) Python libraries. Idefix uses the Kokkos
portability tool (Edwards et al., 2014).

## References

  * Balbus & Hawley (1991) Balbus, S. A. & Hawley, J. F. 1991, ApJ, 376, 214
  * Balman & Revnivtsev (2012) Balman, Ş. & Revnivtsev, M. 2012, A&A, 546, A112
  * Bitner et al. (2007) Bitner, M. A., Robinson, E. L., & Behr, B. B. 2007, ApJ, 662, 564
  * Blondin (2000) Blondin, J. M. 2000, New A, 5, 53
  * Chandrasekhar (1961) Chandrasekhar, S. 1961, Hydrodynamic and hydromagnetic stability
  * Edwards et al. (2014) Edwards, H. C., Trott, C. R., & Sunderland, D. 2014, Journal of Parallel and Distributed Computing, 74, 3202 , domain-Specific Languages and High-Level Frameworks for High-Performance Computing
  * Gammie & Menou (1998) Gammie, C. F. & Menou, K. 1998, ApJ, 492, L75
  * Hameury (2020) Hameury, J. M. 2020, Advances in Space Research, 66, 1004
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357
  * Harten et al. (1983) Harten, A., Lax, P. D., & Van Leer, B. 1983, SIAM Review, 25, 35
  * Heemskerk (1994) Heemskerk, M. H. M. 1994, A&A, 288, 807
  * Hennebelle et al. (2016) Hennebelle, P., Lesur, G., & Fromang, S. 2016, A&A, 590, A22
  * Hirose et al. (2014) Hirose, S., Blaes, O., Krolik, J. H., Coleman, M. S. B., & Sano, T. 2014, ApJ, 787, 1
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90
  * Ju et al. (2016) Ju, W., Stone, J. M., & Zhu, Z. 2016, ApJ, 823, 81
  * Kato (2015) Kato, T. 2015, PASJ, 67, 108
  * Klahr & Hubbard (2014) Klahr, H. & Hubbard, A. 2014, ApJ, 788, 21
  * Kley et al. (2008) Kley, W., Papaloizou, J. C. B., & Ogilvie, G. I. 2008, A&A, 487, 671
  * Kornet & Rozyczka (2000) Kornet, K. & Rozyczka, M. 2000, Acta Astron., 50, 163
  * Lasota (2001) Lasota, J.-P. 2001, New A Rev., 45, 449
  * Latter & Papaloizou (2012) Latter, H. N. & Papaloizou, J. C. B. 2012, MNRAS, 426, 1107
  * Lesur et al. (2022) Lesur, G., Ercolano, B., Flock, M., et al. 2022, arXiv e-prints, arXiv:2203.09821
  * Lesur et al. (2015) Lesur, G., Hennebelle, P., & Fromang, S. 2015, A&A, 582, L9
  * Lesur et al. (2023) Lesur, G. R. J., Baghdadi, S., Wafflard-Fernandez, G., et al. 2023, IDEFIX: Astrophysical fluid dynamics, Astrophysics Source Code Library, record ascl:2306.036
  * Lovelace et al. (1999) Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. F. 1999, ApJ, 513, 805
  * Lubow (1991) Lubow, S. H. 1991, ApJ, 381, 259
  * Makita et al. (2000) Makita, M., Miyawaki, K., & Matsuda, T. 2000, MNRAS, 316, 906
  * Manger et al. (2020) Manger, N., Klahr, H., Kley, W., & Flock, M. 2020, MNRAS, 499, 1841
  * Martin et al. (2019) Martin, R. G., Nixon, C. J., Pringle, J. E., & Livio, M. 2019, New A, 70, 7
  * Masset (2000) Masset, F. 2000, A&AS, 141, 165
  * Matsuda et al. (1990) Matsuda, T., Sekino, N., Shima, E., Sawada, K., & Spruit, H. 1990, A&A, 235, 211
  * Meyer & Meyer-Hofmeister (1984) Meyer, F. & Meyer-Hofmeister, E. 1984, A&A, 132, 143
  * Mignone et al. (2007) Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228
  * Mignone et al. (2012) Mignone, A., Flock, M., Stute, M., Kolb, S. M., & Muscianisi, G. 2012, A&A, 545, A152
  * Mineshige & Osaki (1983) Mineshige, S. & Osaki, Y. 1983, PASJ, 35, 377
  * Nelson et al. (2013) Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610
  * Osaki (1974) Osaki, Y. 1974, PASJ, 26, 429
  * Osaki (1996) Osaki, Y. 1996, PASP, 108, 39
  * Otulakowska-Hypka et al. (2016) Otulakowska-Hypka, M., Olech, A., & Patterson, J. 2016, MNRAS, 460, 2526
  * Paczynski (1977) Paczynski, B. 1977, ApJ, 216, 822
  * Pala et al. (2019) Pala, A. F., Gänsicke, B. T., Marsh, T. R., et al. 2019, MNRAS, 483, 1080
  * Pjanka & Stone (2020) Pjanka, P. & Stone, J. M. 2020, ApJ, 904, 90
  * Reback et al. (2022) Reback, J., jbrockmendel, McKinney, W., et al. 2022, pandas-dev/pandas: Pandas 1.4.2, Zenodo
  * Rozyczka & Spruit (1993) Rozyczka, M. & Spruit, H. C. 1993, ApJ, 417, 677
  * Ruiz-Carmona et al. (2020) Ruiz-Carmona, R., Khangale, Z. N., Woudt, P. A., & Groot, P. J. 2020, MNRAS, 491, 344
  * Savonije & Papaloizou (1983) Savonije, G. J. & Papaloizou, J. C. B. 1983, MNRAS, 203, 581
  * Savonije et al. (1994) Savonije, G. J., Papaloizou, J. C. B., & Lin, D. N. C. 1994, MNRAS, 268, 13
  * Scepi et al. (2018) Scepi, N., Lesur, G., Dubus, G., & Flock, M. 2018, A&A, 609, A77
  * Shakura & Sunyaev (1973) Shakura, N. I. & Sunyaev, R. A. 1973, A&A, 24, 337
  * Smak (1971) Smak, J. 1971, Acta Astron., 21, 15
  * Smak (1984) Smak, J. 1984, Acta Astron., 34, 317
  * Spruit (1987) Spruit, H. C. 1987, A&A, 184, 173
  * Stehle (1999) Stehle, R. 1999, MNRAS, 304, 687
  * Toro et al. (1994) Toro, E. F., Spruce, M., & Speares, W. 1994, Shock Waves, 4, 25
  * Velikhov (1959) Velikhov, E. P. 1959, Zhur. Eksptl’. i Teoret. Fiz., 36
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * Xu & Goodman (2018) Xu, W. & Goodman, J. 2018, MNRAS, 480, 4327

##  Appendix A Linear theory wavelength

The wavelength from linear theory is shortest at the inner edge of the disc
(Eq. 7). Figure 6 shows how this wavelength compares with the grid spacing for
our ‘low-resolution’ runs. Our resolution (with 1081 radial points) is
sufficiently fine up to much higher Mach than explored in this work. In
practice, the $m=1$ mode that develops has a much wider wavelength. This plot
also shows that the width of the wave-killing zones implemented at the inner
boundary is much wider than the linear theory $m=2$ wavelength by several
orders of magnitude (see Sect. 3.2).

Figure 6: Grid spacing versus analytical wavelength for the spiral mode $m=2$
as a function of the radius.

##  Appendix B Reproducibility with Idefix

We verified that we could reproduce the results of Ju et al. (2016) with the
code Idefix. Figures 7, 8, and 9 can be compared with their Figure 12.
Although we only present isothermal simulations in the main text, we also
compared non-isothermal runs with an adiabatic index of $\gamma=1.1$. The
general agreement is very good, with small differences likely due to
variations in implementation between Idefix and Athena++.

As the simulations presented in this paper do not use the same time units (by
a factor of $2\pi$), the units in this section are converted to the same units
as in Ju et al. (2016). We also note that: (1) the $\alpha$ parameter plotted
on Fig. 9 is computed in the same way as in Ju et al. (2016), i.e. from the
radial mass flux; (2) our simulation box extends down to $R=0.02$ but the plot
only extends down to $R=0.04$ to ease comparison with Ju et al. (2016); and
(3) the Mach numbers cannot be directly compared because of differences in the
location of the inner boundary (where it is defined).

Figure 7: Density profile after four binary orbits, averaged over four binary
orbits ($t=25$ to $t=50$ using the units of Ju et al. 2016). Figure 8: Same as
figure 7 but for the accretion rate. Figure 9: Same as figure 7 but for the
effective $\alpha$.

##  Appendix C Definition(s) of $\alpha$

Figure 10: Angular momentum transport parameter $\alpha$ averaged over azimuth
and ten binary orbits for the run $\mathrm{Ma}=250$ and $q=0.3$. The
definition from fluxes (this paper) is plotted with a solid red line and the
definition from radial mass flux (Ju et al. 2016) is plotted as a green dashed
line.

The angular momentum transport parameter can be defined in several ways, and
not all definitions are equivalent, depending on the exact system. Ju et al.
(2016) argue that for binary systems such as the one studied here, this
parameter should be defined from the accretion rate rather than from the
Reynolds (and, in MHD, Maxwell) stress. It has been claimed that this
effective transport parameter $\alpha_{\text{eff}}$ better reflects the
influence of the tidal torque of the companion star. Ju et al. (2016) take

$\alpha_{\text{eff}}=\frac{\dot{M}}{3\pi\Sigma c_{\text{s}}H},$ (11)

where the accretion rate $\dot{M}$ is taken from the inter-cell fluxes
computed by the Riemann solver at each time step. We chose to measure $\alpha$
from the angular momentum flux (Eq. 9). Both definitions can be obtained from
the angular momentum balance equation, using several approximations. In the
case of the definition from the radial mass flux, there is a clear
contribution from the tidal torque. In our simulations, the disc is clearly
not in steady state, and so a definition of the transport parameter based on
the accretion rate $\dot{M}$ cannot easily be decomposed into a Reynolds
stress term (e.g. our $\alpha$) plus a tidal torque term, as this would
require $\dot{M}$ to be constant with respect to the radius. In practice,
there is indeed a slight difference in the values computed using these two
methods; however, these differences do not change the conclusions of this
paper. As shown in Figure 10, regardless of the chosen definition, we obtain
$\alpha\ll 10^{-2}$ at late times.

##  Appendix D Pitch angle measurements

The pitch angle $\theta$ of a spiral $R(\phi)$ is defined as
$\tan\theta=\frac{R^{\prime}}{R}$. From the Fourier decomposition of Eq. (8),
we can define the spiral associated to each mode $m$ as
$\phi_{m}(R)=\arg(\tilde{\Sigma}_{m})$ (up to a constant phase). Therefore,

$\theta_{m}=\arctan\left(\frac{1}{R\phi_{m}^{\prime}(R)}\right).$ (12)

Figure 11 shows the spiral pattern on a deprojected plot. Consistently with
Fig. 2, we see that the one-armed spiral wave dominates the dynamic. The left
panel shows that pitch angle of the spiral wave is consistently $<5^{\circ}$
in the disc. The spiral waves are tightly wound. This explains how, even
though the $m=1$ mode dominates, the disc remains circular in shape overall,
rather than very eccentric.

Figure 11: Density profile properties of the run $\textrm{Ma}=250$ after 38.8
binary orbits. Left panel: Pitch angle of the $m=1$ and $m=2$ spiral modes in
the disc in the density profile. Right panel: De-projected density fluctuation
$\tilde{\Sigma}$, as defined in Eq. (8).

##  Appendix E RWI/KHI-driven accretion

As discussed above, we verified that the mixture of RWI/KHI and centrifugal
instability drives negligible accretion compared to the spiral shocks
triggered by the tidal potential of the secondary star. The setup used to
produce Fig. 12 is the same as for the run with $\mathrm{Ma}=250$ and $q=0.1$.
The mass of the white dwarf is the same but the secondary mass
$M_{\textrm{s}}$ is set to zero in the potential (Eq. 4). This setup is very
similar to the setup of Lesur et al. (2015).

Figure 12: Measured angular momentum transport parameter. Top panel: Measured
angular momentum transport parameter $\alpha$ azimuthally averaged and
smoothed over three binary orbits for a run with $\mathrm{Ma}=250$ and no
companion star. Here we plot the absolute value of the transport parameter.
Bottom panel: Azimuthally and radially (over $0.02\leq R\leq 0.1$) averaged
$\alpha$, smoothed over three binary orbits for three runs. Compared to Fig.
3, the $\alpha$ scale of the bottom panel has been extended to show lower
values of the transport parameter.


# Lie symmetries for the cosmological field equations in brane-world gravity
with bulk scalar field

Andronikos Paliathanasis<EMAIL_ADDRESS>Institute of Systems Science,
Durban University of Technology, Durban 4000, South Africa  Departamento de
Matemáticas, Universidad Católica del Norte, Avda. Angamos 0610, Casilla 1280
Antofagasta, Chile

###### Abstract

We address the group classification problem for gravitational field equations
within the context of brane-world cosmology, considering the presence of a
bulk scalar field. Our investigation revolves around a five-dimensional
spacetime, with the four-dimensional Friedmann–Lemaître–Robertson–Walker
geometry embedded within it. Additionally, we assume that the scalar field
exists in this five-dimensional geometry (bulk) and possesses a nonzero mass.
The resulting field equations constitute a system of nonlinear partial
differential equations. We apply the Lie symmetry condition to identify all
functional forms of the scalar field potential, ensuring that the field
equations remain invariant under one-parameter point transformations.
Consequently, we find that only the exponential potential exhibits Lie
symmetries. Finally, the Lie invariants are used to construct similarity
transformations, which enable us to derive exact solutions for the system.

Brane cosmology; Scalar field; Lie symmetries;

##  I Introduction

The theory of General Relativity is a pure geometrical theory of gravity which
states that physical space is described by a four-dimensional Riemannian
manifold ae1 . General Relativity has been extensively tested and proven
successful in explaining various gravitational phenomena grt1 ; grt2 ; grt3 .
However, in cosmological investigations, the inclusion of the cosmological
constant $\Lambda$ becomes crucial to account for the late-time acceleration
phase of our universe rr1 ; Teg ; Kowal .

The cosmological constant introduces an effective fluid described by an
energy-momentum tensor featuring constant energy density $\rho_{\Lambda}$ and
pressure $p_{\Lambda}$, with a relation $p_{\Lambda}=-\rho_{\Lambda}$.
Although the cosmological constant is a very simple mechanism for the
description of the acceleration of the universe it suffers from two major
problems, the old cosmological problem and the coincidence problem conpr1 ;
conpr2 . The main direction followed by cosmologists to overcome the
aforementioned problems is the introduction of new degrees of freedom in the
cosmological field equations. Scalar fields q4  are particularly relevant in
this context, as they can describe the additional degrees of freedom proposed
in other models q11 ; q1 .

In brane-world gravity, the observable four-dimensional universe is considered
to exist within a higher-dimensional world br1 ; br2 ; br3 ; br4 . This
characteristic offers various mechanisms for addressing the cosmological
constant problems cc1 ; cc2 ; cc3 ; cc4 . Brane-world cosmological scenarios
involving bulk scalar fields have been previously studied bb0 ; bb01 . In bb1
, a massless bulk scalar field was introduced, leading to the derivation of a
new exact solution describing plane waves. Additionally, a scalar field
potential, responsible for imparting mass to the bulk field, was introduced in
bb2 . Specifically, for an exponential scalar field potential, new solutions
were determined, including separable solutions and a solution in which the
scalar field is proportional to the logarithm of the scale factor. Moreover,
in bb3 , brane-world scalar field scenarios with a four-dimensional de Sitter
embedded space were investigated.

In this study, our focus lies on the symmetry analysis of the cosmological
field equations within brane-world gravity, specifically in an extended
spatially flat Friedmann–Lemaître–Robertson–Walker (FLRW) geometry
incorporating a bulk scalar field. The field equations in brane gravity
constitute a system of nonlinear partial differential equations, with an
unknown function representing the scalar field potential.

Symmetry analysis serves as a systematic approach to determine transformations
that preserve the invariance of the given dynamical system Bluman ; Stephani ;
olver . When a symmetry vector exists, we can define a set of new variables to
reduce the given dynamical system through similarity transformations. In the
case of partial differential equations, applying similarity transformations
reduces the number of independent variables, while in the case of ordinary
differential equations, it leads to a dynamical system of lower order.
Additionally, symmetries can be employed to construct conservation laws for
the dynamical system ibra .

The application of Lie symmetries in gravitational physics has found numerous
uses. In the context of the vacuum Einstein field equations within a four-
dimensional manifold, Lie symmetries were explored in ei0 , while a similar
analysis for a higher-dimensional manifold was presented in ei1 . The
inclusion of a perfect fluid was also considered in ei2 .

In cosmology, symmetry analysis has been widely employed to constrain the
unknown functions of the theory and construct conservation laws using
Noether’s second theorem Bluman . The derivation of conservation laws has led
to the discovery of new analytic and exact solutions for various cosmological
models, including those in scalar field theory ns1 ; ns2 , scalar-tensor
theory ns3 ; ns4 ; ns5 , and modified theories of gravity ns6 ; ns7 , as
discussed in the review article ns0 . With this in mind, the structure of the
paper is organized accordingly.

n Section II, we introduce the gravitational model and present the field
equations pertinent to our study, focusing on brane-world cosmology within a
five-dimensional manifold coupled with a bulk scalar field. Moving in Section
III we solve the Lie symmetry classification problem by determining all the
functional forms for the scalar field potential where the field equations
admit Lie symmetries. We demonstrate the application of the Lie symmetries and
we construct similarity solutions for the field equations. Finally, in Section
IV, we summarize our findings and draw conclusions from our analysis.

##  II Brane cosmology with bulk scalar field

We consider a five-dimensional manifold with line element bb1

$ds^{2}=e^{2B\left(t,y\right)}\left(-dt^{2}+dy^{2}\right)+e^{2A\left(t,y\right)}\left(dx_{1}^{2}+dx_{2}^{2}+dx_{3}^{2}\right),$
(1)

in which $A\left(t,y\right)$ and $B\left(t,y\right)$ are the two scalar
factors of the brane. In the line $y=const$, that is, $y=y_{0}$, the line-
element (1) reduces to that of the FLRW, where $A\left(t,y_{0}\right)$ is the
scale factor and $B\left(t,y_{0}\right)$ is the lapse function, and
$H\left(t\right)=e^{-B\left(t,y_{0}\right)}A_{,t}\left(t,y_{0}\right)$ is the
Hubble function. The line element (1) admits a six-dimensional Killing algebra
consisting of the Killing symmetries of the three-dimensional Euclidean space.

Let $\phi$ be a bulk scalar field, we consider the gravitational Action
Integral to be bb2

$S=\int
dx^{5}\sqrt{g}\left(\frac{{}^{\left(5\right)}R}{2}-\frac{1}{2}\nabla_{\kappa}\phi\nabla^{\kappa}\phi-V\left(\phi\right)\right),~{}$
(2)

where $\kappa=1,2,3,4,5$ and ${}^{\left(5\right)}R\,\,\ $is the Ricciscalar
for the five-dimensional metric (1)

${}^{\left(5\right)}R=e^{-2B}\left(B_{,tt}-B_{,yy}+3\left(A_{,tt}+2A_{,t}^{2}-A_{,yy}-2A_{,y}^{2}\right)\right).$
(3)

Thus, the gravitational field equations are

${}^{\left(5\right)}G_{AB}=\left(\nabla_{A}\phi\nabla^{A}\phi-
g_{AB}\left(\frac{1}{2}\nabla_{A}\phi\nabla^{A}\phi+V\left(\phi\right)\right)\right),$
(4)

where ${}^{\left(5\right)}G_{AB}$ is the Einstein tensor for the line element
(1). Furthermore, the equation of motion for the scalar field, that is, the
Klein-Gordon equation reads

$\nabla^{2}\phi+V_{,\phi}=0.$ (5)

For the scalar field $\phi$, we assume that inherits the symmetries of the
background space, that is, $\phi=\phi\left(t,x\right)$.

We follow bb1 ; bb2 , thus we define the variables $\left(u,v\right)$ with
transformation rule

$u=t-y~{},~{}v=t+y,$ (6)

that is

$t=\frac{u+v}{2}~{},~{}y=\frac{v-y}{2}.$ (7)

Hence, in the new variables the line-element (1) becomes

$ds^{2}=e^{2B\left(u,v\right)}\left(-dudv\right)+e^{2A\left(u,v\right)}\left(dx_{1}^{2}+dx_{2}^{2}+dx_{3}^{2}\right),$
(8)

and the field equations (4) are

$\displaystyle A_{,uv}+3A_{,u}A_{,v}$
$\displaystyle=\frac{1}{6}e^{2B}V\left(\phi\right),$ (9)
$\displaystyle-3A_{,uu}+6A_{,u}B_{,u}-3A_{,u}^{2}$
$\displaystyle=\phi_{,u}^{2},$ (10)
$\displaystyle-3A_{,vv}+6A_{,v}B_{,v}-3A_{,v}^{2}$
$\displaystyle=\phi_{,v}^{2},$ (11) $\displaystyle
B_{,uv}+2A_{,uv}+3A_{,u}A_{,v}$
$\displaystyle=-\frac{1}{2}\left(\phi_{,u}\phi_{,v}-\frac{1}{2}e^{2B}V\left(\phi\right)\right),$
(12)

and the Klein-Gordon equation (5)

$4\phi_{,uv}+6\left(A_{,u}\phi_{,v}+A_{,v}\phi_{,u}\right)+e^{2B}V_{,\phi}=0.$
(13)

##  III Lie symmetries for the field equations

In the augmented space $\left\\{u,v,A,B,\Phi\right\\}$ we consider the
infinitesimal transformation

$\displaystyle u^{\prime}$
$\displaystyle=u+\varepsilon\xi^{u}\left(u,v,A,B,\phi\right),$ (14)
$\displaystyle v^{\prime}$
$\displaystyle=v+\varepsilon\xi^{v}\left(u,v,A,B,\phi\right),$ (15)
$\displaystyle A^{\prime}$
$\displaystyle=A+\varepsilon\eta^{A}\left(u,v,A,B,\phi\right),$ (16)
$\displaystyle B^{\prime}$
$\displaystyle=B+\varepsilon\eta^{B}\left(u,v,A,B,\phi\right),$ (17)
$\displaystyle\phi^{\prime}$
$\displaystyle=\phi+\varepsilon\phi\left(u,v,A,B,\phi\right)$ (18)

with generator the vector field

$\mathbf{X}=\xi^{u}\partial_{u}+\xi^{\nu}\partial_{v}+\eta^{A}\partial_{A}+\eta^{B}\partial_{B}+\eta^{\phi}\partial_{\phi}.$
(19)

Let $\mathbf{H}=0,~{}$describe the field equations (9)-(13), the
$\mathbf{H}=0$ is invariant under the action of the latter infinitesimal
transformation if and only if there exist a function $\lambda$ such that the
following condition is satisfied

$\mathbf{X}^{\left[2\right]}\left(\mathbf{H}\right)=\lambda\mathbf{H},$ (20)

or equivalently

$\mathbf{X}^{\left[2\right]}\left(\mathbf{H}\right)=0.$ (21)

$\mathbf{X}^{\left[2\right]}$ is the second extension/prolongation of $X$ in
the jet space
$\left\\{u,v,A,B,A_{,u},A_{,v},A_{uv},A_{uu},A_{,vv},B_{,u},B_{,v},B_{uv},B_{uu},B_{,vv}\right\\}$
defined as

$\displaystyle\mathbf{X}^{\left[2\right]}$
$\displaystyle=\mathbf{X}+\eta_{\left[u\right]}^{A}\partial_{A_{,u}}+\eta_{\left[v\right]}^{A}\partial_{A_{,v}}+\eta_{\left[u\right]}^{B}\partial_{B_{,u}}$
$\displaystyle+\eta_{\left[v\right]}^{B}\partial_{B_{,v}}+\eta_{\left[u\right]}^{\phi}\partial_{\phi_{,u}}+\eta_{\left[v\right]}^{\phi}\partial_{\phi_{,v}}$
(22)
$\displaystyle\eta_{\left[uu\right]}^{A}\partial_{A_{,uu}}+\eta_{\left[vv\right]}^{A}\partial_{A_{,vv}}+\eta_{\left[uv\right]}^{A}\partial_{A_{,uv}}+$
$\displaystyle\eta_{\left[uu\right]}^{B}\partial_{B_{,uu}}+\eta_{\left[vv\right]}^{B}\partial_{B_{,vv}}+\eta_{\left[uv\right]}^{B}\partial_{B_{,uv}}+$
(23)
$\displaystyle\eta_{\left[uu\right]}^{\phi}\partial_{\phi_{,uu}}+\eta_{\left[vv\right]}^{\phi}\partial_{\phi_{,vv}}+\eta_{\left[uv\right]}^{\phi}\partial_{\phi_{,uv}},$
(24)

with

$\eta_{\left[\alpha\right]}^{A}=D_{\alpha}\eta^{A}-A_{,\beta}D_{\alpha}\xi^{\beta},$
(25)
$\eta_{\left[\alpha\right]}^{B}=D_{\alpha}\eta^{B}-B_{,\beta}D_{\alpha}\xi^{\beta},$
(26)
$\eta_{\left[\alpha\right]}^{\phi}=D_{\alpha}\eta^{\phi}-A_{,\beta}D_{\alpha}\xi^{\beta},$
(27)
$\eta_{\left[\alpha\beta\right]}^{A}=D_{\alpha}\eta_{\left[\beta\right]}^{A}-A_{,\alpha\beta}D_{\alpha}\xi^{\beta},$
(28)
$\eta_{\left[\alpha\beta\right]}^{B}=D_{\alpha}\eta_{\left[\beta\right]}^{B}-B_{,\alpha\beta}D_{\alpha}\xi^{\beta},$
(29)
$\eta_{\left[\alpha\beta\right]}^{\phi}=D_{\alpha}\eta_{\left[\beta\right]}^{\phi}-\phi_{,\alpha\beta}D_{\alpha}\xi^{\beta},$
(30)

where $\alpha,\beta=u,v$ and we consider Einstein’s summation convention.

The Lie symmetry condition (20) applied to the field equations (9)-(13) leads
to a system of partial differential equations, from which the coefficients of
the vector field $\mathbf{X}$ are determined. The solution of these symmetry
conditions is contingent upon the functional form of the scalar field
potential $V\left(\phi\right)$, resulting in different allowed Lie symmetries
for various potentials.

From the Lie symmetry conditions it follows that the vector field $X$
components $\xi^{u},~{}\xi^{v},~{}\eta^{A},\eta^{B}$ and$~{}\eta^{\phi}$ are

$\xi^{u}=\xi^{u}\left(u\right)~{},~{}\xi^{v}=\xi^{v}\left(v\right)~{},~{}\eta^{A}=a_{0}$
$\eta^{B}=a_{1}\phi-\xi_{,u}^{u}-\xi_{,v}^{v}+a_{3}~{},~{}\eta^{\phi}=3a_{1}\phi+\alpha_{2},$

with constraint equations for the potential function $V\left(\phi\right)$ and
the free parameters $a_{1}$ and $a_{2}$,

$a_{1}AV_{,\phi}+a_{1}\phi V+a_{2}V_{,\phi}+2a_{3}V=0,$
$7a_{2}V_{,\phi}+21a_{1}AV_{,\phi}+14a_{1}\phi-3a_{1}V_{,\phi}+14a_{3}V=0,$
$3a_{1}AV_{,\phi}+a_{2}V_{,\phi\phi}+2a_{1}V+2a_{1}\phi
V_{,\phi}+2a_{3}V_{,\phi}=0.$

For arbitrary potential function the latter system admits a solution with
$a_{1}=0$, $a_{2}=0,~{}a_{3}=0$; that is, the allowed Lie symmetries are

$X_{1}=F\left(u\right)\partial_{u}-\frac{1}{2}F_{,u}\partial_{B}~{},~{}X_{2}=G\left(v\right)\partial_{v}-\frac{1}{2}G_{,v}\partial_{B}~{},~{}X_{3}=\partial_{A}.$

Vector fields $X_{1}$ and $X_{2}$ describe infinite number of Lie symmetry
vectors. The commutators of the allowed Lie symmetries are zero.

In the case where the potential function is the cosmological constant it
follows that $a_{1}=0$, $a_{3}=0$, i.e. $V\left(\phi\right)=\Lambda$, the
allowed Lie symmetries by the field equations (9)-(13) are

$X_{1}~{},~{}X_{2}~{},~{}X_{3}~{},~{}X_{4}=\partial_{\phi}\text{,}$

with zero commutators.

In the limit where $\Lambda=0$, that is, the scalar field is massless the
admitted symmetry vectors are

$X_{1}~{},~{}X_{2}~{},~{}X_{3}~{},~{}X_{4}~{},~{}X_{5}=\partial_{B}~{},~{}X_{6}=\frac{1}{3}\phi\partial_{B}+A\partial_{\phi}\text{,}$

with non-zero commutator $\left[X_{3},X_{6}\right]=X_{4}$.

Finally, for the exponential potential function
$V\left(\phi\right)=V_{0}e^{\lambda\phi}$, we determine $a_{1}=0$,
$a_{3}=-\frac{\lambda}{2}a_{2}$; from where it follows that the admitted Lie
symmetries for the cosmological field equations (9)-(13) are

$X_{1}~{},~{}X_{2}~{},~{}X_{3}\text{~{}and
}X_{7}=-\frac{\lambda}{2}\partial_{B}+\partial_{\phi}\text{,}$

where all the commutators are zero.

###  III.1 Similarity solutions

We proceed with the application of the Lie symmetries to reduce the field
equations and when it is feasible to determine new solutions.

Let $W=W\left(u,v,A,B\right)$ be a function and $X$ be a symmetry vector for
the function $W$, i.e. $X\left(W\right)=0$. Then, there we can always define a
new set of variables
$\left\\{u,v,A,B\right\\}\rightarrow\left\\{\bar{u},\bar{v},\bar{A},\vec{B}\right\\}$,
where the symmetry vector can be written as $X=\partial_{\bar{u}}$, such that
the symmetry condition becomes $\frac{\partial W}{\partial\bar{u}}=0$, that
is, $W=W\left(\bar{v},\bar{A},\bar{B}\right)$. This transformation is called a
similarity transformation.

In order to demonstrate the application of the Lie symmetries, in the
following lines we consider specific forms of the vector fields $X_{1}$ and
$X_{2}$ where closed-form similarity solutions are recovered.

####  III.1.1 Cosmological solution

Consider now the Lie symmetry vector $X_{1-2}^{0}=\partial_{u}-\partial_{v}$,
which is generated by the $X_{1},~{}X_{2}$ for $F\left(u\right)=1$ and
$G\left(u\right)=1$. The similarity transformation reads
$\bar{t}=u-v~{},~{}\bar{y}=u+v$, $A=A\left(\bar{t}\right)$ and
$B=B\left(\bar{t}\right)$. In the new variables the line element (8) becomes

$ds^{2}=e^{2B\left(\bar{t}\right)}\left(-2d\bar{t}^{2}+2d\bar{y}^{2}\right)+e^{2A\left(\bar{t}\right)}\left(dx_{1}^{2}+dx_{2}^{2}+dx_{3}^{2}\right).$
(31)

The latter spacetime describes a five-dimensional homogeneous and anisotropic
geometry, the extension of the Bianchi I spacetime in the five dimensions with
two scale factors. The process of deriving the analytic solution for the field
equations, whether with or without the exponential potential, follows a
construction method similar to that presented in ns010  for the four-
dimensional Bianchi I spacetime.

Similarly, the vector field $X_{1-2}^{\prime 0}=\partial_{u}+\partial_{v}$
induces a similarity transformation that characterizes a static and
anisotropic five-dimensional geometry with line element

$ds^{2}=e^{2B\left(\bar{y}\right)}\left(-2d\bar{t}^{2}+2d\bar{y}^{2}\right)+e^{2A\left(\bar{y}\right)}\left(dx_{1}^{2}+dx_{2}^{2}+dx_{3}^{2}\right).$
(32)

####  III.1.2 Plane-wave solutions

The application of Lie symmetries $X_{1}$ or $X_{2}$ yields plane-wave
solutions. However, it is worth noting that such solutions do not exist for a
nonzero potential bb2 . Conversely, vacuum plane-wave solutions have been
investigated in bb1 .

####  III.1.3 Travel-wave solution

Vector fields $X_{1}$ and $X_{2}$ correspond to an infinite number of
invariant transformations. Assume now that $F\left(u\right)$ and
$G\left(v\right)$ are quadratic functions, i.e. $F\left(u\right)=u^{2}$ and
$G\left(v\right)=v^{2}$, then $X_{1}^{2}=u^{2}\partial_{u}-u\partial_{B}$ and
$X_{2}^{2}=v^{2}\partial_{v}-v\partial_{B}$. We consider the Lie symmetry
vector $X_{1-2}^{2}=X_{1}-X_{2}$ from where we derive the similarity
transformation

$A\left(u,v\right)=a\left(w\right)~{},~{}B=-\ln\left(uv\right)+b\left(w\right)~{},~{}\phi\left(u,v\right)=\psi\left(w\right),$
(33)

where the new independent variable $w$ is defined as $w=\frac{u+u}{vu}$, or in
the $\left(t,y\right)$ coordinates, $w=\frac{2t}{t^{2}-y^{2}}$.

By replacing (33) in the field equations (9)-(13) we end with a system of
nonlinear ordinary differential equations.

The reduced system is

$\displaystyle 6\left(a_{,ww}+3a_{,w}^{2}\right)-V_{0}e^{2b+\lambda\psi}$
$\displaystyle=0,$ (34) $\displaystyle
3\left(a_{,ww}+a_{,w}^{2}\right)-6a_{,w}b_{,w}+\psi_{,w}^{2}$
$\displaystyle=0,$ (35) $\displaystyle
4\left(\psi_{,ww}+3a_{,w}\psi_{,w}\right)+V_{0}e^{2b+\lambda\psi}$
$\displaystyle=0,$ (36) $\displaystyle
4\left(2a_{,ww}+3a_{,w}^{2}\right)+4b_{,ww}+2\psi_{,w}^{2}+V_{0}e^{2b+\lambda\psi}$
$\displaystyle=0,$ (37)

For zero potential, $V_{0}=0$, the latter dynamical system admits the closed-
form solution

$A\left(u,v\right)=\ln\left(f_{0}w^{p_{1}}\right)~{},~{}B=-\ln\left(uv\right)+\ln\left(b_{0}w^{p_{2}}\right)~{},~{}\phi\left(u,v\right)=\frac{1}{\lambda}\ln\left(\phi_{0}w^{p_{3}}\right),$
(38)

with constraint equations

$p_{1}=\frac{1}{3}~{},~{}p_{2}=\frac{3p_{3}^{2}-2\lambda^{2}}{6\lambda^{2}}\text{.}$

For this exact solution the line element (8) becomes

$ds^{2}=\left(\frac{b_{0}w^{p_{2}}}{uv}\right)^{2}\left(-dudv\right)+\left(f_{0}w^{\frac{1}{3}}\right)^{2}\left(dx_{1}^{2}+dx_{2}^{2}+dx_{3}^{2}\right),$
(39)

or

$ds^{2}=\left(\frac{b_{0}}{\left(t^{2}-y^{2}\right)}\left(\frac{t^{2}-y^{2}}{2t}\right)^{p_{2}}\right)^{2}\left(-dt^{2}+dy^{2}\right)+\left(f_{0}\left(\frac{2t}{t^{2}-y^{2}}\right)^{\frac{1}{3}}\right)^{2}\left(dx_{1}^{2}+dx_{2}^{2}+dx_{3}^{2}\right).$
(40)

In the four-dimensional limit with $y=0$, the cosmological solution reads

$ds^{2}=\left(\bar{b}_{0}\left(t\right)^{p_{2}-1}\right)^{2}\left(-dt^{2}\right)+\left(\bar{f}_{0}\left(\frac{1}{t}\right)^{\frac{1}{3}}\right)^{2}\left(dx_{1}^{2}+dx_{2}^{2}+dx_{3}^{2}\right).$
(41)

The scaling solution described represents a universe in which the cosmological
fluid behaves as an ideal gas with a constant equation of state parameter
$w_{\text{eff}}=2p_{2}-1$. For values of $p_{2}<\frac{1}{2}$, the universe
experiences acceleration.

However, for $V_{0}\neq 0$, it is observed that the dynamical system does not
admit a real solution.

Similarly, we can utilize the Lie symmetry vectors $X_{1}$ and $X_{2}$ to
construct additional solutions in a similar manner.

##  IV Conclusions

In our study, we conducted an extensive classification of Lie symmetries for
the cosmological field equations in a five-dimensional spacetime with brane-
gravity featuring a bulk scalar field. By imposing the condition that the
field equations must possess non-trivial Lie symmetries, we derived the scalar
field potential. Through this approach, we determined the scalar field
potential in a manner that ensures the field equations admit Lie symmetries,
which form various Lie algebras.

For an arbitrary potential function, the field equations exhibit an infinite
number of Lie symmetries, forming the Lie algebra
$G_{\infty}=2A_{\infty}\oplus A_{1}$ pat1 . In the special case of a massless
scalar field, the admitted Lie symmetries give rise to the Lie algebra
$G_{1}=2A_{\infty}\oplus 2A_{1}$. Furthermore, when the scalar field potential
takes on a constant value, representing the cosmological constant, the Lie
symmetries form the $G_{2}=2A_{\infty}\oplus A_{1}\oplus A_{2,1}$ Lie algebra.
Additionally, for the exponential potential, the Lie symmetries yield the
$G_{2}=2A_{\infty}\oplus 2A_{1}$ Lie algebra.

n the last step of our research, we harnessed the Lie symmetries to establish
similarity transformations. These transformations played a crucial role in
reducing the field equations into ordinary differential equations and enabled
us to derive closed-form solutions. This application of Lie symmetries has
proven to be highly beneficial in the study of gravitational theories,
contributing valuable insights to the subject of symmetry analysis in this
field.

Our findings underscore the usefulness and potency of the Lie analysis as a
method for exploring gravitational theories. Building on these results, we
intend to expand our research in the future by investigating the determination
of conservation laws for the gravitational field equations. This forthcoming
work will further enrich our understanding of the dynamics and conservation
principles within the realm of gravitational theories.

Data Availability Statements: Data sharing is not applicable to this article
as no datasets were generated or analyzed during the current study.

###### Acknowledgements.

AP was partially financially supported by the National Research Foundation of
South Africa (Grant Numbers 131604). AP thanks the support of Vicerrectoría de
Investigación y Desarrollo Tecnológico (Vridt) at Universidad Católica del
Norte through Núcleo de Investigación Geometría Diferencial y Aplicaciones,
Resolución Vridt No - 098/2022.

## References

  * (1) A. Einstein, Grundgedanken der allgemeinen Relativit¨atstheorie und Anwendung dieser Theorie in der Astronomie, Preussische Akademie der Wissenschaften, Sitzungsberichte, 315, (1915)
  * (2) B.P. Abbott et al., Observation of Gravitational Waves from a Binary Black Hole Merger, Phys. Rev. Lett. 116, 061102 (2016)
  * (3) M. Kramer et al., Strong-Field Gravity Tests with the Double Pulsar, Phys. Rev. X 11, 041050 (2021)
  * (4) C.M. Will, The Confrontation between General Relativity and Experiment, Living Rev. Relativ. 17, 4 (2014)
  * (5) A.G. Riess et al., Observational evidence from Supernovae for an accelerating universe and cosmological constant, Astron. J. 116, 1009 (1998)
  * (6) M. Tegmark et al., The Three-Dimensional Power Spectrum of Galaxies from the Sloan Digital Sky Survey, Astrophys. J. 606, 702 (2004)
  * (7) M. Kowalski et al., Improved Cosmological Constraints from New, Old and Combined Supernova Datasets, Astrophys. J. 686, 749 (2008)
  * (8) S. Weinberg, The cosmological constant problem, Rev. Mod. Phys. 61, 1 (1989)
  * (9) T. Padmanabhan, Cosmological Constant - the Weight of the Vacuum, Phys. Rept. 380, 235 (2003)
  * (10) P. Ratra and L. Peebles, Cosmological consequences of a rolling homogeneous scalar field, Phys. Rev. D 37, 3406 (1988)
  * (11) T.P. Sotiriou, f(R) gravity and scalar–tensor theory, Class. Quant. Grav. 23, 5117 (2006)
  * (12) J.C. Fabris, T.C. da C. Guio, M.H. Daouda and O.F. Piatella, Scalar models for the generalized Chaplygin gas and the structure formation constraints, Gravitation and Cosmology 17, 259 (2011)
  * (13) N. Arkani-Hamed, S. Dimopoulos and G. Dvali, The Hierarchy Problem and New Dimensions at a Millimeter, Phys. Lett. B 429, 263 (1998)
  * (14) I. Antoniadis, N. Arkani-Hamed, S. Dimopoulos and G. Dvali, New Dimensions at a Millimeter to a Fermi and Superstrings at a TeV, Phys. Lett. B 436, 257 (1998)
  * (15) L. Randall and R. Sundrum, Large Mass Hierarchy from a Small Extra Dimension, Phys. Rev. Lett. 83, 3370 (1999)
  * (16) R. Maartens and K. Koyama, Brane-World Gravity, Living Rev. Relat. 13, 5 (2010)
  * (17) N. Arkani-Hamed, S Dimopoulos, N Kaloper and R. Sundrum, A Small Cosmological Constant from a Large Extra Dimension, Phys. Lett. B 480, 193 (2000)
  * (18) Z. Kakushadze, Bulk Supersymmetry and Brane Cosmological Constant, Phys. Lett. B 489, 207 (2000)
  * (19) S. Kachru, M. Schulz and E. Silverstein, Self-tuning flat domain walls in 5d gravity and string theory, Phys. Rev. D 62, 045021 (2000)
  * (20) S. Forste, Z. Lalak, S. Lavignac and H.P. Nilles, A Comment on Self-Tuning and Vanishing Cosmological Constant in the Brane World, Phys. Lett. B 481, 360 (2000)
  * (21) O. DeWolfe, D.Z. Freedman, S.S. Gubser and A. Karch, Modeling the fifth dimension with scalars and gravity, Phys. Rev. D 62, 046008 (2000)
  * (22) O. Bertolami and C. Carvalho, Spontaneous symmetry breaking in the bulk and the localization of fields on the brane, Phys. Rev. D 76, 104048 (2007)
  * (23) G. T. Horowitz, I. Low and A. Zee, Self-tuning in an outgoing brane wave model, Phys. Rev. D 62, 086005 (2000)
  * (24) D. Langlois and M. Rodriguez-Martinez, Brane cosmology with a bulk scalar field, Phys. Rev. D 64, 123507 (2001)
  * (25) É.É. Flanagan, S.-H. Henry Tye and I. Wasserman, Brane world models with bulk scalar fields, Phys. Lett. B 522, 155 (2001)
  * (26) G.W. Bluman and S. Kumei, Symmetries of Differential Equations, Springer-Verlag, New York (1989)
  * (27) H. Stephani, Differential Equations: Their Solutions Using Symmetry, Cambridge University Press, New York (1989)
  * (28) P.J. Olver, Applications of Lie Groups to Differential Equations, Springer-Verlag, New York (1993)
  * (29) N.H. Ibragimov, CRC Handbook of Lie Group Analysis of Differential Equations, Volume I: Symmetries, Exact Solutions, and Conservation Laws, CRS Press LLC, Florida (2000)
  * (30) C.G. Torre and I.M. Anderson, Phys. Rev. Lett. 70, 23 (1993)
  * (31) L. Marchildon, Lie Symmetries of Einstein’s Vacuum Equations in N Dimensions, J. Nonlin. Math. Phys. 5, 68 (1998)
  * (32) H. Stephani, Symmetries of Einstein’s field equations with a perfect fluid source as examples of Lie–Bäcklund symmetries, J. Math. Phys. 29, 1650 (1988)
  * (33) S. Basilakos, M. Tsamparlis and A. Paliathanasis, Using the Noether symmetry approach to probe the nature of dark energy, Phys. Rev. D 83, 103512 (2011)
  * (34) R. de Ritis, G. Marmo, G. Platania, C. Rubano and P. Scudellaro, New approach to find exact solutions for cosmological models with a scalar field, Phys. Rev. D 42, 1091 (1990)
  * (35) A. Paliathanasis, M. Tsamparlis, G. Leon and A. Paliathanasis, New Conservation Laws and Exact Cosmological Solutions in Brans–Dicke Cosmology with an Extra Scalar Field, Symmetry 13, 1364 (2021)
  * (36) P. A. Terzis, N. Dimakis and T. Christodoulakis, Noether analysis of scalar-tensor cosmology, Phys. Rev. D 90, 123543 (2014)
  * (37) A. Paliathanasis, M. Tsamparlis, S. Basilakos and S. Capozziello, Scalar-Tensor Gravity Cosmology: Noether symmetries and analytical solutions, Phys. Rev. D 89, 063532 (2014)
  * (38) S. Capozziello, K. Dialektopoulos and S.V. Sushkov, Classification of the Horndeski cosmologies via Noether symmetries, Eur. Phys. J. C 78, 447 (2018)
  * (39) N. Sk., Noether symmetry in f(T) teleparallel gravity, Phys. Lett. B 775, 100 (2017)
  * (40) M. Tsamparlis and A. Paliathanasis, Symmetries of Differential Equations in Cosmology, Symmetry 10, 233 (2018)
  * (41) M. Tsamparlis and A. Paliathanasis, The geometric nature of Lie and Noether symmetries, Gen. Relat. Grav. 43, 1861 (2011)
  * (42) J. Patera, R.T. Sharp, P. Wintenritz and H. Zassenhaus, Invariants of real low dimension Lie algebras, J. Math. Phys. 17, 986 (1976)