# The Theory of Generalised Hydrodynamics for the One-dimensional Bose Gas

M. L. Kerr  School of Mathematics and Physics, University of Queensland,
Brisbane, Queensland 4072, Australia  K. V. Kheruntsyan  School of Mathematics
and Physics, University of Queensland, Brisbane, Queensland 4072, Australia

###### Abstract

This article reviews the recent developments in the theory of generalised
hydrodynamics (GHD) with emphasis on the repulsive one-dimensional Bose gas.
We discuss the implications of GHD on the mechanisms of thermalisation in
integrable quantum many-body systems as well as its ability to describe far-
from-equilibrium behaviour of integrable and near integrable systems in a
variety of quantum quench scenarios. We outline the experimental tests of GHD
in cold-atom gases and its benchmarks with other microscopic theoretical
approaches. Finally, we offer some perspectives on the future direction of the
development of GHD.

Keywords: quantum many-body dynamics, integrable systems, generalised
hydrodynamics, ultra-cold quantum gases, Lieb-Liniger model

* Corresponding author<EMAIL_ADDRESS>
##  I Introduction

The study of non-equilibrium phenomena in isolated quantum systems has
received much attention in recent years [1, 2, 3, 4]. A fruitful avenue into
non-equilibrium dynamics has emerged from the study of integrable systems.
These systems are characterised by a large number of conserved quantities and
include several paradigmatic models, such as the Lieb-Liniger model describing
a gas of one-dimensional bosons [5, 6], and the Hubbard model describing
electrons in a solid [7]. Integrability is known to fundamentally affect the
dynamics of isolated quantum systems in contrast to that of generic non-
integrable systems [3, 8, 2, 9]. The most prominent consequence of this is the
preclusion of thermalisation, whereby a system does not relax to a thermal
state. Rather, integrable systems will typically relax to non-thermal,
equilibrium states, described by the so-called generalised Gibbs ensemble [10,
11, 12]. This has been demonstrated in various models and theories, such as
conformal field theories [13, 14], the repulsive Lieb-Liniger gas [15, 16],
the quantum Ising chain in a transverse field [17, 18], and in Heisenberg spin
chains [19].

In classical systems, thermalisation can be well understood in terms of
dynamical chaos and the notions of ergodicity and mixing [20] in which a
system explores its phase space uniformly and densely for almost all initial
conditions. This enables a description of the system via classical statistical
mechanics. However, the mechanisms of thermalisation for isolated quantum
systems are not well-established. Understanding these mechanisms has been the
subject of many theoretical and experimental investigations over the past two
decades.

An indispensable protocol to induce non-equilibrium scenarios in isolated
quantum systems is the so-called quantum quench [13, 11, 21] where an initial
equilibrium state of a many-body Hamiltonian is suddenly made to evolve
unitarily under another Hamiltonian. These quench scenarios constitute a vital
outpost to address the lack of thermalisation in integrable quantum systems.
However, as is the case for any quantum system, a direct theoretical approach
with exact microscopic calculations of physically relevant observables
involves the diagonalisation of the many-body Hamiltonian associated with the
system. The dimension of such a many-body Hilbert space notoriously scales
exponentially with the total number of particles. This often renders direct
computational approaches intractable beyond more than a dozen particles. One
can transcend this exponentially-large space of quantum states by instead
employing an aggregated, coarse-grained description of the system with a
reduced number of degrees of freedom. One of the most prominent and successful
such frameworks is hydrodynamics.

At the heart of any hydrodynamic theory is a set of conservation laws that
govern a specific model’s dynamical degrees of freedom over vast space and
time scales. These quantities typically include energy, momentum, and particle
number density. Hydrodynamics has been particularly successful as an effective
theory for the emergent dynamical behaviour of a variety of interacting many-
body systems. For instance, Landau’s celebrated two-fluid model of superfluid
helium [22], the description of electron currents in graphene [23], and in the
magneto-hydrodynamic description of electrically conducting fluids [24].

Hydrodynamic theories are a valuable resource to address the mechanisms of
thermalisation in isolated integrable quantum systems. Here our interest lies
in the rich many-body physics of one-dimensional integrable systems, which,
besides being a mathematical curiosity, are also relevant to ultracold Bose
and Fermi gases in highly anisotropic traps, and to other low-dimensional
condensed matter systems such as superconducting nanowires [25] and 1D
polariton condensates [26]. In particular, this review focuses on the one-
dimensional Bose gas. This is motivated by the high degree of experimental
control over system parameters and dynamics, owing to unprecedented
advancements over the past twenty years [27, 28, 29].

By utilising the exact integrability of the one-dimensional Bose gas in the
uniform limit, its equilibrium thermodynamic properties can be derived exactly
using Yang-Yang thermodynamics [30]. One can then insert these thermodynamic
quantities into the classical hydrodynamic equations. This approach has been
dubbed “conventional hydrodynamics” (CHD) and provides an excellent
description of the collective excitations of 1D bosons driven out of
equilibrium [31, 32, 33]. The equations of CHD take the form of Euler
hydrodynamic equations, which are valid over large space and time scales.
However, this theory is limited in its applicability by the assumption of
local thermal equilibrium, in which the system is divided into small cells
assumed to thermalise sufficiently fast with their environment. In isolated
integrable quantum systems, – wherein the 1D Bose gas is an example – this
assumption of fast thermalisation times is not a priori justified. Moreover,
simulations of the Euler equations of conventional hydrodynamics are often
plagued with the gradient catastrophe problem, resulting in unphysical shock
waves [34].

In 2016, two papers sparked the development of a hydrodynamic theory
specialised to integrable systems, known as the theory of generalised
hydrodynamics (GHD) [35, 36]. In contrast to CHD, GHD does not rely on the
assumption of local thermal equilibrium in the canonical Gibbs ensemble sense.
Instead, systems are assumed to relax to equilibrium states, described by a
generalised Gibbs ensemble (GGE)111The term “generalised” in GHD is used in
the same sense as in generalised Gibbs ensemble.. GHD avoids the gradient
catastrophe problem typical of scenarios involving, e.g., shock waves, and
thus, can describe systems very far from equilibrium [34]. In particular, GHD
has recently been shown to reproduce the most striking effects observed in the
quantum Newton’s cradle experiment [34, 38]. The applicability of GHD extends
well beyond ultra-cold atomic physics. Indeed, GHD has proven to be a robust
theory applicable in both classical and quantum models over large length and
time scales [36].

GHD has been extended to account for various experimentally relevant physical
effects. This review aims to overview the theory of generalised hydrodynamics
and recent developments since its original formulation in 2016. In doing so,
we restrict ourselves to the applications of GHD to describe the dynamics of
the repulsive one-dimensional (1D) Bose gas described by the Lieb-Liniger
model [5, 6].

The organisation of this review is as follows. In section II, we introduce the
Lieb-Liniger model and the equations of GHD at the Euler scale. We also
describe the quantum Newton’s cradle experiment, its description via GHD, and
the implications on the mechanisms of thermalisation. In section III, we
outline several extensions of the original formulation of GHD. In section IV,
we review the experimental tests of GHD and the benchmarks against other
theoretical approaches. Finally, section V offers some perspectives on the
future direction of GHD and some open problems regarding GHD.

##  II Euler-scale Generalised Hydrodynamics

###  II.1 Theoretical considerations

The Lieb-Liniger Hamiltonian for $N$ bosons in a uniform box of length $L$
(with periodic boundary condition) interacting via the two-body contact
interaction potential $U(x,x^{\prime})=g\delta(x-x^{\prime})$ is given by [5,
*Lieb-Liniger-II]

$H_{\text{LL}}=-\frac{\hbar^{2}}{2m}\sum_{i=1}^{N}\frac{\partial^{2}}{\partial
x_{i}^{2}}+g\sum_{1\leq i<j<N}\delta(x_{i}-x_{j}),$ (1)

where $g$ quantifies the strength of interactions, assumed here to be
repulsive ($g>0$). The Hamiltonian (1) is integrable and exactly solvable
using Bethe ansatz, and as such it admits an infinite number of conservation
laws. In the second quantized form, it can be rewritten as

$\displaystyle\hat{H}_{\text{LL}}=$ $\displaystyle-\frac{\hbar^{2}}{2m}\int
dx\,\hat{\Psi}^{\dagger}(x)\frac{\partial^{2}}{\partial x^{2}}\hat{\Psi}(x)$
$\displaystyle+\frac{g}{2}\int
dx\,\hat{\Psi}^{\dagger}(x)\hat{\Psi}^{\dagger}(x)\hat{\Psi}(x)\hat{\Psi}(x).$
(2)

Such a one-dimensional system can be experimentally realised by confining an
ultracold gas of bosons to, e.g., a highly elongated harmonic trap with
transverse frequency $\omega_{\perp}$ and axial frequency
$\omega\ll\omega_{\perp}$. When the transverse excitation energy is much
larger than all other relevant energies of the problem, such as the average
thermal energy and the chemical potential of the system,
$\hbar\omega_{\perp}\gg\max\\{k_{B}T,\,\mu$}, the transverse excitations are
negligible, and the dynamics take place only along the longitudinal dimension
while being frozen out in the transverse dimension. For sufficiently large
systems, the boundary effects can be neglected, so that the above Hamiltonian
can describe the properties of systems that are not necessarily periodic or
even uniform, wherein the inhomogeneities due to the longitudinal trapping
$V(x)$ can be accounted for within the local density approximation (LDA) [39].

An important parameter encoding the strength of interactions between particles
in a uniform 1D Bose gas is the dimensionless Lieb-Liniger parameter $\gamma$
defined as

$\gamma=\frac{mg}{\hbar^{2}n},$ (3)

where $n=N/L$ is the 1D (linear) density. When $\gamma\ll 1$, the system is
weakly interacting. Conversely, for $\gamma\gg 1$, the interaction energy is
large, and the system is strongly interacting. Note that one can equivalently
enter the strongly interacting regime $\gamma\gg 1$ by either increasing the
interaction strength $g$ or decreasing the density of the system $n$.

In addition to the dimensionless interaction strength $\gamma$, one can also
define a dimensionless temperature parameter, $\mathcal{T}$, by scaling the
temperature of the system $T$ by the temperature of quantum degeneracy
$T_{d}=\hbar^{2}n^{2}/2mk_{B}$,

$\mathcal{T}=\frac{T}{T_{d}}=\frac{2mk_{B}T}{\hbar^{2}n^{2}}.$ (4)

When $T\sim T_{d}$, the thermal de Broglie wavelength of the particles is on
the order of the mean interparticle separation. This represents the
temperature regime below which quantum effects begin to dominate. Overall,
finite temperature equilibrium properties of 1D Bose gas systems can be
studied using the thermodynamics Bethe ansatz first derived and solved by Yang
and Yang [30].

To model the dynamics of a 1D Bose gas in the hydrodynamic sense, we view the
system as a continuum of mesoscopic fluid cells that are thermodynamically
large but small enough to be considered spatially homogeneous. For equilibrium
states, this assumption is equivalent to the LDA. Due to integrability, our
system admits a set of conserved charges $Q_{i}$ ($i=1,2,3,...$), such as
energy, particle number, and momentum, each of which we assume can be written
as an integral of a corresponding charge density, $q_{i}(x,t)$, i.e.,

$Q_{i}(t)=\int dx\,q_{i}(x,t).$ (5)

To each charge density, $q_{i}(x,t)$, we have a corresponding current density
$j_{i}(x,t)$ satisfying a continuity equation,

$\partial_{t}q_{i}+\partial_{x}j_{i}=0.$ (6)

The fundamental assumption of GHD is that after some relaxation time, an
inhomogeneous non-stationary system approaches, within each fluid cell, states
which have maximised entropy with respect to each of the conserved quantities.
These maximal entropy states are described by generalised Gibbs ensembles
(GGE) with density matrices of the form,

$\rho_{\text{GGE}}(x,t)\propto e^{-\sum_{i}\beta^{i}(x,t)Q_{i}},$ (7)

where $\beta^{i}(x,t)$ is the Lagrange multiplier corresponding to the charge
$Q_{i}$. The density matrix $\rho_{\text{GGE}}$ is related to the average of
the charge density $q_{i}(x,t)$ via

$\mathfrak{q}_{i}(x,t):=\braket{q_{i}(x,t)}=\frac{\text{tr}\left(\rho_{\text{GGE}}\,q_{i}(x,t)\right)}{\text{tr}(\rho_{\text{GGE}})}.$
(8)

The set of average conserved densities $\\{\mathfrak{q}_{i}\\}_{i=1}^{\infty}$
can be considered as a set of coordinates for the manifold of maximal entropy
states. In principle, this gives a complete (coarse-grained) description of
our system as the set of average conserved densities specifies a particular
configuration of our system, described as a point in the manifold of maximal
entropy states. An alternative set of coordinates can be obtained by utilising
the quasi-particle formulation of the thermodynamic Bethe ansatz (TBA) [40,
30], in which the maximal entropy states are described via a phase space
density of quasi-particles. As we will see, these coordinates offer a
significant improvement over the set $\\{\mathfrak{q}_{i}\\}$ as the density
of quasi-particles satisfies a single partial differential equation and once
obtained, it is possible to construct each $\mathfrak{q}_{i}$. We note that
whilst the techniques of the thermodynamic Bethe ansatz enable a description
of the equilibrium properties of a uniform system, the applicability of GHD
extends beyond this to non-equilibrium scenarios in spatially inhomogeneous
systems within the so-called local density approximation (LDA), which is
implicitly assumed in any hydrodynamic theory.

The main quantity of GHD is the density of quasi-particles, also known as the
rapidity distribution, denoted $\rho_{p}(\theta,x,t)$222In this review, we
only consider many-body models of identical particles of single species, and
hence, a single quasi-particle distribution. It is straightforward to
generalise the equations of GHD to many species., that carry quasi-momentum
$m\theta$ at the space-time point $(x,t)$. This quantity is analogous to the
equilibrium root density in the thermodynamic Bethe ansatz solution. However,
the key postulate of GHD is that the rapidity distribution in our out-of-
equilibrium problem evolves in time according to a classical Euler-like
hydrodynamic equation. More specifically, the evolution of the rapidity
distribution in GHD is governed by the set of integro-differential equations
[35, 36, 42],

$\partial_{t}\rho_{p}+\partial_{x}\left(v^{\text{eff}}(\theta)\rho_{p}\right)-\frac{1}{m}\left(\partial_{x}V\right)\partial_{\theta}\rho_{p}=0,$
(9)

where $V(x)$ is an external potential, $v^{\text{eff}}(\theta)$ is the
effective velocity, defined via the integral equation,

$v^{\text{eff}}(\theta)=\theta+\int_{-\infty}^{\infty}d\theta^{\prime}\,\Phi(\theta-\theta^{\prime})\left(v^{\text{eff}}(\theta^{\prime})-v^{\text{eff}}(\theta)\right)\rho_{p}(\theta^{\prime}),$
(10)

and $\Phi(\theta-\theta^{\prime})$ is the two-body scattering shift, which in
the Lieb-Liniger model takes the form,

$\Phi(\theta)=\frac{\hbar^{2}}{m}\frac{2g}{g^{2}+\hbar^{2}\theta^{2}}.$ (11)

Note we have suppressed the $x$ and $t$ dependence of $\rho_{p}$ and
$v^{\text{eff}}$ above for simplicity. These are the main equations of GHD at
the Euler scale, i.e., when quantities vary very slowly over space and time.

Figure 1: Evolution of the rapidity distribution over the first ten
oscillation cycles for a strongly-interacting 1D Bose gas in the quantum
Newton’s cradle scenario in a harmonic potential (top row) and with a small
anharmonicity (middle row). Despite the presence of an integrability-breaking
trapping potential and dephasing effects, neither system thermalises; the
dephased states cannot be identified with a thermal state. The bottom row
shows the corresponding density profiles, obtained by integrating the rapidity
distribution $\rho_{p}(\theta,x)$ over all rapidities $\theta$. The blue curve
corresponds to the harmonic potential, whilst the red curve is for the
anharmonic trap. Adapted from [38]

Note that the effective velocity is a functional of the rapidity distribution
and so two equations are coupled. Physically, this effective velocity can be
interpreted as the large-scale, coarse-grained velocity of a quasi-particle as
it travels through the gas, taking into account the scattering shifts it
accumulates at collisions with the other quasi-particles [43]. For the uniform
Lieb-Liniger gas, the initial rapidity distribution $\rho_{p}(\theta,x,0)$
that is supplied to the GHD equations is often taken as the thermal
equilibrium distribution obtained via the Yang-Yang thermodynamic Bethe ansatz
[30]. For the nonuniform case, the initial rapidity distribution can be found
locally for each $x$ using the same thermodynamic Bethe ansatz but within the
LDA [39].

We note that although GHD was first introduced for quantum field theories [35]
and quantum chains [36], equations (9) and (10) have previously been derived
rigorously in the context of the classical hard rod gas [44, 45] and in
soliton gases [46, 47].

Given the rapidity distribution, $\rho_{p}(\theta,x,t)$, at some fixed time,
the expectation value of the conserved charge densities and their
corresponding currents can then be computed via

$\displaystyle\mathfrak{q}_{i}(x,t)$ $\displaystyle=\int
d\theta\,h_{i}(\theta)\,\rho_{p}(\theta),$ (12)
$\displaystyle\mathfrak{j}_{i}(x,t)$ $\displaystyle=\int
d\theta\,h_{i}(\theta)\,v^{\text{eff}}(\theta)\,\rho_{p}(\theta).$ (13)

where $h_{i}(\theta)$ is the one-particle eigenvalue of the $i$-th conserved
charge. For example, the average particle density can be obtained by setting
$h_{i}(\theta)=1$, whilst the average energy density is found by letting
$h_{i}(\theta)=m\theta^{2}/2$. Intuitively, equations (12) and (13) can be
understood as follows: the contribution of a quasi-particle of rapidity
$\theta$ to the $i$-th conserved charge is $h_{i}(\theta)$, moreover, the
number of quasi-particles with rapidity in the interval
$[\theta,\theta+d\theta)$ is given by $\rho_{p}(\theta)d\theta$. Therefore,
the contribution of all quasi-particles with rapidity in the range
$[\theta,\theta+d\theta)$ to the $i$-th charge is given by
$h_{i}(\theta)\rho_{p}(\theta)d\theta$. Summing over all such rapidities
yields Eq. (12). Similarly, one obtains Eq. (13) for the average current
density. A rigorous derivation of Eq. (13) has been derived using relativistic
quantum field theory, and its validity has been verified numerically in both
classical and quantum models [48, 36].

It is often more convenient to encode the thermodynamic properties of a system
via the filling factor,

$\vartheta(\theta,x,t):=\frac{\rho_{p}(\theta)}{\rho_{p}(\theta)+\rho_{h}(\theta)}=\frac{\rho_{p}(\theta)}{\rho_{s}(\theta)}.$
(14)

where $\rho_{h}$ is the density of ‘holes’ and $\rho_{s}(\theta)$ is the
density of states which is related to the density of quasi-particles,
$\rho_{p}$, via the thermodynamic Bethe equation,

$2\pi\rho_{s}(\theta)=1+\int_{-\infty}^{\infty}d\theta^{\prime}\,\Phi(\theta-\theta^{\prime})\rho_{p}(\theta^{\prime}).$
(15)

One can then construct the average charge and current densities using

$\displaystyle\mathfrak{q}_{i}(x,t)$
$\displaystyle=\int\frac{d\theta}{2\pi}\,\vartheta(\theta)\,h^{\text{dr}}_{i}(\theta),$
(16) $\displaystyle\mathfrak{j}_{i}(x,t)$
$\displaystyle=\int\frac{d\theta}{2\pi}\,\vartheta(\theta)\,\theta\,h^{\text{dr}}_{i}(\theta),$
(17)

where the dressing operation $f\mapsto f^{\text{dr}}$ is defined via the
integral equation,

$f^{\text{dr}}(\theta)=f(\theta)+\int\frac{d\theta^{\prime}}{2\pi}\,\Phi(\theta^{\prime}-\theta)\vartheta(\theta^{\prime})f^{\text{dr}}(\theta^{\prime}).$
(18)

In this formulation, Eq. (9) becomes,

$\partial_{t}\vartheta+v^{\text{eff}}\partial_{x}\vartheta-\frac{1}{m}\left(\partial_{x}V\right)\partial_{\theta}\vartheta=0,$
(19)

where the effective velocity is now outside the spatial derivative.
Numerically, Eq. (19) is more convenient to work with than Eq. (9). We also
advertise the open-source Matlab framework “iFluid” that numerically solves
the equations of GHD [49].

###  II.2 Modelling the quantum Newton’s cradle experiment

Arguably, the lack of thermalisation in isolated quantum systems was best
demonstrated in the hallmark quantum Newton’s cradle experiment of Kinoshita
et al. [50]. There, clouds of strongly-interacting rubidium atoms confined to
a one-dimensional harmonic trap undergo repeated collisions without noticeably
thermalising on observable time scales – corresponding to thousands of
collisions. The lack of thermalisation can be attributed to the integrability
of the underlying Lieb-Liniger model in the uniform limit and to weak
integrability in the nonuniform (harmonically trapped) system owing to the
applicability of the local density approximation [39].

In the decade that followed the Kinoshita experiment, a quantitatively
accurate model of the quantum Newton’s cradle with experimentally relevant
parameters remained elusive. GHD has emerged as an ideal tool to model this
experiment. In 2018, Caux et al. used GHD to simulate the dynamics of a
strongly-interacting 1D Bose gas in the quantum Newton’s cradle experiment
[38]. The evolution of the rapidity distribution is shown in Fig. 1 for the
case of a harmonic potential and an anharmonic potential to mimic the trapping
potential in the original quantum Newton’s cradle experiment.

Caux et al. found that even in the presence of a trapping potential that
weakly breaks integrability, the 1D gas does not thermalise. Instead, the
system relaxes to a generalised Gibbs ensemble. The preclusion of
thermalisation here can be attributed to the existence of conserved quantities
that are incompatible with convergence towards thermal equilibrium [38]. These
quantities take the form

$Q[f]=\int d\theta\,dx\,f(\vartheta(\theta,x,t))\rho_{p}(\theta,x,t),$ (20)

where $f$ is an arbitrary function and $\rho_{p}(\theta,x,t)$ is continuous in
$\theta$ and $x$. The fact that these quantities are conserved under GHD
evolution in a trap follows directly from Eq. (9) and Eq. (19). We note that
these quantities are only conserved at the Euler scale.

##  III Beyond Euler-scale Generalised Hydrodynamics

Here we present some of the extensions to the original formulation of GHD.
Namely, the incorporation of diffusive and quantum effects, the description of
dimensional crossover in one-dimensional Bose gases, and the inclusion of
space-time dependent interactions.

###  III.1 Diffusive effects

Despite the excellent predictive power of standard GHD, it is often necessary
to go beyond the, lowest-order, Euler scale. For instance, spin and charge
transport in quantum chains have been shown to exhibit diffusion, and other
non-Eulerian behaviours not captured by Euler scale GHD [51, 52, 53].

Figure 2: Simulation of diffusive GHD for a weakly interacting 1D Bose gas in
the quantum Newton’s cradle setup initiated via a quench of a double-well trap
to a single harmonic well. (a) Evolution of the density profile. (b) The
density profile at three fixed times compared to its distribution in thermal
equilibrium. In contrast to Euler-scale GHD, this system eventually
thermalizes over sufficiently large time scales when diffusive effects are
considered. Adapted from [54].

As a first approximation, we assume that the state of the system and local
observables at the point $(x,t)$ can be described by the averages of conserved
charges in a neighbourhood of this point. We are then permitted to take a
gradient expansion of the local observables around $(x,t)$. In particular,
averages of the current densities, $\mathfrak{j}_{i}(x,t)$ can be written in
the form,

$\mathfrak{j}_{i}(x,t)=\mathcal{F}_{i}(x,t)+\sum_{j}\mathcal{F}_{ij}(x,t)\partial_{x}\mathfrak{q}_{j}(x,t)+\dots,$
(21)

where the space-time dependence of the functions $\mathcal{F}_{i}$,
$\mathcal{F}_{ij}$, etc. is encoded through the charge densities
$\mathfrak{q}_{i}(x,t)$. Retaining only the first order term $\mathcal{F}_{i}$
corresponds to the so-called Euler-scale GHD. The next higher-order correction
is the diffusive Navier-Stokes correction. This term introduces an arrow of
time through the (irreversible) production of entropy. The corresponding
Navier-Stokes diffusive GHD equation was first derived by De Nardis et al.,
and reads

$\partial_{t}\rho_{p}+\partial_{x}\left(v^{\text{eff}}\rho_{p}\right)=\partial_{x}\left(\mathfrak{D}\partial_{x}\rho_{p}\right)+\frac{1}{m}\partial_{x}V\partial_{\theta}\rho_{p},$
(22)

where $\mathfrak{D}$ denotes the integral operator defined by the action,

$\mathfrak{D}f(\theta)=\int
d\theta^{\prime}\,\mathcal{D}(\theta,\theta^{\prime})f(\theta^{\prime}),$ (23)

and $\mathcal{D}(\theta,\theta^{\prime})$ is the diffusion kernel which
satisfies the relation,

$[\mathcal{D}(\theta,\cdot)]^{\text{dr}}(\theta^{\prime})\rho_{s}(\theta^{\prime})=[\rho_{s}(\cdot)\tilde{\mathcal{D}}(\cdot,\theta^{\prime})]^{\text{dr}}(\theta),$
(24)

with

$\displaystyle\tilde{\mathcal{D}}(\theta,\theta^{\prime})=\delta(\theta-\theta^{\prime})\left[\int
d\alpha\,\rho_{p}(\alpha)(1-\vartheta(\alpha))\right.$
$\displaystyle\,\,\,\,\,\,\quad\quad\quad\quad\times\left.\left(\frac{\Phi^{\text{dr}}(\alpha-\theta)}{\rho_{s}(\theta)}\right)^{2}\left|v^{\text{eff}}(\alpha)-v^{\text{eff}}(\theta)\right|\right]$
$\displaystyle\quad-\rho_{p}(\theta)(1-\vartheta(\theta))\left(\frac{\Phi^{\text{dr}}(\theta-\theta^{\prime})}{\rho_{s}(\theta)}\right)^{2}\left|v^{\text{eff}}(\theta)-v^{\text{eff}}(\theta^{\prime})\right|,$
(25)

which are all quantities from Euler-scale GHD [55]. The diffusion kernel
arises due to two-body scattering processes among quasi-particles [56].

Figure 3: Fermi contour of a system at zero temperature. The orange region
$\Gamma_{t}$ corresponds to a unit filling factor $\vartheta(x,\theta)=1$. The
dynamics of the system are reduced to that of the contour
$\partial\Gamma_{t}$. Locally, the region $\Gamma_{t}$ is split into disjoint
Fermi seas defined by the set of Fermi points
$\\{\theta_{1},\dots,\theta_{2n}\\}$. Adapted from [57].

Using diffusive GHD, Bastianello et al. studied the thermalisation of a 1D
Bose gas in the quantum Newton’s cradle experiment, see Fig. 2. The authors
found that diffusion was the leading mechanism that induced thermalisation
towards a stationary state in the presence of an (integrability-breaking)
trapping potential. The system initially relaxes to a pre-thermal state
described by a GGE, and then – at much longer times scales – it eventually
relaxes to the thermal state described by the standard Gibbs ensemble. These
results demonstrate the fundamental role that diffusive effects play in the
late-time dynamics of near-integrable systems.

Figure 4: Evolution of the rapidity distribution for a 1D Bose gas containing
130 atoms at $94$ nK in a quantum Newton’s cradle setup initiated with Bragg
pulses. The top row shows the evolution during the first 100 oscillation
periods according to the equations of standard GHD. The bottom row shows the
evolution according to the extended model of Møller et al that incorporates
population in transversely excited states and hence leads to faster
thermalization. The dashed lines mark the excitation threshold. The small
fraction of transverse excited atoms has a strong influence on the dynamics of
the system; the inclusion of the collision integral enables quasi-particles to
re-distribute across phase space and thermalize. The last two columns show the
best fit for a thermal state at the temperature of the final evolved system.
Adapted from [58].

###  III.2 Quantum Generalised Hydrodynamics

In its original formulation, GHD neglects important quantum effects, such as
quantum fluctuations and entanglement entropy. A fundamental assumption of GHD
is that each fluid cell comprising the system is independent at any fixed
time. This implies that equal-time correlations between fluid cells vanish.
However, as with many quantum systems, these equal-time correlations are
typically non-zero [59]; such effects occur beyond the Euler scale. One may
resolve these shortcomings by re-quantising the theory of GHD. This was first
achieved by Ruggiero et al. in 2020, in which they considered a system
initially at zero temperature, for which quantum fluctuations are most
significant [57]. Since the entropy is initially zero and is conserved under
the evolution of the standard GHD equations, the state remains at zero entropy
at all times.

Zero-entropy states are characterised by the so-called Fermi contour
$\partial\Gamma_{t}$ enclosing the region $\Gamma_{t}$ in $(x,\theta)$ phase
space with unit filling,

$\vartheta(x,\theta,t)=\begin{cases}1,&\text{if }(x,\theta)\in\Gamma_{t},\\\
0,&\text{otherwise}.\end{cases}$ (26)

Locally, a zero-entropy state can be described by a split Fermi sea,
consisting of a collection of disjoint regions in $(x,\theta)$ phase space for
which the filling factor is unity,

$\vartheta(x,\theta,t)=\begin{cases}1,&\text{if
}\theta\in[\theta_{1},\theta_{2}]\cup\dots\cup[\theta_{2p-1},\theta_{2p}]\\\
0,&\text{otherwise},\end{cases}$ (27)

where $\theta_{n}\equiv\theta_{n}(x,t)$, $n\in\\{1,2,\dots,2p\\}$ denote the
Fermi points at position $x$ and time $t$, see Fig. 3. For these zero-entropy
states, it was shown by Doyon et al. that the equations of GHD, Eq. (19),
reduce to a set of equations for the Fermi points [34],

$\partial_{t}\theta_{n}+v^{\text{eff}}\partial_{x}\theta_{n}=-\frac{1}{m}\partial_{x}V.$
(28)

One then considers small fluctuations about the Fermi contour. Locally, the
Fermi points are modified via $\theta_{n}\to\theta_{n}+\delta\theta_{n}$,
where $\delta\theta_{n}$ satisfies the differential equation,

$(\partial_{t}+v^{\text{eff}}\partial_{x})\delta\theta_{n}(x,t)=0,$ (29)

describing the propagation of linear sound waves [60]. To re-quantise this
theory, one promotes the fluctuations about the Fermi contour to linear
operators acting on a Hilbert space, i.e.,
$\delta\theta_{n}\to\delta\hat{\theta}_{n}$. The key insight of Ruggiero et
al. is that the problem of quantising sound waves above the classical GHD
ground state can be recast as a problem of quantising incompressible regions
in phase space, which is well-known in the literature on the quantum Hall
effect [61]. The resulting re-quantised theory is a time-dependent, spatially
inhomogeneous, multi-component Luttinger liquid, and has been dubbed “quantum
generalised hydrodynamics”. It describes quantum fluctuations of non-
equilibrium systems where conventional Luttinger liquid theory fails [57].
Additionally, unlike Luttinger liquid theory, this quantum GHD is not
restricted to low energies.

Figure 5: Measurements of the in situ density profile for a system of
$N=6300\pm 200$ bosons, initially trapped in a 1D double-well potential and
subsequently released to freely expand in 1D. These results (noisy lines) are
compared against the predictions of conventional (dashed lines) and
generalised hydrodynamics (solid lines). Conventional hydrodynamics
incorrectly predicts two distinct density peaks at late times. By contrast,
GHD shows excellent agreement with the experimental data. Adapted from [62].

###  III.3 Dimensional crossover for the 1D Bose Gas

One-dimensional atomic gases are often realised by using highly-elongated
“cigar-shaped” trapping potentials. It is assumed that due to the strong
radial confinement, the transverse modes are effectively frozen out. However,
when two atoms of sufficiently large momenta collide, the collisional energy
may exceed the level spacing of the transverse confinement. This leads to a
non-zero population of transverse excited states, which breaks integrability.
To account for collisions of atoms in transverse excited states, Møller et al.
introduced a phenomenological Boltzmann-type collision integral into the
equations of GHD [58].

Møller et al. studied the dynamics of a quasi-1D Bose gas in the quantum
Newton’s cradle scenario using this modified GHD. The evolution of the
rapidity distributions under the equations of standard GHD and the modified
GHD are shown in Fig. 4. By accounting for collisions with atoms in transverse
excited states, the modified GHD predicts that the system thermalizes. This is
in contrast to standard GHD which does not predict thermalization. The rate of
thermalization of this modified GHD was also found to be consistent with
previous experimental observations [63].

Figure 6: Measurements of the in situ density profile for a system of
$N=3500\pm 140$ atoms, quenched from a 1D double-well potential to a harmonic
potential [62]. These results are compared against the predictions of GHD
shown as smooth lines on top of noisy experimental data. While the standard
GHD reproduces most of the main features of the emergent dynamics, the minor
disagreement with experimental data at later times can be rectified [58] by
accounting for transverse excitations of the gas using the extension of GHD to
dimensional crossover as outlined in Section III.3. Adapted from [62].

###  III.4 Space-time inhomogeneous interactions

Figure 7: Comparison of experimental and GHD rapidity distributions of a
strongly interacting 1D Bose gas following a strong confinement quench. The
rapidity distributions are shown at different times, as indicated in each
panel. The red curves show the experimental data, whilst the blue lines are
the predictions of GHD, demonstrating good overall agreement between the
theory and experiment. Adapted from [64].

Inhomogeneities are ubiquitous in experiments. The effects of a (spatially-
varying) trapping potential were first introduced into GHD by Doyon and
Yoshimura in 2017 [42]. However, standard GHD does not address the effect of
space-time inhomogeneous interactions. Bastianello et al. were the first to
incorporate these effects into GHD. They found that these inhomogeneities
introduce an additional term into the standard GHD equations [65],

$\partial_{t}\rho_{p}+\partial_{x}\left(v^{\text{eff}}\rho_{p}\right)+\partial_{\theta}\left(\frac{f^{\text{dr}}\partial_{t}\alpha+\Lambda^{\text{dr}}\partial_{x}\alpha}{(\partial_{\theta}p)^{\text{dr}}}\right)=0,$
(30)

where $\alpha(x,t)$ denotes the interaction strength, and the forces $f$ and
$\Lambda$ are given by,

$\displaystyle f(\theta)$ $\displaystyle=\int
d\theta^{\prime}\,\frac{\partial\Phi(\theta-\theta^{\prime})}{\partial\alpha}\rho_{p}(\theta^{\prime}),$
(31) $\displaystyle\Lambda(\theta)$ $\displaystyle=\int
d\theta^{\prime}\,\frac{\partial\Phi(\theta-\theta^{\prime})}{\partial\alpha}v^{\text{eff}}(\theta^{\prime})\rho_{p}(\theta^{\prime}),$
(32)

These results exhaust all possible inhomogeneities which can be considered on
a purely hydrodynamical level [65].

##  IV Experimental Tests and Benchmarks of Generalised Hydrodynamics

To date, there have been three major experimental tests of GHD [62, 64, 58].
These have all been within the context of 1D Bose gases confined to highly
anisotropic trapping potentials as to realise the Lieb-Liniger model. In this
section, we outline two of these experiments [62, 64] and their descriptions
via GHD, while the third experiment [58] was already mentioned in Section
III.3 in the context of dimensional crossover in quasi-1D Bose gases.

###  IV.1 Tests of GHD in the weakly interacting regime

In 2019, Schemmer et al. experimentally demonstrated the validity of GHD for a
system of one-dimensional bosons realised on an atom chip [62]. There,
approximately $6300\pm 200$ particles were confined to one dimension via
magnetic trapping techniques developed by [66] and [67]. The particles were
initially confined in a double-well potential at thermal equilibrium and were
subsequently allowed to expand freely. The in situ density profile is shown at
three times during the evolution in Fig. 5. One finds that GHD provides an
excellent description of the evolving density profile, whilst conventional
hydrodynamics incorrectly predicts the formation of two sharp density peaks.

Schemmer et al. also considered the dynamics of a system of 1D bosons induced
via a quench from a double-well potential to a harmonic potential. They found
that GHD reproduces most of the main features of the emergent dynamics, see
Fig. 6. However, GHD predicts that the two density peaks persist much longer
than the experimental data. Schemmer et al. attributed this discrepancy to
atom losses in the experiment [68], such as three-body losses due to three-
body recombination processes [69, 70], that break integrability. It was found
that the total number of atoms in their system was reduced by 15% during the
dynamics. It has been argued by Møller et al. that transverse excitations may
have also influenced the observed dynamics [58]. By using the extension of GHD
to the dimensional crossover regime, as outlined in Section III.3, Møller et
al. found that the density peaks were less pronounced than in standard GHD, in
agreement with the experimental data.

###  IV.2 Tests of GHD in the strongly interacting regime

In the strongly interacting regime – the opposite regime to that of the atom
chip experiment of Schemmer et al. – Malvania et al. used a 2D lattice of
independent 1D gases to test GHD following a strong confinement quench [64].

Figure 8: Comparison of the post-quench rapidity distributions from GHD
(dashed lines) and from exact dynamics (solid lines) in the Tonks-Girardeau
limit of infinite interaction strength for small atom numbers. The red, green,
and blue lines correspond to 1D tubes with $5$, $10$, and $20$ atoms,
respectively. The black lines are averages of the other lines. Adapted from
[64].

In this experiment, bundles of one-dimensional gases of rubidium-87 atoms were
generated using optical potentials. The average number of atoms per gas was on
the order of 10. This ensures that the system is well within the strongly
interacting regime. The axial trapping potential was then dramatically
quenched, and the ensuing dynamics were measured. The evolution of the space-
integrated rapidity distribution is shown in Fig. 7. Note that the reported
rapidity distribution is an average over all 1D tubes in the lattice. It was
found that the dynamics of the gases were well-described by GHD. This is
surprising as the low number of particles challenges the fundamental
hydrodynamic assumption of GHD, which assumes that each fluid cell is
thermodynamically large.

To address the unreasonable predictive power of GHD at low particle numbers,
Malvania et al. numerically simulated the experiment for a 1D Bose gas with a
significantly stronger interaction strength, for which a comparison with exact
results is possible. They found that the dynamics predicted by GHD were
identical to that obtained exactly, up to small ripples (Friedel oscillations)
in the rapidity distribution, which arise at low particle numbers, see Fig. 8.

Figure 9: Dynamics of dispersive quantum shock waves forming in evolution of
an initial localized density bump in an otherwise uniform 1D Bose gas. Shown
are snapshots of the evolving density profile at different dimensionless times
$\tau=t\hbar/mL^{2}$, where $L$ is the length of the system, for different
total atom numbers $N$ and dimensionless interaction strength
$\gamma_{\mathrm{bg}}$ in the background. Due to the reflectional symmetry
about the origin, we only show the densities for $x>0$. See text for further
discussion. Adapted from [71].

###  IV.3 Benchmarks of generalised hydrodynamics

Recently, Watson et al. [71] benchmarked the performance of GHD against a
variety of theoretical approaches for an interacting 1D Bose gas in four
dynamical scenarios. Namely, the free propagation of a localised density bump
and dip induced via a confinement quench and known to lead to dispersive
quantum shock waves [72], and the quantum Newton’s cradle setup via a quench
of a double-well trap to a single-well and via a Bragg pulse.

Examples of simulations in the dispersive quantum shock wave scenario
initiated from a localised density bump are shown in Figs. 9 (a)–(d) for a
weakly interacting regime of the 1D Bose gas, $\gamma_{\mathrm{bg}}\ll 1$, and
in Fig. 9 (e) in the Tonks-Girardeau regime of infinitely strong interactions
($\gamma_{\mathrm{bg}}\to\infty$). Here,
$\gamma_{\mathrm{bg}}=mg/\hbar^{2}\rho_{\mathrm{bg}}$ is the dimensionless
interaction strength at background density $\rho_{\mathrm{bg}}$. In (a) and
(b), the predictions of GHD are compared with the mean-field Gross-Pitaevskii
equation (GPE), and the truncated Wigner approach (TWA) which takes into
account the effect of quantum fluctuations on top of the mean-filed
description. In (c) the prediction of GHD at time $\tau=0.0007$ is compared
with coarse-grained averages of the GPE and TWA results and show good
agreement. Such coarse-graining mimics finite imaging resolution in ultracold
atom experiments and may explain the success of GHD when compared to
experiments, even though its predictions may depart from other theoretical
approaches that are valid at short wavelengths.

In Fig. 9 (d), the GHD results at $\tau=0.0007$ are compared with the
predictions of stochastic projected Gross-Pitaevskii equation (SPGPE) for a
finite-temperature initial system, for two different initial dimensionless
temperatures $\overline{\mathcal{T}}$ (where
$\overline{\mathcal{T}}\\!=\\!T/T_{d}$, with
$T_{d}\\!=\\!\hbar^{2}\rho_{\mathrm{bg}}^{2}/2mk_{B}$ being the temperature of
quantum degeneracy), again showing excellent agreement. Finally, in (d) we
compare a snapshot of the density profile at $\tau=0.00004$ (evolved from the
same initial density bump) in the Tonks-Girardeau regime calculated using GHD
and exact diagonalization (ED) of an equivalent free-fermion problem to which
a Tonks-Girardeau gas can be mapped. Apart from missing density oscillations
on short (microscopic) length scales, GHD agrees very well with the ED result
in this strongly interacting regime as well.

For intermediate interactions strengths, Watson et al. [71] have also compared
the GHD results for the evolution of an initial bump with the results of the
exact matrix product state (MPS) method for $N=50$ particles. They found that
the performance of GHD generally improves with stronger interactions and
higher temperatures as short-wavelength interference phenomena become
increasingly suppressed due to the reduced phase coherence in the system. On
the other hand, GHD fails to capture interference phenomena at short-length
scales for very low temperatures and very weak interaction strengths, where
the local density approximation, intrinsic to any hydrodynamic approach, also
fails.

Similar conclusions have been obtained from simulations of an initial density
dip (see Fig. 10), known to shed a train of gray solitons in the weakly
interacting regime in the mean-field GPE approximation. Here, the situation is
similar to that of Figs. 9 (b) and (c): while GHD fails to capture individual
solitons, whose characteristic width (on the order of the microscopic healing
length) lies beyond the intended range of applicability of any hydrodynamic
theory, it adequately captures the coarse-grained average density over the
soliton train.

Figure 10: Same as in Fig. 9, except for the evolution of an initial density
dip. Panel (b), which is snapshot of density profile at time $\tau=0.002$, a
train of three gray solitons can be clearly identified in the mean-field GPE
results; their visibility diminishes once quantum fluctuations are taken into
account through the TWA approach. While GHD fails to capture individual
soliton profiles or short-wavelength structures, it agrees well with the
coarse-grained averages of the GPE and TWA results, shown in panels (c) and
(d) at $\tau=0.0005$ and $\tau=0.002$, respectively. Panel (e) is a comparison
of GHD and ED results in the Tonks-Girardeau regime, showing excellent
agreement. Finally, panel (f) shows comparison of GHD with exact MPS
calculations for $N=42$ in the nearly ideal Bose gas regime, where LDA is not
applicable, hence explaining the disagreement of GHD with MPS results. Adapted
from [71].

In the benchmarks of GHD in the quantum Newton’s cradle setup, performed for a
double- well to single-well trap quench of a weakly interacting
quasicondensate, Watson et al. [71] observed excellent agreement with the
SPGPE results in both the transient dynamics and the final relaxed state, as
well as in the overall characteristic thermalization timescale. However, in
the comparison for the quantum Newton’s cradle initiated by Bragg pulses, they
observed disagreement with SPGPE in the characteristic thermalization
timescale [73], with the GHD generally predicting slower thermalization. This
discrepancy, however, can be attributed to qualitatively very different ways
that the Bragg pulses are implemented in GHD and the SPGPE (see Ref. [71] for
further details).

##  V Conclusion and Perspectives

Generalised hydrodynamics has recently emerged as a broadly applicable
hydrodynamic theory for modelling the quantum many-body dynamics of integrable
and near-integrable systems on large scale. Since its inception, GHD has been
extended to account for various experimentally relevant effects such as
diffusion, dimensional crossover in the 1D Bose gas, inhomogeneous
interactions, and quantum effects. Experimentally, tests of GHD are still in
their infancy. However, the first few experimental investigations show
excellent agreement between the observed results and GHD predictions.

The unreasonable effectiveness of GHD at low particle numbers in [64] warrants
further investigation, as do tunnelling effects between 1D tubes. Tunnelling
between tubes breaks integrability and is therefore expected to lead to
complete thermalisation over sufficiently large time scales. Future
experimental tests of GHD would benefit by exploring more challenging
scenarios that push the assumptions of GHD, such as by straying further from
integrability via long-range dipolar interactions. Experimental verification
of diffusive GHD and hence characteristic thermalization rates [74, 73, 71] in
different regimes of the 1D Bose gas would also be an important achievement.

In the context of verifying the required conditions for GHD, we also mention
the recent work by Le et al. [75], which probed the rapid onset of
hydrodynamics – referred to as hydrodynamization – in an array of 1D Bose
gases in a strongly interacting regime. Hydrodynamization precedes local
prethermalization, and it was conjectured that GHD can be applied immediately
after hydrodynamization, i.e., during local prethermalization, even though the
local GGE is not yet established. If confirmed, this conjecture would imply a
further relaxation of the conditions required for applicability of GHD, which
in turn would explain the success of GHD under a broader range of experimental
conditions.

While we have reviewed GHD in the context of a repulsive 1D Bose gas, GHD is
also applicable to the attractive 1D Bose gas [76] and to other integrable
models such as the Hubbard model [77], the XXZ chain [36], and the Yang-Gaudin
model, describing the spin–$1/2$ Fermi gas [78, 79, 80]. As such, it would be
beneficial to see experimental verification of GHD in these and other
condensed matter systems. Finally, it would be interesting to use GHD to model
the non-equilibrium dynamics of more general integrable systems such as the
multi-component Bose gas [81] and Bose-Fermi mixtures [82].

###### Acknowledgements.

K. V. K. acknowledges support by the Australian Research Council Discovery
Project Grant DP190101515.

## References

  * Cazalilla and Rigol [2010] M. A. Cazalilla and M. Rigol, Focus on dynamics and thermalization in isolated quantum many-body systems, New Journal of Physics 12, 055006 (2010).
  * Polkovnikov  _et al._ [2011] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Colloquium: Nonequilibrium dynamics of closed interacting quantum systems, Rev. Mod. Phys. 83, 863 (2011).
  * Eisert  _et al._ [2015] J. Eisert, M. Friesdorf, and C. Gogolin, Quantum many-body systems out of equilibrium, Nature Physics 11, 124 (2015).
  * Durnin  _et al._ [2021] J. Durnin, M. J. Bhaseen, and B. Doyon, Nonequilibrium dynamics and weakly broken integrability, Phys. Rev. Lett. 127, 130601 (2021).
  * Lieb and Liniger [1963] E. H. Lieb and W. Liniger, Exact analysis of an interacting Bose gas. I. the general solution and the ground state, Phys. Rev. 130, 1605 (1963).
  * Lieb [1963] E. H. Lieb, Exact analysis of an interacting Bose gas. ii. the excitation spectrum, Phys. Rev. 130, 1616 (1963).
  * Hubbard [1963] J. Hubbard, Electron correlations in narrow energy bands, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 276, 238 (1963).
  * Gogolin and Eisert [2016] C. Gogolin and J. Eisert, Equilibration, thermalisation, and the emergence of statistical mechanics in closed quantum systems, Reports on Progress in Physics 79, 056001 (2016).
  * Rigol [2009] M. Rigol, Breakdown of thermalization in finite one-dimensional systems, Phys. Rev. Lett. 103, 100403 (2009).
  * Rigol  _et al._ [2007] M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii, Relaxation in a completely integrable many-body quantum system: An ab initio study of the dynamics of the highly excited states of 1D lattice hard-core bosons, Phys. Rev. Lett. 98, 050405 (2007).
  * Rigol  _et al._ [2008] M. Rigol, V. Dunjko, and M. Olshanii, Thermalization and its mechanism for generic isolated quantum systems, Nature 452, 854 (2008).
  * Langen  _et al._ [2015] T. Langen, S. Erne, R. Geiger, B. Rauer, T. Schweigler, M. Kuhnert, W. Rohringer, I. E. Mazets, T. Gasenzer, and J. Schmiedmayer, Experimental observation of a generalized Gibbs ensemble, Science 348, 207 (2015).
  * Calabrese and Cardy [2006] P. Calabrese and J. Cardy, Time dependence of correlation functions following a quantum quench, Phys. Rev. Lett. 96, 136801 (2006).
  * Calabrese and Cardy [2016] P. Calabrese and J. Cardy, Quantum quenches in $1+1$ dimensional conformal field theories, Journal of Statistical Mechanics: Theory and Experiment 2016, 064003 (2016).
  * Caux [2016] J.-S. Caux, The quench action, Journal of Statistical Mechanics: Theory and Experiment 2016, 064006 (2016).
  * De Nardis  _et al._ [2014] J. De Nardis, B. Wouters, M. Brockmann, and J.-S. Caux, Solution for an interaction quench in the Lieb-Liniger Bose gas, Phys. Rev. A 89, 033601 (2014).
  * Essler and Fagotti [2016] F. H. Essler and M. Fagotti, Quench dynamics and relaxation in isolated integrable quantum spin chains, Journal of Statistical Mechanics: Theory and Experiment 2016, 064002 (2016).
  * Cazalilla  _et al._ [2012] M. A. Cazalilla, A. Iucci, and M.-C. Chung, Thermalization and quantum correlations in exactly solvable models, Phys. Rev. E 85, 011133 (2012).
  * Ilievski  _et al._ [2016] E. Ilievski, M. Medenjak, T. Prosen, and L. Zadnik, Quasilocal charges in integrable lattice systems, Journal of Statistical Mechanics: Theory and Experiment 2016, 064008 (2016).
  * Loskutov [2007] A. Y. Loskutov, Dynamical chaos: systems of classical mechanics, Physics-Uspekhi 50, 939 (2007).
  * Mitra [2018] A. Mitra, Quantum quench dynamics, Annual Review of Condensed Matter Physics 9, 245 (2018).
  * Landau [1941] L. Landau, Theory of the superfluidity of helium II, Phys. Rev. 60, 356 (1941).
  * Levitov and Falkovich [2016] L. Levitov and G. Falkovich, Electron viscosity, current vortices and negative nonlocal resistance in graphene, Nature Physics 12, 672 (2016).
  * Moreau [1990] R. J. Moreau,  _Magnetohydrodynamics_ , Vol. 3 (Springer Science & Business Media, Berlin, 1990).
  * Kumanek and Janas [2019] B. Kumanek and D. Janas, Thermal conductivity of carbon nanotube networks: a review, Journal of materials science 54, 7397 (2019).
  * Pigeon  _et al._ [2011] S. Pigeon, I. Carusotto, and C. Ciuti, Hydrodynamic nucleation of vortices and solitons in a resonantly excited polariton superfluid, Phys. Rev. B 83, 144513 (2011).
  * Reichel and Vuletic [2011] J. Reichel and V. Vuletic,  _Atom chips_ (John Wiley & Sons, Berlin, 2011).
  * Bloch  _et al._ [2008] I. Bloch, J. Dalibard, and W. Zwerger, Many-body physics with ultracold gases, Rev. Mod. Phys. 80, 885 (2008).
  * Bouchoule  _et al._ [2011] I. Bouchoule, N. Van Druten, and C. I. Westbrook, Atom chips and one-dimensional Bose gases, Atom chips , 331 (2011).
  * Yang and Yang [1969] C.-N. Yang and C. P. Yang, Thermodynamics of a one-dimensional system of bosons with repulsive delta-function interaction, Journal of Mathematical Physics 10, 1115 (1969).
  * Bouchoule  _et al._ [2016] I. Bouchoule, S. S. Szigeti, M. J. Davis, and K. V. Kheruntsyan, Finite-temperature hydrodynamics for one-dimensional Bose gases: Breathing-mode oscillations as a case study, Phys. Rev. A 94, 051602 (2016).
  * Menotti and Stringari [2002] C. Menotti and S. Stringari, Collective oscillations of a one-dimensional trapped Bose-Einstein gas, Phys. Rev. A 66, 043610 (2002).
  * Stringari [1996] S. Stringari, Collective excitations of a trapped Bose-condensed gas, Phys. Rev. Lett. 77, 2360 (1996).
  * Doyon  _et al._ [2017] B. Doyon, J. Dubail, R. Konik, and T. Yoshimura, Large-scale description of interacting one-dimensional Bose gases: Generalized hydrodynamics supersedes conventional hydrodynamics, Phys. Rev. Lett. 119, 195301 (2017).
  * Castro-Alvaredo  _et al._ [2016] O. A. Castro-Alvaredo, B. Doyon, and T. Yoshimura, Emergent hydrodynamics in integrable quantum systems out of equilibrium, Phys. Rev. X 6, 041065 (2016).
  * Bertini  _et al._ [2016] B. Bertini, M. Collura, J. De Nardis, and M. Fagotti, Transport in out-of-equilibrium ${XXZ}$ chains: Exact profiles of charges and currents, Phys. Rev. Lett. 117, 207201 (2016).
  * Note [1] The term “generalised” in GHD is used in the same sense as in generalised Gibbs ensemble.
  * Caux  _et al._ [2019] J.-S. Caux, B. Doyon, J. Dubail, R. Konik, and T. Yoshimura, Hydrodynamics of the interacting Bose gas in the quantum Newton cradle setup, SciPost Phys. 6, 70 (2019).
  * Kheruntsyan  _et al._ [2005] K. V. Kheruntsyan, D. M. Gangardt, P. D. Drummond, and G. V. Shlyapnikov, Finite-temperature correlations and density profiles of an inhomogeneous interacting one-dimensional Bose gas, Phys. Rev. A 71, 053615 (2005).
  * van Tongeren [2016] S. J. van Tongeren, Introduction to the thermodynamic Bethe ansatz, Journal of Physics A: Mathematical and Theoretical 49, 323005 (2016).
  * Note [2] In this review, we only consider many-body models of identical particles of single species, and hence, a single quasi-particle distribution. It is straightforward to generalise the equations of GHD to many species.
  * Doyon and Yoshimura [2017] B. Doyon and T. Yoshimura, A note on generalized hydrodynamics: inhomogeneous fields and other concepts, SciPost Phys. 2, 014 (2017).
  * Doyon [2020] B. Doyon, Lecture notes on generalised hydrodynamics, SciPost Physics Lecture Notes (2020).
  * Percus [1976] J. Percus, Equilibrium state of a classical fluid of hard rods in an external field, Journal of Statistical Physics 15, 505 (1976).
  * Boldrighini  _et al._ [1983] C. Boldrighini, R. Dobrushin, and Y. M. Sukhov, One-dimensional hard rod caricature of hydrodynamics, Journal of Statistical Physics 31, 577 (1983).
  * El and Kamchatnov [2005] G. El and A. M. Kamchatnov, Kinetic equation for a dense soliton gas, Phys. Rev. Lett. 95, 204101 (2005).
  * El [2003] G. El, The thermodynamic limit of the Whitham equations, Physics Letters A 311, 374 (2003).
  * Cao  _et al._ [2019] X. Cao, V. B. Bulchandani, and H. Spohn, The GGE averaged currents of the classical Toda chain, Journal of Physics A: Mathematical and Theoretical 52, 495003 (2019).
  * Møller and Schmiedmayer [2020] F. S. Møller and J. Schmiedmayer, Introducing iFluid: a numerical framework for solving hydrodynamical equations in integrable models, SciPost Physics 8, 10.21468/scipostphys.8.3.041 (2020).
  * Kinoshita  _et al._ [2006] T. Kinoshita, T. Wenger, and D. S. Weiss, A quantum Newton’s cradle, Nature 440, 900 (2006).
  * Steinigeweg  _et al._ [2017] R. Steinigeweg, F. Jin, H. De Raedt, K. Michielsen, and J. Gemmer, Charge diffusion in the one-dimensional Hubbard model, Phys. Rev. E 96, 020105 (2017).
  * Sirker  _et al._ [2009] J. Sirker, R. G. Pereira, and I. Affleck, Diffusion and ballistic transport in one-dimensional quantum systems, Phys. Rev. Lett. 103, 216602 (2009).
  * Ljubotina  _et al._ [2017] M. Ljubotina, M. Žnidarič, and T. Prosen, Spin diffusion from an inhomogeneous quench in an integrable system, Nature communications 8, 1 (2017).
  * Bastianello  _et al._ [2020a] A. Bastianello, A. De Luca, B. Doyon, and J. De Nardis, Thermalization of a trapped one-dimensional Bose gas via diffusion, Phys. Rev. Lett. 125, 240604 (2020a).
  * De Nardis  _et al._ [2018] J. De Nardis, D. Bernard, and B. Doyon, Hydrodynamic diffusion in integrable systems, Phys. Rev. Lett. 121, 160603 (2018).
  * Nardis  _et al._ [2019] J. D. Nardis, D. Bernard, and B. Doyon, Diffusion in generalized hydrodynamics and quasiparticle scattering, SciPost Phys. 6, 49 (2019).
  * Ruggiero  _et al._ [2020] P. Ruggiero, P. Calabrese, B. Doyon, and J. Dubail, Quantum generalized hydrodynamics, Phys. Rev. Lett. 124, 140603 (2020).
  * Møller  _et al._ [2021] F. Møller, C. Li, I. Mazets, H.-P. Stimming, T. Zhou, Z. Zhu, X. Chen, and J. Schmiedmayer, Extension of the generalized hydrodynamics to the dimensional crossover regime, Phys. Rev. Lett. 126, 090602 (2021).
  * Cazalilla  _et al._ [2011] M. A. Cazalilla, R. Citro, T. Giamarchi, E. Orignac, and M. Rigol, One dimensional bosons: From condensed matter systems to ultracold gases, Rev. Mod. Phys. 83, 1405 (2011).
  * Ruggiero  _et al._ [2021] P. Ruggiero, P. Calabrese, B. Doyon, and J. Dubail, Quantum generalized hydrodynamics of the Tonks–Girardeau gas: density fluctuations and entanglement entropy, Journal of Physics A: Mathematical and Theoretical 55, 024003 (2021).
  * Wen [1990] X. G. Wen, Chiral Luttinger liquid and the edge excitations in the fractional quantum Hall states, Phys. Rev. B 41, 12838 (1990).
  * Schemmer  _et al._ [2019] M. Schemmer, I. Bouchoule, B. Doyon, and J. Dubail, Generalized hydrodynamics on an atom chip, Phys. Rev. Lett. 122, 090601 (2019).
  * Li  _et al._ [2020] C. Li, T. Zhou, I. Mazets, H.-P. Stimming, F. S. Møller, Z. Zhu, Y. Zhai, W. Xiong, X. Zhou, X. Chen, and J. Schmiedmayer, Relaxation of bosons in one dimension and the onset of dimensional crossover, SciPost Phys. 9, 058 (2020).
  * Malvania  _et al._ [2021] N. Malvania, Y. Zhang, Y. Le, J. Dubail, M. Rigol, and D. S. Weiss, Generalized hydrodynamics in strongly interacting 1D Bose gases, Science 373, 1129 (2021).
  * Bastianello  _et al._ [2019] A. Bastianello, V. Alba, and J.-S. Caux, Generalized hydrodynamics with space-time inhomogeneous interactions, Phys. Rev. Lett. 123, 130602 (2019).
  * Reichel  _et al._ [1999] J. Reichel, W. Hänsel, and T. W. Hänsch, Atomic micromanipulation with magnetic surface traps, Phys. Rev. Lett. 83, 3398 (1999).
  * Folman  _et al._ [2000] R. Folman, P. Krüger, D. Cassettari, B. Hessmo, T. Maier, and J. Schmiedmayer, Controlling cold atoms using nanofabricated surfaces: Atom chips, Phys. Rev. Lett. 84, 4749 (2000).
  * Bouchoule and Dubail [2022] I. Bouchoule and J. Dubail, Generalized hydrodynamics in the one-dimensional Bose gas: theory and experiments, Journal of Statistical Mechanics: Theory and Experiment 2022, 014003 (2022).
  * Söding  _et al._ [1999] J. Söding, D. Guéry-Odelin, P. Desbiolles, F. Chevy, H. Inamori, and J. Dalibard, Three-body decay of a rubidium Bose–Einstein condensate, Applied physics B 69, 257 (1999).
  * Tolra  _et al._ [2004] B. L. Tolra, K. M. O’Hara, J. H. Huckans, W. D. Phillips, S. L. Rolston, and J. V. Porto, Observation of reduced three-body recombination in a correlated 1D degenerate Bose gas, Phys. Rev. Lett. 92, 190401 (2004).
  * Watson  _et al._ [2022] R. Watson, S. Simmons, and K. Kheruntsyan, Benchmarks of generalized hydrodynamics for 1D Bose gases, arXiv preprint arXiv:2208.06614 10.48550/arXiv.2208.06614 (2022).
  * Simmons  _et al._ [2020] S. A. Simmons, F. A. Bayocboc, J. C. Pillay, D. Colas, I. P. McCulloch, and K. V. Kheruntsyan, What is a quantum shock wave?, Phys. Rev. Lett. 125, 180401 (2020).
  * Thomas  _et al._ [2021] K. F. Thomas, M. J. Davis, and K. V. Kheruntsyan, Thermalization of a quantum Newton’s cradle in a one-dimensional quasicondensate, Phys. Rev. A 103, 023315 (2021).
  * Bastianello  _et al._ [2020b] A. Bastianello, A. De Luca, B. Doyon, and J. De Nardis, Thermalization of a trapped one-dimensional Bose gas via diffusion, Phys. Rev. Lett. 125, 240604 (2020b).
  * Le  _et al._ [2022] Y. Le, Y. Zhang, S. Gopalakrishnan, M. Rigol, and D. S. Weiss, Direct observation of hydrodynamization and local prethermalization, arXiv preprint arXiv:2210.07318 10.48550/arXiv.2210.07318 (2022).
  * Koch  _et al._ [2022] R. Koch, J.-S. Caux, and A. Bastianello, Generalized hydrodynamics of the attractive non-linear Schrödinger equation, Journal of Physics A: Mathematical and Theoretical 55, 134001 (2022).
  * Nozawa and Tsunetsugu [2020] Y. Nozawa and H. Tsunetsugu, Generalized hydrodynamic approach to charge and energy currents in the one-dimensional Hubbard model, Phys. Rev. B 101, 035121 (2020).
  * Gaudin [1967] M. Gaudin, Un systeme a une dimension de fermions en interaction, Physics Letters A 24, 55 (1967).
  * Scopa  _et al._ [2022] S. Scopa, P. Calabrese, and L. Piroli, Generalized hydrodynamics of the repulsive spin-$\frac{1}{2}$ fermi gas, Phys. Rev. B 106, 134314 (2022).
  * Yang [1967] C. N. Yang, Some exact results for the many-body problem in one dimension with repulsive delta-function interaction, Phys. Rev. Lett. 19, 1312 (1967).
  * Fuchs  _et al._ [2005] J. N. Fuchs, D. M. Gangardt, T. Keilmann, and G. V. Shlyapnikov, Spin waves in a one-dimensional spinor Bose gas, Phys. Rev. Lett. 95, 150402 (2005).
  * Pâţu and Klümper [2019] O. I. Pâţu and A. Klümper, Momentum reconstruction and contact of the one-dimensional Bose-Fermi mixture, Phys. Rev. A 99, 013628 (2019).


MI-HET-782

# Search Prospects for Axion-like Particles at
Rare Nuclear Isotope Accelerator Facilities

Wooyoung Jang  Department of Physics, University of Texas, Arlington, TX
76019, USA  Doojin Kim  Mitchell Institute for Fundamental Physics and
Astronomy, Department of Physics and Astronomy, Texas A&M University, College
Station, TX 77843, USA  Kyoungchul Kong  Department of Physics and Astronomy,
University of Kansas, Lawrence, KS 66045, USA  Youngjoon Kwon  Department of
Physics & IPAP, Yonsei University, Seoul 03722, Republic of Korea  Jong-Chul
Park  Department of Physics and IQS, Chungnam National University, Daejeon
34134, Republic of Korea  Min Sang Ryu  Department of Physics, University of
Seoul, Seoul 02504, Republic of Korea  Department of Physics, Kyungpook
National University, Daegu 41566, Republic of Korea  Seodong Shin  Department
of Physics, Jeonbuk National University, Jeonju, Jeonbuk 54896, Republic of
Korea  Richard G. Van de Water  Los Alamos National Laboratory, Los Alamos, NM
87545, USA  Un-Ki Yang  Department of Physics, Seoul National University,
Seoul 08826, Republic of Korea  Jaehoon Yu  Department of Physics, University
of Texas, Arlington, TX 76019, USA

###### Abstract

We propose a novel experimental scheme, called DAMSA (Dump-produced Aboriginal
Matter Searches at an Accelerator), for searching for dark-sector particles,
using rare nuclear isotope accelerator facilities that provide high-flux
proton beams to produce a large number of rare nuclear isotopes. The high-
intensity nature of their beams enables the investigation of dark-sector
particles, including axion-like particles (ALPs) and dark photons. By
contrast, their typical beam energies are not large enough to produce the
backgrounds such as neutrinos resulting from secondary charged particles. The
detector of DAMSA is then placed immediate downstream of the proton beam dump
to maximize the prompt decay signals of dark-sector particles, which are often
challenging to probe in other beam-dump-type experiments featuring a longer
baseline, at the expense of an enormous amount of the beam-related neutron
(BRN) background. We demonstrate that BRN can be significantly suppressed if
the signal accompanies multiple, correlated visible particles in the final
state. As an example physics case, we consider ALPs interacting with the
Standard Model photon and their diphoton decay signal at DAMSA implemented at
a rare nuclear isotope facility similar to the Rare isotope Accelerator
complex for ON-line experiment under construction in South Korea. We show that
the close proximity of the detector to the ALP production dump makes it
possible to probe a high-mass region of ALP parameter space that the existing
experiments have never explored.

Introduction. Dark matter or more generally dark-sector physics is a well-
motivated paradigm for extensions of the Standard Model (SM), and a wide range
of models have been proposed, predicting various new particles. Of them, new
particles including axion-like particles (ALPs) and dark photon are receiving
particular attention, as they often mediate the interactions between SM and
dark-sector particles. Their feebly interacting nature, however, makes it
challenging for energy-frontier experiments to discover such mediators; it is
more probable that the first discovery of them would be made in an experiment
at the facilities with high-intensity beams, such as a beam-dump type
experiment.

A number of ongoing beam-dump type experiments mostly focus on the decay of
mediators to SM particles and have set stringent limits in their parameter
space. By construction, once a mediator is produced in a target, it should be
sufficiently long-lived to enter the detector before decaying to lighter
particles. It, however, should then decay inside the detector to be detected.
The balance between its survival and decay governs the shape of the
sensitivity reaches at a given experiment, highly depending on the details of
the experimental configuration such as the detector geometric coverage, the
distance between the production location and the detector, etc. While higher
beam intensity enables expanding the coverage of parameter space, the region
associated with high-mass mediators with large coupling strengths is
challenging to explore due to the prompt decay of the mediators; the access to
the high-mass range of the parameter space is determined by the distance
between the mediator production vertex and the detector.

In this regard, we point out that high-intensity proton beam facilities for
rare nuclear isotope production (e.g., FRIB Wei:2019oew , ISOLDE
Catherall:2017kbr , and RAON Jeon:2014oja ) possess the possibility of placing
a precision large-scale detector immediate downstream of the dump, greatly
reducing the distance between the source and the detector. Since the available
center-of-mass frame energy (e.g., $\sim 280$ MeV from a 600-MeV beam
collision on a target) is barely above twice the pion mass, the flux of
charged pions is significantly reduced and so is that of neutrinos from their
decay, leaving the beam-related neutrons (BRN) to be the major background
which can be significantly suppressed for the signals with features.

Exploiting these observations, we propose an experimental scheme, named Dump-
produced Aboriginal Matter Searches at an Accelerator or DAMSA,111 In Korean
characters mj담사, /dms/: deep thought, rumination. for the dark-sector
investigation. For illustration purposes, we consider a RAON-like facility
Jeon:2014oja  and study the detection prospect of ALPs as a benchmark physics
case for demonstrating capability of DAMSA in the search for dark-sector
physics and how the BRN can be evaded.

Indeed, the low-energy, high-intensity, beam facilities are particularly
suited for probing ALPs interacting with SM photons. First of all, a highly
intensified beam can produce photons copiously in the beam dump through
various production mechanisms including bremsstrahlung of secondary charged
particles Dutta:2020vop ; Brdar:2020dpr . These photons can be converted to
ALPs via the Primakoff effect, and thus an enormous flux of ALPs is generally
expected. Second, since ALPs can decay into a pair of photons which are
massless, the relevant experiment can probe even low-mass ALPs (say lower than
twice the electron mass) in the decay channel, fully exploiting the photon
flux as long as kinematics is allowed. Third, as discussed earlier, the
expected (neutrino-induced) background rates are small, compared to those
utilizing a higher-energy beam. Therefore, we expect that the DAMSA experiment
at RAON-like facilities can provide great opportunities in terms of ALP
searches.

Figure 1: A conceptual layout of an example DAMSA experiment implemented at a
RAON-like facility.

Experimental Concept. We envision the situation that a 660 $\mu$A, 600-MeV
proton beam provided by a RAON-like facility impinges on a 1-meter long
cylindrical tungsten target, immediately followed by a polyurethane neutron
moderator and a vacuum decay chamber whose volume is defined by a 10-meter
long cylinder and whose wall is framed by 0.6-cm thick steel plates (see Fig.
1). Electromagnetic particles are detected for identification at a high-
granularity, high resolution ECAL located immediate downstream behind the
decay chamber. We assume the ECAL angular coverage $\theta_{\rm det}$ to be as
large as 0.5 radians to improve the ALP signal sensitivity, i.e., when a line
emanating from the center of the top base of the target intersects the beam
axis by $\theta_{\rm det}=0.5$ rad. We further assume that the energy
threshold for photons $E_{\rm th}$ and the angular resolution $\theta_{\rm
res}$ should be as low as 15 MeV and $1^{\circ}$, respectively, not only to
maximize the ALP signal sensitivity but to achieve significant background
rejection.

Given the models of ALP interacting with photons, it is crucial to estimate
spectral information of photons emerging in the target as precisely as
possible for precise estimates of the ALP-signal sensitivity. To this end, we
employ a dedicated detector-level simulation code package, GEANT4, with the
QGSP_BIC_AllHP library Agostinelli:2002hh . We deliver $10^{6}$ 600-MeV
protons to the target and then find $\sim 1.9\times 10^{8}$ photons inside the
target. For simplicity, we considered the target geometry and no other
components which would be part of the target complex, expecting that such
details would make an insignificant impact on our results.

Figure 2: Two-dimensional correlation between energy and $\cos\theta_{\gamma}$
of photons produced in the target of the proposed experiment at a RAON-like
facility. The angle is measured with respect to the beam axis.

We display the correlation between the energy and the $\cos\theta_{\gamma}$ of
photons produced in the target in Fig. 2. First, soft photons dominate the
energy spectrum. Nevertheless, a non-negligible fraction of photons with
energy greater than 100 MeV can be produced, allowing DAMSA to be sensitive to
$\mathcal{O}(100)$ MeV ALPs and set the new limits for them. Second, the
angular spectrum suggests that such energetic photons move preferentially in
the forward direction. Given the fact that our detector system covers up to
$\theta_{\rm det}=0.5$ rad in the forward region, we find that about 7.5% of
photons will move toward the detector system.

ALP Signal. We consider models of ALPs interacting with the SM photon for
which the Lagrangian contains $\mathcal{L}_{\rm int}\supset-
g_{a\gamma\gamma}\,aF_{\mu\nu}\tilde{F}^{\mu\nu}/4$ where $g_{a\gamma\gamma}$
parameterizes the effective coupling strength between the ALP and the SM
photon and $F_{\mu\nu}$ ($\tilde{F}_{\mu\nu}$) denotes the (dual) field
strength tensor of the SM photon. The presence of the above interaction allows
a photon to convert to an ALP via the Primakoff process Primakoff:1951iae ,
$\gamma+A\to a+A$ with $A$ being the atomic system of interest. Once produced,
an ALP flies to the detector and can decay into two photons inside the decay
chamber, $a\to 2\gamma$. The decay width of ALP to a photon pair is well known
and expressed as $\Gamma_{a}=g_{a\gamma\gamma}^{2}m_{a}^{3}/(64\pi)$ with
$m_{a}$ being the ALP mass. Obviously, an ALP decays rather promptly if
$m_{a}$ and/or $g_{a\gamma\gamma}$ are large. Therefore, it is highly
beneficial to place the detector as close to the target as possible to achieve
improved sensitivities to decay signals of (relatively) short-lived ALPs. In
addition, an expected high flux of ALPs at DAMSA enables us to probe new
regions of parameter space.

The data analysis strategy for the ALP signal of this sort is well known and
standardized in literature (see, for example, Refs. Dent:2019ueq ;
Brdar:2020dpr  and a review article Fortin:2021cog ). We follow the same
strategy, modifying detector-specific parameters accordingly (see also
Appendix A for details).

Finally, we discuss kinematic distributions of ALP events at the detector, as
they are related to the detector requirements and the cuts to discriminate the
signal from potential backgrounds that we discuss shortly. Since the signal of
our interest involves a pair of photons in the final state, we display spectra
of $E_{\gamma,{\rm det}}$, the energy of individual photons, and spectra of
$\theta_{\gamma\gamma,{\rm det}}$, the angular separation between the two
photons in the left and right panels of Fig. 3.

Figure 3: Left: Energy spectra of the photons from the ALP decay for three ALP
mass values. The black dashed vertical line marks $E_{\rm th}=15$ MeV. Right:
Spectra in the angular separation between the two final-state photons for the
same three mass points.

We remark that the energy spectra develop a log-symmetric structure with
respect to the energy value being half the mass of ALP, since the ALP is a
spinless particle and decays to two massless photons Agashe:2012bn ;
Boddy:2016hbp . This implies that lower-mass ALP signals are more affected by
the energy threshold or cut in their acceptance, since one of the two photons
is less likely to pass the energy criterion. We show an example cut of 15 MeV
by a dashed vertical line in the left panel of Fig. 3. One can see that most
of the photon pairs in the $m_{a}=100$ MeV case satisfy the cut, while large
(some) fraction in the $m_{a}=1$ (10) MeV case fails in passing the cut. In
principle, softer energy cuts imposed on the final-state photons would be
advantageous for increasing the signal acceptance in a wider range of
parameter space, but the detailed choice also depends on the rejection
capability of backgrounds.

Regarding $\theta_{\gamma\gamma,{\rm det}}$, most of the two final-state
photons are separated enough to be resolved with $\theta_{\rm res}=1^{\circ}$,
as shown in the right panel of Fig. 3. Angular separations are likely to be
larger with increasing $m_{a}$ since ALPs are less boosted.
$\theta_{\gamma\gamma,{\rm det}}$ can be an additional handle to enhance the
signal-over-background depending on $m_{a}$.

Background Considerations. We now discuss potential backgrounds to the ALP
signal, considering non-neutron backgrounds followed by neutron-induced
backgrounds and their rejection strategies.

i) Non-neutron backgrounds: Given the diphoton signature, possible SM
backgrounds would be neutrino neutral-current single $\pi^{0}$ events
occurring at the steel vacuum vessel, e.g., $\nu+{\rm Fe}\to\nu+\pi^{0}+{\rm
Fe}$, since the $\pi^{0}$ decay gives rise to the same signature. However,
RAON-like facilities are essentially stopped-pion sources and the resulting
neutrinos enter the detector with insufficient energy to form a $\pi^{0}$.
Therefore, we expect that this potential backgrounds will not arise at the
proposed detector.

The second class of backgrounds related to neutrinos is the scattering of
neutrino off an electron in the vacuum vessel together with a hard photon
emission, i.e., $\nu+e^{-}\to\nu+e^{-}+\gamma$, in which the final-state
electron is misidentified as a photon. Our simulation suggests that for
$E_{\nu}>15$ MeV, the scattering cross-section of the above process with
$E_{\gamma}>15~{}{\rm MeV}$ is about $10^{-49}~{}{\rm cm}^{2}$ and $\sim 0.11$
neutrinos are produced per proton collision. Considering a $\mathcal{O}(1)$-cm
thickness of the vacuum vessel and the number density of target electrons of
steel, $2.3\times 10^{24}$ cm-3, we find that the number of this type of
neutrino events created in a year is much less than 1.

When it comes to non-neutrino-related backgrounds, particles coming out of the
target end would induce signal-faking events. Although the rate is expected to
be small, the DAMSA detector system in immediate contact with the target may
take a significant flux of such particles, making them non-negligible. Given
the beam energy, however, a 1-m long tungsten target also plays a role of
shielding material. Therefore, the (charged) particles will be almost
completely absorbed before escaping from the target. Our GEANT4 simulation
confirms this expectation and only $\mathcal{O}(2-3)$ photons and electrons
may exit the target for $10^{5}$ proton collisions. If any pairs of such
photons and/or photon-faking electrons are accidentally reconstructed as if
they came from the same vertex point inside the vacuum chamber within the
timing and position resolutions of the ECAL, they will be mis-identified as
signal events. A similar type of potential background events can arise from
neutron-induced events with a much higher rate, which will be discussed in the
next part.

A similar issue would arise with cosmic particles because the detector will be
essentially placed near the surface. While beam triggering would greatly
suppress such backgrounds, a comic tagger or active veto regions surrounding
the detector fiducial volume may be needed.

ii) Neutron-induced backgrounds: The most important feature of the DAMSA
concept is the proximity of the detector to the target, potentially resulting
in a formidable level of neutron flux entering the detector system. The
experimental feasibility relies on the ability to distinguish the neutron-
induced photons mimicking the signal. We here propose a viable strategy to
significantly reduce the neutron-induced backgrounds and discuss its
performance.

First of all, as mentioned earlier, we employ a moderator to suppress the
number of neutrons potentially coming near the decay chamber or the detector.
Our simulation with a 20-cm polyurethane moderator suggests that about 97% of
neutrons be shielded.222The shielding rate, of course, depends on the energy
of the neutrons, and in general, low-energy neutrons can be effectively
shielded by this method. Second, we find that a 15-MeV energy threshold can
reject about 99.99% of photons induced by the remaining neutrons without
deteriorating the signal acceptance. Finally, we devise a few selection
criteria applicable to any photon pairs, making use of the diphotonic nature
of our ALP signal.

  * C1.

Distance-of-closest-approach (DCA): We introduce DCA ($d_{\rm DCA}$) which
defines the closest distance between given two photon tracks. If any two
photons' extrapolated tracks cross over each other at a common point inside
the decay chamber, such a photon pair is potentially mis-tagged as a signal.
DCA allows us to exclude the unrelated photon pairs from signal candidates.
The black line in the left panel of Fig. 4 shows the fractional number of
photon pairs in the $d_{\rm DCA}$ cut.

  * C2.

Arrival time difference: The temporal correlation of any two photons can be
examined, as signal photons are produced simultaneously. We estimate the time
difference ($\Delta t_{\rm arr}$) of two photons at the DCA position and
discard any photon pairs showing a large time difference. The red line in the
left panel of Fig. 4 shows the fractional number of photon pairs in the
$\Delta t_{\rm arr}$ cut.

  * C3.

Fiducial volume: Neutron-induced photons are produced mostly at the steel
frame of the vacuum decay chamber, whereas signal photons emerge inside the
decay volume. We reject the photon pairs satisfying C1 and C2, if their DCA
position is located inside the steel frame.

  * C4.

Trace-back analysis: The vectorial sum of the two photon three momenta, which
would be considered as the ALP momentum, should be traced back to the target
to be misidentified as a signal. We extrapolate back the momentum vector from
the DCA position and examine whether or not this extended line stems from the
target.

  * C5.

Invariant mass window: Finally, the reconstructed invariant mass of the two
photons ($M_{\gamma\gamma}$) should agree with the $m_{a}$ value under
consideration within detector mass resolution $\sigma_{m}$ for which we assume
1 MeV. The right panel of Fig. 4 shows the cut efficiency as a function of the
invariant mass window.

Figure 4: Left: Fractional number of photon pairs vs. $d_{\rm DCA}$ cut (black
solid line) [$\Delta t_{\rm arr}$ cut (red dotted line)]. Right: Fractional
number of photon pairs against invariant mass.

As above, five different independent selection criteria are used to mitigate
the neutron-induced background, and their individual performance is summarized
in Table 1.

Description (per pulse) | Number
---|---
Protons | $4.8\times 10^{7}$
Beam-induced neutrons | $1.29\times 10^{5}$
Neutron-induced photons | $2.74\times 10^{5}$
Photons with $E_{\gamma}>15$ MeV | $25.1$
Photons hitting the detector | $2.94$
$<40$ ns arrival time cut | $1.47$
Photon pairs | $1.08$
Yearly number of pulses | $2.93\times 10^{15}$
Photon pair selection criteria | Efficiency
$d_{\rm DCA}<1$ cm | $4.23\times 10^{-3}$
$\Delta t_{\rm arr}<0.1$ ns | $2.01\times 10^{-1}$
Fiducial volume cut | $6.13\times 10^{-1}$
Back-tracing | $4.16\times 10^{-2}$
$M_{\gamma\gamma}\in 50\pm 1$ MeV | $8.34\times 10^{-3}$
$M_{\gamma\gamma}\in 100\pm 1$ MeV | $1.25\times 10^{-3}$
$M_{\gamma\gamma}\in 200\pm 1$ MeV | $1.91\times 10^{-4}$
Table 1: A summary of the numbers related to the neutron-induced background
mitigation and individual performance of the cuts.

Results. In our analysis, we assume that $\sim 1.4\times 10^{23}$ protons are
annually delivered to the target, and the results are reported for one-year
and 10-year exposures. Our GEANT4 simulation shows that most of the photons
are created at the position about 20 cm away from the facade of the target, so
we assume that all of the photons emerge at this position for simplicity.
Furthermore, since a 20-cm long polyurethane moderator is placed between the
target and the decay chamber to mitigate the neutron flux out of the target,
$\ell$, the distance between the ALP production vertex and the decay chamber
is $\sim 1$ m.

However, depending on the incident angle of ALP, $\ell$ varies. In our
analysis, we set $\ell$ to be the average value, 1.04 m. Likewise, the cord
length that ALPs would travel along is about 10 m, and $L=10.4$ m is used as a
baseline value. The production cross-section of ALPs is computed in the
collinear limit, i.e., $E_{\gamma}\approx E_{a}$, as the momentum transfer to
the target nucleus is much smaller than $E_{\gamma}$. $\sigma_{\rm SM}$ is
generally a function of $E_{\gamma}$ xcom , and we consider photons of
$E_{\gamma}>15$ MeV. If the photon energy is large enough, the change in
$\sigma_{\rm SM}$ is not large. However, we will utilize $\sigma_{\rm
SM}=\sigma_{\rm SM}(E_{\gamma})$ for a more precise estimate, according to the
measurement data archived in xcom . Finally, we envision a detector with
$E_{\rm th}=15~{}{\rm MeV}$ and $\theta_{\rm res}=1^{\circ}$ as mentioned
earlier, and require two final-state photons satisfying the following baseline
criteria:

  * •

their reconstructed tracks intersect at a point inside the decay chamber
within the timing and position resolutions of ECAL, $t_{\rm res}=0.1$ ns and
$\vec{x}_{\rm res}=1$ cm,

  * •

the diphoton invariant mass is consistent with $m_{a}$ under consideration
within $\pm 1$ MeV window, and

  * •

their vectorial 3-momentum sum, which is considered to be the same as the ALP
momentum, traces back toward the target and the vertex position should be
located inside the fiducial volume of the decay chamber.

As we discussed in the previous section, these criteria allow us to suppress
backgrounds including the neutron-induced by more than 7 orders of magnitude.
As conventional electromagnetic calorimetric systems can attain the above-
described level of detector performance, we expect that the resultant
sensitivity reaches can be considered rather conservative.

Figure 5: Expected 90% C.L. sensitivity reaches of the ALP decay signal in the
DAMSA experiment at a RAON-like facility (based on expected statistical error
only) with 1-year (solid) and 10-year (dotted) exposures, under the assumption
of $1.4\times 10^{23}$ POTs per year.

We now report our results of the 90% C.L. sensitivity reaches expected at
DAMSA, in the $(m_{a},g_{a\gamma\gamma})$ plane in Fig. 5. The limits are
evaluated in the basis of expected statistical error only, for both one-year
(solid) and 10-year (dotted) exposures. We also show the existing excluded
regions based on the limits compiled in e.g., Refs. Bauer:2018uxu ;
Lanfranchi:2020crw . The gray-colored covers the regions excluded by the
laboratory-produced ALP searches, and the boundaries are set by various lepton
collider and beam-dump experiments. By contrast, the regions constrained by
astrophysical ALP searches (e.g., HB stars, supernova, etc) are shown in
yellow.

Our findings suggest that DAMSA can explore the regions beyond the current
beam-dump limits, i.e., $m_{a}\sim 100$ MeV in-between $g_{a\gamma\gamma}\sim
10^{-6}~{}{\rm GeV}^{-1}$ and $\sim 10^{-4}~{}{\rm GeV}^{-1}$, which none of
the existing (in)direct experiments have ever probed. The somewhat limited
beam energy relative to those in ongoing and upcoming beam-dump type
experiments would not allow for copious production of $\mathcal{O}(100)$ MeV
ALPs. However, close proximity of the DAMSA detector system to the target can
offset such lack of statistics by capturing a sizable fraction of rather
promptly decaying such ALPs.

Conclusions. We have proposed an experiment that would exploit the high-
intensity proton beams at rare isotope nuclear physics facilities to search
for ALPs, envisioning a facility similar to RAON in South Korea for
illustration. The results show a great potential to cover a wide range of
unexplored parameter space even with one year of data taking, thanks not only
to the high flux proton beams but also to the lower proton energy and the
proximity of the DAMSA detector system to the dump. These advantages could
also open a great potential for expanding physics opportunities beyond the
ALP.

###### Acknowledgements.

We thank Yeonsei Chung and Taeksoo Shin in the RAON Construction Team at the
Institute for Basic Science, Deuksoon Ahn at Korea Basic Science Institute,
and Hyunjun Shin at POSTech for insightful discussions. The work of WYJ and JY
is supported by the U.S. Department of Energy (DOE), HEP Award No. DE-
SC0011686. The work of DK is supported by DOE under Grant No. DE-
FG02-13ER41976, DE-SC0009913, DE-SC0010813. The work of KK is supported by the
DOE No. DE-SC0021447. The work is supported by the National Research
Foundation of Korea (NRF) [NRF-2021R1A4A2001897 (JCP, YK),
NRF-2019R1C1C1005073 (JCP),
NRF-2018R1A6A1A06024977/NRF-2020K1A3A2A01084612/NRF-2018R1A6A1A06024970/NRF-2019R1A2C1008544
(MSR), and NRF-2020R1I1A3072747/NRF-2022R1A4A5030362 (SS)].

##  Appendix A Data Analysis

We focus on the calculation of the expected signal rate at the detector for a
given single photon of energy $E_{\gamma}$. The production cross-section
$\sigma_{P}$ as a function of the ALP angular variables is given by

$\frac{d^{2}\sigma_{P}}{d\theta^{\prime}_{a}d\phi^{\prime}_{a}}=\frac{1}{8\pi}g_{a\gamma\gamma}^{2}\alpha
Z^{2}\left[F(t)\right]^{2}\frac{p_{a}^{4}\sin^{3}\theta^{\prime}_{a}}{t^{2}}\,,$
(1)

where $\alpha$, $Z$, and $p_{a}$ are the electromagnetic fine structure
constant, the atomic number of the target material, and the momentum of the
outgoing ALP, and momentum transfer
$t=m_{a}^{2}-2E_{\gamma}(E_{a}-p_{a}\cos\theta^{\prime}_{a})$. In the
collinear limit where the momentum transfer to the atomic system is
negligible, $E_{a}\approx E_{\gamma}$, and we work on in this limit throughout
this paper. Here $\theta^{\prime}_{a}$ and $\phi^{\prime}_{a}$ are the polar
and the azimuthal angles of ALP measured from the incoming photon direction.
See also Fig. 6 for more detailed configuration. Finally, $F$ describes a form
factor for which we take the Helm parameterization,

$F(t)=\frac{3j_{1}(\sqrt{|t|}R_{1})}{\sqrt{|t|}R_{1}}\exp\left(-\frac{|t|s^{2}}{2}\right)\,,$
(2)

where $j_{1}$ denotes the spherical Bessel function. Here $s=0.9$ fm and
$R_{1}=\sqrt{(1.23A^{1/3}-0.6)^{2}+2.18}$ fm with $A$ being the atomic mass
number of the target material as per Ref. Lewin:1995rx .

Figure 6: An example configuration of an incoming photon (black solid arrow
coming from the origin), an outgoing ALP (red dashed arrow coming from the
origin), and a detector. The photon is assumed to emerge and convert to an ALP
at $(0,0,0)$. The unprimed angles are defined with respect to $\hat{z}$ (i.e.,
the beam axis), while the primed angles are defined with respect to
$\hat{p}_{\gamma}$ (i.e., the momentum direction of the initial-state photon).

Inside the target material, the Primakoff process in (1) competes with the SM
processes such as pair production and photoelectric absorption that the
photons usually go through. Therefore, the probability of ALP production, say
$P_{\rm prod}$, is given by the ratio of the total ALP production cross-
section to the total cross-section of SM interactions $\sigma_{\rm SM}$:

$P_{\rm prod}=\frac{1}{\sigma_{\rm SM}+\sigma_{P}}\int
d\theta^{\prime}_{a}d\phi^{\prime}_{a}\frac{d^{2}\sigma_{P}}{d\theta^{\prime}_{a}d\phi^{\prime}_{a}}\,.$
(3)

To determine the integration limits, we consider the coordinates and angles
defined in Fig. 6. Assuming that a photon emerging at the origin intersects
with the beam axis by an angle of $\theta_{\gamma}$, we see that the unit
vector of the photon direction $\hat{p}_{\gamma}$ is given by
$\hat{p}_{\gamma}=(\sin\theta_{\gamma},0,\cos\theta_{\gamma})$. Then
$\theta^{\prime}_{a}$ and $\phi^{\prime}_{a}$ mentioned above define the unit
vector of the outgoing ALP $\hat{p}^{\prime}_{a}$ about $\hat{p}_{\gamma}$:
$\hat{p}^{\prime}_{a}=(\cos\phi^{\prime}_{a}\sin\theta^{\prime}_{a},\sin\phi^{\prime}_{a}\sin\theta^{\prime}_{a},\cos\theta^{\prime}_{a})$.
A simple algebra enables us to find that the polar angle of ALP $\theta_{a}$
measured from the beam axis is

$\theta_{a}=\cos^{-1}(\cos\theta_{a}^{\prime}\cos\theta_{\gamma}+\cos\phi_{a}^{\prime}\sin\theta_{a}^{\prime}\sin\theta_{\gamma})\,.$
(4)

If $\theta_{a}\leq\theta_{\rm det}$, then the associated ALP can enter the
detector. One can find that the maximum and minimum $\theta^{\prime}_{a}$
values obeying this requirement are

$\displaystyle\theta^{\prime}_{a,+}$ $\displaystyle=$
$\displaystyle\min(\pi,\theta_{\gamma}+\theta_{\rm det})\,,$ (5)
$\displaystyle\theta^{\prime}_{a,-}$ $\displaystyle=$
$\displaystyle\min\left[\max(0,\theta_{\gamma}-\theta_{\rm
det}),2\pi-\theta_{\gamma}-\theta_{\rm det}\right]\,.$ (6)

When it comes to the limits of $\phi^{\prime}_{a}$, we again rely on the
$\theta_{a}\leq\theta_{\rm det}$ condition from which we obtain

$\phi^{\prime}_{a,\pm}=\pm\cos^{-1}\left(\max\left[-1,\frac{\cos\theta_{\rm
det}-\cos\theta^{\prime}_{a}\cos\theta_{\gamma}}{\sin\theta^{\prime}_{a}\sin\theta_{\gamma}}\right]\right)\,.$
(7)

Before moving onto the next step, it is instructive to see the (area-
normalized) ALP spectra in $\theta_{a}$ predicted by the formalism thus far in
combination with the photon flux generated, since they allow for quantitative
estimates on what fraction of the produced ALP flux would traverse the
detector. We show the expected angular spectra in Fig. 7 for three
representative mass values: $m_{a}=1$ MeV (red), 10 MeV (green), and 100 MeV
(blue). For reference purposes, we add a black, dashed, vertical line at
$\theta_{a}=0.5$ rad below which the ALPs would enter the detector. We see
that the spectrum for $m_{a}=1$ MeV does not show a substantial difference
from the corresponding photon angular spectrum, but those for larger $m_{a}$
are more forward-directed. The reason is as follows. Basically, most of the
ALPs are created in the forward region with respect to $\hat{p}_{\gamma}$, as
also suggested in (1). Energetic photons responsible for production of heavy
ALPs, however, are populated relatively more along the beam direction, while
soft photons are more broadly distributed. We find that about 10% of the ALP
flux (i.e., below the black, dashed, vertical line) would move toward the
detector for $m_{a}=100$ MeV.

Figure 7: Angular spectra of ALPs produced in the target for three
representative masses: $m_{a}=1$ MeV, 10 MeV, and 100 MeV. The black dashed
vertical line indicates $\theta_{a}=0.5$ rad by which the coverage of the
proposed detector will be defined.

Once a produced ALP lies in a right direction, it should survive, not decay
before reaching the facade of the detector fiducial region. It is simply
described by the usual decay law. Suppose that the distance between the ALP
production vertex and the detector is $\ell$. The survival probability $P_{\rm
surv}$ is then given by

$P_{\rm surv}=\exp\left(-\frac{\ell}{\ell_{a}^{\rm lab}}\right)\,,$ (8)

where $\ell_{a}^{\rm lab}$ denotes the laboratory-frame mean decay length of
ALP which is a well-known function of the decay width and the boost factor of
ALP $\gamma_{a}$:

$\ell_{a}^{\rm lab}=\frac{\sqrt{\gamma_{a}^{2}-1}}{\Gamma_{a}}\,.$ (9)

On the other hand, an ALP reaching the detector should decay to a photon pair
before escaping from the detector fiducial region to leave an experimental
signature. It is again described by the decay law. Assuming that the cord
length along which the ALP would sweep is denoted by $L$, we write the decay
probability $P_{\rm dec}$ as follows:

$P_{\rm dec}=1-\exp\left(-\frac{L}{\ell_{a}^{\rm lab}}\right)\,.$ (10)

Collecting all the probabilities together, we find that the expected number of
ALP events for any $i$th injection photon (denoted by $N_{a,i}$) can be
written as

$N_{a,i}=P_{\rm prod}\times P_{\rm surv}\times P_{\rm dec}\times 1\,.$ (11)

Note here that the last factor ``1'' in the right-hand side is added
explicitly since we consider a single injection photon. Therefore, if
$N_{\gamma}$ photons are produced in the target for a given exposure, the
total number of ALP events $N_{a,{\rm tot}}$ equals to the summation over
$N_{\gamma}$:

$N_{a,{\rm tot}}=\sum_{i=1}^{N_{\gamma}}N_{a,i}\,.$ (12)

One can calculate an average of $N_{a,i}$ in (11) over a sufficiently large
sample of photons, say $\langle N_{a,i}\rangle$. This can be interpreted as
acceptance of the ALP signal, and $N_{a,{\rm tot}}$ in (12) can be calculated
alternatively by

$N_{a,{\rm tot}}=N_{\gamma}\langle N_{a,i}\rangle=N_{\rm POT}R_{\gamma}\langle
N_{a,i}\rangle\,,$ (13)

where in the second equality we factorize $N_{\gamma}$ into POT for a given
exposure $N_{\rm POT}$ and the average production rate of photons per incident
proton $R_{\gamma}$.

## References

  * (1) J. Wei et al., _Advances of the FRIB project_ , _Int. J. Mod. Phys. E_ 28 (2019) 1930003.
  * (2) R. Catherall et al., _The ISOLDE facility_ , _J. Phys. G_ 44 (2017) 094002.
  * (3) D. Jeon et al., _Design of the RAON accelerator systems_ , _J. Korean Phys. Soc._ 65 (2014) 1010–1019.
  * (4) B. Dutta, D. Kim, S. Liao, J.-C. Park, S. Shin, L. E. Strigari et al., _Searching for dark matter signals in timing spectra at neutrino experiments_ , _JHEP_ 01 (2022) 144, [2006.09386].
  * (5) V. Brdar, B. Dutta, W. Jang, D. Kim, I. M. Shoemaker, Z. Tabrizi et al., _Axionlike Particles at Future Neutrino Experiments: Closing the Cosmological Triangle_ , _Phys. Rev. Lett._ 126 (2021) 201801, [2011.07054].
  * (6) GEANT4 collaboration, S. Agostinelli et al., _GEANT4: A Simulation toolkit_ , _Nucl. Instrum. Meth._ A506 (2003) 250–303.
  * (7) H. Primakoff, _Photoproduction of neutral mesons in nuclear electric fields and the mean life of the neutral meson_ , _Phys. Rev._ 81 (1951) 899.
  * (8) J. B. Dent, B. Dutta, D. Kim, S. Liao, R. Mahapatra, K. Sinha et al., _New Directions for Axion Searches via Scattering at Reactor Neutrino Experiments_ , _Phys. Rev. Lett._ 124 (2020) 211804, [1912.05733].
  * (9) J.-F. Fortin, H.-K. Guo, S. P. Harris, D. Kim, K. Sinha and C. Sun, _Axions: From magnetars and neutron star mergers to beam dumps and BECs_ , _Int. J. Mod. Phys. D_ 30 (2021) 2130002, [2102.12503].
  * (10) K. Agashe, R. Franceschini and D. Kim, _Simple “invariance” of two-body decay kinematics_ , _Phys. Rev. D_ 88 (2013) 057701, [1209.0772].
  * (11) K. K. Boddy, K. R. Dienes, D. Kim, J. Kumar, J.-C. Park and B. Thomas, _Boxes, Boosts, and Energy Duality: Understanding the Galactic-Center Gamma-Ray Excess through Dynamical Dark Matter_ , _Phys. Rev. D_ 95 (2017) 055024, [1609.09104].
  * (12) M. Berger, J. Hubbell, S. Seltzer, J. Chang, J. Coursey, R. Sukumar et al., _Xcom: Photon cross section database (version 1.5)_ , 2010.
  * (13) M. Bauer, M. Heiles, M. Neubert and A. Thamm, _Axion-Like Particles at Future Colliders_ , _Eur. Phys. J. C_ 79 (2019) 74, [1808.10323].
  * (14) G. Lanfranchi, M. Pospelov and P. Schuster, _The Search for Feebly-Interacting Particles_ , 2011.02157.
  * (15) J. Lewin and P. Smith, _Review of mathematics, numerical factors, and corrections for dark matter experiments based on elastic nuclear recoil_ , _Astropart. Phys._ 6 (1996) 87–112.