#  $\mathcal{N}$-field cosmology in hyperbolic field space: stability and
general solutions

Perseas Christodoulidis  Andronikos Paliathanasis

###### Abstract

We study the dynamics of a cosmological model with a perfect fluid and
$\mathcal{N}$ fields on a hyperbolic field space interacting via a symmetric
potential. We list all late-time solutions, investigate their stability and
briefly discuss predictions of the theory. Moreover, for the case of two
scalar fields and an exponential potential we prove that the field equations
are Liouville integrable and we provide for the first time the general
solution for a region of the parameter space.

98.80.-k, 95.35.+d, 95.36.+x

##  1 Introduction

Over the last years, two-field models with a symmetric potential and
hyperbolic field space have been extensively studied in the literature in the
context of multi-field inflation [1, 2] or late-time universe [3, 4]. These
models have displayed interesting phenomenology while remaining
observationally viable [2, 5, 6, 7, 8]. 111Here we refer to the predictions of
this theory during the inflationary era, whereas the viability of
quintessence-like models during late-time cosmology has recently been
challenged in the literature (see e.g. [9, 10]). Though most works have
focused on the two-field regime, certain many-field constructions have also
been proposed as in e.g. [11, 12]. Similarly, some progress has been made in
the derivation of general solutions in scalar-field cosmology. On the
contrary, multi-field generalizations have been proved more challenging and up
to date only a few solutions are known for arbitrary number of fields [13, 14,
15, 16, 17].

The existence of exact and analytic solutions is an essential property for the
mathematical description of a physical theory. Although a dynamical system can
be solved by using numerical techniques we do not know that the numerical
trajectories correspond always to real solution of the problem, thus we should
investigate if the dynamical system posses the integrability property. There
are various techniques for the study of the integrability in the literature.
In cosmological studies, due to the fact that the gravitational field
equations for scalar field theories admit a minisuperspace description,
techniques from analytic mechanics can be applied. The theory of similarity
transformations for the derivation of conservation laws has been applied in
[18, 19, 20, 21, 22] while some other approaches can be found in [23, 24, 25].
Another important approach for the study of a cosmological model is the
determination of the stationary points. The latter points can be used for the
determination of the asymptotic behaviour of a specific theory and to extract
important information and criteria for the cosmological evolution of the
specific model [26, 27, 28, 29, 30].

In this work we will investigate the $\mathcal{N}$-field generalization of the
two-field hyperbolic problem in the presence of a perfect fluid. We will first
list all critical-point solutions and investigate their stability. Next, we
will apply the Noether method to derive two-field general solutions in some
cases and then extend them to $\mathcal{N}$-fields.

The paper is organized as follows: in Sec. 3 we revisit the stability analysis
for the two-field hyperbolic problem in the presence of a fluid. Next, we
generalize the discussion for $\mathcal{N}$ fields in Sec. 4 and list all new
solutions as well as their stability properties. In Sec. 5 we calculate the
general solution for a subset of the parameter space using the Noether method.
Finally, in Sec. 6 we offer our conclusions.

##  2 Chiral (multi-field) cosmology

The Chiral cosmological model belongs to the family of the Einstein-nonlinear
$\sigma$-model where in the Einstein-Hilbert action two scalar fields
minimally coupled to gravity are introduced such that the gravitational action
integral can be written as follows [31, 32, 33, 34]

$S=\int\sqrt{-g}\,\mathrm{d}^{4}X\left(R-\frac{1}{2}g^{\mu\nu}\nabla_{\mu}\phi\nabla_{\nu}\phi-\frac{1}{2}g^{\mu\nu}F\left(\phi\right)\nabla_{\mu}\psi\nabla_{\nu}\psi-V\left(\phi\right)\right)+S_{m}\,,$
(2.1)

where there is a coupling in the kinetic term between the two scalar fields.
When the coupling function $F\left(\phi\right)$ is constant the Action
Integral (2.1) describes two quintessence fields; however, in Chiral cosmology
the two dynamics of the two fields evolve in a space of constant non-zero
curvature, that is, $F\left(\phi\right)=e^{\kappa\phi}$. At this point it is
important to mention that we refer to a two-dimensional space defined by the
kinetic terms of the scalar fields and not in the background space with metric
$g_{\mu\nu}$ and Ricci scalar $R$.

According to the cosmological principle for the background space we assume
that of spatially flat Friedmann–Lemaître–Robertson–Walker (FLRW) described by
the line element

$\,\mathrm{d}s^{2}=-N_{l}\left(t\right)^{2}\,\mathrm{d}t^{2}+a\left(t\right)^{2}\left(\,\mathrm{d}X^{2}+\,\mathrm{d}Y^{2}+\,\mathrm{d}Z^{2}\right),$
(2.2)

where $a\left(t\right)$ is the scale factor and $N_{l}\left(t\right)$ the
lapse function. Moreover, we assume that the scalar fields inherit the
symmetries of the background space, that is,
$\phi\left(x^{\mu}\right)=\phi\left(t\right)$,
$\psi\left(x^{\mu}\right)=\psi\left(t\right)$, while the Action Integral
$S_{m}$ describes an ideal gas with energy density $\rho$, pressure $p$ and
constant equation of state parameter $p=w\rho$.

Hence, the gravitational field equations are [34]

$\displaystyle 3H^{2}$ $\displaystyle=\rho_{f}+\rho,$ (2.3)
$\displaystyle-\left(2\dot{H}+3H^{2}\right)$ $\displaystyle=p_{f}+P,$ (2.4)

in which $\rho_{f},$ $p_{f}$ are the energy density and pressure components of
the two scalar fields, that is,

$\displaystyle\rho_{f}$
$\displaystyle=\frac{1}{2}\dot{\phi}^{2}+\frac{1}{2}e^{\kappa\phi}\dot{\psi}^{2}+V\left(\phi\right),$
(2.5) $\displaystyle p_{f}$
$\displaystyle=\frac{1}{2}\dot{\phi}^{2}+\frac{1}{2}e^{\kappa\phi}\dot{\psi}^{2}-V\left(\phi\right).$
(2.6)

Finally, the equations of motion for the two fields are

$\displaystyle\ddot{\phi}+3H\dot{\phi}-\kappa\frac{1}{2}e^{\kappa\phi}\dot{\psi}^{2}+V_{,\phi}\left(\phi\right)=0\,,$
(2.7) $\displaystyle\ddot{\psi}+3H\dot{\psi}+\kappa\dot{\phi}\dot{\psi}=0~{},$
(2.8)

while the ideal gas satisfies the continuity equation

$\dot{\rho}+3H(1+w)\rho=0\,,$ (2.9)

from which it follows $\rho=\rho_{m0}a^{-3\left(1+w\right)}$.

At this point it is important to mention that the gravitational field
equations can be derived by the variation of the point-like Lagrangian

$L\left(a,\dot{a},\phi,\dot{\phi},\psi,\dot{\psi}\right)=\frac{1}{2N_{l}\left(t\right)}\left(-6a\dot{a}^{2}+a^{3}\left(\dot{\phi}^{2}+e^{\kappa\phi}\dot{\psi}^{2}\right)\right)-N\left(t\right)a^{3}V\left(\phi\right)+N_{l}\left(t\right)\rho_{m0}a^{-3w_{m}}.$
(2.10)

That it is an important observation, because we can apply techniques from
analytic mechanics for the determination of exact solutions. As far as the
scalar field potential $V\left(\phi\right)$ is concerned, we will assume the
exponential function $V\left(\phi\right)=V_{0}e^{\lambda\phi}$ for the most
part of this paper which has been shown to provide interesting physical
results [35, 36].

In order to study the asymptotic behaviour of the previous model it is better
to switch from cosmic time to the e-folding number defined from
$\,\mathrm{d}N=H\,\mathrm{d}t$. In this way the set of evolution equations for
the scalar fields and the fluid can be written as an autonomous dynamical
system. Introducing a new variable $z=\rho/H^{2}$, describing the fluid energy
density, the Friedman constraint becomes

$3={1\over 2}v^{i}v_{i}+z+V\,,$ (2.11)

which implies that the allowed values for the field velocities and $z$ should
satisfy

${1\over 2}v^{i}v_{i}+z\leq 3\,.$ (2.12)

The slow-roll parameter becomes

$\epsilon={1\over 2}v^{i}v_{i}+{1\over 2}(1+w)z\,,$ (2.13)

while the potential satisfies

${V\over H^{2}}=3-\epsilon+{1\over 2}(w-1)z\,.$ (2.14)

The evolution equations for the fields and the fluid are

$\displaystyle(v^{i})^{\prime}+\Gamma^{i}_{jk}v^{j}v^{k}+(3-\epsilon)(v^{i}+\lambda^{i})+{1\over
2}(w-1)\lambda^{i}z=0\,,$ (2.15) $\displaystyle
z^{\prime}+(3+3w-2\epsilon)z=0\,.$ (2.16)

However, critical points for the velocities are not expected to be found when
a generic field metric is considered because Christoffel symbols are field
dependent. Instead, it is better to study equations for the normalized
velocities $\sqrt{G_{ii}}v^{i}$ (no sum is assumed in $i$) that enter the
definition of $\epsilon$. For the field metric with an isometry these
velocities are defined as

$y\equiv\phi^{\prime}\,,\qquad x\equiv\sqrt{F(\phi)}\psi^{\prime}\,.$ (2.17)

The variables of the two-field dynamical system in first order form are
$\phi,\psi,y,x,z$.

##  3 Recap of the two-field problem with a perfect fluid

###  3.1 Hyperbolic field metric

Specializing to the hyperbolic field metric

$\,\mathrm{d}s^{2}=\,\mathrm{d}\phi^{2}+e^{\kappa\phi}\,\mathrm{d}\psi^{2}\,,$
(3.1)

the system in first order form is

$\displaystyle\phi^{\prime}=y\,,$ (3.2)
$\displaystyle\psi^{\prime}=xe^{-\kappa/2\phi}\,,$ (3.3) $\displaystyle
y^{\prime}+(3-\epsilon)(y+\lambda)-{\kappa\over 2}x^{2}+{1\over 2}(w-1)\lambda
z=0\,,$ (3.4) $\displaystyle x^{\prime}+\left(3-\epsilon+{\kappa\over
2}y\right)x=0\,,$ (3.5) $\displaystyle z^{\prime}+(3+3w-2\epsilon)z=0\,.$
(3.6)

The second equation can be discarded because $\psi$ is a cyclic variable and
does not affect dynamics and, moreover, we observe that the last three
equations do not depend on $\phi$ and so the first equation can also be
omitted; the reduced $3\times 3$ system is sufficient to extract information
regarding the solution. Although the stability analysis for this model has
been presented recently in Ref. [3], in this section we will revisit it in
order to facilitate the transition to more fields. In addition, we will
mention the stability criteria for generic $w$ that were missing for some
solutions of the aforementioned work.

For our set of variables the eigenvalues of the Jacobian matrix follow
straightforwardly from our analysis because that matrix can always be written
in block-diagonal form. We have the following critical points:

  1. 1.

First, we have scalar field domination solutions, which generalize the
solutions presented in [36, 37] with the addition of $z=0$. Their stability
properties remain unchanged with the additional requirement
$\epsilon<3/2(1+w)$. This happens because the stability matrix acquires an
upper triagonal form with zeros below the main diagonal and so the presence of
the fluid does not affect eigenvalues (for this type of solutions). There are
three types of scalar-dominated solutions:

    1. (a)

The scalar-field gradient solution,

$(y,x,z)_{\rm gr}=\left(-\lambda,0,0\right)\,.$ (3.7)

Motion is aligned with the potential gradient flow. It is stable provided

$-\sqrt{6+{\kappa^{2}\over 4}}-{\kappa\over
2}<\lambda<\sqrt{6+{\kappa^{2}\over 4}}-{\kappa\over
2}\,,\qquad|\lambda|<\sqrt{6}\,,\qquad|\lambda|<\sqrt{3(1+w)}\,.$ (3.8)

The solution is depicted at the left panel of Fig. 1.

Figure 1:  The numerical solutions for a wide range of initial conditions
(drawn uniformly from the surface defined from $\epsilon\approx 3$). Left: For
$\lambda=\kappa=1$ and $w=0$ the solution asymptotes to the gradient critical
point. Right: For a fluid with $w=2$, $\lambda=3$ and $\kappa=-1$ the solution
asymptotes to kinetic domination. Blue dots correspond to the respective
critical points and the semi-transparent blue surface denotes the region of
definition for $x,y,z$.

    2. (b)

The scalar-field kinetic domination solution given as

$(y,x,z)_{\rm kin}=\left(\pm\sqrt{6},0,0\right)\,,$ (3.9)

and is stable for

$|\lambda|>\sqrt{6}\,,\qquad\lambda\cdot\kappa<0\,,\qquad 1+w>2\,,$ (3.10)

(see left panel of Fig. 1).

    3. (c)

The hyperbolic solution

$(y,x,z)_{\rm
hyper}=\left(-{6\over\kappa+\lambda},\pm{\sqrt{6}\sqrt{\lambda^{2}+\kappa\lambda-6}\over\kappa+\lambda},0\right)\,.$
(3.11)

This solution exists provided $k\neq-\lambda$, $\lambda^{2}+\kappa\lambda-6>0$
or

$\lambda>\sqrt{6+{\kappa^{2}\over 4}}-{\kappa\over
2}\,,\qquad\lambda<-\sqrt{6+{\kappa^{2}\over 4}}-{\kappa\over 2}\,,$ (3.12)

and it is stable whenever the previous two solutions are unstable, namely when
the following conditions are satisfied

$\displaystyle\kappa>0\,,\qquad\lambda>\sqrt{6+{\kappa^{2}\over
4}}-{\kappa\over 2}\,,\qquad{\lambda\over\kappa+\lambda}<{1\over 2}(1+w)$
(3.13) $\displaystyle\kappa<0\,,\qquad\lambda<-\sqrt{6+{\kappa^{2}\over
4}}-{\kappa\over 2}\qquad{\lambda\over\kappa+\lambda}<{1\over 2}(1+w)\,.$
(3.14)

The solution is depicted at the left panel of Fig. 2.

Figure 2:  The numerical solutions assuming a fluid with $w=0$ for a wide
range of initial conditions (drawn uniformly from the surface defined from
$\epsilon\approx 3$). Left: For $\lambda=1$, $\kappa=15$ the solution
asymptotes to one of the two hyperbolic critical points (depending on the
initial $x$). Right: For $\lambda=2.2$, $\kappa=1$ the solution asymptotes to
the scaling critical point. Blue dots correspond to the respective critical
points and the semi-transparent blue surface denotes the region of definition
for $x,y,z$.

  2. 2.

The second type of solutions describe fluid domination $z=3$, where fields
have zero values. In order for this to be a solution of the dynamical system
the parenthesis of Eq. (3.6) should be zero and this gives the value of
$\epsilon$

$\epsilon={3\over 2}(1+w)\,,$ (3.15)

(this value of $\epsilon$ is compatible with Eq. (3.4)). These solutions are
stable for $w<-1$, which yields $\epsilon<0$, and thus describe contracting
universes.

  3. 3.

Finally, we find the scaling solution with the fluid and at least one of the
fields non-zero. Again we require the parenthesis of Eq. (3.6) to vanish and
so $\epsilon$ has the same value as in fluid domination. The solution is

$(y,x,z)_{\rm
scal}=\left(-{3(w+1)\over\lambda},0,3-{9(1+w)\over\lambda^{2}}\right)\,,$
(3.16)

which exists provided $z_{\rm scal}\geq 0$ or

$\lambda^{2}\geq 3(1+w)\,.$ (3.17)

The eigevalues of the Jacobian matrix are

$\displaystyle m_{1}$ $\displaystyle={3\over
2\lambda}\left[(w+1)\kappa+(w-1)\lambda\right]\,,$ (3.18) $\displaystyle
m_{\pm}$ $\displaystyle={3\over
4}\left(w-1\pm\sqrt{(w-1)^{2}+8(w^{2}-1)-24(w-1)(w+1)^{2}\lambda^{-2}}\right)\,,$
(3.19)

and they are non-positive for $w>-1$ with the following additional
restrictions on $w,\lambda,\kappa$

$\displaystyle-1<w<1\,,\qquad{\kappa\over\lambda}<{1-w\over 1+w}\,.$ (3.20)

For $w=0$ we recover the relations mentioned in Ref. [3]. The solution is
illustrated at the right panel of Fig. 2. Note that there are no real
solutions with $y,x,z\neq 0$.

###  3.2 Generic field metric with isometry

We will briefly comment on the case of a general field metric with isometry in
$\psi$

$\,\mathrm{d}s^{2}=\,\mathrm{d}\phi^{2}+F(\phi)\,\mathrm{d}\psi^{2}\,.$ (3.21)

In this case $\kappa\equiv F_{,\phi}/F$ is field dependent and Eq. (3.2) can
not be omitted. Choosing an exponential potential, or a potential that
asymptotes to an exponential at e.g. $-\infty$, all previous solutions (except
for the hyperbolic one) may exist for appropriate choices of the metric
function (see Fig. 3 for a model with $F=e^{\phi^{2}}$), while the hyperbolic
solution is replaced by a de Sitter asymptotic state $y,x\rightarrow 0$.

Figure 3:  The numerical solutions assuming a fluid with $w=0$ and
$\phi_{0}=2$ for a wide range of initial conditions (drawn uniformly from the
surface defined from $\epsilon\approx 3$). Left: For $\lambda=1$ the solution
asymptotes to the gradient critical point. Right: For $\lambda=2.2$ the
solution asymptotes to the scaling critical point. Blue dots correspond to the
respective critical points and the semi-transparent blue surface denotes the
region of definition for $x,y,z$.

This can be shown as follows: if the gradient or kinetic solutions are
unstable then a solution, that resembles the hyperbolic one, can be obtained
if $\kappa$ diverges to plus/minus infinity at the boundary of the space. In
this case the combination $\kappa y$ is required to be constant and so the
parenthesis of Eq. (3.5) can vanish. Plugging back into Eq. (3.4) gives the
asymptotic solution for $y$ and $x$ which has exactly the same form as in the
hyperbolic case, albeit $\kappa$ is field dependent and growing in norm [36,
4]. Even though a proper stability analysis for a field-dependent $\kappa$
requires study of the dynamical system at infinity, we can use a simpler
argument to understand the behaviour of these solutions. For the $4\times 4$
dynamical system calculating Lyapunov exponents, results to one zero
eigenvalues which is associated with $\phi$. Since $\phi$ will eventually roll
towards decreasing values of the potential (otherwise the Friedman constraint
would be violated) no instability related to this marginal direction will be
present. Therefore, one can apply the previous formulae for solutions and
stability criteria after taking the limit $\phi\rightarrow\infty$. Note though
that convergence towards the de Sitter critical point is much slower compared
to other critical points, while the system may pass through other critical
points first (see Fig. 4 for a model with $F=e^{\phi^{3}}$ as well as the
discussion in Ref. [4]).

Figure 4:  Numerical solutions that asymptote to the de Sitter critical point
with $\phi_{0}=2$, $\lambda=1$.

##  4 The $\mathcal{N}$\- field hyperbolic solution

To properly generalize the model to $\mathcal{N}$ fields we choose the
following form of the field metric

$\,\mathrm{d}s^{2}=\,\mathrm{d}\phi^{2}+\sum
e^{\kappa_{i}\phi}\,\mathrm{d}\psi_{i}^{2}\,,$ (4.1)

with Ricci scalar

$R=-{1\over 4}\left(\sum\kappa_{i}\right)^{2}-{1\over 2}\sum\kappa_{i}^{2}\,.$
(4.2)

Equations of motion are

$\displaystyle\phi^{\prime}=y\,,$ (4.3)
$\displaystyle\psi^{\prime}_{i}=x_{i}e^{-\kappa_{i}\phi}\,,$ (4.4)
$\displaystyle y^{\prime}+(3-\epsilon)(y+\lambda)-\sum{\kappa_{i}\over
2}x_{i}^{2}+{1\over 2}(w-1)\lambda z=0\,,$ (4.5) $\displaystyle
x_{i}^{\prime}+\left(3-\epsilon+{\kappa_{i}\over 2}y\right)x_{i}=0\,,$ (4.6)
$\displaystyle z^{\prime}+(3+3w-2\epsilon)z=0\,.$ (4.7)

To analyse the problem we will distinguish between two cases.

###  4.1 $\kappa_{i}$ are all equal

In the symmetric case $\kappa_{i}=\kappa$ the solutions presented in Sec. 3
carry over with the substitution $x^{2}\rightarrow\sum x_{i}^{2}$. This
becomes apparent if we write the differential equation of the slow-roll
parameter $\epsilon$ which is found by contracting Eq. (2.15) with $v_{i}$:

$\epsilon^{\prime}+(3-\epsilon)(2\epsilon+\lambda y)+{1\over 2}(w-1)\lambda
yz=0\,,$ (4.8)

and substituting $\sum x_{i}^{2}=2\epsilon-y^{2}$ in Eq. (4.5)

$y^{\prime}+(3-\epsilon)(y+\lambda)-{\kappa\over 2}(2\epsilon-y^{2})+{1\over
2}(w-1)\lambda z=0\,.$ (4.9)

This shows that the set of Eqs. (4.6) can be replaced with Eq. (4.8) and
$x_{i}$ are left undetermined (a similar argument was used in Ref. [38]).

###  4.2 $\kappa_{i}$ are different

The situation is drastically different when $\kappa_{i}$ are different. In the
next we list all critical-point solutions and their stability properties.

  1. 1.

The analogue of the hyperbolic solution with $y$ and all $x_{i}$ different
than zero is inconsistent as it requires

$(3-\epsilon)+{\kappa_{i}\over 2}y=0\,,$ (4.10)

to hold for every $\kappa_{i}$. Therefore, we conclude that only one $x_{j}$
can be non-zero and $x_{i}=0$ for $i\neq j$. To study the stability we
calculate the Jacobian matrix evaluated on the hyperbolic solution and we
observe that it always acquires a block diagonal form (some permutations of
rows and columns may be necessary)

$\begin{pmatrix}A_{2\times 2}&0_{2\times\mathcal{N}-1}\\\
0_{\mathcal{N}-1\times
2}&B_{\mathcal{N}-1\times\mathcal{N}-1}\end{pmatrix}\,,$ (4.11)

where

$A_{2\times
2}={1\over(\kappa_{j}+\lambda)^{2}}\begin{pmatrix}36-3(\kappa_{j}-\lambda)(\kappa_{j}+2\lambda)&\sqrt{6}\sqrt{\lambda^{2}+\kappa_{j}\lambda-6}[(\kappa_{j}-\lambda)^{2}-6]\\\
-\sqrt{{3\over
2}}(12+\kappa_{j}^{2}+\lambda\kappa_{j})\sqrt{\lambda^{2}+\kappa_{j}\lambda-6}&6[(\kappa_{j}-\lambda)^{2}-6]\end{pmatrix}\,,$
(4.12)

is exactly the stability matrix of the reduced two-field problem, while the
other matrix is diagonal

$B_{\mathcal{N}-1\times\mathcal{N}-1}=\text{diag}\left({-3(\kappa_{j}-\kappa_{i})\over\kappa_{j}+\lambda},\cdots,3+3w-2\epsilon\right)\,,$
(4.13)

for $i=1,\cdots,\mathcal{N}-2\neq j$. The eigenvalues of $A$ need to satisfy
the inequalities (3.13)-(3.14) (with the substitution
$\kappa\rightarrow\kappa_{j}$) while the rest $\mathcal{N}-1$ eigenvalues are
the diagonal elements of $B$ and so a stable solution requires

$\displaystyle\kappa_{j}>\kappa_{i}~{}~{}\text{for}~{}~{}\kappa_{j}>0\,,\qquad\kappa_{j}<\kappa_{i}~{}~{}\text{for}~{}~{}\kappa_{j}<0\,,$
(4.14)

for $i\neq j$.

Figure 5:  The numerical solution for three fields and a wide range of initial
conditions (drawn uniformly from the surface of sphere with $\epsilon\approx
3$) with $\lambda=3$, $\kappa_{1}=1$, $z=0$ and $\kappa_{2}=-1$ (left) and
$\kappa_{2}=10$ (right). The blue dots correspond to the two hyperbolic
solutions.

  2. 2.

The fluid domination and the scaling solution with non-zero $y,z$ and
$x_{i}=0$

$(y,x_{i},z)_{\rm
scal}=\left(-{3(w+1)\over\lambda},0,\cdots,3-{9(1+w)\over\lambda^{2}}\right)\,,$
(4.15)

both exist with the same stability properties as previously.

  3. 3.

Finally, in addition to the usual kinetic-domination solution

$(y,x_{i},z)_{\rm kin}=\left(\pm\sqrt{6},0,\cdots,0\right)\,,$ (4.16)

which is stable for $|\lambda|>\sqrt{6}$ and $\lambda\kappa_{i}<0$, new
solutions ($\epsilon=3$) exist with $y=0$ and $x_{i}\neq 0$ which follow from

$\kappa_{i}x_{i}^{2}=0\,,\qquad\sum x_{i}^{2}=6.$ (4.17)

For three fields the solution is trivially found to be

$(y,x_{1},x_{2},z)_{\rm
kin}=\left(0,\pm\sqrt{6}\sqrt{{\kappa_{1}\over\kappa_{1}-\kappa_{2}}},\pm\sqrt{6}\sqrt{-{\kappa_{2}\over\kappa_{1}-\kappa_{2}}},0\right)\,,$
(4.18)

which exists provided $\kappa_{1}\kappa_{2}<0$. It can be shown that the
stability matrix evaluated on the solution contains at least one positive
eigenvalue and thus this kinetic domination is unstable. For $\mathcal{N}>4$
fields this type of kinetic solution is defined on the $\mathcal{N}-2$
hypersurface containing points that satisfy the relations of (4.17) and it is
again unstable.

Having thoroughly examined the background in the next subsection we move to
the study of quantum fluctuations neglecting the fluid’s energy density and
fluctuations.

###  4.3 Observables

To extract physical quantities it is necessary to express gauge-invariant
perturbations in the orthonormal basis (local Frenet system)
$Q^{i}=E^{~{}i}_{A}q^{A}$ where the matrix $E^{~{}i}_{A}$ has as columns the
components of the orthonormal vectors

$E^{~{}i}_{A}=\left(t^{i},n^{i},b^{i},\cdots\right)\,.$ (4.19)

Here, $t^{i}$ is the tangent unit vector, $n^{i}$ the normal vector, $b^{i}$
the binormal vector and so on. The projection along the tangent vector is
related the curvature perturbation, while projections along the orthogonal
directions are related to isocurvature perturbations. The covariant time
derivatives of the orthonormal vectors satisfy the Frenet-Serret equations

$\,\text{D}_{t}E^{~{}i}_{A}=C^{~{}B}_{A}E^{~{}i}_{B}\,,$ (4.20)

where the matrix $C$ is antisymmetric with non-zero elements in the upper and
lower diagonals.

It is well known that the curvature perturbation ($Q_{\sigma}\equiv
Q_{i}E^{~{}i}_{1}=q_{1}$) is sourced by the first isocurvature perturbation
($Q_{s}\equiv Q_{i}E^{~{}i}_{2}=q_{2}$) [39, 40], while the rest orthogonal
perturbations interact via a “mass matrix” as well as through the curvatures
that appear in the Frenet-Serret equations (see e.g. Refs. [41, 42, 43] for
examples including up to three fields and Refs. [44, 45] for a formal
discussion including an arbitrary number of fields). More specifically,
orthogonal fields are coupled through the following terms in the second order
action

$2\dot{q}^{T}\cdot\Omega\cdot q+q^{T}\cdot\Omega^{T}\cdot\Omega\cdot
q-q^{T}\cdot M^{2}\cdot q\,,$ (4.21)

where $\Omega$ is the truncated matrix obtained from $C$ after removing
elements of the first row and column, and the “mass matrix” is defined from

$M^{2}_{AB}\equiv
E^{~{}i}_{A}E^{~{}j}_{B}\left(V_{;ij}+\dot{\sigma}^{2}t^{k}t^{l}R_{kilj}+3\omega^{2}\delta_{A2}\delta_{B2}\right)\,,$
(4.22)

for $A,B>1$. Extra orthogonal fields decouple from $Q_{\sigma}$ and $Q_{s}$
when the matrices $\Omega$ and $M$ are block diagonal. A necessary condition
for this decoupling is the vanishing of the torsion of the
$\mathcal{N}$-dimensional field-space trajectory. Recall that the torsion of a
curve is found by calculating the rate of change of the normal unit vector.
Using the Frenet-Serret equations the torsion ($\tau$) is defined from

$\,\text{D}_{t}n^{i}=-\omega t^{i}+\tau b^{i}\,,$ (4.23)

where $t^{i},n^{i},b^{i}$ are the first three unit vectors of the orthonormal
frame at some point of the curve and $\omega$ is the turn rate. The case of
zero torsion is reminiscent to geodesic motion where the curvature
perturbation decouples from isocurvature perturbations.

For the hyperbolic solution the vectors $t^{i},n^{i}$ have the first two
components non-zero and the rest zero, which forces the next vectors in the
series to have the first two components zero and the rest non-zero:

$\displaystyle t^{i}$ $\displaystyle=(t^{\phi},t^{\chi},0,\cdots)\,,$ (4.24)
$\displaystyle n^{i}$ $\displaystyle=(n^{\phi},n^{\chi},0,\cdots)\,,$ (4.25)
$\displaystyle b^{i}$ $\displaystyle=(0,0,b^{\chi_{2}},\cdots)\,,$ (4.26)
$\displaystyle\cdots$ (4.27)

This means that the matrix $E^{~{}i}_{A}$ is block diagonal. Using the Frenet-
Serret equation for $n^{i}$ we find $\text{D}_{t}n^{i}=-\omega t^{i}$ and,
hence, the torsion is zero, which implies that $\Omega$ is block diagonal.
Moreover, $M^{2}$ turns out to be block diagonal and so we conclude that
$Q_{\sigma}$ and $Q_{s}$ evolve independently from the rest fields. Therefore,
the basic predictions for this model, namely the spectral index and the
tensor-to-scalar ratio, are identical to the two-field case.

###  4.4 Field metric with $\mathcal{N}-1$ isometries

Similar to Sec. 3.2 we can consider metrics that can support $\mathcal{N}-1$
integrals of motion, and hence they should have $\mathcal{N}-1$ isometries.
This type of metrics naturally admit a $1\times\mathcal{N}-1$ decomposition,
with the $\phi$ field is canonically normalized and the rest
$\mathcal{N}-1\times\mathcal{N}-1$ matrix depending only on $\phi$. To
simplify calculations and retain some analytical control we will consider the
case of a diagonal metric

$G_{ij}=\text{diag}(1,F_{1}(\phi),F_{2}(\phi),\cdots)\,,$ (4.28)

where the metric functions $F_{1},F_{2},\cdots$ can be different. Following
the two-field discussion, in addition to the kinetic and gradient solutions we
can find de Sitter asymptotic solutions in the region where at least one of
$\kappa_{i}\equiv(\ln F_{i})_{,\phi}$ diverges and the inequalities (4.14) are
satisfied.

##  5 The general solution

###  5.1 Without matter source

For the case without matter source we find that there are general solutions
for the case where $\kappa$ and the lapse function are given by

$\kappa=-\left(\lambda+\sqrt{6}\right)\,,\qquad
N_{l}\left(t\right)=a^{-3-\sqrt{6}\lambda}\,.$ (5.1)

From the asymptotic analysis of Sec. 3 we know that this combination of
$\kappa$ and $\lambda$ will give rise to either kinetic-
($|\lambda|>\sqrt{6}$) or gradient- ($|\lambda|<\sqrt{6}$) type solutions.
Note also that for this choice of the lapse function $t$ is not the cosmic
time, but rather another time variable. For the latter selection the field
equations admit the Noetherian conservation laws

$I_{1}\left(a,\dot{a},\phi,\dot{\phi},\psi,\dot{\psi}\right)=\frac{\,\mathrm{d}}{\,\mathrm{d}t}\left(a^{2+\frac{\sqrt{6}}{2}\lambda}e^{-{1\over
2}(\lambda+\sqrt{6})\phi}\right)\,,$ (5.2)
$I_{2}\left(a,\dot{a},\phi,\dot{\phi},\psi,\dot{\psi}\right)=a^{6+\sqrt{6}\lambda}e^{-(\lambda+\sqrt{6})\phi}\dot{\psi}\,.$
(5.3)

When $\lambda=-\sqrt{6}$ the field-space curvature is zero and the solution
for this case can be found in Refs. [13, 14]. In the following we will
consider $\lambda\neq-\sqrt{6}$.

Applying the coordinate transformation

$\displaystyle
a\left(\chi\left(t\right),\xi\left(t\right),\zeta\left(t\right)\right)$
$\displaystyle=a_{0}\left(\left(\lambda+\sqrt{6}\right)^{2}\left(\chi\left(t\right)\zeta\left(t\right)-\xi\left(t\right)^{2}\right)\right)^{\left(6+\sqrt{6}\lambda\right)^{-1}}\,,$
(5.4)
$\displaystyle\phi\left(\chi\left(t\right),\xi\left(t\right),\zeta\left(t\right)\right)$
$\displaystyle=\frac{2}{\lambda+\sqrt{6}}\ln\left(\frac{\sqrt{\left(\lambda+\sqrt{6}\right)^{2}\left(\chi\left(t\right)\zeta\left(t\right)-\xi\left(t\right)^{2}\right)}}{2\chi\left(t\right)}\right)\,,$
$\displaystyle\psi\left(\chi\left(t\right),\xi\left(t\right),\zeta\left(t\right)\right)$
$\displaystyle=\frac{\xi\left(t\right)}{\chi\left(t\right)}\,,$

the point-like Lagrangian is expressed in the new coordinates as

$L\left(\chi,\dot{\chi},\xi,\dot{\xi},\zeta,\dot{\zeta}\right)=2\left(\dot{\xi}^{2}-\dot{\chi}\dot{\zeta}\right)-\bar{V}_{0}\chi^{-\bar{\lambda}}\,,$
(5.5)

where we defined for simplicity

$\bar{V}_{0}=4^{-{\lambda\over\sqrt{6}+\lambda}}V_{0}\,,\qquad\bar{\lambda}={2\lambda\over\sqrt{6}+\lambda}\,.$
(5.6)

The field equations are written as follows

$\ddot{\xi}=0\,,\qquad\ddot{\chi}=0\,,\qquad\ddot{\zeta}-\frac{\bar{\lambda}\bar{V}_{0}}{2}\chi^{\bar{\lambda}-1}=0\,,$
(5.7)

along with the constraint

$2\left(\dot{\xi}^{2}-\dot{\chi}\dot{\zeta}\right)+\bar{V}_{0}\chi^{-\bar{\lambda}}=0\,.$
(5.8)

Consequently, the analytic solution of the field equations is

$\chi\left(t\right)=\chi_{0}\left(t-t_{0}\right)\,,\qquad\xi\left(t\right)=\xi_{0}\left(t-t_{1}\right)\,,$

and

$\zeta\left(t\right)=\frac{\bar{V}_{0}\chi_{0}^{\bar{\lambda}-1}}{2\left(\bar{\lambda}+1\right)}\left(t-t_{0}\right)^{\bar{\lambda}+1}+\zeta_{0}\left(t-t_{2}\right)\,,\qquad\text{
for }\bar{\lambda}\neq-1,0\,,$ (5.9)

or

$\zeta\left(t\right)=\frac{\bar{V}_{0}}{2\chi_{0}^{2}}\ln\left(t-t_{0}\right)+\zeta_{0}\left(t-t_{2}\right)\,,\qquad\text{
for }\bar{\lambda}=-1\,,$ (5.10)

with the constraint $\xi_{0}^{2}-\chi_{0}\zeta_{0}=0$. For the special case of
$\bar{\lambda}=0$, that is the potential is a cosmological constant with
$\lambda=0$, the exact solution is

$\chi\left(t\right)=\chi_{0}\left(t-t_{0}\right)\,,\qquad\xi\left(t\right)=\xi_{0}\left(t-t_{1}\right)\,,\qquad\zeta\left(t\right)=\zeta_{0}\left(t-t_{2}\right)\,,$
(5.11)

with the constraint equation $\xi_{0}^{2}-\chi_{0}\zeta_{0}=0$. With
$\chi,\xi$ and $\zeta$ known the solution in terms of the original variables
of the problem can be found using the transformation (5.4).

In Fig. 6 we plot the qualitative evolution of the effective equation of state
parameter $w_{\rm eff}=w_{\rm eff}\left(a\right)$, which is defined as $w_{\rm
eff}\left(a\left(t\right)\right)=-1-\frac{2}{3}\frac{\dot{H}}{H^{2}}$, for two
different values of the parameter $\lambda$. We observe that the $w_{\rm
eff}\left(a\right)$ is that of the hyperbolic expansion.

Figure 6: Qualitative evolution of the effective equation of state parameter
$w_{\rm eff}=-1-\frac{2}{3}\frac{\dot{H}}{H^{2}}$ for the analytic solution of
our consideration for
$\left(\lambda,\chi_{0},\zeta_{0},t_{0},t_{1},t_{2},\bar{V}_{0}\right)$ equal
to $\left(-1,0.1,0.1,0,0.1,10,1\right)$ (left) and
$\left(\sqrt{6},1,\zeta_{0},0,0,0,1\right)$ (right).

It was pointed out recently in Ref. [46] that when the second field is phantom
we can reconstruct the analytic solution for the two-field model in a same way
as before by setting $\tilde{\psi}=i\psi$. Therefore, the analytic solution
for the second case for $\kappa=\left(\lambda+\sqrt{6}\right)$ is determined
in a similar way and we omit the presentation. Indeed by considering
$\xi(t)=i\xi(t)$ we obtain the coordinate transformation

$\displaystyle
a\left(\chi\left(t\right),\xi\left(t\right),\zeta\left(t\right)\right)$
$\displaystyle=a_{0}\left(\left(\lambda+\sqrt{6}\right)^{2}\left(\chi\left(t\right)\zeta\left(t\right)+\xi\left(t\right)^{2}\right)\right)^{\left(6+\sqrt{6}\lambda\right)^{-1}}\,,$
(5.12)
$\displaystyle\phi\left(\chi\left(t\right),\xi\left(t\right),\zeta\left(t\right)\right)$
$\displaystyle=\frac{2}{\lambda+\sqrt{6}}\ln\left(\frac{\sqrt{\left(\left(\lambda+\sqrt{6}\right)^{2}\left(\chi\left(t\right)\zeta\left(t\right)+\xi\left(t\right)^{2}\right)\right)}}{2\chi\left(t\right)}\right)\,,$
(5.13)
$\displaystyle\psi\left(\chi\left(t\right),\xi\left(t\right),\zeta\left(t\right)\right)$
$\displaystyle=\frac{\xi\left(t\right)}{\chi\left(t\right)}\,,$ (5.14)

which produces the same second-order differential equations as before.
Consequently, the solution for the variables
$\chi\left(t\right),~{}\xi\left(t\right)$ and $\zeta\left(t\right)$ is the
same as before, while the constraint equation is now
$\xi_{0}^{2}+\chi_{0}\zeta_{0}=0$.

###  5.2 In the presence of matter source

Consider now the existence of an additional matter field. We conclude that
when the equation of state parameter is
$w\left(\lambda\right)=-1+\frac{\sqrt{6}}{3}\lambda$, then the solutions found
in the previous section hold for this cosmological model as well, with the
only difference that the constraint equations for the integration constants
$\xi_{0},\chi_{0},\zeta_{0}$ are now

$\rho_{m0}=2\xi_{0}^{2}-\chi_{0}\zeta_{0}\,,$ (5.15)

for the Chiral model and

$\rho_{m0}=\xi_{0}^{2}+\chi_{0}\zeta_{0}\,,$ (5.16)

for the Chiral-quintom model.

We observe that for $w_{m}=-\frac{1}{3}$ the point-like Lagrangian (2.10)
describes the gravitational field equations for the case of a FLRW spacetime
with spatially curvature $k=\rho_{m0}$, that is, for the line element

$\,\mathrm{d}s^{2}=-N^{2}\left(t\right)\,\mathrm{d}t^{2}+\frac{a^{2}\left(t\right)}{1-\frac{k}{4}\left(X^{2}+Y^{2}+Z^{2}\right)}\left(\,\mathrm{d}X^{2}+\,\mathrm{d}Y^{2}+\,\mathrm{d}Z^{2}\right)\,.$
(5.17)

Hence, we also presented for the first time an analytic solution for a two-
field model in a non-flat FLRW background space.

##  6 Discussion

In this work we studied the dynamics and the existence of analytic solutions
for a multi-field cosmological model in a spatially flat FLRW background
space. In particular, we considered a cosmological model consisting of
$\mathcal{N}$ scalar fields minimally coupled to gravity with a hyperbolic
interaction in the kinetic terms in the presence of an ideal gas. For this
gravitational model we studied the asymptotic behaviour of the field
equations. We recovered previous results for the two-field model, that is all
critical points for the quintessence, i.e. of single-field theory, as well as
a pair of points which describe a hyperbolic solution where both fields
contribute.

In the multi-field scenario with $\mathcal{N}>2$ fields we showed that the
number of dynamical fields for any late-time solution, and most importantly
for the hyperbolic one, remains two. Therefore, at the background level the
late-time behaviour of this problem is identical to the two-field case.
Moreover, in Sec. 4.3 we showed that the same is true for first order
perturbations. This is an important observation because it is clear that by
adding additional scalar fields in this specific theory we do not get new
physical results (regarding late-time solutions), while no information can be
extracted from current observations to support this $\mathcal{N}$-field theory
with $\mathcal{N}>2$. It would be interesting to investigate whether this
holds for non-Gaussianities and other higher-order correlators as well.

Finally, for the two-field model with an exponential potential we proved for
the first time the Liouville integrability of the field equations while we
derive the analytic solution of the model. We applied the theory of similarity
transformations to construct conservation laws. We found that for a specific
combination of the exponents $\lambda$ and $\kappa$ (associated with the
potential gradient and the curvature of the hyperbolic space respectively) new
conservation laws exist which facilitate the derivation of a closed-form
solution of the dynamical system. We demonstrated that the solution holds for
the case of the presence of additional matter, while it can be used to
construct the analytic solution of the multi-field model in a non-spatially
flat FLRW background space.

## Acknowledgments

PC acknowledges financial support from the Dutch Organisation for Scientific
Research (NWO).

## References

  * [1] A. R. Brown, “Hyperbolic Inflation,” Phys. Rev. Lett. 121 (2018) no.25, 251601 [arXiv:1705.03023 [hep-th]].
  * [2] S. Mizuno and S. Mukohyama, “Primordial perturbations from inflation with a hyperbolic field-space,” Phys. Rev. D 96 (2017) no.10, 103533 [arXiv:1707.05125 [hep-th]].
  * [3] M. Cicoli, G. Dibitetto and F. G. Pedro, “New accelerating solutions in late-time cosmology,” Phys. Rev. D 101 (2020) no.10, 103524 [arXiv:2002.02695 [gr-qc]].
  * [4] M. Cicoli, G. Dibitetto and F. G. Pedro, “Out of the Swampland with Multifield Quintessence?,” JHEP 10 (2020), 035 [arXiv:2007.11011 [hep-th]].
  * [5] J. Fumagalli, S. Garcia-Saenz, L. Pinol, S. Renaux-Petel and J. Ronayne, “Hyper-Non-Gaussianities in Inflation with Strongly Nongeodesic Motion,” Phys. Rev. Lett. 123 (2019) no.20, 201302 [arXiv:1902.03221 [hep-th]].
  * [6] T. Bjorkmo, R. Z. Ferreira and M. C. D. Marsh, “Mild Non-Gaussianities under Perturbative Control from Rapid-Turn Inflation Models,” JCAP 12 (2019), 036 [arXiv:1908.11316 [hep-th]].
  * [7] R. Z. Ferreira, “Non-Gaussianities in models of inflation with large and negative entropic masses,” JCAP 08 (2020), 034 [arXiv:2003.13410 [astro-ph.CO]].
  * [8] M. Bounakis, I. G. Moss and G. Rigopoulos, “Observational constraints on Hyperinflation,” [arXiv:2010.06461 [gr-qc]].
  * [9] S. Basilakos, G. Leon, G. Papagiannopoulos and E. N. Saridakis, “Dynamical system analysis at background and perturbation levels: Quintessence in severe disadvantage comparing to $\Lambda$CDM,” Phys. Rev. D 100 (2019) no.4, 043524 [arXiv:1904.01563 [gr-qc]].
  * [10] A. Banerjee, H. Cai, L. Heisenberg, E. Ó. Colgáin, M. M. Sheikh-Jabbari and T. Yang, [arXiv:2006.00244 [astro-ph.CO]].
  * [11] T. Bjorkmo and M. C. D. Marsh, “Hyperinflation generalised: from its attractor mechanism to its tension with the ‘swampland conditions’,” JHEP 04 (2019), 172 [arXiv:1901.08603 [hep-th]].
  * [12] V. Aragam, S. Paban and R. Rosati, “The Multi-Field, Rapid-Turn Inflationary Solution,” [arXiv:2010.15933 [hep-th]].
  * [13] L. P. Chimento, “General solution to two-scalar field cosmologies with exponential potentials,” Class. Quant. Grav. 15 (1998), 965-974
  * [14] P. Christodoulidis, “Probing the inflationary evolution using analytical solutions,” [arXiv:1811.06456 [astro-ph.CO]].
  * [15] A. Paliathanasis, G. Leon and S. Pan, “Exact Solutions in Chiral Cosmology,” Gen. Rel. Grav. 51 (2019) no.9, 106 [arXiv:1811.10038 [gr-qc]].
  * [16] J. Socorro, S. Pérez-Payán, R. Hernández, A. Espinoza-García and L. R. Díaz-Barrón, “Classical and quantum exact solutions for a FRW in chiral like cosmology,” [arXiv:2012.11108 [gr-qc]].
  * [17] J. Socorro, S. Pérez-Payán, A. Espinoza-García and L. R. Díaz-Barrón, [arXiv:2101.05973 [gr-qc]].
  * [18] S. Basilakos, M. Tsamparlis and A. Paliathanasis, “Using the Noether symmetry approach to probe the nature of dark energy,” Phys. Rev. D 83 (2011), 103512 [arXiv:1104.2980 [astro-ph.CO]].
  * [19] A. Paliathanasis, M. Tsamparlis and S. Basilakos, “Constraints and analytical solutions of $f(R)$ theories of gravity using Noether symmetries,” Phys. Rev. D 84 (2011), 123514 [arXiv:1111.4547 [astro-ph.CO]].
  * [20] J. A. Belinchón, T. Harko and M. K. Mak, Astrophys. Space Sci. 361 (2016) no.2, 52 doi:10.1007/s10509-015-2642-7 [arXiv:1512.08054 [gr-qc]].
  * [21] M. Demianski, E. Piedipalumbo, C. Rubano and C. Tortora, “Accelerating universe in scalar tensor models: Confrontation of theoretical predictions with observations,” Astron. Astrophys. 454 (2006), 55-66 [arXiv:astro-ph/0604026 [astro-ph]].
  * [22] P. A. Terzis, N. Dimakis and T. Christodoulakis, “Noether analysis of Scalar-Tensor Cosmology,” Phys. Rev. D 90 (2014) no.12, 123543[arXiv:1410.0802 [gr-qc]].
  * [23] A. Paliathanasis, J. D. Barrow and P. G. L. Leach, “Cosmological Solutions of $f(T)$ Gravity,” Phys. Rev. D 94 (2016) no.2, 023525 [arXiv:1606.00659 [gr-qc]].
  * [24] J. D. Barrow, S. Cotsakis and A. Tsokaros, “Series expansions and sudden singularities,” d[arXiv:1301.6523 [gr-qc]].
  * [25] S. Cotsakis, “Asymptotic Poincaré compactification and finite-time singularities,” Grav. Cosmol. 19 (2013), 240-245 [arXiv:1301.4778 [gr-qc]].
  * [26] L. Amendola, R. Gannouji, D. Polarski and S. Tsujikawa, “Conditions for the cosmological viability of f(R) dark energy models,” Phys. Rev. D 75 (2007), 083504 [arXiv:gr-qc/0612180 [gr-qc]].
  * [27] A. Coley and G. Leon, “Static Spherically Symmetric Einstein-aether models I: Perfect fluids with a linear equation of state and scalar fields with an exponential self-interacting potential,” Gen. Rel. Grav. 51 (2019) no.9, 115 [arXiv:1905.02003 [gr-qc]].
  * [28] G. Leon, Y. Leyva and J. Socorro, “Quintom phase-space: beyond the exponential potential,” Phys. Lett. B 732 (2014), 285-297 [arXiv:1208.0061 [gr-qc]].
  * [29] R. Lazkoz, G. Leon and I. Quiros, “Quintom cosmologies with arbitrary potentials,” Phys. Lett. B 649 (2007), 103-110 [arXiv:astro-ph/0701353 [astro-ph]].
  * [30] Y.F. Cai, E.N. Saridakis, M.R. Setare and J.-Q. Xia, “Quintom cosmology: theoretical implications and observations,” Phys. Rep. 493, 1 (2010)
  * [31] S.V. Chervon, “On the chiral model of cosmological inflation,” Russ. Phys. J. 38, 539 (1995)
  * [32] S. V. Ketov, “Quantum Non-linear Sigma Models,”, Springer-Verlag, Berlin, (2000).
  * [33] J. Lee, T.H. Lee, T. Moon and P. Oh, “de Sitter nonlinear sigma model and accelerating universe,” Phys. Rev. D 80, 065016 (2009)
  * [34] S.V. Chervon, “Chiral Cosmological Models: Dark Sector Fields Description,” Quantum Matter 2, 71 (2013)
  * [35] A. Paliathanasis, “ Dynamics of chiral cosmology”, Class. Quantum Grav. 37, 195014 (2020)
  * [36] P. Christodoulidis, D. Roest and E. I. Sfakianakis, “Scaling attractors in multi-field inflation,” JCAP 12, 059 (2019) [arXiv:1903.06116 [hep-th]].
  * [37] N. Dimakis, A. Paliathanasis, P. A. Terzis and T. Christodoulakis, Eur. Phys. J. C 79, no.7, 618 (2019) [arXiv:1904.09713 [gr-qc]].
  * [38] A. Paliathanasis and G. Leon, “Asymptotic behavior of N-fields Chiral Cosmology.” Eur. Phys. J. C 80 (2020) no.9, 847 [arXiv:2007.13223 [gr-qc]].
  * [39] C. Gordon, D. Wands, B. A. Bassett and R. Maartens, “Adiabatic and entropy perturbations from inflation,” Phys. Rev. D 63 (2000), 023506 [arXiv:astro-ph/0009131 [astro-ph]].
  * [40] K. A. Malik and D. Wands, “Cosmological perturbations,” Phys. Rept. 475 (2009), 1-51 [arXiv:0809.4944 [astro-ph]].
  * [41] D. I. Kaiser, E. A. Mazenc and E. I. Sfakianakis, “Primordial Bispectrum from Multifield Inflation with Nonminimal Couplings,” Phys. Rev. D 87 (2013), 064004 [arXiv:1210.7487 [astro-ph.CO]].
  * [42] D. I. Kaiser and E. I. Sfakianakis, “Multifield Inflation after Planck: The Case for Nonminimal Couplings,” Phys. Rev. Lett. 112 (2014) no.1, 011302 [arXiv:1304.0363 [astro-ph.CO]].
  * [43] V. Aragam, S. Paban and R. Rosati, “Multi-field Inflation in High-Slope Potentials,” JCAP 04 (2020), 022 [arXiv:1905.07495 [hep-th]].
  * [44] A. Achúcarro, S. Céspedes, A. C. Davis and G. A. Palma, “Constraints on Holographic Multifield Inflation and Models Based on the Hamilton-Jacobi Formalism,” Phys. Rev. Lett. 122 (2019) no.19, 191301 [arXiv:1809.05341 [hep-th]].
  * [45] L. Pinol, “Multifield inflation beyond $N_{\mathrm{field}}=2$: non-Gaussianities and single-field effective theory,” [arXiv:2011.05930 [astro-ph.CO]].
  * [46] A. Paliathanasis and G. Leon, Dynamics of a two scalar field cosmological model with phantom terms, [arXiv:2009.12874]


* Cleveland et al. (1992) Cleveland, W. S., Grosse, E., & Shyu, W. M. 1992, in Statistical Models in S, ed. J. Chambers & T. Hastie (Wadsworth & Brooks/Cole).  https://www.bibsonomy.org/bibtex/249dc786266d86d0cc54ed6ccb2fe9894/gron
  * Coffaro et al. (2022) Coffaro, M., Stelzer, B., & Orlando, S. 2022, A&A, 661, A79, doi: 10.1051/0004-6361/202142298
  * Coffaro et al. (2020) Coffaro, M., Stelzer, B., Orlando, S., et al. 2020, A&A, 636, A49, doi: 10.1051/0004-6361/201936479
  * Cohen et al. (2017) Cohen, O., Yadav, R., Garraffo, C., et al. 2017, ApJ, 834, 14, doi: 10.3847/1538-4357/834/1/14
  * Donati et al. (1997) Donati, J. F., Semel, M., Carter, B. D., Rees, D. E., & Collier Cameron, A. 1997, MNRAS, 291, 658, doi: 10.1093/mnras/291.4.658
  * Donati et al. (2007) Donati, J. F., Jardine, M. M., Gregory, S. G., et al. 2007, MNRAS, 380, 1297, doi: 10.1111/j.1365-2966.2007.12194.x
  * Drake et al. (2014) Drake, J. J., Braithwaite, J., Kashyap, V., Günther, H. M., & Wright, N. J. 2014, ApJ, 786, 136, doi: 10.1088/0004-637X/786/2/136
  * Drake et al. (2009) Drake, J. J., Ercolano, B., Flaccomio, E., & Micela, G. 2009, ApJ, 699, L35, doi: 10.1088/0004-637X/699/1/L35
  * Dupuy et al. (2018) Dupuy, R., Bertin, M., Féraud, G., et al. 2018, Nature Astronomy, 2, 796, doi: 10.1038/s41550-018-0532-y
  * Durney et al. (1993) Durney, B. R., De Young, D. S., & Roxburgh, I. W. 1993, Sol. Phys., 145, 207, doi: 10.1007/BF00690652
  * Ercolano & Pascucci (2017) Ercolano, B., & Pascucci, I. 2017, Royal Society Open Science, 4, 170114, doi: 10.1098/rsos.170114
  * Evans et al. (2018) Evans, D. W., Riello, M., De Angeli, F., et al. 2018, A&A, 616, A4, doi: 10.1051/0004-6361/201832756
  * Favata et al. (2005) Favata, F., Flaccomio, E., Reale, F., et al. 2005, ApJS, 160, 469, doi: 10.1086/432542
  * Feigelson (1996) Feigelson, E. D. 1996, ApJ, 468, 306, doi: 10.1086/177691
  * Feigelson (2010) —. 2010, Proceedings of the National Academy of Science, 107, 7153, doi: 10.1073/pnas.0913952107
  * Feigelson et al. (2002) Feigelson, E. D., Broos, P., Gaffney, James A., I., et al. 2002, ApJ, 574, 258, doi: 10.1086/340936
  * Feigelson & Montmerle (1999) Feigelson, E. D., & Montmerle, T. 1999, ARA&A, 37, 363, doi: 10.1146/annurev.astro.37.1.363
  * Feigelson et al. (2004) Feigelson, E. D., Hornschemeier, A. E., Micela, G., et al. 2004, ApJ, 611, 1107, doi: 10.1086/422248
  * Feigelson et al. (2013) Feigelson, E. D., Townsley, L. K., Broos, P. S., et al. 2013, ApJS, 209, 26, doi: 10.1088/0067-0049/209/2/26
  * Flaccomio et al. (2018) Flaccomio, E., Micela, G., Sciortino, S., et al. 2018, A&A, 620, A55, doi: 10.1051/0004-6361/201833308
  * Flaischlen et al. (2021) Flaischlen, S., Preibisch, T., Manara, C. F., & Ercolano, B. 2021, arXiv e-prints, arXiv:2103.03039.  https://arxiv.org/abs/2103.03039
  * Fromang et al. (2002) Fromang, S., Terquem, C., & Balbus, S. A. 2002, MNRAS, 329, 18, doi: 10.1046/j.1365-8711.2002.04940.x
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1, doi: 10.1051/0004-6361/201629272
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1, doi: 10.1051/0004-6361/202039657
  * Gallet & Bouvier (2015) Gallet, F., & Bouvier, J. 2015, A&A, 577, A98, doi: 10.1051/0004-6361/201525660
  * García Muñoz et al. (2021) García Muñoz, A., Fossati, L., Youngblood, A., et al. 2021, ApJ, 907, L36, doi: 10.3847/2041-8213/abd9b8
  * Garmire et al. (2003) Garmire, G. P., Bautz, M. W., Ford, P. G., Nousek, J. A., & Ricker, George R., J. 2003, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 4851, X-Ray and Gamma-Ray Telescopes and Instruments for Astronomy., ed. J. E. Truemper & H. D. Tananbaum, 28–44, doi: 10.1117/12.461599
  * Gehrels (1986) Gehrels, N. 1986, ApJ, 303, 336, doi: 10.1086/164079
  * Getman et al. (2017) Getman, K. V., Broos, P. S., Kuhn, M. A., et al. 2017, ApJS, 229, 28, doi: 10.3847/1538-4365/229/2/28
  * Getman & Feigelson (2021) Getman, K. V., & Feigelson, E. D. 2021, ApJ, 916, 32, doi: 10.3847/1538-4357/ac00be
  * Getman et al. (2010) Getman, K. V., Feigelson, E. D., Broos, P. S., Townsley, L. K., & Garmire, G. P. 2010, ApJ, 708, 1760, doi: 10.1088/0004-637X/708/2/1760
  * Getman et al. (2021) Getman, K. V., Feigelson, E. D., & Garmire, G. P. 2021, ApJ, 920, 154, doi: 10.3847/1538-4357/ac1746
  * Getman et al. (2014a) Getman, K. V., Feigelson, E. D., & Kuhn, M. A. 2014a, ApJ, 787, 109, doi: 10.1088/0004-637X/787/2/109
  * Getman et al. (2018a) Getman, K. V., Feigelson, E. D., Kuhn, M. A., et al. 2018a, MNRAS, 476, 1213, doi: 10.1093/mnras/sty302
  * Getman et al. (2019) Getman, K. V., Feigelson, E. D., Kuhn, M. A., & Garmire, G. P. 2019, MNRAS, 487, 2977, doi: 10.1093/mnras/stz1457
  * Getman et al. (2008) Getman, K. V., Feigelson, E. D., Micela, G., et al. 2008, ApJ, 688, 437, doi: 10.1086/592034
  * Getman et al. (2018b) Getman, K. V., Kuhn, M. A., Feigelson, E. D., et al. 2018b, MNRAS, 477, 298, doi: 10.1093/mnras/sty473
  * Getman et al. (2005) Getman, K. V., Flaccomio, E., Broos, P. S., et al. 2005, ApJS, 160, 319, doi: 10.1086/432092
  * Getman et al. (2014b) Getman, K. V., Feigelson, E. D., Kuhn, M. A., et al. 2014b, ApJ, 787, 108, doi: 10.1088/0004-637X/787/2/108
  * Glassgold et al. (2000) Glassgold, A. E., Feigelson, E. D., & Montmerle, T. 2000, in Protostars and Planets IV, ed. V. Mannings, A. P. Boss, & S. S. Russell, 429
  * Glocer et al. (2012) Glocer, A., Kitamura, N., Toth, G., & Gombosi, T. 2012, Journal of Geophysical Research (Space Physics), 117, A04318, doi: 10.1029/2011JA017136
  * Gregory et al. (2016) Gregory, S. G., Adams, F. C., & Davies, C. L. 2016, MNRAS, 457, 3836, doi: 10.1093/mnras/stw259
  * Gregory et al. (2014) Gregory, S. G., Donati, J. F., Morin, J., et al. 2014, in Magnetic Fields throughout Stellar Evolution, ed. P. Petit, M. Jardine, & H. C. Spruit, Vol. 302, 40–43, doi: 10.1017/S1743921314001677
  * Gregory et al. (2010) Gregory, S. G., Jardine, M., Gray, C. G., & Donati, J. F. 2010, Reports on Progress in Physics, 73, 126901, doi: 10.1088/0034-4885/73/12/126901
  * Gressel et al. (2013) Gressel, O., Nelson, R. P., Turner, N. J., & Ziegler, U. 2013, ApJ, 779, 59, doi: 10.1088/0004-637X/779/1/59
  * Gronoff et al. (2020) Gronoff, G., Arras, P., Baraka, S., et al. 2020, Journal of Geophysical Research (Space Physics), 125, e27639, doi: 10.1029/2019JA027639
  * Güdel (2004) Güdel, M. 2004, A&A Rev., 12, 71, doi: 10.1007/s00159-004-0023-2
  * Güdel et al. (1997) Güdel, M., Guinan, E. F., & Skinner, S. L. 1997, ApJ, 483, 947, doi: 10.1086/304264
  * Güdel & Nazé (2009) Güdel, M., & Nazé, Y. 2009, A&A Rev., 17, 309, doi: 10.1007/s00159-009-0022-4
  * Güdel et al. (2007) Güdel, M., Briggs, K. R., Arzner, K., et al. 2007, A&A, 468, 353, doi: 10.1051/0004-6361:20065724
  * Henderson & Stassun (2012) Henderson, C. B., & Stassun, K. G. 2012, ApJ, 747, 51, doi: 10.1088/0004-637X/747/1/51
  * Hu et al. (2022) Hu, J., Airapetian, V. S., Li, G., Zank, G., & Jin, M. 2022, Science Advances in press.
  * Iben (1965) Iben, Icko, J. 1965, ApJ, 141, 993, doi: 10.1086/148193
  * Ilgner & Nelson (2006) Ilgner, M., & Nelson, R. P. 2006, A&A, 455, 731, doi: 10.1051/0004-6361:20065308
  * Jardine & Unruh (1999) Jardine, M., & Unruh, Y. C. 1999, A&A, 346, 883
  * Johnstone (2020) Johnstone, C. P. 2020, ApJ, 890, 79, doi: 10.3847/1538-4357/ab6224
  * Johnstone et al. (2021a) Johnstone, C. P., Bartel, M., & Güdel, M. 2021a, A&A, 649, A96, doi: 10.1051/0004-6361/202038407
  * Johnstone et al. (2021b) Johnstone, C. P., Lammer, H., Kislyakova, K. G., Scherf, M., & Güdel, M. 2021b, Earth and Planetary Science Letters, 576, 117197, doi: 10.1016/j.epsl.2021.117197
  * Kaplan & Meier (1958) Kaplan, E. L., & Meier, P. 1958, Journal of the American Statistical Association, 53, 457, doi: 10.1080/01621459.1958.10501452
  * Khodachenko et al. (2021) Khodachenko, M. L., Shaikhislamov, I. F., Lammer, H., et al. 2021, MNRAS, 507, 3626, doi: 10.1093/mnras/stab2366
  * Kite & Barnett (2020) Kite, E. S., & Barnett, M. N. 2020, Proceedings of the National Academy of Science, 117, 18264.  https://arxiv.org/abs/2006.02589
  * Kuhn et al. (2013) Kuhn, M. A., Getman, K. V., Broos, P. S., Townsley, L. K., & Feigelson, E. D. 2013, ApJS, 209, 27, doi: 10.1088/0067-0049/209/2/27
  * Kuhn et al. (2015) Kuhn, M. A., Getman, K. V., & Feigelson, E. D. 2015, ApJ, 802, 60, doi: 10.1088/0004-637X/802/1/60
  * Kuhn et al. (2019) Kuhn, M. A., Hillenbrand, L. A., Sills, A., Feigelson, E. D., & Getman, K. V. 2019, ApJ, 870, 32, doi: 10.3847/1538-4357/aaef8c
  * Lammer et al. (2003) Lammer, H., Selsis, F., Ribas, I., et al. 2003, ApJ, 598, L121, doi: 10.1086/380815
  * Lammer et al. (2014) Lammer, H., Stökl, A., Erkaev, N. V., et al. 2014, MNRAS, 439, 3225, doi: 10.1093/mnras/stu085
  * Loader (1999) Loader, C. 1999, Local Regression and Likelihood (Springer-Verlag New York), doi: 10.1007/b98858
  * Loader (2020) Loader, C. 2020, locfit: Local Regression, Likelihood and Density Estimation, codes.  https://cran.r-project.org/web/packages/locfit/index.html
  * Luhman (2018) Luhman, K. L. 2018, AJ, 156, 271, doi: 10.3847/1538-3881/aae831
  * Luhman & Esplin (2020) Luhman, K. L., & Esplin, T. L. 2020, AJ, 160, 44, doi: 10.3847/1538-3881/ab9599
  * Magaudda et al. (2021) Magaudda, E., Stelzer, B., Raetz, S., & A, K. 2021, arXiv e-prints, arXiv:2106.14548.  https://arxiv.org/abs/2106.14548
  * Maggio et al. (1987) Maggio, A., Sciortino, S., Vaiana, G. S., et al. 1987, ApJ, 315, 687, doi: 10.1086/165170
  * Maschberger (2013) Maschberger, T. 2013, MNRAS, 429, 1725, doi: 10.1093/mnras/sts479
  * Massol et al. (2016) Massol, H., Hamano, K., Tian, F., et al. 2016, Space Sci. Rev., 205, 153, doi: 10.1007/s11214-016-0280-1
  * Ng & Maechler (2007) Ng, P., & Maechler, M. 2007, Statistical Modelling, 7, 315.  http://smj.sagepub.com/content/7/4/315.abstract
  * Ng & Maechler (2020) Ng, P. T., & Maechler, M. 2020, COBS – Constrained B-splines (Sparse matrix based).  https://CRAN.R-project.org/package=cobs
  * Nuñez et al. (2021) Nuñez, E. H., Povich, M. S., Binder, B. A., Townsley, L. K., & Broos, P. S. 2021, AJ, 162, 153, doi: 10.3847/1538-3881/ac0af8
  * Ortenzi et al. (2020) Ortenzi, G., Noack, L., Sohl, F., et al. 2020, Scientific Reports, 10, 10907, doi: 10.1038/s41598-020-67751-7
  * Owen (2019) Owen, J. E. 2019, Annual Review of Earth and Planetary Sciences, 47, 67, doi: 10.1146/annurev-earth-053018-060246
  * Owen et al. (2012) Owen, J. E., Clarke, C. J., & Ercolano, B. 2012, MNRAS, 422, 1880, doi: 10.1111/j.1365-2966.2011.20337.x
  * Owen et al. (2013) Owen, J. E., Hudoba de Badyn, M., Clarke, C. J., & Robins, L. 2013, MNRAS, 436, 1430, doi: 10.1093/mnras/stt1663
  * Picogna et al. (2019) Picogna, G., Ercolano, B., Owen, J. E., & Weber, M. L. 2019, MNRAS, 487, 691, doi: 10.1093/mnras/stz1166
  * Pizzolato et al. (2003) Pizzolato, N., Maggio, A., Micela, G., Sciortino, S., & Ventura, P. 2003, A&A, 397, 147, doi: 10.1051/0004-6361:20021560
  * Preibisch (2012) Preibisch, T. 2012, Research in Astronomy and Astrophysics, 12, 1, doi: 10.1088/1674-4527/12/1/001
  * Preibisch & Feigelson (2005) Preibisch, T., & Feigelson, E. D. 2005, ApJS, 160, 390, doi: 10.1086/432094
  * Preibisch et al. (2017) Preibisch, T., Flaischlen, S., Gaczkowski, B., Townsley, L., & Broos, P. 2017, A&A, 605, A85, doi: 10.1051/0004-6361/201730874
  * Preibisch et al. (2005) Preibisch, T., Kim, Y.-C., Favata, F., et al. 2005, ApJS, 160, 401, doi: 10.1086/432891
  * R Core Team (2020) R Core Team. 2020, R: A language and environment for statistical computing, R Foundation for Statistical Computing, Vienna, Austria.  https://www.R-project.org
  * Rab et al. (2017) Rab, C., Güdel, M., Padovani, M., et al. 2017, A&A, 603, A96, doi: 10.1051/0004-6361/201630241
  * Rab et al. (2020) Rab, C., Padovani, M., Güdel, M., et al. 2020, in Origins: From the Protosun to the First Steps of Life, ed. B. G. Elmegreen, L. V. Tóth, & M. Güdel, Vol. 345, 310–311, doi: 10.1017/S174392131900156X
  * Reggiani et al. (2011) Reggiani, M., Robberto, M., Da Rio, N., et al. 2011, A&A, 534, A83, doi: 10.1051/0004-6361/201116946
  * Richert et al. (2018) Richert, A. J. W., Getman, K. V., Feigelson, E. D., et al. 2018, MNRAS, 477, 5191, doi: 10.1093/mnras/sty949
  * Riello et al. (2021) Riello, M., De Angeli, F., Evans, D. W., et al. 2021, A&A, 649, A3, doi: 10.1051/0004-6361/202039587
  * Sanz-Forcada et al. (2003) Sanz-Forcada, J., Brickhouse, N. S., & Dupree, A. K. 2003, ApJS, 145, 147, doi: 10.1086/345815
  * Shaikhislamov et al. (2020) Shaikhislamov, I. F., Fossati, L., Khodachenko, M. L., et al. 2020, A&A, 639, A109, doi: 10.1051/0004-6361/202038363
  * Shang et al. (2002) Shang, H., Glassgold, A. E., Shu, F. H., & Lizano, S. 2002, ApJ, 564, 853, doi: 10.1086/324197
  * Sheather (2009) Sheather, S. 2009, A Modern Approach to Regression with R (Springer New York), doi: 10.1007/978-0-387-09608-7
  * Siess et al. (2000) Siess, L., Dufour, E., & Forestini, M. 2000, A&A, 358, 593.  https://arxiv.org/abs/astro-ph/0003477
  * Skumanich (1972a) Skumanich, A. 1972a, ApJ, 171, 565, doi: 10.1086/151310
  * Skumanich (1972b) —. 1972b, ApJ, 171, 565, doi: 10.1086/151310
  * Smith et al. (2001) Smith, R. K., Brickhouse, N. S., Liedahl, D. A., & Raymond, J. C. 2001, ApJ, 556, L91, doi: 10.1086/322992
  * Stelzer (2017) Stelzer, B. 2017, Astronomische Nachrichten, 338, 195, doi: 10.1002/asna.201713330
  * Stelzer et al. (2005) Stelzer, B., Flaccomio, E., Montmerle, T., et al. 2005, ApJS, 160, 557, doi: 10.1086/432375
  * Stelzer et al. (2009) Stelzer, B., Robrade, J., Schmitt, J. H. M. M., & Bouvier, J. 2009, A&A, 493, 1109, doi: 10.1051/0004-6361:200810540
  * Telleschi et al. (2007) Telleschi, A., Güdel, M., Briggs, K. R., Audard, M., & Palla, F. 2007, A&A, 468, 425, doi: 10.1051/0004-6361:20066565
  * Therneau (2020) Therneau, T. M. 2020, A Package for Survival Analysis in R.  https://CRAN.R-project.org/package=survival
  * Townsley et al. (2019) Townsley, L. K., Broos, P. S., Garmire, G. P., & Povich, M. S. 2019, ApJS, 244, 28, doi: 10.3847/1538-4365/ab345b
  * Tu et al. (2015) Tu, L., Johnstone, C. P., Güdel, M., & Lammer, H. 2015, A&A, 577, L3, doi: 10.1051/0004-6361/201526146
  * Vidotto et al. (2014) Vidotto, A. A., Gregory, S. G., Jardine, M., et al. 2014, MNRAS, 441, 2361, doi: 10.1093/mnras/stu728
  * Vilhu (1984) Vilhu, O. 1984, A&A, 133, 117
  * Waggoner & Cleeves (2019) Waggoner, A. R., & Cleeves, L. I. 2019, ApJ, 883, 197, doi: 10.3847/1538-4357/ab3d38
  * Waggoner & Cleeves (2022) —. 2022, arXiv e-prints, arXiv:2202.06962.  https://arxiv.org/abs/2202.06962
  * Warnecke & Käpylä (2020) Warnecke, J., & Käpylä, M. J. 2020, A&A, 642, A66, doi: 10.1051/0004-6361/201936922
  * Williams & Cieza (2011) Williams, J. P., & Cieza, L. A. 2011, ARA&A, 49, 67, doi: 10.1146/annurev-astro-081710-102548
  * Wilms et al. (2000) Wilms, J., Allen, A., & McCray, R. 2000, ApJ, 542, 914, doi: 10.1086/317016
  * Wright et al. (2011a) Wright, N. J., Drake, J. J., Mamajek, E. E., & Henry, G. W. 2011a, ApJ, 743, 48, doi: 10.1088/0004-637X/743/1/48
  * Wright et al. (2011b) —. 2011b, ApJ, 743, 48, doi: 10.1088/0004-637X/743/1/48
  * Yadav et al. (2015) Yadav, R. K., Christensen, U. R., Morin, J., et al. 2015, ApJ, 813, L31, doi: 10.1088/2041-8205/813/2/L31
  * Zhang et al. (2022) Zhang, M., Knutson, H. A., Wang, L., et al. 2022, AJ, 163, 68, doi: 10.3847/1538-3881/ac3f3b
  * Zhu et al. (2017) Zhu, H., Tian, W., Li, A., & Zhang, M. 2017, MNRAS, 471, 3494, doi: 10.1093/mnras/stx1580


# High-dimensional and Permutation Invariant Anomaly Detection

Vinicius Mikuni<EMAIL_ADDRESS>National Energy Research Scientific Computing
Center, Berkeley Lab, Berkeley, CA 94720, USA  Benjamin Nachman
<EMAIL_ADDRESS>Physics Division, Lawrence Berkeley National Laboratory,
Berkeley, CA 94720, USA  Berkeley Institute for Data Science, University of
California, Berkeley, CA 94720, USA

###### Abstract

Methods for anomaly detection of new physics processes are often limited to
low-dimensional spaces due to the difficulty of learning high-dimensional
probability densities. Particularly at the constituent level, incorporating
desirable properties such as permutation invariance and variable-length inputs
becomes difficult within popular density estimation methods. In this work, we
introduce a permutation-invariant density estimator for particle physics data
based on diffusion models, specifically designed to handle variable-length
inputs. We demonstrate the efficacy of our methodology by utilizing the
learned density as a permutation-invariant anomaly detection score,
effectively identifying jets with low likelihood under the background-only
hypothesis. To validate our density estimation method, we investigate the
ratio of learned densities and compare to those obtained by a supervised
classification algorithm.

##  I Introduction

Anomaly detection (AD) has emerged as a complementary strategy to classical
model-dependent searches for new particles at the Large Hadron Collider and
elsewhere. These tools are motivated by the current lack of excesses and the
vast parameter space of possibilities [1, 2]. Machine learning (ML) techniques
are addressing these motivations and also allowing for complex particle
physics data to be probed holistically in their natural high dimensionality
[3].

Nearly all searches for new particles begin by positing a particular signal
model, simulating the signal and relevant Standard Model (SM) backgrounds, and
then training (with or without ML) a classifier to distinguish the signal and
background simulations. Machine learning–based AD tries to assume as little as
possible about the signal while also maintaining the ability to estimate the
SM background. Two main classes of ML approaches are unsupervised and
weakly/semi-supervised. Unsupervised methods use ‘no’ information about the
signal in training while weakly/semi-supervised methods use limited or noisy
labels. The ‘no’ is in quotes because there is often implicit signal
information used through event and feature selection.

At their core, unsupervised methods select events that are rare, while
weakly/semi supervised methods focus on events that have a high likelihood
ratio with respect to some reference(s). The first ML-based AD proposals in
high energy physics explored both weakly/semi-supervised classifiers [4, 5, 6]
as well as unsupervised learning via a type of ML tool called an autoencoder
[7, 8, 9]. Since that time, there have been many proposals in the literature
(see e.g. Ref. [10]), community challenges comparing a large number of
approaches [11, 12], and first physics results using a variety of methods [13,
14, 15, 16, 17, 18]. Even though a number of weakly supervised methods have
statistical guarantees of optimality that unsupervised method lack [19, 20],
there has been significant interest in unsupervised AD because of its
flexibility.

The flexibility of unsupervised learning leads to a number of challenges.
There is no unique way to estimate the probability density of a given dataset,
with some methods offering only an implicit approximation through proxy
quantities like the reconstruction fidelity of compression algorithms. The
probability density itself is not invariant under coordinate transformations,
so the selected rare events will depend on the feature selection [21]. Even
though particle physics data are often described by high- (and
variable-)dimensional, permutation-invariant sets (‘point clouds’), there has
not yet been a proposal to use explicit density estimation techniques for AD
that account for all of these properties. Implicit density estimation has been
studied with a variety of high-dimensional, but mostly fixed-length
representations, such as (variational) autoencoders and related approaches [8,
9, 22, 23, 24, 25, 26, 27]. Since our validation protocol requires access to
the density, we focus only on explicit methods. So far, the only111Except for
Ref. [28, 29], which discretize the phase space and turn the problem into a
multi-class classification task. high-dimensional explicit density estimators
in particle physics [30, 31, 32, 33, 34, 35] have been based on normalizing
flows [36, 37]. These works process fixed-length and ordered inputs, but
recent work has shown with higher-level observables how to accommodate
variable-length and permtuation invariance with normalizing flows [38].

However, variable-length is not a natural property for normalizing flows which
are built on bijective maps from the data space to a fixed-length latent
space. In contrast, a newer class of methods called score-matching or
diffusion models do not have this restriction. These techniques estimate the
gradient of the density instead of the density itself, and therefore have
fewer restrictions than normalizing flows. Diffusion models have been shown to
accurately model both high- [39] and/or variable- [40, 41, 42, 43] dimensional
feature spaces. Despite these early successes, such models have not been used
yet for explicit density estimation in particle physics.

We propose to use point cloud diffusion models combined with explicit density
estimation for AD. Our approach is based on Ref. [41], and inherits the
ability to process variable-length and permutation-invariant sets. From the
learned score function, we estimate the data density and provide results for
two different diffusion models; one trained with standard score-matching
objective and one trained using maximum likelihood estimation. Since the true
density is not known, we quantify the performance of the density estimation
with likelihood ratios. Finally, we demonstrate the performance of the density
as an anomaly score for top quark jets as well as jets produced from dark
showers in a hidden valley model. Other tasks that require access to the data
density could also benefit from our method.

This paper is organized as follows. Section II introduces the methodology of
maximum likelihood-based diffusion modeling for permutation-invariant density
estimation. The datasets used for our numerical examples are presented in Sec.
III and the results themselves appear in Sec. IV. The paper ends with
conclusions and outlook in Sec. V.

##  II Score Matching and Maximum Likelihood Training of Diffusion Models

Score-based generative models are a class of generative algorithms that aim to
generate data by learning the score function, or gradients of the logarithm of
the probability density of the data. The training strategy presented in Ref.
[44] introduces the idea of denoising score-matching, where data can be
perturbed by a smearing function and matching the score of the smeared data is
equivalent to matching the score of the smearing function Ref. [45]. Given
some high-dimensional distribution $\mathbf{x}\in\mathbb{R}^{\mathrm{D}}$, the
score function we want to approximate, $\mathbf{\nabla_{x}}\log
p_{\text{data}}$, with $\mathbf{x}\sim p_{\text{data}}$, is obtained by
minimizing the following quantity

$\frac{1}{2}\mathbb{E}_{t}\mathbb{E}_{p_{t}(\mathbf{x})}\left[\lambda(t)\left\|\mathbf{s_{\theta}(\mathbf{x}_{t}},t\mathbf{)}-\mathbf{\nabla_{\mathbf{x}_{t}}}\log
p_{t}(\mathbf{x}_{t}|x_{0})\right\|^{2}_{2}\right].$ (1)

The goal of a neural network $\mathbf{s_{\theta}(\mathbf{x}_{t}},t\mathbf{)}$
with trainable parameters $\theta$ and evaluated with data $\mathbf{x}_{t}$
that have been perturbed at time $t$ is to give a time-dependent approximation
of the score function. The time dependence of the score function is introduced
to address the different levels of perturbation used in each time step. At
times near $0$, at the beginning of the diffusion process
($\mathbf{x}(0):=\mathbf{x}_{0}:=\mathbf{x}$), the smearing applied to data is
small, gradually increasing as time increases and ensures that at longer time
scales the distribution is completely overridden by noise. Similarly, the
positive weighing function $\lambda(t)$ can be chosen independently and
determines the relative importance of the score-matching loss at different
time scales.

The score function of the perturbed data is calculated by using a Gaussian
perturbation kernel
$p_{\sigma}(\tilde{x}|x):=\mathcal{N}(\mathbf{x},\sigma^{2})$ and
$p_{\sigma}(\mathbf{\tilde{\mathbf{x}}}):=\int
p_{\mathrm{data}}(\mathbf{x})p_{\sigma}(\mathbf{\tilde{x}|x})\mathbf{\mathrm{d}x}$,
simplifying the last term of Eq. 1

$\nabla_{\tilde{\mathbf{x}}}\log
p_{\sigma}(\tilde{\mathbf{x}}|\mathbf{x})=\frac{\mathbf{x}-\tilde{\mathbf{x}}}{\sigma^{2}}\sim\frac{\mathcal{N}(0,1)}{\sigma}.$
(2)

The learned approximation to the score function can then be used to recover
the data probability density by solving the following equation:

$\log p_{0}(\mathbf{x}_{0})=\log
p_{T}(\mathbf{x}_{T})+\int_{0}^{T}\nabla\cdot\tilde{\mathbf{f}}_{\theta}(\mathbf{x}_{t},t)\mathrm{d}t,$
(3)

with

$\tilde{\mathbf{f}}_{\theta}(\mathbf{x}_{t},t)=[f(t)\mathbf{x}_{t}-\frac{1}{2}g(t)^{2}s_{\theta}(\mathbf{x}_{t},t)].$
(4)

The drift ($f$) and diffusion ($g$) coefficients are associated with the
parameters of the Gaussian perturbation kernel. In our studies, we use the
VPSDE [46] framework that defines $f(t)=-\frac{1}{2}\beta(t)$ and
$g(t)=\sqrt{\beta(t)}$ with the parameter
$\beta(t)=\beta_{\textrm{min}}+t\left(\beta_{\textrm{max}}-\beta_{\textrm{min}}\right)$
with $\beta_{\textrm{min}}=0.1$ and $\beta_{\textrm{max}}=20$, resulting in
the induced perturbation kernel
$\mathcal{N}(\mathbf{x}(0)e^{-\frac{1}{2}\int_{0}^{t}\beta(s)\mathrm{d}s},1-e^{-\int_{0}^{t}\beta(s)\mathrm{d}s})$
with $\mathbf{x}(T):=\mathbf{x}(1)\sim\mathcal{N}(0,1)$.

To solve Eq. 3, we need to calculate the divergence of
$\tilde{\mathbf{f}}_{\theta}(\mathbf{x}_{t},t)$, which is computationally
expensive in high dimensions. Instead, we use the Skilling-Hutchinson trace
estimation method [47, 48], reducing the overall computation to near linear
complexity.

While the estimation of the data probability density is independent from the
choice of the weighing function $\lambda(t)$ described in Eq. 1, different
choices can enforce different properties to the learned score function. For
example, a common choice is to set $\lambda(t)=\sigma(t)^{2}$, which avoids
the last ratio in Eq. 2 that diverges as $\sigma(t)\rightarrow 0$ at times
near 0. On the other hand, Ref. [49] shows that choosing $\lambda(t)=g(t)^{2}$
turns the training objective in Eq. 1 into an upper bound to the negative log-
likelihood of the data, effectively allowing the maximum likelihood training
of diffusion models and leading to more precise estimates of the data
probability density. The negative aspect of this choice is that the lack of
the multiplicative $\sigma^{2}$ term can lead to unstable training. This issue
can be mitigated by using an importance sampling scheme that reduces the
variance of the loss function. During the training of the likelihood weighted
objective we implement the same importance sampling scheme as used in Ref.
[49] while in the standard implementation the time component is sampled from
an uniform distribution.

##  III Top Quark Tagging Dataset and Semi-Visible Jets

The top quark tagging dataset is the widely-used community standard benchmark
from Ref. [50, 51]. Events are simulated with Pythia 8 [52, 53] and Delphes
[54, 55] (ATLAS card). The background consists of dijets produced via Quantum
Chromodynamics (QCD) and the signal is top quark pair production with all-
hadronic decays. The default energy flow algorithm in Delphes is used to
create jet consituents, which are clustered using the anti-$k_{T}$ algorithm
with $R=0.8$. [56, 57, 58]. All jets in the range $550$ GeV $<p_{T}<650$ GeV
and $|\eta|<2$ are saved for processing. Each jet is represented with up to
200 constituents (zero-padded if fewer; truncated if more).

In practice, supervised learning should be used to look for top quark
jets222Top quark jet modeling has known inaccuracies, so there still may be
utility in training directly with (unlabeled) data, but since it is possible
to isolate relatively pure samples of top quark jets in data, this is far from
‘anomaly detection’.. To illustrate the anomaly detection abilities of our
approach, we also simulate jets produced from a dark shower within a hidden
valley model [59, 60, 61, 62]. Our dark showers are motivated by333In contrast
to Ref. [63], our mesons have much higher masses, which makes the substructure
more non-trivial. Ref. [63], and consist of a $Z^{\prime}$ with a mass of 1.4
TeV that decays to two dark fermions charged under a strongly coupled U(1)’.
These fermions have a mass of 75 GeV and hadronize into dark pion and $\rho$
mesons, each of which can decay back to the Standard Model. The meson masses
are 150 GeV, resulting in two-prong jet substructure.

##  IV Results

The network implementation and training scheme used to train the diffusion
model are the same ones introduced in Ref. [41], based on the DeepSets [64]
architecture with Transformer layers [65]. This model is trained to learn the
score function of the jet constituents in
$(\Delta\eta,\Delta\phi,\log\mathrm{p_{T}})$ coordinates, while conditioned to
the overall jet kinematics described by
$(\eta_{\text{jet}},\text{mass},N_{\text{part}})$, with the relative particle
coordinates $\Delta\eta=\eta_{\text{part}}-\eta_{\text{jet}}$ and
$\Delta\phi=\phi_{\text{part}}-\phi_{\text{jet}}$, calculated with respect to
the jet axis. The overall jet kinematic information is learned
(simultaneously) by a second diffusion model as done in Ref. [41] using a
model based on the ResNet [66] architecture.

All features are first normalized to the range [-1,1] before training.
Additional pre-processing choices, such as normalizing the data to the range
[0,1], or standardization, by changing the input distributions to have mean 0
and standard deviation of 1, were tested, but did not yield a better density
estimation performance. The probability density is calculated with Eq. 3 using
the Skilling-Hutchinson approximate trace estimation with random samples drawn
from a normal distribution. The integral is solved using SciPy [67] with
explicit Runge-Kutta method of order 5(4) [68, 69] with absolute and relative
tolerances of $10^{-5}$ and $10^{-4}$, respectively. Lower values of the
absolute and relative tolerances were tested with results remaining unchanged.

First, we demonstrate the permutation invariance of the probability density by
evaluating the estimated negative log-likelihood (nll) of the data, trained
using exclusively QCD jets. We show a single jet using different permutations
of the input particles. These results are presented in Figure. 1.
Uncertainties are derived from the standard deviation of 10 independent
estimations of the nll.

Figure 1: Estimated negative log-likelihood in the model trained exclusively
on QCD jets, evaluated on a single jet under multiple permutations of the
input particles. Uncertainties represent the standard deviation from 10
independent evaluations for each permutation.

Since the model was trained only on QCD jet events, the estimated nll tends to
be lower for QCD jets compared to the other classes. We can confirm this
observation by looking at the nll distribution for 10,000 jets shown in Fig.
2.

Figure 2: Negative log-likelihood for QCD, top quark, and $Z^{\prime}$ jets
evaluated on the model trained exclusively on QCD jet events.

We observe lower values of the nll for QCD jet events compared to other jet
categories. The diffusion model training using maximum likelihood
($\lambda(t)=g(t)^{2}$) also presents, on average, lower nll compared to the
standard diffusion approach ($\lambda(t)=\sigma(t)^{2}$). Since the nll,
estimated from the sample containing QCD jets, shows a different distribution
when evaluated against other jet categories, we use the nll as the anomaly
score and we show in Fig. 3 the significance improvement characteristic curve
(SIC).

Figure 3: Significance improvement characteristic curve using the negative-
log-likelihood as the anomaly score.

The nll as an anomaly score results in maximum values for the SIC curve above
1, supporting the choice of metric for anomaly detection. Conversely, the
maximum-likelihood training results in slightly lower SIC curve for anomalous
jets containing the decay products of top quarks. Even though the nll is lower
in the maximum-likelihood training for QCD jets, this value also tends to be
lower for the other jet classes. Similarly, we can train the diffusion model
on a dataset containing only top quark initiated jets and evaluate the
estimated data probability density using different jet categories. The result
is shown in Fig. 4.

Figure 4: Negative log-likelihood for QCD, top quark, and $Z^{\prime}$ jets
evaluated on the model trained exclusively on top quark jet events.

In this case, the nll values for top quark initiated jets are lower on average
compared to the other categories.

A key challenge with unsupervised AD is how to compare different methods.
Weakly supervised methods based on likelihood ratios can be compared with an
optimal classifier using the same noisy label inputs [70, 71] and they
converge to the fully supervised classifier in the limit of large signal,
independent of the signal properties. Unfortunately, there is no analog for
this in the unsupervised case. The existing papers on unsupervised AD compare
methods by demonstrating the performance on a (usually small) set of benchmark
signal models, as we have also done in Fig. 3. However, this is a model-
dependent comparison whose conclusions can easily be altered by simply
changing the physics model(s) considered [20]. As the unsupervised AD
hypothesis is that the new physics, if it exists, is rare given some set of
coordinates, then one could instead directly compare the fidelity of the
density estimation in the background-only case. Since the true probability
density is unknown, this can be achieved using likelihood-ratio methods.

Recent studies have used classifier-based likelihood ratio estimation to
assess and/or improve deep generative models [72, 30, 31, 32, 33, 35, 39, 73,
74, 75]. These classifiers are trained using samples drawn from the generative
model and from the target dataset. As with the training of a Generative
Adversarial Network (GAN) [76], when the classifier is not able to distinguish
the generative model and the target, then the generative model is accurate.
Density estimators are a subclass of generative models and could be evaluated
in this way. However, being able to effectively produce samples and being able
to estimate the probability density are often at odds with each other and so
it is desirable to have a comparison method that uses the probability density
without relying on sampling.

Following Ref. [29], we use another approach that directly assesses the
quality of explicit density-based AD methods. Given two samples (e.g. top
quark and QCD jets), we take the ratio of learned densities (see also Ref.
[70]) and compare the resulting score to a fully supervised classifier trained
to distinguish the same two samples. The likelihood ratio is the optimal
classifier [77] and if the density estimation is exactly correct and the
classifier is optimal, then these two approaches should have the same
performance. Training a supervised classifier is an easier problem (Ref. [70]
versus Ref. [71]), so a reasonable approximation is that the classifier can be
made nearly optimal. For the top-tagging dataset, this problem has already
been studied extensively (see e.g. Ref. [50] and papers that cite it). This
approach does depend on the samples used to estimate the likelihood ratio, but
it is still a sensitive test to the density across the phase space. In Fig. 5,
we calculate the receiver operating characteristic (ROC) curves obtained in
the anomaly detection task using the nll as the anomaly score metric. We also
provide the ROC curves obtained using the log-likelihood ratio between two
dedicated diffusion models, trained exclusively on QCD or top quark jets, and
the one obtained from the outputs of a classifier. The classification network
is trained using the same network architecture as the diffusion model for
particle generation, with additional pooling operation after the last
transformer layer, followed by a fully connected layer with LeakyRelu
activation function and 128 hidden nodes. The output of the classifier is a
single number with a Sigmoid activation function.

Figure 5: Receiver operating characteristic curve obtained from the
unsupervised anomaly detection strategy (left), direct density estimation, and
a supervised classifier (right). The density ratio uses the estimated
densities from individual diffusion models trained either with QCD or top
quark jets. Similarly, the classifier is trained to separate QCD from top
quark jets.

The ROC curve obtained using the log-likelihood ratio has similar area under
the curve (AUC) as the dedicated classifier, even though the performance still
differs significantly in the whole true positive range. Similar results are
found in Ref. [29]. This suggests that even though we are using a state-of-
the-art density estimation strategy, there is still plenty of room to innovate
in order to close the performing gap. Additionally, this illustrates the
danger of relying only on AUC, since it may not be sensitive to tails of phase
space relevant for AD. Similarly to the previous study, we only observe
marginal differences between the results obtained from the different
strategies used to train the diffusion model.

In Table 1, we present a summary of the results consisting of the average nll,
maximum SIC value, and AUC for the anomaly detection task and supervised
density estimation.

Table 1: Comparison of different quality metrics for the anomaly detection task using different datasets. Results are reported using the standard diffusion training with maximum-likelihood training results in parenthesis. For comparison, we present the AUC obtained from the classification of top quarks from QCD jets using the ratio of the estimated densities, or directly training a classifier on the same dataset. Bold quantities represent the best model for a given metric. Dataset |  | Max.
---
SIC
| Unsupervised
---
AUC
| Density
---
Ratio AUC
| Sup.
---
AUC
Top |  1.77 (1.71) |  0.852 (0.847) | 0.954 (0.954) | 0.982
$Z^{\prime}$ |  2.72 (3.99) |  0.908 (0.905) | - | -

##  V Conclusions and Outlook

In this work we presented an unsupervised anomaly detection methodology based
on diffusion models to perform density estimation. Our method approximates the
score function to estimate the probability density of the data. The diffusion
model is trained directly on low-level objects, represented by particles
clustered inside jets. The model for the score function is equivariant with
respect to permutations between particles, leading to a permutation invariant
density estimation. We test different strategies to train the diffusion model,
including a standard implementation and a maximum-likelihood training of the
score model. The maximum-likelihood training presents a lower negative-log-
likelihood, indicating improved probability density estimation. However, using
the negative log-likelihood as the anomaly score does not yield notable
improvements in the anomaly detection task.

Additionally, we evaluate the density estimation performance by studying the
log-likelihood ratio for two density estimators; one trained on QCD jet events
and the other exclusively on top quark jet events. The dedicated classifier
shows a better performance compared to the individual estimation of the log-
likelihood ratio, indicating room for improvement.

For future studies, we plan to investigate alternative diffusion strategies
beyond our implementation to improve the density estimation. Those include
high-order denoising score-matching [78] or using a learnable reweighing
scheme presented in Ref. [79], both showing promising density estimation
performance. There may also be additional applications of high-dimensional,
permutation-invariant density estimation beyond anomaly detection.

## Code Availability

The code for this paper can be found at
https://github.com/ViniciusMikuni/PermutationInvariantAD.git.

## Acknowledgments

We thank Julia Gonski and Manuel Sommerhalder for feedback on the manuscript.
VM and BN are supported by the U.S. Department of Energy (DOE), Office of
Science under contract DE-AC02-05CH11231. This research used resources of the
National Energy Research Scientific Computing Center, a DOE Office of Science
User Facility supported by the Office of Science of the U.S. Department of
Energy under Contract No. DE-AC02-05CH11231 using NERSC award HEP-
ERCAP0021099.

## References

  * Craig  _et al._ [2019] N. Craig, P. Draper, K. Kong, Y. Ng, and D. Whiteson, Acta Phys. Polon. B 50, 837 (2019), arXiv:1610.09392 [hep-ph] .
  * Kim  _et al._ [2020] J. H. Kim, K. Kong, B. Nachman, and D. Whiteson, JHEP 04, 030 (2020), arXiv:1907.06659 [hep-ph] .
  * Karagiorgi  _et al._ [2021] G. Karagiorgi, G. Kasieczka, S. Kravitz, B. Nachman, and D. Shih, (2021), arXiv:2112.03769 [hep-ph] .
  * Collins  _et al._ [2019] J. H. Collins, K. Howe, and B. Nachman, Phys. Rev. D99, 014038 (2019), arXiv:1902.02634 [hep-ph] .
  * Collins  _et al._ [2018] J. H. Collins, K. Howe, and B. Nachman, Phys. Rev. Lett. 121, 241803 (2018), arXiv:1805.02664 [hep-ph] .
  * D’Agnolo and Wulzer [2019] R. T. D’Agnolo and A. Wulzer, Phys. Rev. D99, 015014 (2019), arXiv:1806.02350 [hep-ph] .
  * Hajer  _et al._ [2018] J. Hajer, Y.-Y. Li, T. Liu, and H. Wang,  (2018), 10.1103/PhysRevD.101.076015, arXiv:1807.10261 [hep-ph] .
  * Farina  _et al._ [2018] M. Farina, Y. Nakai, and D. Shih,  (2018), 10.1103/PhysRevD.101.075021, arXiv:1808.08992 [hep-ph] .
  * Heimel  _et al._ [2019] T. Heimel, G. Kasieczka, T. Plehn, and J. M. Thompson, SciPost Phys. 6, 030 (2019), arXiv:1808.08979 [hep-ph] .
  * [10] HEP ML Community, “A Living Review of Machine Learning for Particle Physics,” .
  * Aarrestad  _et al._ [2021] T. Aarrestad _et al._ , (2021), arXiv:2105.14027 [hep-ph] .
  * Kasieczka  _et al._ [2021] G. Kasieczka _et al._ , (2021), arXiv:2101.08320 [hep-ph] .
  * Aad  _et al._ [2020] G. Aad _et al._ (ATLAS), Phys. Rev. Lett. 125, 131801 (2020), arXiv:2005.02983 [hep-ex] .
  * ATL [2022] _Anomaly detection search for new resonances decaying into a Higgs boson and a generic new particle $X$ in hadronic final states using $\sqrt{s}$ = 13 TeV $pp$ collisions with the ATLAS detector_, Tech. Rep. (CERN, Geneva, 2022).
  * Shih  _et al._ [2021] D. Shih, M. R. Buckley, L. Necib, and J. Tamanas, (2021), arXiv:2104.12789 [astro-ph.GA] .
  * Shih  _et al._ [2023] D. Shih, M. R. Buckley, and L. Necib, (2023), arXiv:2303.01529 [astro-ph.GA] .
  * Pettee  _et al._ [2023] M. Pettee, S. Thanvantri, B. Nachman, D. Shih, M. R. Buckley, and J. H. Collins, (2023), arXiv:2305.03761 [astro-ph.GA] .
  * ATL [2023] (2023).
  * Metodiev  _et al._ [2017] E. M. Metodiev, B. Nachman, and J. Thaler, JHEP 10, 174 (2017), arXiv:1708.02949 [hep-ph] .
  * Collins  _et al._ [2021] J. H. Collins, P. Martín-Ramiro, B. Nachman, and D. Shih, (2021), arXiv:2104.02092 [hep-ph] .
  * Kasieczka  _et al._ [2023] G. Kasieczka, R. Mastandrea, V. Mikuni, B. Nachman, M. Pettee, and D. Shih, Phys. Rev. D 107, 015009 (2023), arXiv:2209.06225 [hep-ph] .
  * Roy and Vijay [2019] T. S. Roy and A. H. Vijay, (2019), arXiv:1903.02032 [hep-ph] .
  * Amram and Suarez [2020] O. Amram and C. M. Suarez,  (2020), 10.1007/JHEP01(2021)153, arXiv:2002.12376 [hep-ph] .
  * Kahn  _et al._ [2021] A. Kahn, J. Gonski, I. Ochoa, D. Williams, and G. Brooijmans, (2021), arXiv:2105.09274 [hep-ph] .
  * Atkinson  _et al._ [2021] O. Atkinson, A. Bhardwaj, C. Englert, V. S. Ngairangbam, and M. Spannowsky, JHEP 08, 080 (2021), arXiv:2105.07988 [hep-ph] .
  * Dillon  _et al._ [2022] B. M. Dillon, L. Favaro, T. Plehn, P. Sorrenson, and M. Krämer, (2022), arXiv:2206.14225 [hep-ph] .
  * Atkinson  _et al._ [2022] O. Atkinson, A. Bhardwaj, C. Englert, P. Konar, V. S. Ngairangbam, and M. Spannowsky, Front. Artif. Intell. 5, 943135 (2022), arXiv:2204.12231 [hep-ph] .
  * Andreassen  _et al._ [2018] A. Andreassen, I. Feige, C. Frye, and M. D. Schwartz, (2018), arXiv:1804.09720 [hep-ph] .
  * Finke  _et al._ [2023] T. Finke, M. Krämer, A. Mück, and J. Tönshoff, (2023), arXiv:2303.07364 [hep-ph] .
  * Krause and Shih [2021a] C. Krause and D. Shih, (2021a), arXiv:2106.05285 [physics.ins-det] .
  * Krause and Shih [2021b] C. Krause and D. Shih, (2021b), arXiv:2110.11377 [physics.ins-det] .
  * Diefenbacher  _et al._ [2023] S. Diefenbacher, E. Eren, F. Gaede, G. Kasieczka, C. Krause, I. Shekhzadeh, and D. Shih, (2023), arXiv:2302.11594 [physics.ins-det] .
  * Krause  _et al._ [2022] C. Krause, I. Pang, and D. Shih, (2022), arXiv:2210.14245 [physics.ins-det] .
  * Käch  _et al._ [2022] B. Käch, D. Krücker, I. Melzer-Pellmann, M. Scham, S. Schnake, and A. Verney-Provatas, (2022), arXiv:2211.13630 [hep-ex] .
  * Buckley  _et al._ [2023] M. R. Buckley, C. Krause, I. Pang, and D. Shih, (2023), arXiv:2305.11934 [physics.ins-det] .
  * Rezende and Mohamed [2015] D. J. Rezende and S. Mohamed, International Conference on Machine Learning 37, 1530 (2015).
  * Kobyzev  _et al._ [2020] I. Kobyzev, S. Prince, and M. Brubaker, IEEE Transactions on Pattern Analysis and Machine Intelligence , 1 (2020).
  * Verheyen [2022] R. Verheyen, SciPost Phys. 13, 047 (2022), arXiv:2205.01697 [hep-ph] .
  * Mikuni and Nachman [2022] V. Mikuni and B. Nachman, Phys. Rev. D 106, 092009 (2022), arXiv:2206.11898 [hep-ph] .
  * Leigh  _et al._ [2023] M. Leigh, D. Sengupta, G. Quétant, J. A. Raine, K. Zoch, and T. Golling, (2023), arXiv:2303.05376 [hep-ph] .
  * Mikuni  _et al._ [2023] V. Mikuni, B. Nachman, and M. Pettee, (2023), arXiv:2304.01266 [hep-ph] .
  * Buhmann  _et al._ [2023] E. Buhmann, S. Diefenbacher, E. Eren, F. Gaede, G. Kasieczka, A. Korol, W. Korcari, K. Krüger, and P. McKeown, (2023), arXiv:2305.04847 [physics.ins-det] .
  * Butter  _et al._ [2023] A. Butter, N. Huetsch, S. P. Schweitzer, T. Plehn, P. Sorrenson, and J. Spinner, (2023), arXiv:2305.10475 [hep-ph] .
  * Song  _et al._ [2021a] Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon, and B. Poole, ArXiv abs/2011.13456 (2021a).
  * Vincent [2011] P. Vincent, Neural Computation 23, 1661 (2011).
  * Ho  _et al._ [2020] J. Ho, A. Jain, and P. Abbeel, Advances in Neural Information Processing Systems 33, 6840 (2020).
  * Skilling [1989] J. Skilling, “The eigenvalues of mega-dimensional matrices,” in  _Maximum Entropy and Bayesian Methods: Cambridge, England, 1988_, edited by J. Skilling (Springer Netherlands, Dordrecht, 1989) pp. 455–466.
  * Hutchinson [1990] M. Hutchinson, Communications in Statistics - Simulation and Computation 19, 433 (1990), https://doi.org/10.1080/03610919008812866 .
  * Song  _et al._ [2021b] Y. Song, C. Durkan, I. Murray, and S. Ermon, Advances in Neural Information Processing Systems 34, 1415 (2021b).
  * Butter  _et al._ [2019] A. Butter _et al._ , SciPost Phys. 7, 014 (2019), arXiv:1902.09914 [hep-ph] .
  * Kasieczka  _et al._ [2019] G. Kasieczka, T. Plehn, J. Thompson, and M. Russel, “Top quark tagging reference dataset,”  (2019).
  * Sjostrand  _et al._ [2006] T. Sjostrand, S. Mrenna, and P. Z. Skands, JHEP 05, 026 (2006), arXiv:hep-ph/0603175 .
  * Sjöstrand  _et al._ [2015] T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten, S. Mrenna, S. Prestel, C. O. Rasmussen, and P. Z. Skands, Comput. Phys. Commun. 191, 159 (2015), arXiv:1410.3012 [hep-ph] .
  * Mertens [2015] A. Mertens, J. Phys. Conf. Ser. 608, 012045 (2015).
  * Selvaggi [2014] M. Selvaggi, J. Phys. Conf. Ser. 523, 012033 (2014).
  * Cacciari and Salam [2006] M. Cacciari and G. P. Salam, Phys. Lett. B 641, 57 (2006), arXiv:hep-ph/0512210 .
  * Cacciari  _et al._ [2012] M. Cacciari, G. P. Salam, and G. Soyez, Eur. Phys. J. C 72, 1896 (2012), arXiv:1111.6097 [hep-ph] .
  * Cacciari  _et al._ [2008] M. Cacciari, G. P. Salam, and G. Soyez, JHEP 04, 063 (2008), arXiv:0802.1189 [hep-ph] .
  * Strassler and Zurek [2007] M. J. Strassler and K. M. Zurek, Phys. Lett. B 651, 374 (2007), arXiv:hep-ph/0604261 .
  * Carloni and Sjostrand [2010] L. Carloni and T. Sjostrand, JHEP 09, 105 (2010), arXiv:1006.2911 [hep-ph] .
  * Carloni  _et al._ [2011] L. Carloni, J. Rathsman, and T. Sjostrand, JHEP 04, 091 (2011), arXiv:1102.3795 [hep-ph] .
  * Knapen  _et al._ [2021] S. Knapen, J. Shelton, and D. Xu, Phys. Rev. D 103, 115013 (2021), arXiv:2103.01238 [hep-ph] .
  * Buss  _et al._ [2022] T. Buss, B. M. Dillon, T. Finke, M. Krämer, A. Morandini, A. Mück, I. Oleksiyuk, and T. Plehn, (2022), arXiv:2202.00686 [hep-ph] .
  * Zaheer  _et al._ [2017] M. Zaheer, S. Kottur, S. Ravanbakhsh, B. Póczos, R. Salakhutdinov, and A. J. Smola, CoRR abs/1703.06114 (2017), 1703.06114 .
  * Vaswani  _et al._ [2017] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. Kaiser, and I. Polosukhin, CoRR abs/1706.03762 (2017), 1706.03762 .
  * He  _et al._ [2015] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” (2015), arXiv:1512.03385 [cs.CV] .
  * Virtanen  _et al._ [2020] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, Nature Methods 17, 261 (2020).
  * Dormand and Prince [1980] J. R. Dormand and P. J. Prince, Journal of computational and applied mathematics 6, 19 (1980).
  * Shampine [1986] L. F. Shampine, Mathematics of computation 46, 135 (1986).
  * Nachman and Shih [2020] B. Nachman and D. Shih, Phys. Rev. D 101, 075042 (2020), arXiv:2001.04990 [hep-ph] .
  * Hallin  _et al._ [2021] A. Hallin, J. Isaacson, G. Kasieczka, C. Krause, B. Nachman, T. Quadfasel, M. Schlaffer, D. Shih, and M. Sommerhalder, (2021), arXiv:2109.00546 [hep-ph] .
  * Diefenbacher  _et al._ [2020] S. Diefenbacher, E. Eren, G. Kasieczka, A. Korol, B. Nachman, and D. Shih, Journal of Instrumentation 15, P11004 (2020), arXiv:2009.03796 [hep-ph] .
  * Winterhalder  _et al._ [2021] R. Winterhalder, M. Bellagente, and B. Nachman, (2021), arXiv:2106.00792 [stat.ML] .
  * Nachman and Winterhalder [2023] B. Nachman and R. Winterhalder, (2023), arXiv:2305.07696 [hep-ph] .
  * Das  _et al._ [2023] R. Das, L. Favaro, T. Heimel, C. Krause, T. Plehn, and D. Shih, (2023), arXiv:2305.16774 [hep-ph] .
  * Goodfellow  _et al._ [2014] I. J. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio, (2014), arXiv:1406.2661 [stat.ML] .
  * Neyman and Pearson [1933] J. Neyman and E. S. Pearson, Phil. Trans. R. Soc. Lond. A 231, 289 (1933).
  * Lu  _et al._ [2022] C. Lu, K. Zheng, F. Bao, J. Chen, C. Li, and J. Zhu, in  _International Conference on Machine Learning_ (PMLR, 2022) pp. 14429–14460.
  * Kingma  _et al._ [2021] D. Kingma, T. Salimans, B. Poole, and J. Ho, in  _Advances in Neural Information Processing Systems_, Vol. 34, edited by M. Ranzato, A. Beygelzimer, Y. Dauphin, P. Liang, and J. W. Vaughan (Curran Associates, Inc., 2021) pp. 21696–21707.