RBI-ThPhys-2023-11

#  Optimizing the Evolution of Perturbations in the $\Lambda$CDM Universe

Nicholas Choustikov<EMAIL_ADDRESS>Sub-department of
Astrophysics, University of Oxford, Keble Road, Oxford OX1 3RH, UK  Institute
of Astronomy, Madingley Road, Cambridge, CB3 0HA, UK  Zvonimir Vlah
<EMAIL_ADDRESS>Division of Theoretical Physics, Ruđer Bošković Institute, 10000
Zagreb, Croatia  Kavli Institute for Cosmology Cambridge, Madingley Road,
Cambridge, CB3 0HA, UK  DAMTP, Centre for Mathematical Sciences, Wilberforce
Road, Cambridge CB3 0WA, UK  Anthony Challinor<EMAIL_ADDRESS>Institute of Astronomy, Madingley Road, Cambridge, CB3 0HA, UK  Kavli
Institute for Cosmology Cambridge, Madingley Road, Cambridge, CB3 0HA, UK
DAMTP, Centre for Mathematical Sciences, Wilberforce Road, Cambridge CB3 0WA,
UK

###### Abstract

Perturbation theory is a powerful tool for studying large-scale structure
formation in the universe and calculating observables such as the power
spectrum or bispectrum. However, beyond linear order, typically this is done
by assuming a simplification in the time-dependence of gravitational-coupling
kernels between the matter and velocity fluctuations. Though the true
dependencies are known for Lambda cold dark matter cosmologies, they are
ignored due to the computational costs associated with considering them in
full and, instead, are replaced by simpler dependencies valid for an
Einstein–de-Sitter cosmology. Here we develop, implement and demonstrate the
effectiveness of a new numerical method for finding the full dynamical
evolution of these kernels to all perturbative orders based upon spectral
methods using Chebyshev polynomials. This method is found to be orders of
magnitude more efficient than direct numerical solvers while still producing
highly accurate and reliable results. A code implementation of the Chebyshev
spectral method is then presented and characterised. The code has been made
publicly available alongside this paper. We expect our method to be of use for
interpretation of upcoming galaxy clustering measurements.

##  I Introduction

Studying the history of large-scale structure (LSS) formation in our Universe
is crucial to modern cosmology. The highly structured cosmic web of galaxy
clusters, sheets, walls, filaments and voids present in the late-time Universe
arises from evolution under gravity of small, primordial fluctuations in the
density and velocity of matter, believed to have been sourced during a period
of cosmic inflation. LSS therefore encodes information on the primordial
fluctuations as well as the expansion history, geometry and matter content of
the Universe, which affect the subsequent evolution of the perturbations.

Theoretical studies of LSS aim to predict the statistical properties of the
clustering of the matter density (and associated velocity), for example, the
power spectrum and higher-point correlation functions. These can be compared
with the observed statistics of the clustering of galaxies since the galaxy
over-density on large scales traces the matter over-density. Such comparisons
are complicated by two issues, however. Firstly, galaxies need not necessarily
follow the underlying matter distribution exactly, although on the largest
cosmological scales the relation for galaxies is essentially linear with a
constant of proportionality known as galaxy bias, while on smaller, mildly
nonlinear scales, corrections can be treated perturbatively [1]. Secondly, we
have the issue that each galaxy’s redshift depends not only on distance but
also on its peculiar velocity via the Doppler effect [2, 3]. Furthermore,
these velocities are not random but instead correlate with the matter density
field itself. This alters galaxy statistics by producing redshift-space
distortions (RSDs), also introducing features such as the well-known “Fingers
of God” [3]. These issues are compounded by the fact that we are now entering
the era of high-precision cosmology with large galaxy surveys such as Euclid
[4], DESI [5], Rubin [6] and Roman [7]. Here, it becomes necessary to quantify
these effects fully, lest they introduce greater uncertainties or issues of
biased parameter estimation.

Given the small amplitude of the primordial fluctuations, perturbation theory
is a very powerful tool for calculating the evolution of fluctuations in
cosmological fields. For early-time observables, such as the cosmic microwave
background (CMB), linear perturbation theory is an excellent description.
Accurate predictions for the CMB power spectra can be calculated in this way
(e.g., [8, 9]) and their comparison with measurements of the CMB temperature
and polarization anisotropies has been instrumental in establishing the
standard model of cosmology, the Lambda cold dark matter ($\Lambda$CDM) model,
and determining its parameters to high precision (e.g., [10]). However, linear
perturbation theory fails to describe the late-time Universe, where,
particularly at smaller scales, complex non-linearities take hold and baryon
physics begins to play a role [11]. Here, higher-order perturbation theory is
crucial to building analytic models. The standard procedure in these
perturbative approaches is to use the quasi-static Einstein–de-Sitter (EdS)
approximation, henceforth referred to as the quasi-EdS approximation (qEDS).
Here, a blend of models is used, with non-linear interaction kernels
calculated in the dramatically simpler EdS (completely matter-dominated)
universe. These are then combined with linear growth-rate results for the full
$\Lambda$CDM dynamical Universe in order to extract observables [12]. In
essence, this approximation assumes that any non-linearities are constant in
time, to then be scaled up to their present-day values.

The qEDS approach has been fairly successful thus far, yielding percent-level
accuracy for all necessary two-point statistics close to the linear regime
[13]. Such models have already yielded concrete predictions, in particular for
the one-loop power spectrum and tree-level bispectrum, and these have been
used to extract relevant cosmological information from galaxy surveys [14, 15,
16]. However, for future surveys, this accuracy will not be enough.
Particularly in the context of RSDs, the full dynamics of $\Lambda$CDM must be
taken into account. A large body of work has already explored the extent of
these effects at the level of one-loop results [17, 18, 19, 20, 21, 22, 23,
24], and recently to two-loop order[25, 26, 27]. Unfortunately, accounting for
the full $\Lambda$CDM dynamics introduces several issues, including that some
of these features are slow to calculate using brute-force numerical methods
(especially at higher perturbative order). This makes exploration of parameter
space computationally challenging when interpreting survey data. Therefore,
this paper aims to optimize the process by which the full $\Lambda$CDM
dynamics can be calculated at each perturbative order, using previous results
for time-dependent coefficients as a starting point [19, 26].

We tackle this problem with a novel numerical method, expanding the time-
dependent interaction kernels in shifted Chebyshev polynomials [28] as in the
Chebyshev Spectral Method (CSM; see, e.g., [29] for the explicit, matrix-based
approach that we follow). Doing so, we convert the process of calculating
these dynamical coefficients into a linear-algebra problem that can readily be
solved. We have implemented this method into a new Python library that we
present below. The code is freely available at
https://github.com/Chousti/CSMethod.git, with the hope that it can be
integrated into likelihood analyses of forthcoming survey data, allowing more
accurate parameter estimates to be derived from these data. While we focus on
$\Lambda$CDM cosmologies here, it is expected that the method can be
straightforwardly extended to more generalised cosmologies, including features
such as clustered quintessence [30, 31].

This paper is arranged as follows. In Sec. II we review certain theoretical
results underpinning the problem. Section III discusses the implementation of
our spectral method with explanations of the code philosophy, and Sec. IV
presents results and code tests. Finally, we discuss our results in Sec.V.
More extensive auxiliary functions are given in Appendix A.

All numerical results (unless stated otherwise) are calculated using the
Planck best-fit $\Lambda$CDM cosmology, with parameters as follows: the
present-day matter density parameter $\Omega_{m_{0}}=0.315$; present-day dark-
energy density parameter $\Omega_{\Lambda_{0}}=0.685$, giving a flat universe;
and Hubble constant
$H_{0}=$67.74\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{Mpc}}^{-1}$$
[10].

##  II Dynamics in the $\Lambda$CDM Universe

###  II.1 Deriving the Equations of Motion

On the large scales for which a perturbative treatment is valid, we may
approximate the cosmic matter density field as a single, self-gravitating,
cold, pressureless fluid. Baryons follow the dark matter on these scales and
pressure and other baryonic effects are negligible compared to gravitational
interactions. The matter density may therefore be characterised by its density
contrast $\delta(\bm{x},a)$ and peculiar velocity $\bm{v}(\bm{x},a)$, where
$\bm{x}$ is comoving position and $a$ is the scale factor (we use $a$ and
conformal time $\tau$, interchangeably, as time variables in the following).
It is convenient to decompose the peculiar velocity into its divergence
$\theta=\bm{\nabla}\cdot\bm{v}(\bm{x},a)$ and vorticity
$\bm{\omega}=\bm{\nabla}\times\bm{v}(\bm{x},a)$. These quantities obey the
following equations of motion [13, 1]:

$\displaystyle\frac{\partial\delta_{\bm{k}}}{\partial\tau}+\theta_{\bm{k}}$
$\displaystyle=-\int_{\bm{q}_{1},\bm{q}_{2}}(2\pi)^{3}\delta^{D}_{\bm{k}-\bm{q}_{12}}\alpha(\bm{q}_{1},\bm{q}_{2})\theta_{\bm{q}_{1}}\delta_{\bm{q}_{2}},$
(1)
$\displaystyle\frac{\partial\theta_{\bm{k}}}{\partial\tau}+\mathcal{H}\theta_{\bm{k}}+\frac{3}{2}\Omega_{m}\mathcal{H}^{2}\delta_{\bm{k}}$
$\displaystyle=-\int_{\bm{q}_{1},\bm{q}_{2}}(2\pi)^{3}\delta^{D}_{\bm{k}-\bm{q}_{12}}\beta(\bm{q}_{1},\bm{q}_{2})\theta_{\bm{q}_{1}}\theta_{\bm{q}_{2}},$
(2)
$\displaystyle\frac{\partial\bm{\omega}_{\bm{k}}}{\partial\tau}+\mathcal{H}\bm{\omega}_{\bm{k}}$
$\displaystyle=i\bm{k}\times\int_{\bm{q}_{1},\bm{q}_{2}}(2\pi)^{3}\delta^{D}_{\bm{k}-\bm{q}_{12}}\bm{v}_{\bm{q}_{1}}\bm{\omega}_{\bm{q}_{2}},$
(3)

where $\delta_{\bm{q}}^{D}$ is the Dirac delta function,
$\bm{k}_{12}=\bm{k}_{1}+\bm{k}_{2}$,
$\int_{\bm{q}_{1},\ldots,\bm{q}_{n}}\equiv(2\pi)^{-3n}\int
d^{3}\bm{q}_{1}\cdots d^{3}\bm{q}_{n}$, $\Omega_{m}(a)$ is the matter-density
parameter and $\mathcal{H}=d\ln a/d\ln\tau$ is the conformal Hubble parameter.
The kernels are defined by
$\alpha(\bm{q}_{1},\bm{q}_{2})=1+(\bm{q}_{1}\cdot\bm{q}_{2})/\bm{q}_{1}^{2}$
and
$\beta(\bm{q}_{1},\bm{q}_{2})=(\bm{q}_{12})^{2}(\bm{q}_{1}\cdot\bm{q}_{2})/(2\bm{q}_{1}^{2}\bm{q}_{2}^{2})$.
Equation (3) implies that $\bm{\omega}$ remains zero for all time in the case
of vanishing primordial vorticity, while it decays as $1/a$ at linear order if
vorticity were present in the early universe. As a result, any small initial
vorticity rapidly decays as the universe expands. Therefore, we proceed to
assume the velocity field to be irrotational – an assumption that holds well
up to shell crossing (and the formation of shocks) [12].

The remaining system of Eqs. (1) and (2) are now closed and well-suited to be
solved perturbatively, with an ansatz of

$\displaystyle\delta(\bm{k},a)$
$\displaystyle=\sum_{n=1}^{\infty}\int_{\bm{q}_{1},\ldots,\bm{q}_{n}}(2\pi)^{3}\delta^{D}_{\bm{k}-\bm{q}_{12...n}}F_{n}(\bm{q}_{1},\dots,\bm{q}_{n};a)D_{+}^{n}(a)\delta^{\text{in}}_{\bm{q}_{1}}\dots\delta^{\text{in}}_{\bm{q}_{n}},$
$\displaystyle\theta(\bm{k},a)$
$\displaystyle=-\mathcal{H}f_{+}(a)\sum_{n=1}^{\infty}\int_{\bm{q}_{1},\ldots,\bm{q}_{n}}(2\pi)^{3}\delta^{D}_{\bm{k}-\bm{q}_{12...n}}G_{n}(\bm{q}_{1},\dots,\bm{q}_{n};a)D_{+}^{n}(a)\delta^{\text{in}}_{\bm{q}_{1}}\dots\delta^{\text{in}}_{\bm{q}_{n}}.$
(4)

Here, $F_{n}$ and $G_{n}$ are our solution kernels (defined to be fully
symmetrized with respect to momenta), $D_{\pm}$ are the growing and decaying
linear growth factors satisfying

$\displaystyle
a^{2}\frac{d^{2}D_{\pm}}{da^{2}}+a\left(2+\frac{d\ln\mathcal{H}}{d\ln
a}\right)\frac{dD_{\pm}}{da}-\frac{3}{2}\Omega_{m}(a)D_{\pm}=0,$ (5)

in $\Lambda$CDM cosmologies, $f_{\pm}=d\ln D_{\pm}/d\ln a$ are the associated
linear growth rates and $\delta^{\mathrm{in}}_{\bm{q}}$ describe the initial
density contrast. At linear order, we find trivial results of $F_{1}=G_{1}=1$,
giving

$\delta_{\bm{k}}^{(1)}(a)=D_{+}(a)\delta^{\mathrm{in}}_{\bm{k}},\qquad\theta_{\bm{k}}^{(1)}=-\mathcal{H}(a)f_{+}(a)D_{+}(a)\delta^{\mathrm{in}}_{\bm{k}}.$
(6)

Figure 1: Comparison of the growth factors for the growing mode (left), along
with the logarithmic growth rates of the growing (middle) and decaying (right)
modes, between the $\Lambda$CDM and quasi-EdS universes.

The linear growth factor for the growing mode, and its associated growth rate,
are given by the following:

$\displaystyle D_{+}(a)$
$\displaystyle=\frac{5}{2}H_{0}^{2}\Omega_{m_{0}}H(a)\int_{0}^{a}\frac{dx}{[xH(x)]^{3}},$
(7) $\displaystyle f_{+}(a)$
$\displaystyle=\frac{\Omega_{m_{0}}}{(1-\Omega_{m_{0}})a^{3}+\Omega_{m_{0}}}\left(\frac{5a}{2D_{+}(a)}-\frac{3}{2}\right),$
(8)

where $H(a)=H_{0}[\Omega_{m_{0}}a^{-3}+(1-\Omega_{m_{0}})]^{1/2}$. We have
normalised $D_{+}(a)$ so that it approaches $a$ as $a\rightarrow 0$. The
linear growth factor for the decaying mode is $D_{-}(a)\propto H(a)$ and the
associated growth rate is

$\displaystyle
f_{-}(a)=-\frac{3}{2}\frac{\Omega_{m_{0}}}{(1-\Omega_{m_{0}})a^{3}+\Omega_{m_{0}}}=-\frac{3}{2}\Omega_{m}(a).$
(9)

Although we only consider growing-mode solutions, $f_{-}(a)$ appears below in
the equation of motion for the $F_{n}$ and $G_{n}$ kernels through $d\ln
H/d\ln a$, or, equivalently, $\Omega_{m}(a)$. Figure 1 shows the growth factor
$D_{+}(a)$ along with the logarithmic growth rates $f_{\pm}$ in both the
$\Lambda$CDM and EdS universes. In the latter, $D_{+}(a)=a$, $f_{+}(a)=1$ and
$f_{-}(a)=-3/2$. The behaviour of these functions in $\Lambda$CDM differs from
the EdS limits at late times and these deviations eventually compound to
affect observables. The constancy of $f_{\pm}$ in EdS helps to explain how
this approximation is viewed as ‘static’, as these growth rates do not
introduce any dynamics into the system. Typically, the quasi-EdS approximation
is made, whereby the kernels $F_{n}$ and $G_{n}$ are assumed to be time-
independent and are replaced by their counterparts in EdS, while the full
$\Lambda$CDM growth factor ($D_{+}$) is used in Eq. (4). This is done due to
the fact that such kernels are much easier to compute in the simpler EdS case
[12, 13, 19, 26].

Dropping the qEdS approximation, these kernels obey the following equations of
motion:

$\displaystyle\frac{a}{f_{+}}\frac{dF_{n}}{da}+nF_{n}-G_{n}=h_{\alpha}^{(n)}(\bm{q}_{1},...,\bm{q}_{n},a),$
$\displaystyle\frac{a}{f_{+}}\frac{dG_{n}}{da}+(n-1)G_{n}-\frac{f_{-}}{f_{+}^{2}}(G_{n}-F_{n})=h_{\beta}^{(n)}(\bm{q}_{1},...,\bm{q}_{n},a),$
(10)

with source terms given by Eq. (8) in [19]. In order to proceed, we assume a
separable solution at each order, with ansatz [19, 26]

$\displaystyle F_{n}(\bm{q}_{1},...,\bm{q}_{n};a)$
$\displaystyle=\sum_{l=1}^{N(n)}\lambda_{n}^{(l)}(a)\;H_{n}^{(l)}(\bm{q}_{1},...,\bm{q}_{n}),$
$\displaystyle G_{n}(\bm{q}_{1},...,\bm{q}_{n};a)$
$\displaystyle=\sum_{l=1}^{N(n)}\kappa_{n}^{(l)}(a)\;H_{n}^{(l)}(\bm{q}_{1},...,\bm{q}_{n}),$
(11)

where $H_{n}^{(l)}$ are the momentum operators (see [26]) and
$\lambda_{n}^{(l)}$ and $\kappa_{n}^{(l)}$ are the time-dependent
coefficients, which are the focus of this paper. The numbering function $N(n)$
gives us the simplest way to ensure the full dimension of each kernel is
accounted for, albeit allowing for some redundancies. These are imposed
primarily through physical constraints such as conservation of mass and
momentum, along with the equivalence principle — though these are not explored
here. The $N(n)$ are found recursively, with explicit form given in Appendix
A; the first few terms of which are $N(1)=1$, $N(2)=2$, $N(3)=6$, $N(4)=25$
and $N(5)=111$. For the purpose of this paper, we focus on the time-dependent
coefficients, which can be written as [26]

$\displaystyle\lambda_{n}^{(l)}(a)$
$\displaystyle=\delta^{K}_{\frac{n}{2},\lfloor\frac{n}{2}\rfloor}\sum_{i=1}^{N(n/2)}\bigg{[}\sum_{j=1}^{N(n/2)}[W_{\alpha;\frac{n}{2},\frac{n}{2}}^{n(i,j)}]\delta^{K}_{l,\phi_{1}}+\sum_{j=i}^{N(n/2)}[W_{\beta;\frac{n}{2},\frac{n}{2}}^{n(i,j)}]\delta^{K}_{l,\phi_{2}}\bigg{]}\qquad\qquad\qquad\qquad\qquad$
$\displaystyle\qquad+\sum_{m=1}^{\lfloor(n-1)/2\rfloor}\sum_{i=1}^{N(m)}\sum_{j=1}^{N(n-m)}\bigg{[}[W_{\alpha;m,n-m}^{n(i,j)}]\delta^{K}_{l,\phi_{3}}+[W_{\alpha;n-m,m}^{n(j,i)}]\delta^{K}_{l,\phi_{4}}+[W_{\beta;m,n-m}^{n(i,j)}]\delta^{K}_{l,\phi_{5}}\bigg{]},$
(12) $\displaystyle\kappa_{n}^{(l)}(a)$
$\displaystyle=\delta^{K}_{\frac{n}{2},\lfloor\frac{n}{2}\rfloor}\sum_{i=1}^{N(n/2)}\bigg{[}\sum_{j=1}^{N(n/2)}[U_{\alpha;\frac{n}{2},\frac{n}{2}}^{n(i,j)}]\delta^{K}_{l,\phi_{1}}+\sum_{j=i}^{N(n/2)}[U_{\beta;\frac{n}{2},\frac{n}{2}}^{n(i,j)}]\delta^{K}_{l,\phi_{2}}\bigg{]}\qquad\qquad\qquad\qquad\qquad$
$\displaystyle\qquad+\sum_{m=1}^{\lfloor(n-1)/2\rfloor}\sum_{i=1}^{N(m)}\sum_{j=1}^{N(n-m)}\bigg{[}[U_{\alpha;m,n-m}^{n(i,j)}]\delta^{K}_{l,\phi_{3}}+[U_{\alpha;n-m,m}^{n(j,i)}]\delta^{K}_{l,\phi_{4}}+[U_{\beta;m,n-m}^{n(i,j)}]\delta^{K}_{l,\phi_{5}}\bigg{]},$
(13)

for $n>1$. For the case $n=1$, we have $\lambda^{(1)}_{1}=\kappa^{(1)}_{1}=1$.
In these expressions, $\delta^{K}$ is the Kronecker delta and $\phi_{k}$ are
bijective maps of their arguments ($n$ and all of the integers being summed
over) that serve to identify the correct coefficient for each momentum
operator. Their form is given explicitly in Appendix A. For
$\lambda_{n}^{(l)}$, exactly one of the $W_{\alpha}(a)$ or $W_{\beta}(a)$
coefficients are selected, while for $\kappa_{n}^{(l)}$ it is one of the
$U_{\alpha}(a)$ or $U_{\beta}(a)$. Though this formulation appears terse, we
reiterate that this represents an algorithm that is capable of readily
producing all information necessary to each perturbative order. For
pedagogical reasons, we summarise the amount of information necessary to
predict typical observables in Tab. 1 [30, 32].

Required order ($n$) | Observable | Number of dynamical coefficients
---|---|---
1 | Linear power spectrum | 2 (trivial)
2 | Tree-level bispectrum | 6
3 | One-loop power spectrum | 18
4 | One-loop bispectrum | 68
5 | Two-loop power spectrum | 290
Table 1: Cumulative information required to calculate typically used
observables of the matter and velocity fluctuations. The number of
coefficients required is calculated as $2\sum_{i=1}^{n}N(i)$.

The $W$ and $U$ functions satisfy coupled differential equations, with the
following structure:

$\displaystyle\dot{W}_{\alpha;m_{1},m_{2}}^{n(i,j)}+nW_{\alpha;m_{1},m_{2}}^{n(i,j)}-U_{\alpha;m_{1},m_{2}}^{n(i,j)}$
$\displaystyle=\kappa_{m_{1}}^{(i)}\lambda_{m_{2}}^{(j)},$
$\displaystyle\dot{U}_{\alpha;m_{1},m_{2}}^{n(i,j)}+(n-1)U_{\alpha;m_{1},m_{2}}^{n(i,j)}-\frac{f_{-}}{f_{+}^{2}}\bigg{[}U_{\alpha;m_{1},m_{2}}^{n(i,j)}-W_{\alpha;m_{1},m_{2}}^{n(i,j)}\bigg{]}$
$\displaystyle=0,$ (14)
$\displaystyle\dot{W}_{\beta;m_{1},m_{2}}^{n(i,j)}+nW_{\beta;m_{1},m_{2}}^{n(i,j)}-U_{\beta;m_{1},m_{2}}^{n(i,j)}$
$\displaystyle=0,$
$\displaystyle\dot{U}_{\beta;m_{1},m_{2}}^{n(i,j)}+(n-1)U_{\beta;m_{1},m_{2}}^{n(i,j)}-\frac{f_{-}}{f_{+}^{2}}\bigg{[}U_{\beta;m_{1},m_{2}}^{n(i,j)}-W_{\beta;m_{1},m_{2}}^{n(i,j)}\bigg{]}$
$\displaystyle=\kappa_{m_{1}}^{(i)}\kappa_{m_{2}}^{(j)},$ (15)

where $n=m_{1}+m_{2}$, and overdots denote differentiation with respect to
$\eta$, where $\eta=\ln D_{+}$ is used as a reduced time coordinate, so that
$d/d\eta=(a/f_{+})d/da$. Equations (14) and (15) have initial conditions as
$a\rightarrow 0$ given by their time-independent EdS values, found using the
recursion relations

$\displaystyle[W_{\alpha;m_{1},m_{2}}^{n(i,j)}]^{\mathrm{EdS}}$
$\displaystyle=\bigg{(}\frac{2n+1}{2n^{2}+n-3}\bigg{)}[\kappa_{m_{1}}^{(i)}]^{\mathrm{EdS}}[\lambda_{m_{2}}^{(j)}]^{\mathrm{EdS}},$
$\displaystyle[U_{\alpha;m_{1},m_{2}}^{n(i,j)}]^{\mathrm{EdS}}$
$\displaystyle=\bigg{(}\frac{3}{2n^{2}+n-3}\bigg{)}[\kappa_{m_{1}}^{(i)}]^{\mathrm{EdS}}[\lambda_{m_{2}}^{(j)}]^{\mathrm{EdS}},$
(16) $\displaystyle[W_{\beta;m_{1},m_{2}}^{n(i,j)}]^{\mathrm{EdS}}$
$\displaystyle=\bigg{(}\frac{2}{2n^{2}+n-3}\bigg{)}[\kappa_{m_{1}}^{(i)}]^{\mathrm{EdS}}[\kappa_{m_{2}}^{(j)}]^{\mathrm{EdS}},$
$\displaystyle[U_{\beta;m_{1},m_{2}}^{n(i,j)}]^{\mathrm{EdS}}$
$\displaystyle=\bigg{(}\frac{2n}{2n^{2}+n-3}\bigg{)}[\kappa_{m_{1}}^{(i)}]^{\mathrm{EdS}}[\kappa_{m_{2}}^{(j)}]^{\mathrm{EdS}}.$
(17)

The efficient solution of the equations of motion (14) and (15) is the main
goal of this paper, as the iterative source terms make their numerical
evaluation relatively computationally intensive. This is accentuated by the
information in Tab. 1, which shows just how many of these functions must be
calculated to predict each observable. As a result, this work aims to find an
efficient method to determine the solution of non-linear, coupled ordinary
differential equations of this form. Doing so will allow us to take into
account efficiently the full time-dependence of the non-linear density and
velocity fields in perturbation-theory calculations.

###  II.2 Direct Numerical Solution for the Equations of Motion

To establish a benchmark, we first solve Eqs. (12) and (13) directly with a
basic LSODA solver for first-order ordinary differential equations (ODEs) [33,
34]. We use initial conditions defined in Eqs. (16) and (17), giving for
instance

$\displaystyle[W_{\alpha;1,1}^{2(1,1)}]^{\mathrm{EdS}}=\frac{5}{7},\;\;\;[U_{\alpha;1,1}^{2(1,1)}]^{\mathrm{EdS}}=\frac{3}{7},\qquad\qquad\;\;\;\;$
$\displaystyle[W_{\beta;1,1}^{2(1,1)}]^{\mathrm{EdS}}=\frac{2}{7},\;\;\;[U_{\beta;1,1}^{2(1,1)}]^{\mathrm{EdS}}=\frac{4}{7},\;\;\;\mathrm{for}\;n=2.$

This process was completed up to third-order, with second-order solutions
being used in turn to generate higher-order results. Figure 2 shows the
solutions for $\lambda_{3}^{(l)}$ and $\kappa_{3}^{(l)}$ for $l=1$–$6$ (since
$N(3)=6$), found using this method. These have been normalised by their
respective EdS values in order to make visual comparison easier. We find that
in the majority of cases these coefficients tend to depart from EdS with
increasing $a$. It is important to note that the iterative nature of Eqs. (14)
and (15) makes this method progressively more time-consuming as we aim to
solve for higher $n$, corresponding to higher-order dynamics. Therefore, a
more efficient method of solution is clearly desirable.

(a) $\lambda_{3}^{(l)}$

(b) $\kappa_{3}^{(l)}$

Figure 2: Numerical solutions for $\lambda_{n}^{(l)}$ (left) and
$\kappa_{n}^{(l)}$ (right) for $n=3$ obtained from Eqs. (12) and (13) with the
ODEs (14) and (15) solved by direct numerical integration. The solutions are
normalised by their respective EdS values. These coefficients start to depart
from their EdS values with increasing $a$, reaching around 1% and 4%
differences for $\lambda_{n}^{(l)}$ and $\kappa_{n}^{(l)}$, respectively.

##  III The Chebyshev Spectral Method

As stated above, we aim to find an efficient method to solving Eqs. (12–15) in
a full $\Lambda$CDM universe. We do so using the Chebyshev spectral method
(CSM), based around an expansion in Chebyshev polynomials. A code
implementation of the CSM in Python is available at
https://github.com/Chousti/CSMethod.git, with its structure presented in Sec.
III.2.

###  III.1 Shifted Chebyshev Polynomials

Before we proceed to describe the CSM, we introduce the Chebyshev polynomials
themselves [35]. These are two families of polynomials defined in relation to
trigonometric functions:

$\displaystyle T_{n}(\cos\theta)$ $\displaystyle=\cos(n\theta),$
$\displaystyle U_{n}(\cos\theta)\sin\theta$ $\displaystyle=\sin(n+1)\theta,$
(18)

typically defined in the domain $\theta\in[-\pi/2,\pi/2]$. These are referred
to as the Chebyshev polynomials of the first and second kind, respectively. We
use the first kind, which can also be generated from the recurrence relation

$T_{0}=1,\qquad T_{1}=x,\qquad T_{n+1}(x)=2xT_{n}(x)-T_{n-1}(x),$ (19)

where $x\in[-1,1]$.

These polynomials are particularly useful as they are orthogonal within their
domain with respect to the weight function $w(x)=(1-x^{2})^{-1/2}$. This is
because they are solutions to the Chebyshev differential equations, which are
of Sturm–Liouville form.

In order to make the Chebyshev polynomials suitable for our problem, we must
rescale them, such that we can use the cosmological scale factor $a\in[0,1]$
as the argument. This is done by defining the shifted Chebyshev polynomials,

$\tilde{T}_{n}(x)=T_{n}\left(\frac{2x-c-b}{c-b}\right),$ (20)

where, for generality, we have considered the domain $[b,c]$. We now proceed
to outline several key properties (for this arbitrary shift) which will be
particularly useful. These include their orthogonality relations:

$\int_{b}^{c}\tilde{T}_{n}(x)\tilde{T}_{m}(x)\frac{dx}{\sqrt{x(b+c-x)-bc}}=\begin{cases}\mbox{0}&\mbox{if
}n\neq m,\\\ \pi&\mbox{if }n=m=0,\\\ \frac{\pi}{2}&\mbox{if }n=m\neq
0;\end{cases}$ (21)

product relations:

$2\tilde{T}_{n}(x)\tilde{T}_{m}(x)=\tilde{T}_{n+m}(x)+\tilde{T}_{|n-m|}(x);$
(22)

and integral relation:

$\frac{4}{c-b}\int\tilde{T}_{n}(x)\,dx=\frac{\tilde{T}_{n+1}(x)}{n+1}-\frac{\tilde{T}_{n-1}(x)}{n-1}\qquad(n\geq
2)\,.$ (23)

These properties were derived based on known results for the original
polynomials. For the full derivations of the original properties, along with a
myriad of others, the reader may consult [28].

It is now possible to approximate an arbitrary, smooth function $y(x)$, valid
in the range $[b,c]$, with a truncated sum of shifted polynomials:

$y(x)=\sum_{i=0}^{N}a_{i}\tilde{T}_{i}(x)=\bm{a}\cdot\bm{t},$ (24)

where the expansion coefficients are given by

$a_{i}=\frac{c_{i}}{\pi}\int_{b}^{c}y(u)\tilde{T}_{i}(u)\frac{du}{\sqrt{u(b+c-u)-bc}},\;\;\;c_{i}=\begin{cases}2&\mbox{if
}i\neq 0,\\\ 1&\mbox{if }i=0.\end{cases}$ (25)

Henceforth, we will refer to $\bm{a}$ as the vector of components of $y(x)$,
of length $N$. In this vectorial language, it is possible to express
operations such as derivatives and products of functions as matrix operations.
Namely, we have the derivative

$\frac{dy}{dx}=\bm{a}^{\prime}\cdot\bm{t}=(\bm{D}\cdot\bm{a})\cdot\bm{t}\;\rightarrow\;\bm{a}^{\prime}=\bm{D}\cdot\bm{a},$
(26)

where we have defined

$\bm{D}=\frac{4}{c-b}\begin{bmatrix}0&\frac{1}{2}&0&\frac{3}{2}&0&\frac{5}{2}&\cdots\\\
0&0&2&0&4&0\\\ 0&0&0&3&0&5\\\ 0&0&0&0&4&0\\\ 0&0&0&0&0&5\\\
\vdots&&&&&&\ddots\end{bmatrix}.$ (27)

This follows from integrating $dy/dx=\bm{a}^{\prime}\cdot\bm{t}$ and using Eq.
(23). Similarly, it is possible to expand the product of two functions using
Eq. (22) as

$y(x)=g(x)h(x)=(\bm{g}\cdot\bm{t})(\bm{h}\cdot\bm{t})=(\bm{P}(\bm{g})\cdot\bm{h})\cdot\bm{t}\;\rightarrow\;\bm{a}=\bm{P}(\bm{g})\cdot\bm{h},$
(28)

where we have defined

$\bm{P}(\bm{g})=\frac{1}{2}\begin{bmatrix}2g_{0}&g_{1}&g_{2}&g_{3}&\cdots\\\
2g_{1}&2g_{0}+g_{2}&g_{1}+g_{3}&g_{2}+g_{4}&\\\
2g_{2}&g_{1}+g_{3}&2g_{0}+g_{4}&g_{1}+g_{5}&\\\
2g_{3}&g_{2}+g_{4}&g_{1}+g_{5}&2g_{0}+g_{6}&\\\
\vdots&&&&\ddots\end{bmatrix}.$ (29)

Here, we note that two different $\bm{P}$ matrices commute between each other
(i.e., $\bm{P}(\bm{g})\cdot\bm{h}=\bm{P}(\bm{h})\cdot\bm{g}$). Furthermore,
$\bm{P}$ and $\bm{D}$ matrices will not commute, providing a useful sanity
check. This formalism will prove particularly useful and will be implemented
fully in Section III.2. Though these results are given for an arbitrary
rescaling of the argument, henceforth we shall proceed with $\tilde{T}$
denoting a shift to the domain $[0,1]$.

Given the ODEs (14) and (15) that we need to solve, when expressed in terms of
derivatives with respect to $a$ we shall require decompositions of the
functions $f_{-}/f_{+}^{2}$ and $1/f_{+}$ into the polynomial basis. To allow
for future generalisability, this is done numerically using Eqs. (24) and (25)
with the accuracy of the recomposition for $f_{-}/f_{+}^{2}$ shown in Fig. 3
for $N\in[2,4,6]$ Chebyshev components. Here we find that using $N=4$ yields
sub-percent-level accuracy across the whole function. The same is true for the
decomposition of $1/f_{+}$, but this is not shown for brevity.

Figure 3: Numerical decomposition of $\delta f=f_{-}/f_{+}^{2}+3/2$ into the
shifted Chebyshev polynomial basis. Both the residue (top) and absolute values
of the relative error (bottom) are shown. The residue is given by $\Delta
X=X-X_{\mathrm{true}}$.

###  III.2 Implementation of CSM

We now discuss the implementation and structure of the CSM. Generally,
spectral methods work by numerically determining the coefficients in a given
basis of a differential equation subject to boundary conditions [36].
Particularly, the CSM makes use of Chebyshev polynomials because they are easy
to compute and rapidly convergent as compared to other basis functions such as
the Legendre Polynomials [36, 37, 29].

On the surface, this method works by using the polynomials’ properties to
convert differential equations into matrix equations, corresponding to a
system of linear equations for a set of unknown components, with the number of
components ($N$) corresponding to the user’s desired accuracy. Therefore, the
problem’s complexity has been reduced to one of simple linear algebra, for
which a variety of optimised codes exist. As a result, the method will naively
always be faster than one using numerical integration routines.

We proceed by defining the following expansions, using the formalism outlined
in Sec. III.1 of shifted Chebyshev polynomials:

$\displaystyle\lambda(a)$
$\displaystyle=\bm{L}\cdot\bm{t}(a),\;\;\kappa(a)=\bm{K}\cdot\bm{t}(a),$ (30)
$\displaystyle W(a)$
$\displaystyle=\bm{x}\cdot\bm{t}(a),\;\;U(a)=\bm{y}\cdot\bm{t}(a),$ (31)
$\displaystyle\frac{f_{-}}{f_{+}^{2}}(\Omega_{m_{0}},a)$
$\displaystyle=\bm{d}(\Omega_{m_{0}})\cdot\bm{t}(a),\;\;f_{+}^{-1}(\Omega_{m_{0}},a)=\bm{c}(\Omega_{m_{0}})\cdot\bm{t}(a),\;\;a=\bm{e}\cdot\bm{t}(a),$
(32)

where $\bm{e}=(1,0,0,\ldots)$ and all subscripts and superscripts have been
implied but omitted for brevity. The components $c$ and $d$ are calculated as
shown in Fig. 3 on the first iteration of each run.

Next, we substitute these relations into Eqs. (12) and (13) for a given
cosmology and utilise properties of the polynomials to find:

$\displaystyle\bm{L}_{n}^{(l)}=\delta^{K}_{\frac{n}{2},\lfloor\frac{n}{2}\rfloor}\sum_{i=1}^{N(n/2)}\bigg{[}\sum_{j=1}^{N(n/2)}[\bm{x}_{\alpha}]^{n(i,j)}_{\frac{n}{2},\frac{n}{2}}\delta^{K}_{l,\phi_{1}}+\sum_{j=i}^{N(n/2)}[\bm{x}_{\beta}]^{n(i,j)}_{\frac{n}{2},\frac{n}{2}}\delta^{K}_{l,\phi_{2}}\bigg{]}\qquad\qquad\qquad\qquad\qquad$
$\displaystyle+\sum_{m=1}^{\lfloor(n-1)/2\rfloor}\sum_{i=1}^{N(m)}\sum_{j=1}^{N(n-m)}\bigg{[}[\bm{x}_{\alpha}]^{n(i,j)}_{m,n-m}\delta^{K}_{l,\phi_{3}}+[\bm{x}_{\alpha}]^{(ji)}_{n-m,m}\delta^{K}_{l,\phi_{4}}+[\bm{x}_{\beta}]^{n(i,j)}_{m,n-m}\delta^{K}_{l,\phi_{5}}\bigg{]},$
(33)
$\displaystyle\bm{K}_{n}^{(l)}=\delta^{K}_{\frac{n}{2},\lfloor\frac{n}{2}\rfloor}\sum_{i=1}^{N(n/2)}\bigg{[}\sum_{j=1}^{N(n/2)}[\bm{y}_{\alpha}]^{n(i,j)}_{\frac{n}{2},\frac{n}{2}}\delta^{K}_{l,\phi_{1}}+\sum_{j=i}^{N(n/2)}[\bm{y}_{\beta}]^{n(i,j)}_{\frac{n}{2},\frac{n}{2}}\delta^{K}_{l,\phi_{2}}\bigg{]}\qquad\qquad\qquad\qquad\qquad$
$\displaystyle+\sum_{m=1}^{\lfloor(n-1)/2\rfloor}\sum_{i=1}^{N(m)}\sum_{j=1}^{N(n-m)}\bigg{[}[\bm{y}_{\alpha}]^{n(i,j)}_{m,n-m}\delta^{K}_{l,\phi_{3}}+[\bm{y}_{\alpha}]^{(ji)}_{n-m,m}\delta^{K}_{l,\phi_{4}}+[\bm{y}_{\beta}]^{n(i,j)}_{m,n-m}\delta^{K}_{l,\phi_{5}}\bigg{]},$
(34)

where the coupled system of ODEs (14) and (15) becomes a system of algebraic
equations:

$\displaystyle\bigg{[}\bm{P}(\bm{c})\cdot\bm{P}(\bm{e})\cdot\bm{D}+n\bm{I}\bigg{]}\cdot[\bm{x}_{\alpha}]_{m_{1},m_{2}}^{n(i,j)}-\bm{I}\cdot[\bm{y}_{\alpha}]_{m_{1},m_{2}}^{n(i,j)}$
$\displaystyle=\bm{P}(\bm{K}_{m_{1}}^{(i)})\cdot\bm{L}_{m_{2}}^{(j)},$
$\displaystyle\bm{P}(\bm{d})\cdot[\bm{x}_{\alpha}]_{m_{1},m_{2}}^{n(i,j)}+\bigg{[}\bm{P}(\bm{c})\cdot\bm{P}(\bm{e})\cdot\bm{D}+(n-1)\bm{I}-\bm{P}(\bm{d})\bigg{]}\cdot[\bm{y}_{\alpha}]_{m_{1},m_{2}}^{n(i,j)}$
$\displaystyle=0,$ (35)

and

$\displaystyle\bigg{[}\bm{P}(\bm{c})\cdot\bm{P}(\bm{e})\cdot\bm{D}+n\bm{I}\bigg{]}\cdot[\bm{x}_{\beta}]_{m_{1},m_{2}}^{n(i,j)}-\bm{I}\cdot[\bm{y}_{\beta}]_{m_{1},m_{2}}^{n(i,j)}$
$\displaystyle=0,$
$\displaystyle\bm{P}(\bm{d})\cdot[\bm{x}_{\beta}]_{m_{1},m_{2}}^{n(i,j)}+\bigg{[}\bm{P}(\bm{c})\cdot\bm{P}(\bm{e})\cdot\bm{D}+(n-1)\bm{I}-\bm{P}(\bm{d})\bigg{]}\cdot[\bm{y}_{\beta}]_{m_{1},m_{2}}^{n(i,j)}$
$\displaystyle=\bm{P}(\bm{K}_{m_{1}}^{(i)})\cdot\bm{K}_{m_{2}}^{(j)}$ (36)

for the $\alpha$ and $\beta$ systems, respectively. It is important to note
several aspects. First, the length ($N+1$) of the unknown vectors $\bm{x}$ and
$\bm{y}$ corresponds to the greatest order of Chebyshev polynomial used ($N$)
and is set by the user, with the implications discussed in Section IV.1.
Secondly, due to the iterative nature of these equations, in principle, the
source terms on the right-hand sides of Eqs. (35) and (36) are just vectors of
constants, which we shall henceforth define as $\bm{\sigma}_{\alpha/\beta}$
and $\bm{\tau}_{\alpha/\beta}$, respectively. We note that
$\bm{\tau}_{\alpha}=\bm{0}$ and $\bm{\sigma}_{\beta}=\bm{0}$. These source
terms are calculated recursively, making use of the code’s inherent structure
to ensure this is done efficiently. Finally, one small approximation is made
to make the code even more efficient. In particular, it was noted that for all
functions decomposed in this way, the components $a_{i}$ were
$\mathcal{O}(10^{-i})$. Therefore, in Eq. (24), the code forces all components
of orders greater than $N$ to zero.

Schematically, we can simplify Eqs. (35) and (36) into matrix equations of
dimension $2N+2$:

$\displaystyle\left(\begin{array}[]{@{}c|c@{}}\bm{E}&\bm{F}\\\
\hline\cr\bm{G}&\bm{H}\\\
\end{array}\right)\left(\begin{array}[]{@{}c@{}}\bm{x}_{\mu;m_{1},m_{2}}^{n(i,j)}\\\
\bm{y}_{\mu;m_{1},m_{2}}^{n(i,j)}\end{array}\right)=\left(\begin{array}[]{@{}c@{}}\bm{\sigma}_{\mu;m_{1},m_{2}}^{n(i,j)}\\\
\bm{\tau}_{\mu;m_{1},m_{2}}^{n(i,j)}\end{array}\right),$ (43)

where $\mu=\alpha,\beta$ denotes which system is being computed. We must also
implement our boundary conditions, namely that our dynamical coefficients
reduce to their EdS counterparts as $a\rightarrow 0$. Component-wise, this
becomes for $W$:

$\displaystyle[W_{\mu;m_{1},m_{2}}^{n(i,j)}]^{\mathrm{EdS}}=\sum_{k=0}^{N}[x_{\mu;m_{1},m_{2}}^{n(i,j)}]_{k}\tilde{T}_{k}(0)=\sum_{k=0}^{N}(-1)^{k}[x_{\mu;m_{1},m_{2}}^{n(i,j)}]_{k},$
(44)

with a similar result for $U$. We then force this constraint on the system by
replacing the bottom row of $\bm{E}$ with $[1,-1,...,(-1)^{N+1}]$, the bottom
row of $\bm{F}$ with an ($N+1$)-tuple of zeros, and finally setting
$[\sigma_{\mu;m_{1},m_{2}}^{n(i,j)}]_{N}=[W_{\mu;m_{1},m_{2}}^{n(i,j)}]^{\mathrm{EdS}}$.
We then repeat this process for $\bm{H}$, $\bm{G}$ and
$\bm{\tau}_{\mu;m_{1},m_{2}}^{n(i,j)}$, respectively.

Therefore, we have successfully reduced our system of coupled ODEs to a basic
matrix multiplication problem, for which there exist a variety of efficient
linear algebra methods. For the purposes of the CSM, a well-established LAPACK
routine was used [38] to invert Eq. (43). Once these unknown components have
been found, we can recompose them to find $W$ and $U$. In turn, Eqs. (33) and
(34) can then be used to find the dynamical solutions for Eqs. (12) and (13),
therefore solving our problem.

##  IV Results and comparison of solutions

###  IV.1 Testing the Chebyshev Spectral Method

The method described in Section III has been implemented, and is here tested
by means of solving Eqs. (12) and (13) for the case of a $\Lambda$CDM universe
with initial conditions given by the EdS limit (Eqs. 16 and 17) and with
$\Omega_{m_{0}}=0.315$. This case was chosen for illustrative purposes only,
as the CSM is capable of efficiently producing results for any $\Lambda$CDM
universe, allowing it to be used to scan across the full parameter space of
$\Omega_{m_{0}}$.

Figure 4 shows all solutions to Eq. (12) for the $\lambda_{n}^{(l)}$ (red) and
Eq. (13) for the $\kappa^{(l)}_{n}$ (blue) for $n\in[2,3,4,5]$ and for all
associated values of $l$ as generated by the CSM with $N=5$. This demonstrates
both the increasing number of dynamical functions and their greater deviation
from the EdS limit with increased perturbative order. These functions
represent all those necessary in order to derive the matter-density and
associated velocity power spectra to two-loop-order, as given in [26].

(a) $n=2$

(b) $n=3$

(c) $n=4$

(d) $n=5$

Figure 4: Normalised solutions to Eqs. (12) and (13) up to fifth perturbative
order in a $\Lambda$CDM universe, calculated using the Chebyshev spectral
method with $N=6$. The functions are shown divided by their EdS counterparts
(which are constant in time).

Figure 5 shows a direct comparison between results for $\lambda_{3}^{(1)}$ as
calculated by direct numerical integration of the differential equations and
by the CSM with $N\in[2,4,6]$ components. Results obtained using the EdS
approximation are also shown for illustrative purposes. Both the residues are
shown (top) as well as the absolute value of the relative error (bottom),
displaying several clear facts. First, the CSM successfully manages to
reproduce the dynamics of this third-order coefficient, with an expansion
truncated at $N=2$ already achieving an accuracy of greater than 0.03%.
Furthermore, we find the expected result that increasing $N$ yields more
accurate results. This suggests that the user is effectively able to control
the output accuracy of results, though with slight sacrifices in efficiency
(as discussed in Section IV.1). Finally, we also find that the relative
accuracy of the CSM tends to improve at late times ($a\rightarrow 1$) away
from zero crossings for all component numbers and at all orders.

Figure 5: Comparison of approximate results for $\lambda_{3}^{(1)}$ from the
Chebyshev spectral method and direct numerical integration. Differences with
respect to the direct solution are shown in the top panel keeping different
numbers of terms in the expansions in Chebyshev polynomials, while the
(absolute values) of the fractional differences are shown in the bottom panel.
The Chebyshev spectral method is very accurate for all $a$ even with few terms
in the expansion. On the other hand, the EdS approximation (dotted lines)
diverges rapidly away from the true solution.

Next, Fig. 6 shows the normalised dynamical dependence of
$\lambda_{5}^{(1)}(a)$ as a one-parameter-family of $\Omega_{m_{0}}$,
calculated by the CSM with $N=5$. Here, examples of dark-energy dominated
(red), $\Lambda$CDM (black) and matter-dominated (blue) cosmologies are shown,
along with the shaded full parameter space. Here, it is found that
$\lambda_{5}^{(1)}$ shows a non-trivial $\Omega_{m_{0}}$ dependence, further
implying the importance of fully taking these dynamics into account. Finally,
the CSM successfully reproduces the expected EdS result of constant
coefficients as $\Omega_{m_{0}}\rightarrow 1$.

Figure 6: Dynamical results for $\lambda_{5}^{(1)}$ calculated using the CSM
as a one-parameter-family of $\Omega_{m_{0}}$. Specifically, matter-dominated
(blue), $\Lambda$CDM (black) and dark-energy-dominated (red) results are
highlighted, with matter density parameters labelled. The boundary of the
shaded region is for $\Omega_{m_{0}}\rightarrow 0$.

The final comparison we make is to an alternative solution to the problem
considered in this paper, which we shall introduce here following results in
[26]. The key observation there is that the parametric dependence of the
equations of motion (14) and (15) on $\Omega_{m_{0}}$ can be absorbed into a
new independent variable $q$, where

$q\equiv\left(\frac{1-\Omega_{m_{0}}}{\Omega_{m_{0}}}\right)a^{3}.$ (45)

This follows by noting that $D_{+}(a;\Omega_{m_{0}})$ may be written as

$D_{+}(a;\Omega_{m_{0}})=\left(\frac{\Omega_{m_{0}}}{1-\Omega_{m_{0}}}\right)^{1/3}\hat{D}_{+}(q),$
(46)

where the rescaled growth function $\hat{D}_{+}$ depends only on $q$. It
follows that $f_{\pm}(a;\Omega_{m_{0}})$ also depend only on $q$ and, since
$d\ln D_{+}=d\ln\hat{D}_{+}$, so do the solutions $U_{\alpha/\beta}$ and
$W_{\alpha/\beta}$ of Eqs (14) and (15). Expanding $f_{\pm}(q)$ as power
series in $q$, or, equivalently, in $\zeta\equiv\hat{D}_{+}^{3}$ on noting
that around $\Omega_{m_{0}}=1$ (i.e., $q=0$)
$\hat{D}_{+}(q)=q^{1/3}\left[1+\mathcal{O}(q)\right]$, we have

$\frac{f_{-}}{f_{+}^{2}}(a;\;\Omega_{m_{0}})=-\frac{3}{2}+\sum_{i=1}^{\infty}c_{i}\zeta^{i}.$
(47)

Truncating this expansion at low order works well, with [26] suggesting
dropping $c_{4}$ and higher. Expanding the $\lambda_{n}^{(l)}$ and
$\kappa_{n}^{(l)}$ similarly,

$\displaystyle\lambda_{n}^{(l)}(a;\Omega_{m_{0}})\approx\sum_{i=0}^{3}[g_{n}^{l}]_{i}c_{i}\zeta^{i},\qquad\kappa_{n}^{(l)}(a;\Omega_{m_{0}})\approx\sum_{i=0}^{3}[h_{n}^{l}]_{i}c_{i}\zeta^{i},$
(48)

the coefficients may be determined analytically from the equations of motion.

Figure 7 compares this truncated power-series expansion in $\zeta$ and the
CSM, providing plots of the relative errors for all $\lambda_{n}^{(l)}$ and
$\kappa_{n}^{(l)}$ for perturbative orders $n\in[2,3,4]$ compared to the
direct numerical solutions. Though all $\zeta$-expansion (blue) and CSM with
$N=4$ (red) and $N=6$ (orange) curves are given, the averages have been
extracted and shown in bold for clarity. Here, it is found that although both
methods are very successful at reproducing all solutions, we find that the CSM
consistently produces greater accuracy as $a\rightarrow 1$, making it more
useful for computation of late-time observations. The CSM is also helped by
greater efficiency in calculating components and significantly greater
tunability in both accuracy and parameter space, as discussed later in Sec. V.

Figure 7: Comparisons between the relative errors in solutions found for the
dynamical coefficients $\lambda_{n}^{(l)}$ and $\kappa_{n}^{(l)}$ for the
perturbative orders $n=2,3,4$ in the $\Lambda$CDM universe. Results are shown
for the truncated power-series expansion in $\zeta$ given by Eq. (48) (blue;
labelled Perturbative) and the CSM with $N=4$ (red) and $N=6$ (orange). Curves
are shown for all $l$ modes, with averages given in bold.

Finally, Fig. 8 shows the accuracy of the CSM as applied to calculating the
present-day one- (left) and two-loop (right) matter–matter (black),
matter–velocity (red) and velocity–velocity (blue) power spectra compared to
numerically evaluated spectra. These were produced by combining dynamical
coefficients produced by the CSM with momentum kernels calculated as described
in [26], omitted here for brevity. Here, we show both the loop-order
contributions to the $\Lambda$CDM power spectrum (top) and absolute-value
relative error plots (${\rm Er}(x)=\log_{10}(|\Delta x|/x)$; bottom) for the
overall $\Lambda$CDM power spectrum up to each loop order using $N=2$
(dotted), $N=4$ (dot-dashed) and $N=6$ (solid) Chebyshev components. Clearly,
in the case of the one-loop power spectrum, using $N=2$ components is
sufficient to produce results to sub-percent-level accuracy. In the case of
the two-loop power spectrum contributions, it was found that such gains were
dwarfed by numerical uncertainties produced by the numerical integration of
the momentum kernels. In principle, however, this demonstrates the
effectiveness of the CSM in predicting actual observables.

Figure 8: Present-day one-loop (left) and two-loop (right) contributions to
the matter–matter (black), matter–velocity (red) and velocity–velocity (blue)
power spectra, computed using the CSM with $N=2$ (dotted), $N=4$ (dot-dashed)
and $N=6$ (solid) Chebyshev polynomials. We show the loop-order contributions
to the $\Lambda$CDM power spectrum (top) and relative errors in the
$\Lambda$CDM power spectra as compared to the numerical values (${\rm
Er}(x)=\log_{10}(|\Delta x|/x)$; bottom) for each case. Errors in the two-loop
contributions are dominated by numerical uncertainties in calculations of the
momentum kernels.

###  IV.2 Benchmarking the Code

In Fig. 9 we compare the average computational time required to solve Eqs.
(12) and (13) with direct numerical integration (as in Sec. II.2) and our
implementation of the CSM. Specifically, we show the time taken to calculate
all dynamical components for a given perturbative order $n$ for the Planck
best-fit $\Lambda$CDM universe. For the direct method (red) we use 50
subdivisions of the range $a\in[0,1]$, while for the CSM (blue) we use a
variety of components in the Chebyshev expansion ($N\in[3,4,10]$). Of these,
$N=4$ would be sufficiently accurate for most uses (and is thus discussed
below), while $N=10$ is far more accurate than necessary for the time taken.
We find that the terms needed for the one-loop power spectrum ($n=2$ and $3$)
are calculated two orders of magnitude faster than by direct numerical
integration; those for the two-loop power spectrum ($n=2,3,4,5$) are
calculated more than four orders of magnitude faster. We also comment that the
CSM code used here is inherently iterative, therefore returning the results
for all previous-order coefficients while calculating the dynamics of a target
$n$. This further exemplifies the efficiency of this method, particularly
within the context of its flexibility. Finally, due to this iterative nature,
our implementation of the CSM is able to calculate the dynamical coefficients
to any perturbative order. This implies that this code will remain useful as
momentum operators at increasingly greater orders are computed in the future.

Figure 9: Comparison of the time taken to compute all $\lambda_{n}^{(l)}$ and
$\kappa_{n}^{(l)}$ for varying perturbative order $n$ in the Planck best-fit
$\Lambda$CDM universe using both direct numerical integration (red) and
Chebyshev spectral methods (blue) with varying degrees of accuracy. Direct
numerical solutions for higher $n$ are not included as their computation is
prohibitively slow.

##  V Discussion and Conclusions

The results in Section IV, including, most notably, Figs. 5, 7 and 8
demonstrate the accuracy and effectiveness of the Chebyshev spectral method
for solving the full dynamical evolution of the dark matter density and
velocity fluctuation fields in a $\Lambda$CDM universe. This is coupled with
the fact that the CSM is, in practice, at least an order of magnitude faster
than direct numerical ODE solution methods. The true utility of this method is
felt, however, at higher perturbative orders $n$, with Fig. 9 demonstrating
how much faster the CSM is at such orders.

From a practical standpoint, the code implementation produced as part of this
work has the following features. First, the method can be used with any number
of Chebyshev components $N$. Doing so increases the accuracy of the resultant
solution somewhat but incurs an extra computational time cost, as seen in Fig.
9. As such, the user is able to make this decision actively and tailor the
method to their particular situation and use case. Next, it is found that the
magnitude of the calculated components drops almost exponentially with $N$,
further showing the futility of finding these high-order components. As a
result, in testing it was found that using $N=4$ components is more than
suitable, with $N=2$ yielding sufficient results in most cases, as can be seen
explicitly in Figs. 5 and 8.

Furthermore, the code has been generalised to work for arbitrary
$\Omega_{m_{0}}$, extending the parameter space of this method to encompass
$\Lambda$CDM. The effect of this has been explored in Fig. 6. Along these
lines, we have also compared the effectiveness of the CSM to an alternative
perturbative method introduced in [26], which is also valid for $\Lambda$CDM
only. It was found that the CSM was able to achieve a similar or better degree
of accuracy, particularly for late times. The CSM also grants the user a
choice between accuracy and efficiency, which is important when calculating
model predictions for cosmological likelihood analyses. Furthermore, while the
perturbative method is on the surface faster to evaluate for a given
$\Omega_{m_{0}}$, the CSM allows the user to generate dynamical-coefficient
data to any perturbative order efficiently. Next, the framework discussed in
Section II can be extended to a more general set of cosmologies, including
clustered quintessence [30, 31]. On the other hand, it is thought that the
perturbative method outlined above in Eq. (48) will only work in the case of
$\Lambda$CDM. Indeed, such an extension of the CSM has already been done to
first-loop-order for the power spectrum, with work underway to expand this
further.111This refers to work completed in a masters thesis by one of the
authors, along with an upcoming paper.

Another interesting direction to explore would be to consider the implications
of other models of dark matter and energy. These might include decaying dark
matter [39] or a dynamical dark energy equation of state [40, 41]. Finally, it
would be particularly interesting to explore the impact of various dark energy
effective field theories (EFTs), of which clustered quintessence is an example
[42, 43, 44, 45, 30, 31, 46]. However, the most natural extension would be to
implement similar methods to study other observables to the same loop order,
such as the bispectrum.

To conclude, we have shown how the CSM can optimize the problem of solving the
full, perturbative dynamics of perturbations in the $\Lambda$CDM Universe,
reducing computation time by several orders of magnitude, with a view of being
potentially extended to a more general set of cosmologies. Our implementation
of the CSM as a Python library is presented alongside this paper
(https://github.com/Chousti/CSMethod.git). It should be straightforward to
integrate this library into galaxy-clustering likelihood frameworks, allowing
accurate predictions of observables in the perturbative regime fully
accounting for the $\Lambda$CDM dynamics.

###### Acknowledgements.

We thank Tomohiro Fujita and Matteo Fasiello for useful discussions. N.C. and
Z.V. acknowledge the support of the Kavli Foundation. N.C. acknowledges
support from the Science and Technology Facilities Council (STFC) for a Ph.D.
studentship. A.C. acknowledges support from the STFC (grant numbers
ST/N000927/1 and ST/S000623/1).

##  Appendix A Auxiliary Functions

Here we review all auxiliary functions necessary to calculate Eqs. (12) and
(13). These are derived in [26] and represent a recursive method to establish
the maximum number of dynamical coefficients. In principle, this number can be
reduced by utilising physical constraints such as conservation of mass and
momentum or the equivalence principle — though these are not discussed here.
Specifically, we have the numbering function,

$N(n)=\frac{1}{2}\delta^{K}_{\frac{n}{2},\lfloor\frac{n}{2}\rfloor}N(n/2)[1+3N(n/2)]+3\sum_{m=1}^{\lfloor(n-1)/2\rfloor}N(m)N(n-m),$
(49)

and bijective maps,

$\displaystyle\phi_{1}(n,i,j)$ $\displaystyle=N(n/2)(i-j)+i,$ (50)
$\displaystyle\phi_{2}(n,i,j)$
$\displaystyle=[N(n/2)]^{2}-\frac{1}{2}i(i-1)+\phi_{1}(n,i,j),$ (51)
$\displaystyle\phi_{3}(n,m,i,j)$
$\displaystyle=\frac{1}{2}\delta^{K}_{\frac{n}{2},\lfloor\frac{n}{2}\rfloor}N(n/2)[1+3N(n/2)]+\sum_{k=1}^{m-1}N(k)N(n-k)+(i-1)N(n-m)+j,$
(52) $\displaystyle\phi_{4}(n,m,i,j)$
$\displaystyle=\sum_{k=1}^{\lfloor(n-1)/2\rfloor}N(k)N(n-k)+\phi_{3}(n,m,i,j),$
(53) $\displaystyle\phi_{5}(n,m,i,j)$
$\displaystyle=2\sum_{k=1}^{\lfloor(n-1)/2\rfloor}N(k)N(n-k)+\phi_{3}(n,m,i,j).$
(54)

The CSM makes use of these by looping over all maps while evaluating Eqs. (14)
and (15), hence algorithmically calculating all necessary functions.

## References

  * Desjacques  _et al._ [2018] V. Desjacques, D. Jeong, and F. Schmidt, Large-Scale Galaxy Bias, Phys. Rept. 733, 1 (2018), arXiv:1611.09787 [astro-ph.CO] .
  * Kaiser [1992] N. Kaiser, Weak gravitational lensing of distant galaxies, Astrophys. J. 388, 272 (1992).
  * Hamilton [1997] A. J. S. Hamilton, Linear redshift distortions: A Review, in  _Ringberg Workshop on Large Scale Structure_ (1997) arXiv:astro-ph/9708102 .
  * Amendola  _et al._ [2018] L. Amendola _et al._ , Cosmology and fundamental physics with the Euclid satellite, Living Rev. Rel. 21, 2 (2018), arXiv:1606.00180 [astro-ph.CO] .
  * DESI Collaboration [2016] DESI Collaboration, The DESI Experiment Part I: Science,Targeting, and Survey Design, arXiv e-prints , arXiv:1611.00036 (2016), arXiv:1611.00036 [astro-ph.IM] .
  * LSST Dark Energy Science Collaboration [2012] LSST Dark Energy Science Collaboration, Large Synoptic Survey Telescope: Dark Energy Science Collaboration, arXiv e-prints , arXiv:1211.0310 (2012), arXiv:1211.0310 [astro-ph.CO] .
  * Spergel  _et al._ [2015] D. Spergel _et al._ , Wide-Field InfrarRed Survey Telescope-Astrophysics Focused Telescope Assets WFIRST-AFTA 2015 Report, arXiv e-prints , arXiv:1503.03757 (2015), arXiv:1503.03757 [astro-ph.IM] .
  * Peebles and Yu [1970] P. J. E. Peebles and J. T. Yu, Primeval adiabatic perturbation in an expanding universe, Astrophys. J. 162, 815 (1970).
  * Bond and Efstathiou [1984] J. R. Bond and G. Efstathiou, Cosmic background radiation anisotropies in universes dominated by nonbaryonic dark matter, Astrophys. J. Lett. 285, L45 (1984).
  * Aghanim  _et al._ [2020] N. Aghanim _et al._ (Planck), Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Chisari  _et al._ [2018] N. E. Chisari, M. L. A. Richardson, J. Devriendt, Y. Dubois, A. Schneider, A. L. Brun, M. C., R. S. Beckmann, S. Peirani, A. Slyz, and C. Pichon, The impact of baryons on the matter power spectrum from the Horizon-AGN cosmological hydrodynamical simulation, Mon. Not. Roy. Astron. Soc. 480, 3962 (2018), arXiv:1801.08559 [astro-ph.CO] .
  * Bernardeau  _et al._ [2002] F. Bernardeau, S. Colombi, E. Gaztanaga, and R. Scoccimarro, Large scale structure of the universe and cosmological perturbation theory, Phys. Rept. 367, 1 (2002), arXiv:astro-ph/0112551 .
  * Fasiello and Vlah [2016a] M. Fasiello and Z. Vlah, Nonlinear fields in generalized cosmologies, Phys. Rev. D 94, 063516 (2016a), arXiv:1604.04612 [astro-ph.CO] .
  * D’Amico  _et al._ [2020] G. D’Amico, J. Gleyzes, N. Kokron, K. Markovic, L. Senatore, P. Zhang, F. Beutler, and H. Gil-Marín, The Cosmological Analysis of the SDSS/BOSS data from the Effective Field Theory of Large-Scale Structure, JCAP 05, 005, arXiv:1909.05271 [astro-ph.CO] .
  * Ivanov  _et al._ [2020] M. M. Ivanov, M. Simonović, and M. Zaldarriaga, Cosmological Parameters from the BOSS Galaxy Power Spectrum, JCAP 05, 042, arXiv:1909.05277 [astro-ph.CO] .
  * Chen  _et al._ [2022] S.-F. Chen, Z. Vlah, and M. White, A new analysis of galaxy 2-point functions in the BOSS survey, including full-shape information and post-reconstruction BAO, JCAP 02 (02), 008, arXiv:2110.05530 [astro-ph.CO] .
  * Bernardeau [1994] F. Bernardeau, Skewness and Kurtosis in large scale cosmic fields, Astrophys. J. 433, 1 (1994), arXiv:astro-ph/9312026 .
  * Takahashi [2008] R. Takahashi, Third-order density perturbation and one-loop power spectrum in dark-energy-dominated universe, Progress of Theoretical Physics 120, 549 (2008).
  * Fasiello and Vlah [2016b] M. Fasiello and Z. Vlah, Nonlinear fields in generalized cosmologies, Physical Review D 94, 10.1103/physrevd.94.063516 (2016b).
  * Lewandowski  _et al._ [2017] M. Lewandowski, A. Maleknejad, and L. Senatore, An effective description of dark matter and dark energy in the mildly non-linear regime, JCAP 05, 038, arXiv:1611.07966 [astro-ph.CO] .
  * Fujita and Vlah [2020] T. Fujita and Z. Vlah, Perturbative description of biased tracers using consistency relations of LSS, JCAP 10, 059, arXiv:2003.10114 [astro-ph.CO] .
  * Donath and Senatore [2020] Y. Donath and L. Senatore, Biased Tracers in Redshift Space in the EFTofLSS with exact time dependence, JCAP 10, 039, arXiv:2005.04805 [astro-ph.CO] .
  * Schmidt [2021] F. Schmidt, An $n$-th order Lagrangian Forward Model for Large-Scale Structure, JCAP 04, 033, arXiv:2012.09837 [astro-ph.CO] .
  * Rampf  _et al._ [2022] C. Rampf, S. O. Schobesberger, and O. Hahn, Analytical growth functions for cosmic structures in a $\Lambda$CDM Universe, Mon. Not. Roy. Astron. Soc. 516, 2840 (2022), arXiv:2205.11347 [astro-ph.CO] .
  * Garny and Taule [2021] M. Garny and P. Taule, Loop corrections to the power spectrum for massive neutrino cosmologies with full time- and scale-dependence, JCAP 01, 020, arXiv:2008.00013 [astro-ph.CO] .
  * Fasiello  _et al._ [2022] M. Fasiello, T. Fujita, and Z. Vlah, Perturbation theory of large scale structure in the $\Lambda$ CDM Universe: Exact time evolution and the two-loop power spectrum, Phys. Rev. D 106, 123504 (2022), arXiv:2205.10026 [astro-ph.CO] .
  * Garny and Taule [2022] M. Garny and P. Taule, Two-loop power spectrum with full time- and scale-dependence and EFT corrections: impact of massive neutrinos and going beyond EdS, JCAP 09, 054, arXiv:2205.11533 [astro-ph.CO] .
  * Karjanto [2020] N. Karjanto, Properties of chebyshev polynomials, arXiv: History and Overview (2020).
  * Sezer and Kaynak [1996] M. Sezer and M. Kaynak, Chebyshev polynomial solutions of linear differential equations, International Journal of Mathematical Education in Science and Technology 27, 607 (1996), https://doi.org/10.1080/0020739960270414 .
  * Sefusatti and Vernizzi [2011] E. Sefusatti and F. Vernizzi, Cosmological structure formation with clustering quintessence, JCAP 03, 047, arXiv:1101.1026 [astro-ph.CO] .
  * Fasiello and Vlah [2019] M. Fasiello and Z. Vlah, On observables in a dark matter-clustering quintessence system, Journal of Cosmology and Astroparticle Physics 2019 (05), 033.
  * Bose and Taruya [2018] B. Bose and A. Taruya, The one-loop matter bispectrum as a probe of gravity and dark energy, JCAP 10, 019, arXiv:1808.01120 [astro-ph.CO] .
  * Hindermarsh [1983] A. C. Hindermarsh, Odepack, a systemized collection of ode solvers, Scientific Computing, R. S. Stepleman et al. 1 of IMACS transcations on Scientific Computation, 55 (1983).
  * Petzold [1983] L. Petzold, Automatic selection of methods for solving stiff and nonstiff systems of ordinary differential equations, SIAM Journal on Scientific and Statistical Computing 4 (1983).
  * Chebyshev [1854] P. L. Chebyshev, Theorie des mecanismes connus sous le nom de parallelogrammes (1854).
  * Boyd [2000] J. P. Boyd,  _Chebyshev and Fourier Spectral Methods_ (DOVER Publications, Inc., 2000).
  * Clenshaw [1957] C. W. Clenshaw, The numerical solution of linear differential equations in chebyshev series, Mathematical Proceedings of the Cambridge Philosophical Society 53, 134–149 (1957).
  * Anderson  _et al._ [1999] E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. D. Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen,  _LAPACK Users’ Guide_ (Society for Industrial and Applied, 1999).
  * Bell  _et al._ [2011] N. F. Bell, A. J. Galea, and R. R. Volkas, Model for late dark matter decay, Physical Review. D, Particles Fields 83, 10.1103/PHYSREVD.83.063504 (2011).
  * Chevallier and Polarski [2001] M. Chevallier and D. Polarski, Accelerating universes with scaling dark matter, International Journal of Modern Physics D 10, 213 (2001).
  * Linder [2003] E. V. Linder, Exploring the expansion history of the universe, Physical Review Letters 90, 10.1103/physrevlett.90.091301 (2003).
  * Ferreira and Joyce [1997] P. G. Ferreira and M. Joyce, Structure formation with a self-tuning scalar field, Physical Review Letters 79, 4740 (1997).
  * Zlatev  _et al._ [1999] I. Zlatev, L. Wang, and P. J. Steinhardt, Quintessence, cosmic coincidence, and the cosmological constant, Physical Review Letters 82, 896 (1999).
  * Frusciante and Perenon [2020] N. Frusciante and L. Perenon, Effective field theory of dark energy: A review, Physics Reports 857, 1 (2020).
  * Aviles  _et al._ [2018] A. Aviles, M. A. Rodriguez-Meza, J. De-Santiago, and J. L. Cervantes-Cota, Nonlinear evolution of initially biased tracers in modified gravity, Journal of Cosmology and Astroparticle Physics 2018 (11), 013.
  * Wang and Steinhardt [1998] L. Wang and P. J. Steinhardt, Cluster abundance constraints for cosmological models with a time-varying, spatially inhomogeneous energy component with negative pressure, The Astrophysical Journal 508, 483 (1998).


# Computational design of optimal heterostructures for $\beta$-Ga2O3

Sierra Seacat  Department of Physics and Astronomy, University of Kansas,
Lawrence, KS 66045, USA<EMAIL_ADDRESS>John L. Lyons  Center for
Computational Materials Science, US Naval Research Laboratory, Washington, DC
20375, USA  Hartwin Peelaers  Department of Physics and Astronomy, University
of Kansas, Lawrence, KS 66045, USA<EMAIL_ADDRESS>
###### Abstract

Ga2O3 is a wide-bandgap material of interest for a wide variety of devices,
many of these requiring heterostructures, for instance to achieve carrier
confinement. A common method to create such heterostructures is to alloy with
In2O3 or Al2O3. However, the lattice constants of these materials are
significantly different from those of Ga2O3, leading to large amounts of
strain in the resulting heterostructure. By considering alloys of In2O3 and
Al2O3, the lattice constants can be tailored to those of Ga2O3, while still
keeping a sizeable conduction-band offset. We use density functional theory
with hybrid functionals to investigate the structural and electronic
properties of In2O3 and Al2O3 alloys in the bixbyite, corundum, and monoclinic
structures. We find that the lattice constants increase with In incorporation.
Bandgaps decrease with increasing In concentration, demonstrating bowing, with
$b=1.98$ eV for monoclinic, $b=4.27$ eV for corundum, and $b=4.9$ eV for
bixbyite. We find the (In0.25Al0.75)2O3 monoclinic structure to be of
particular interest, as it closely matches the Ga2O3 lattice constants while
providing a conduction-band offset of 1 eV.

As a wide-bandgap semiconductor (4.7 eV) Tippins (1965); Matsumoto  _et al._
(1974a), Ga2O3 is a promising material for solar blind photodetectors, which
have a variety of applications across different fields including military
surveillance, medical imaging, chemical and biological analysis, flame
detection, and secure communications. Suzuki  _et al._ (2009); Hou  _et al._
(2021); Guo  _et al._ (2019); Zou  _et al._ (2014); Zhang  _et al._ (2020)
Additionally, Ga2O3 possesses a large breakdown field (5-9 MVcm-1) Pearton
_et al._ (2018a); Higashiwaki  _et al._ (2017), which makes it a suitable
material for high-power devices such as field effect transistors and Schottky
diodes. Higashiwaki  _et al._ (2017); Pearton  _et al._ (2018b)

Fabricating such devices commonly requires creating heterostructures with a
conduction-band offset to confine charge carriers. One way to achieve this
without changing the crystal structure is by alloying Ga2O3 with oxides such
as In2O3 or Al2O3. At the same time, alloying enables the tailoring of
properties such as lattice constants and bandgaps.

Figure 1:  a) The primitive unit cell of corundum Al2O3. b) The primitive unit
cell of bixbyite In2O3 with the distorted octahedra shown in dark gray/purple
and the regular octahedra sites shown in light pink/purple. c) The primitive
unit cell of monoclinic Ga2O3 with the cations colored by site coordination.
The octahedral sites are shown in orange, while the tetrahedral sites are
shown in green.

In and Al are common candidates for alloying with Ga2O3, therefore we have
previously studied their alloys with Ga2O3 Peelaers  _et al._ (2018, 2015);
Seacat, Lyons, and Peelaers (2020, 2021). The ground state structures of
Ga2O3, In2O3, and Al2O3 are all different: Ga2O3 occurs in a monoclinic
crystal structure [referred to as $\beta$-Ga2O3, and shown in Fig. 1(a)], with
C2/m symmetry, In2O3 in a cubic bixbyite structure with Ia$\bar{3}$ symmetry
[Fig. 1(b)], and Al2O3 in a corundum structure with R3c symmetry (known as
$\alpha$-Al2O3) [Fig. 1(c)]. The lattice constants of these three materials
are all different, with those of In2O3 being the largest, followed by Ga2O3,
and those of Al2O3 being the smallest. This ordering of lattice constance
lengths occurs for both the ground state structures and for the (meta)stable
polymorphs. Peelaers  _et al._ (2018, 2015)

This mismatch in lattice constants is undesirable, as it creates strain in the
resulting alloy. Here we propose a solution to this strain by considering
alloys of In2O3 with Al2O3, which we can tailor to lattice-match monoclinic
Ga2O3. Furthermore, the design space for devices increases, allowing for
different bandgap and lattice constant values than those that are possible
when only looking at (AlxGa${}_{1{\rm-}x})_{2}$O3 and
(InxGa${}_{1{\rm-}x})_{2}$O3 alloys.

We use density functional theory (DFT) at the hybrid functional level to
investigate the structural and electronic properties of
(InxAl${}_{1{\rm-}x})_{2}$O3 alloys. We determine the lowest energy structures
for each In concentration, considering the ground state structures (bixbyite
and corundum), as well as the monoclinic structure (the desired structure as
it is the ground state structure of Ga2O3). The lowest energy configurations
for each of the structures are subsequently used to find the changes in
lattice constants, bandgaps, and valence- and conduction-band offsets as a
function of In concentration. As expected from Vegard’s law, the lattice
constants increase linearly with In content, as In2O3 has larger lattice
constants compared to Al2O3. In contrast, the bandgap decreases non-linearly
with increasing In content, since the bandgap of In2O3 is 3.14 eV, and thus
smaller than the 8.81 eV bandgap of Al2O3. We find that an alloy with 25% In
and 75% Al forms a perfect lattice match with Ga2O3, while still having a
large conduction-band offset.

All DFT calculations were performed using projector augmented wave (PAW)
potentials Blöchl (1994) as implemented in the Vienna Ab-initio Simulation
Package (VASP) Kresse and Hafner (1993); Kresse and Furthmüller (1996). In
order to obtain the necessary accuracy for the structural and electronic
properties, the HSE06 hybrid functional with the mixing parameter set to 32%
was used Heyd, Scuseria, and Ernzerhof (2003); *Heyd2006. To ensure
consistency with previous results and guarantee accurate energetic ordering,
the In pseudopotential with $d$ electrons as part of the valence was used
Peelaers  _et al._ (2018); Mu, Peelaers, and Van de Walle (2019); Sabino, de
Oliveira, and Da Silva (2014). All calculations were performed using a plane-
wave expansion cutoff of 500 eV with all structures relaxed until the forces
were smaller than 10 meV/Å and stresses smaller than 0.5 meV/Å3. A Monkhorst-
Pack 4$\times$4$\times$4 k-point grid was used for the corundum and primitive
monoclinic calculations, a 6$\times$4$\times$4 k-point grid for the
calculations using the conventional monoclinic cell, and a 2$\times$2$\times$2
k-point grid for the bixbyite calculations. Oxide surfaces are generated from
primitive alloy cells. Non-polar (110) surfaces are used, with at least 25 Å
of oxide and 25 Å of vacuum.

Table 1:  Calculated and experimental lattice constants and band gaps for the
ground state structures of Ga2O3, Al2O3, and In2O3.

|  $\beta$-Ga2O3 |  $\alpha$-Al2O3 | In2O3
---|---|---|---
| Calc. | Exp. | Calc. | Exp. | Calc. | Exp.
Symmetry | monoclinic | corundum | bixbyite
Spacegroup | C2/m | R3c | Ia$\bar{3}$
$a$ (Å) | 12.25 | 12.21a | 4.74 | 4.76c | 10.11 | 10.12e
$b$ (Å) | 3.04 | 3.04a | 4.74 | 4.76c |  |
$c$ (Å) | 5.79 | 5.80a | 12.94 | 12.99c |  |
$\beta$ | 103.82∘ | 103.83∘a |  |  |  |
$E^{\rm direct}_{\rm gap}$ (eV) | 4.83 | 4.76b | 8.81 | 8.8d | 3.14 | 2.9f

aRef. Åhman, Svensson, and Albertsson, 1996
bRef. Matsumoto  _et al._ , 1974b
cRef. Haan, 1962
dRef. French, 1990
eRef. Marezio, 1966
fRef. Irmscher  _et al._ , 2014

We consider alloys with monoclinic, corundum, and bixbyite structures. For the
monoclinic alloys, both the 10-atom primitive unit cell and the 20-atom
conventional unit cell are used to perform calculations. We use the 10-atom
primitive unit cell for the corundum alloys, and bixbyite calculations are
performed using the 40-atom primitive unit cell. Our calculated lattice
constants and bandgaps are listed in Table 1 with a comparison to experimental
values.

The different ground state structures, shown in Fig. 1, have different cation
bonding environments. The corundum and bixbyite cells only have octahedrally
coordinated cations, while the primitive monoclinic cell consists of two
octahedrally coordinated and two tetrahedrally coordinated cation sites. All
the cation sites in the corundum cell are identical, while there are two
different types of octahedral sites in the bixbyite structure: regular
octahedra with longer bond lengths and distorted octahedra with shorter bond
lengths. To determine the lowest energy structures for each composition, these
differences in site coordination are be considered.

We use the enthalpy of formation ($\Delta H$) to characterize the relative
stability of each possible alloy, which is defined as:

$\displaystyle\Delta
H[(\text{In}_{x}\text{Al}_{1-x})_{2}\text{O}_{3}]=E[(\text{In}_{x}\text{Al}_{1-x})_{2}\text{O}_{3}]$
$\displaystyle\qquad-(1-x)\,E[\text{Al}_{2}\text{O}_{3}]-x\,E[\text{In}_{2}\text{O}_{3}]\
,$ (1)

where $E[\text{Al}_{2}\text{O}_{3}]$ is the energy per f.u. of corundum Al2O3
and $E[\text{In}_{2}\text{O}_{3}]$ is the energy per f.u. of bixbyite In2O3.

We first determine the most favorable cation coordination environments for the
monoclinic and bixbyite cells. Starting with monoclinic Al2O3, we replace one
Al atom with an In atom for both an octahedral and a tetrahedral site. The
enthalpy of formation is then calculated for each structure to determine the
most likely configuration. The results for a single In atom replacement are
then used to repeat the analysis for larger concentrations of In. A similar
procedure is followed for the bixbyite structure, considering the distorted
and undistorted octahedral sites. In the monoclinic structure, In prefers to
occupy octahedral sites. For the bixbyite structure, single Al or In atoms
prefer to occupy the regular octahedra, but this preference is lost at higher
concentrations.

Figure 2:  Enthalpy of formation per formula unit (f.u.) as function of In
content for the most favorable configurations for each crystal structure. The
colors and symbols represent each of the structures: blue circles represent
corundum, red triangles represent bixbyite, and green squares represent
monoclinic. The dashed lines connect each of the calculated points with
straight lines and are shown as a visual aid.

Fig. 2 shows the formation enthalpies of the lowest energy configuration as
function of In concentration for the three different crystal structures, with
the corundum structures represented by blue circles, the bixbyite structures
by red triangles, and the monoclinic structures by green squares. Consistent
with the ground state structures of Al2O3 and In2O3, corundum structures have
the lowest formation enthalpy for In concentrations of $x=0.25$ and smaller.
For larger In concentrations ($x=0.75$ and greater), bixbyite structures are
energetically favored. For $x=0.5$ the monoclinic structure is the most
energetically favored. This indicates that 50-50 InAlO alloys will prefer
monoclinic crystal structures, and not the crystal structures of their ground
states (bixbyite or corundum).

For each of the lowest enthalpy configurations, we plot in Fig. 3 the
pseudocubic lattice constant (the cube root of the volume) as a function of In
content ($x$) for each of the crystal structures. Independent of the crystal
structure, the pseudocubic lattice constant increases with the In
concentration, consistent with the oxygen bond lengths in the ground state
structures: the In-O bixbyite In2O3 is longer than the Al-O bond length in
$\alpha$-Al2O3. The increase in the pseudocubic lattice constant is linear as
In concentration increases, following Vegard’s law. However, the monoclinic
structure shows a kink at $x=0.5$, with a different linear slope before and
after that concentration. The reason is that a different cation sublattice is
being occupied with In: for In concentrations smaller than $x=0.5$ the
octahedral sites are occupied with In atoms. Since there are only 50%
octahedral sites in the monoclinic structure [see Fig. 1(c)], the tetrahedral
sites starts to be filled after $x=0.5$, giving rise to a different slope.

Of particular interest is the 25% In monoclinic structure. It has monoclinic
lattice constants of $a$=3.03 Å, $b$=12.22 Å, $c$=5.79 Å, and $\beta$=102.88∘.
These lattice constants are quite close to those of monoclinic Ga2O3 ($a$=3.04
Å, $b$=12.25 Å, $c$=5.79 Å, and $\beta$=103.82∘). Creating heterostructures of
monoclinic Ga2O3 and (In0.25Al0.75)2O3 would therefore introduce minimal
strain.

Figure 3:  Pseudocubic lattice constant as function of In content for the
corundum (blue), bixbyite (red), and monoclinic (green) structures given by
Fig. 2. The solid lines are linear extrapolations.

Next we investigate how the bandgap changes with increasing In concentration,
as shown in Fig. 4. For each of the crystal structures, the bandgap decreases
as In content increases, which is commensurate with In2O3 possessing a smaller
bandgap than Al2O3. This decrease is not linear, which we can characterize by
the bowing parameter $b$, by fitting to

$\displaystyle E_{g}(x)=(1{\rm-}x)E_{g}[{\rm Al_{2}O_{3}}]+xE_{g}[{\rm
In_{2}O_{3}}]-bx(1{\rm-}x)\,,$ (2)

where $E_{g}[{\rm Al_{2}O_{3}}]$ and $E_{g}[{\rm In_{2}O_{3}}]$ are the
bandgaps of the 100% Al structure and the 100% In structure respectively. The
bowing parameter for each of the structures, obtained from Eq. 2 are $b=1.98$
eV for monoclinic, $b=4.27$ eV for corundum, and $b=4.9$ eV for bixbyite. This
indicates that the monoclinic structures exhibit less bowing (so are more
linear) compared to the bixbyite and corundum structures. The fitted curves
are also shown in Fig. 4.

While the values shown by Fig. 4 are for the direct bandgap at $\Gamma$, many
of the structures have an indirect bandgap. For the monoclinic structures, the
bandgap is indirect for all concentrations of In, with monoclinic In2O3
possessing an indirect gap of 3.04 eV (0.02 eV smaller than the direct gap)
and monoclinic Al2O3 an indirect gap of 7.21 eV (0.31 eV difference). The
bandgap is direct for the corundum In2O3 and Al2O3 endpoints, but becomes
indirect for the intermediate concentrations. Fitting Eq.(2) to the indirect
band gaps, gives slightly different values for the bowing, with $b=1.51$ eV
for monoclinic, $b=4.91$ eV for corundum, and $b=5.6$ eV for bixbyite.

For the 25% In alloy of interest, the direct bandgap is 5.94 eV and the
indirect bandgap is 5.70 eV. This is larger than the bandgap of Ga2O3,
indicating that confinement of carriers in Ga2O3 might be possible.

Figure 4:  Direct bandgap at $\Gamma$ as function of In content in each of the
crystal structures considered. The corundum structures are represented by the
blue circles, the bixbyite structures by the red triangles, and the monoclinic
structures by the green squares. The solid lines are parabolic fits to Eq.(2)
with $b=4.27$ for corundum, $b$=4.9 for bixbyite, and $b=1.98$ for monoclinic.

To ascertain whether the difference in bandgaps leads to a beneficial
conduction-band offset, we calculated the alignment of the valence- and
conduction-band edges by relating them to the vacuum level (which forms a
common reference energy). To do so we constructed different slabs, with at
least 25 Å of oxide and 25 Å of vacuum. We then use the planar-averaged
electrostatic potential in the direction perpendicular to the surfaceFranciosi
and Van de Walle (1996) to obtain the alignment. The results are shown in Fig.
5, with all energies relative to the VBM of monoclinic Al2O3, which is set to
zero energy.

We find that most of the band-gap change in the (In0.25Al0.75)2O3 alloy is
reflected in the conduction band, which steadily decreases as a function of
indium content. At 50% alloy content, the band positions are quite similar to
those of monoclinic Ga2O3. We also note that the band alignment between
monoclinic and corundum Al2O3 is within 0.1 eV of prior studies.Peelaers  _et
al._ (2018) Most importantly, we find that the (In0.25Al0.75)2O3 monoclinic
alloy has a conduction-band offset of 1 eV, so that any charge carriers in a
heterostructure with Ga2O3 should localize in the conduction band of Ga2O3.
Combined with the lattice matching, this makes (In0.25Al0.75)2O3 the ideal
alloy to use in device applications.

Figure 5:  The absolute energies of the valence- and conduction-band edges (in
eV, relative to vacuum), as determined with slab calculations for alloys of
corundum (blue dashed line) and monoclinic (green dotted line) (In1-xAlx)2O3.

In summary, we used hybrid density functional theory to design a monoclinic
alloy that lattice matches monoclinic Ga2O3, while providing a large
conduction-band offset. To do so, we investigated the structural and
electronic properties of (InxAl${}_{1{\rm-}x})_{2}$O3 alloys. We considered
bixbyite, corundum, and monoclinic crystal structures. Lattice constants
change linearly as a function of In concentration, except for monoclinic,
which exhibits a change in slope at 50% In alloys. Bandgaps vary non-linearly
with In concentration with bixbyite and corundum structures exhibiting more
non-linearity (bowing parameters $b$ of 4.9 and 4.27 eV) compared to
monoclinic alloys ($b$=1.98 eV). The 25% In-content monoclinic alloy
(In0.25Al0.75O3) is a perfect lattice match to monoclinic Ga2O3 and also has a
1 eV conduction-band offset, making in the ideal alloy to use for
heterostructures with Ga2O3.

##  I Acknowledgments*

We thank Hari Nair for helpful discussions. JLL was supported by the Office of
Naval Research (ONR) Basic Research Program. SC was supported by the
University of Kansas General Research Fund (2151089). Computational resources
were provided by Bridges-2 at Pittsburgh Supercomputing Center through
allocation DMR200079 from the Advanced Cyberinfrastructure Coordination
Ecosystem: Services & Support (ACCESS) program, which is supported by National
Science Foundation grants #2138259, #2138286, #2138307, #2137603, and
#2138296, by the DoD Major Shared Resource Center at AFRL, and by the
University of Kansas Center for Research Computing (CRC), including the BigJay
Cluster resource funded through NSF Grant MRI-2117449.

##  II Data Availability Statement*

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Tippins (1965) H. H. Tippins, Phys. Rev. 140, A316 (1965).
  * Matsumoto  _et al._ (1974a) T. Matsumoto, M. Aoki, A. Kinoshita, and T. Aono, Jpn. J. Appl. Phys. 13, 737 (1974a).
  * Suzuki  _et al._ (2009) R. Suzuki, S. Nakagomi, Y. Kokubun, N. Arai, and S. Ohira, Appl. Phys. Lett. 94, 222102 (2009).
  * Hou  _et al._ (2021) X. Hou, Y. Zou, M. Ding, Y. Qin, Z. Zhang, X. Ma, P. Tan, S. Yu, X. Zhou, X. Zhao, G. Xu, H. Sun, and S. Long, J. Phys. D: Appl. Phys. 54, 043001 (2021).
  * Guo  _et al._ (2019) D. Guo, Q. Guo, Z. Chen, Z. Wu, P. Li, and W. Tang, Mater. Today Phys. 11, 100157 (2019).
  * Zou  _et al._ (2014) R. Zou, Z. Zhang, Q. Liu, J. Hu, L. Sang, M. Liao, and W. Zhang, Small 10, 1848 (2014).
  * Zhang  _et al._ (2020) J. Zhang, J. Shi, D.-C. Qi, L. Chen, and K. H. L. Zhang, APL Mater. 8, 020906 (2020).
  * Pearton  _et al._ (2018a) S. J. Pearton, F. Ren, M. Tadjer, and J. Kim, J. Appl. Phys. 124, 220901 (2018a).
  * Higashiwaki  _et al._ (2017) M. Higashiwaki, A. Kuramata, H. Murakami, and Y. Kumagai, J. Phys. D: Appl. Phys. 50, 333002 (2017).
  * Pearton  _et al._ (2018b) S. J. Pearton, J. Yang, P. H. Cary, F. Ren, J. Kim, M. J. Tadjer, and M. A. Mastro, Appl. Phys. Rev. 5, 011301 (2018b).
  * Peelaers  _et al._ (2018) H. Peelaers, J. B. Varley, J. S. Speck, and C. G. Van de Walle, Appl. Phys. Lett. 112, 242101 (2018).
  * Peelaers  _et al._ (2015) H. Peelaers, D. Steiauf, J. B. Varley, A. Janotti, and C. G. Van de Walle, Phys. Rev. B 92, 085206 (2015).
  * Seacat, Lyons, and Peelaers (2020) S. Seacat, J. L. Lyons, and H. Peelaers, Appl. Phys. Lett. 116, 232102 (2020).
  * Seacat, Lyons, and Peelaers (2021) S. Seacat, J. L. Lyons, and H. Peelaers, Appl. Phys. Lett. 119, 042104 (2021).
  * Blöchl (1994) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Kresse and Hafner (1993) G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  * Heyd, Scuseria, and Ernzerhof (2003) J. Heyd, G. E. Scuseria, and M. Ernzerhof, J. Chem. Phys. 118, 8207 (2003).
  * Heyd, Scuseria, and Ernzerhof (2006) J. Heyd, G. E. Scuseria, and M. Ernzerhof, J. Chem. Phys. 124, 219906 (2006).
  * Mu, Peelaers, and Van de Walle (2019) S. Mu, H. Peelaers, and C. G. Van de Walle, Appl. Phys. Lett. 115, 242103 (2019).
  * Sabino, de Oliveira, and Da Silva (2014) F. P. Sabino, L. N. de Oliveira, and J. L. F. Da Silva, Phys. Rev. B 90, 155206 (2014).
  * Åhman, Svensson, and Albertsson (1996) J. Åhman, G. Svensson, and J. Albertsson, Acta Crystallogr. C 52, 1336 (1996).
  * Matsumoto  _et al._ (1974b) T. Matsumoto, M. Aoki, A. Kinoshita, and T. Aono, Jpn. J. Appl. Phys. 13, 1578 (1974b).
  * Haan (1962) Y. D. Haan, Z. Krist.-Cryst. Mate. 117, 235 (1962).
  * French (1990) R. H. French, J. Am. Ceram. Soc. 73, 477 (1990).
  * Marezio (1966) M. Marezio, Acta Crystallogr. 20, 723 (1966).
  * Irmscher  _et al._ (2014) K. Irmscher, M. Naumann, M. Pietsch, Z. Galazka, R. Uecker, T. Schulz, R. Schewski, M. Albrecht, and R. Fornari, Phys. Status Solidi A 211, 54 (2014).
  * Franciosi and Van de Walle (1996) A. Franciosi and C. G. Van de Walle, Surface Science Reports 25, 1 (1996).