# Integrable hydrodynamics of Toda chain: case of small systems

Aritra Kundu  Department of Physics and Materials Science, University of
Luxembourg, L-1511 Luxembourg  SISSA, via Bonomea, 265 - 34136 Trieste Italy

###### Abstract

Passing from a microscopic discrete lattice system with many degrees of
freedom to a mesoscopic continuum system described by a few coarse-grained
equations is challenging. The common folklore is to take the thermodynamic
limit so that the physics of the discrete lattice describes the continuum
results. The analytical procedure to do so relies on defining a small length
scale (typically the lattice spacing) to coarse-grain the microscopic
evolution equations. Moving from the microscopic scale to the mesoscopic scale
then requires careful approximations. In this work, we numerically test the
coarsening in a Toda chain, which is an interacting integrable system, i.e. a
system with a macroscopic number of conserved charges. Specifically, we study
the spreading of fluctuations by computing the spatiotemporal thermal
correlations with three different methods: (a) Using microscopic molecular
dynamics simulation with a large number of particles; (b) solving the
generalized hydrodynamics equation; (c) solving the linear Euler scale
equations for each conserved quantities. Surprisingly, the results for the
small systems (c) match the thermodynamic results in (a) and (b) for
macroscopic systems. This reiterates the importance and validity of integrable
hydrodynamics in describing experiments in the laboratory, where we typically
have microscopic systems.

The hydrodynamic theory is a cornerstone for understanding coarse-grained
phenomena in many-body systems. The local evolution of slowly changing fields
of conservation laws manifests the large-scale space-time dynamics of the
system. This old theory has found new interest to predict transport in low-
dimensional systems relevant to the current technologies. In nonintegrable
systems, with few relevant conservation laws, the non-linear fluctuating
hydrodynamic theory helps to understand diffusive and non-diffusive transport
spohn2016 . In integrable systems with many conservation laws, transport is
historically pictured as an underlying quasiparticle undergoing deterministic
scattering. Recent experiments have made novel classical and quantum dynamics
in these systems kinoshita2006 . However, such a formalism can only sometimes
allow for the analytical tractability of physically relevant quantities
measured in the lab or computer experiments. Integrable models such as
harmonic chain and hard point particle gas (or models reducible to them) have
been studied within the standard kinetic theory formalism to connect this
quasiparticle picture and to understand coarse-grained transport properties
bulchandani2019 . The assumed simplicity of a many-body integrable system only
sometimes implies the analytical tractability of physically relevant
quantities. The complexity arises for interacting systems such as the XXZ
chain or the Toda chain spohn2018 . The solvability in these systems are
addressed using Lax pairs in classical Arutyunov2019  and using the Bethe
ansatz yang1969  in quantum setup but characterizing ,dynamical properties
remain challenging using these formalisms. It is typically expected that,
because of the lack of chaos in an integrable system, the transport is
ballistic in nature. This conjecture is verified in the context of energy
diffusion in integrable systems, for example, in Toda chain, Lieb-Linegar gas,
sinh-gordan model etc. theodorakopoulos1984 . However, exceptions can be found
for example, in the integrable XXZ model which is integrable but the spin
transport is non-ballistic das2019 .

Early attempts to describe large-scale dynamics included merging soliton
theory with kinetic theory el2005 ; bulchandani2017 . The general picture of
quasi-particles has recently been incorporated in the development of
integrable hydrodynamics that address the full spatio-temporal structure of
fluctuations and beyond bertini2016 . There is a serious effort to extend
hydrodynamic theory to the quantum domain with growing interest in both
experimental and theoretical perspectives castro-alvaredo2016 ; ruggiero2020 .
Curiously, experiments which contain few atoms are well described by the
hydrodynamic theory initially developed in the thermodynamic limit. This
simple question has no obvious answer and has deep connections to coarse-
graining and taking the hydrodynamic limit varadhan1990 ; spohn2012 . For
example, the 1D Bose gas in the experiment kinoshita2006  composed of $40-250$
atoms may not seem enough for the thermodynamic limit, but the experimental
results can be explained using integrable hydrodynamics.

In this work, we theoretically explore the question: can a small system be
described by hydrodynamic theory? Our analysis is based on the classical
integrable hydrodynamics of the classical Toda chain as developed in spohn2021
; doyon2019 ; nardis2018 ; nardis2019 . The framework can be generalized to
other integrable systems, classical or quantum schemmer2019 . In the quantum
domain, numerical simulations with time-dependent density matrix
renormalization group of a thermal expansion in the XXZ model lead to a
similar conclusion: hydrodynamic predictions with smooth initial conditions
are accurate, even for small system sizes bulchandani2017 .

Here, we construct the theoretical framework to study spatio-temporal
correlation in the Toda chain at equilibrium with three different but related
approaches. The first uses a microscopic molecular dynamics (MD) simulation
similar to kundu2016 , the second is by solving the generalized hydrodynamics
equation for the Toda chain spohn2021 ; doyon2019  and the third by solving
the linear Euler scale equations considering all the conserved quantities. The
three methods have different levels of difficulty / accessibility using finite
computational resources. The MD simulations can be performed for a relatively
large chain (say $N=2048$ particles) for $\sim 3000$ computational hours,
while solving the integrodifferential GHD equation in the thermodynamic limit
is almost instantaneous, but requires a large time to setup the framework.
Solving the Euler equation is limited to a few atoms only and is of medium
complexity. En route to testing the accuracy of the three methods, we find
that a small system with $N=8$ particles has a spatio-temporal correlation at
equilibrium almost indistinguishable from that obtained using the GHD equation
in the thermodynamic limit.

##  I Euler scale hydrodynamics

We discuss the classical many-body interacting integrable system in a general
setting, formulate the conserved quantities, and give the prescription to
construct the Euler scale hydrodynamics using their respective currents.

Integrable ingredients: Consider a system of $N$ particles with position
$q_{i}$ and momentum $p_{i}$ described by the Hamiltonian
$H(\\{p_{i},q_{i}\\})$. The system is integrable if the equation of motion
generated by the Poisson flow $\dot{A}=\\{A,H\\}$ can be written in terms of
the Lax equation, where the dots refer to the time derivative. The common form
of the Lax equation with the Lax pairs denoted by the $N\times N$ matrices
$L^{\pm}\equiv L^{\pm}(\\{p_{i},q_{i}\\})$ and $M\equiv M(\\{p_{i},q_{i}\\})$,
is

$\displaystyle\dot{L}^{\pm}=[M,L^{\pm}]+\lambda L^{\pm}.$ (1)

Using the cyclic property of the trace, the global conserved quantities of the
system are $I^{k}=\mathrm{Tr}(L^{+}L^{-})^{k}$. Note that these conserved
quantities are not unique and are defined to a constant prefactor. These are
chosen such that the equilibrium correlation matrices for energy have an order
of magnitude governed by the temperature $T$ of the system. $M$ is typically a
skew-symmetric square matrix, i.e. a matrix whose transpose equals its
negative. We are interested in determining the local variations of these
globally conserved quantities that describe the current flowing in and out of
a mesoscopic unit cell in the hydrodynamic regime.

Conserved quantities and currents: The local evolution of the conserved
quantity ($I^{k}=\sum_{i}I^{k}_{i}$) can be written in the continuity form by
observing that they follow the same evolution under the Lax flow. Associated
with the evolution, locally conserved quantities $I^{k}_{i}$, is given as,

$\displaystyle\dot{I}^{k}_{i}=(MI^{k})_{ii}-(I^{k}M)_{ii},$ (2)

where we define the local net currents flowing in and out of the site $i$ as
$(MI^{k})_{ii}=j^{k}_{i}$ and $(I^{k}M)_{ii}=j^{k}_{i+1}$ respectively.
Typically for a $H$ with only nearest-neighbour interaction, $M$ is a two-
sided off-diagonal matrix, resulting in the current flowing in and out of
neighboring sites. However, for a long-range interacting system, the current
flows in from all bonds. Another interesting observation is that the
$I^{k}_{i}$ are increasingly non-local in real space in a short-ranged system,
i.e. $I^{k}_{i}$ depend on the product of lattice variables in between $i\pm
k$, pictorially shown in Fig. 1. This will be important to understand our
results later. The eigenvalues of the matrix $L^{+}L^{-}$ is are also
conserved quantities, but, however, they are not local in real space.

Figure 1: Schematic representation of increasing non-locality of conserved
quantities of interacting integrable models with local interaction. The orange
lines show the support of two sites for a local energy, while the grey line
represents a larger support for higher order conserved quantity with support
of four sites.

However, the eigenvalues of the $L^{+}L^{-}$ matrix only capture some of the
conserved quantities for translation-invariant Hamiltonian systems. There is
an additional conserved quantity called the stretch variable, which is related
to the inverse of the system density and plays a major role in transport. For
the nearest-neighbour interaction, if the translation operator is defined as
$\\{T,H\\}=0$, then $I_{0}=\sum_{i}(T-1)q_{i}$ and the total momentum
$I_{1}=\sum_{i}p_{i}$ are the additional conserved quantity for systems with
periodic boundary conditions. It is important to exercise some caution:
Although these look like all the conserved quantities of the system, there can
be others that go unnoticed, i.e. any function of the conserved quantity is
also conserved. See, for example, davies2011  for the role in the non-linear
Schrodinger equation also in the context of Toda lattice, which was
investigated in dhar2021 .

General Gibbs ensemble (GGE): The integrable system at equilibrium is
described formally by GGE,

$\displaystyle\centering
P(\\{p_{i},q_{i}\\})=Z^{-1}\text{exp}[-\sum_{k=0}^{N}\mu_{k}I^{k}(\\{p_{i},q_{i}\\})],\@add@centering$
(3)

where $\mu_{k}$ are the dual thermodynamic variables corresponding to each
conserved quantity. For example, a system at Gibbs equilibrium is described by
dual variable temperature ($\mu_{2}=\beta$, corresponding to Hamiltonian $H$ )
and pressure ($\mu_{0}=P$, corresponding to stretch or the inverse density).
In the rest of the text, any average $\left\langle.\right\rangle$ is
understood to be taken with respect to the GGE characterized by $\beta,P$.

Euler scale (fluctuating) hydrodynamics at equilibrium: The coarse-grained
large-scale space-time behaviour of the system describing the macroscopic
evolution is governed by local fluctuations of the globally conserved
quantities. We are interested in the system prepared in a GGE characterized by
canonically conjugate thermodynamic variables such as temperature, pressure,
etc. In equilibrium, each conserved quantity has a well-defined average that
does not change with time, i.e. $\partial_{t}\langle I^{k}_{i}\rangle=0$. This
contributes to the trivial shift of the conserved quantity of the system and
is typically subtracted. The local fluctuation of the total conserved quantity
$I^{k}$ is denoted $u^{k}_{i}(t)=I^{k}_{i}(t)-\left\langle
I^{k}_{i}\right\rangle$. The linearized hydrodynamic Euler equations for the
evolution of $u^{k}(x,t)$ in the continuum limit are

$\displaystyle\centering\partial_{t}{u}^{k}(x,t)=\partial_{x}j^{k}(x,t),\@add@centering$
(4)

where we have switched from $i\to x$ in going to the continuum and
$j^{k}(x,t)=j^{k}(x,t)-\left\langle j^{k}_{x}\right\rangle$ is the local
current fluctuation associated with the $k$th conserved quantity. Following
the standard prescription of hydrodynamics to the first order in small
fluctuations, we express the current in terms of the local fluctuations of the
conserved fields i.e. $j^{k}(x,t)=\int
A_{kl}(x,t,x^{\prime},t^{\prime})u^{l}(x^{\prime},t^{\prime})dx^{\prime}dt^{\prime}$,
where the kernel $A_{kl}(x,t,x^{\prime},t^{\prime})=\frac{\delta
j^{k}(x,t)}{\delta u^{l}(x^{\prime},t^{\prime})}$ encapsulates all information
of the evolution. We have neglected the higher-order expansions which are not
important in the ballistic scaling limit addressed in the current manuscript
but become important for super-diffusive transport spohn2014 ; spohn2016 . The
kernel is generally hard to compute and much simplification is obtained in the
averaged current at equilibrium, where the time dependence can be dropped due
to time translation invariance at equilibrium. The dependence on positions is
dropped using spatial invariance and local interactions.

We switch to vector notation, with $\mathbf{u}$ the row vector consisting of
the local variation of conserved quantities and $\mathbf{A}$ a square matrix
with dimensions of the number of conserved quantities. The average current
over the GGE at equilibrium is

$\displaystyle\centering{\mathbf{j}(x,t)}=\left\langle\mathbf{A}\right\rangle\mathbf{u}(x,t)+\partial_{x}\mathbf{D}(x,t)\mathbf{u}(x,t)+\mathbf{B}(x,t)\mathbf{\zeta}(x,t).\@add@centering$
(5)

Here $\left\langle\mathbf{A}\right\rangle=\left\langle\mathbf{\frac{\delta
j}{\delta u}}\right\rangle$ can be interpreted as the projection of the
current to the conserved quantity and $\zeta$ a Gaussian white noise with
$\left\langle\zeta(x,t)\zeta(x^{\prime},t^{\prime})\right\rangle=\delta(x-x^{\prime})\delta(t-t^{\prime})$.
The diffusion matrix $\mathbf{D}$ and the noise matrix $\mathbf{B}$ are added
phenomenologically and cannot be easily derived from a microscopical picture.
It may seem surprising that these terms are not zero here, as the chaos
required for the origin of diffusion and noise is absent in an integrable
system. Typically, the origin of these terms is argued as follows: The system
is as a fluid cell where each cell is a ’mesoscopic’ region, i.e. a region of
finite extent, which is large compared to the distance between particles, but
small compared to the macroscopic spatial variation scales $L$. The dynamics
inside the fluid cell can be split into two parts: the first being the net
current flowing through the fluid cell. Each of these currents are projected
on the basis of conserved quantities as discussed before. The second part
comes from the dynamics in the rapid variation in the fluid cell due to
interactions with neighbour cells. This cannot be taken into account through
the projection to the conserved quantities and is modelled as dissipation and
noise.

This gives the Euler scale hydrodynamic equation

$\displaystyle\centering\dot{\mathbf{u}}(x,t)\@add@centering$
$\displaystyle=\partial_{x}\left[\left\langle\mathbf{A}\right\rangle\mathbf{u}(x,t)+\frac{1}{2}\partial_{x}\mathbf{D}(x,t)\mathbf{u}(x,t)+\mathbf{B}(x,t)\mathbf{\zeta}(x,t)\right],$
(6)

The above equation maintains conservation laws as in the microscopic system,
keeping the total fluctuation across the system constant in time. Although
this looks formal, a more rigorous derivation using large deviation theory can
be obtained along the lines of saito2021 . Further, it can be shown that (with
abuse of notation, we drop the average) the matrix
$\mathbf{A}=\left\langle\mathbf{\frac{\delta j}{\delta
u}}\right\rangle\equiv\langle\mathbf{ju}\rangle\langle\mathbf{u^{2}}\rangle^{-1}$,
where we take time independent averages at GGE 111A naive way to see this is
to multiply the numerator and denominator by same fluctuation $u$ and taking a
average. The matrix
$\left\langle\mathbf{ju}\right\rangle_{\alpha\beta}=\left\langle
j_{\alpha}u_{\beta}\right\rangle-\left\langle
j_{\alpha}\right\rangle\left\langle u_{\beta}\right\rangle$, and
$\left\langle\mathbf{u^{2}}\right\rangle_{\alpha\beta}=\left\langle
u_{\alpha}u_{\beta}\right\rangle-\left\langle
u_{\alpha}\right\rangle\left\langle u_{\beta}\right\rangle$ does not have any
time dependence and is computed at equilibrium. This has a nice geometrical
interpretation since the change of curvature of the differential equation is
interpreted as a projection of the current to the conserved quantities in the
system. A more rigorous procedure will be to follow the methods presented in
saito2021 ; spohn2021 . Before going to the next section, where we discuss the
ballistic correlation, a small remark is that for a non-integrable system like
the Fermi-Pasta-Ulam (FPU) chain with three conserved quantities, Eq. (6)
consisting of $3\times 3$ matrix corresponding to the stretch, density and
energy conservation. In this case, truncating the local current to linear
order is not enough and one needs to keep higher-order terms. This gives rise
to evolution with non-ballistic scaling and is described by the non-linear
fluctuating hydrodynamic theory.

Ballistic spread of correlations:  We are interested in spatiotemporal
correlation functions defined as
$\mathbf{C}(x,t)=\langle\mathbf{u(x,t)}\mathbf{u^{T}(0,0)}\rangle$ which
follow the evolution as Eq. (6). It is easy to check that the total
correlations dictated by this evolution is constant, i.e.
$\partial_{t}\sum_{x}\mathbf{C}(x,t)=0$ as the system undergoes conservative
dynamics. In general, the late-time evolution of the correlation functions
follows a self-similar scaling from
$\mathbf{C}(x,t)=\frac{1}{t^{\gamma}}f(\frac{x}{t^{\gamma}})$. It is expected
[although not obvious] that since our system is integrable, the transport will
be ballistic, that is, $\gamma=1$. Using this in Eq. (6), we see that the
contribution of the diffusive term is subleading for ballistic scaling, and we
can safely ignore the last term and are left with

$\displaystyle\centering\dot{\mathbf{C}}(x,t)\@add@centering$
$\displaystyle=\partial_{x}{\mathbf{A}}\mathbf{C}(x,t).$ (7)

Solving the correlation equations: The evolution of the correlations is now a
straightforward linear equation and using the Fourier transformation
$f(x)=\frac{1}{L}\sum_{-L/2}^{L/2}f(k)e^{-i2\pi k/L},~{}x\in-L/2,-L/2+1\dots
L/2$ and $f(k)=\sum_{-L/2}^{L/2}f(x)e^{i2\pi xk/L},~{}k\in[-L/2,-L/2+1\dots
L/2]$, Eq. 7 becomes $\partial_{t}\mathbf{C}(k,t)=-\frac{2\pi
ik}{L}\mathbf{A}\mathbf{C}(k,t)$. The solution is then written as

$\displaystyle\centering\mathbf{C}(k,t)=e^{-2i\pi\frac{k}{L}\mathbf{A}t}\mathbf{C}(k,0),\@add@centering$
(8)

where $\mathbf{C}(k,0)$ is the Fourier transform of the equal time
correlations between the conserved quantities prepared in GGE. The matrix
$\mathbf{A}$ is not symmetric, however the eigenvalues are always real as
expected physically from the spatiotemporal symmetries of the system. This can
be undertood as follows, when the equilibrium system is flipped spatially and
its currents are time-reversed, it behaves identically to the original system.
This requires that the eigenvalues of the matrix $\mathbf{A}$ are real. The
matrix $\mathbf{C}$ is a symmetric matrix. The inverse Fourier transform then
gives the time-dependent solution

$\displaystyle\centering\mathbf{C}(x,t)=\@add@centering$
$\displaystyle\frac{1}{L}\sum_{-L/2}^{L/2}\mathbf{C}(k,t)e^{-i2\pi kx/L}.$ (9)

Taking the continuum limit and scaling the correlation functions, this can be
rewritten in scale-free form as

$\displaystyle\centering f(z)=tC\left(z=\frac{x}{t}\right)=\@add@centering$
$\displaystyle\frac{1}{2\pi}\int_{-\pi}^{\pi}d{k}~{}\mathbf{C}(0,0)e^{-i(\mathbf{A}{k}+{k}{z})}.$
(10)

It is convenient to go to a diagonal basis of $\mathbf{A}$, with
$\mathbf{R}\mathbf{A}\mathbf{R}^{-1}=\text{diag}[\lambda]$, where $\mathbf{R}$
diagonalize the matrix. The eigenvalues of the matrix $\mathbf{A}$ are
physically relevant quantities which can be interpreted as generalized
velocities of ballistic propagating modes. In FPU-like non-integrable systems,
with three conserved quantities, the matrix $\mathbf{A}$ is for example the
$3\times 3$ matrix: the eigenvalues give the sound speed of the system along
with a zero eigenvalue, which corresponds to the stationary heat peak of the
system. In this case, the Euler approximation does not reveal the structure of
the peaks itself, and one must go to higher order expansions in the current to
reveal the structure of the peaks. In an integrable system with ballistic
scaling, the situation is interesting, as there is only one timescale. The
system of equations gives rise to $N$ nonzero velocities corresponding to the
number of conserved quantities in the system. The eigenvalues are ordered in
the sense that $\lambda_{1}\leq\lambda_{2}\dots\leq\lambda_{N}$ and form a
linear space. The spectral gaps between the eigenvalues are smaller toward the
edges and larger in the bulk. Taking these into account renormalizes the
amplitude of the scale-free function at the $k^{th}$ site by
$\lambda_{k+1}-\lambda_{k}$. The scale-free equation Eq. (10) can be written
in terms of these velocities as,

$\displaystyle\centering f_{ij}(z)=\@add@centering$
$\displaystyle\sum_{k}R_{i,k}(R^{-1}C)_{k,j}\frac{\delta({z}-\lambda_{k})}{\lambda_{k+1}-\lambda_{k}}.$
(11)

The normalization can also be understood from a more physical point of view:
Since the integrable dynamics of the system is conservative, normalization
eliminates the gauge degree of freedom from the eigenvectors, and we have the
sum rule discussed above, $\sum_{z}\mathbf{f}(z)dz=\mathbf{C}(0,0)$, where the
$dz=\lambda_{k+1}-\lambda_{k}$ are the discrete inverse density of the ordered
eigenvalues and the right-hand side is the value of correlations computed at
equilibrium. The density fluctuations are studied under the “unfolded”
spectral gap of the eigenvalues of the integrable system akin to that used in
random matrix theory. A discussion of this formula for the Toda lattice is
shown in Sec. III. As we shall see, the above formula for a small system size,
i.e. when $N$ is of the order of $O(1)$ reproduces the thermodynamic limit
spatio-temporal correlations. Note that we have not taken any hydrodynamic or
large system size limit in the above formula. This makes us wonder if the
dynamics of small systems at equilibrium are sufficient to be described by
hydrodynamics.

Figure 2: Comparison of scaled equilibrium correlations in a small finite
sized Toda lattice with thermodynamic limit results. The circle and the stars
are data from numerically solving the Euler scale hydrodynamic equation at
Gibbs equilibrium as in Eq. 11 for small system sizes with $N=8$ and $N=16$.
The noisy grey line is microscopic molecular dynamics simulation of the
macroscopic number of $N=2048$ particles with Toda interaction with an initial
condition average of $10^{7}$ initial conditions from the Gibbs distribution
with $T=1,P=1$. The black dashed line is obtained by numerically solving the
thermodynamic GHD equation. We see that the results of the finite short system
already have an excellent overlap with the exact results in the thermodynamic
limit.

##  II Generalized hydrodynamics: GHD

Generalized hydrodynamics assumes that the dynamics of an integrable system
can be described by a quasi-particle which undergoes ballistic diffusion and a
drift with effective velocity. At thermal equilibrium, this is obtained by
averaging the extensive number of conserved charges and their dynamics in a
local GGE state. The natural hydrodynamic fields are the density of states
(DOS) of the Lax matrix and the stretch, $r$, which we anticipated earlier.

For completeness, we first briefly review the basic concepts of the GHD in
classical systems following the notation introduced in mendl2022 ; spohn2021 .
The DOS is represented as $r\rho_{p}(r)$, where the distribution $\rho_{p}$
encodes all information of the conserved quantities. Then the generalized
hydrodynamic equations are

$\displaystyle\centering\partial_{t}r+\partial_{x}q_{1}=0,\quad\partial_{t}(r\rho_{p})+\partial_{x}\left((v^{eff}-q_{1})\rho_{p}\right)=0,\@add@centering$
(12)

where $q_{1}=r\int dw~{}w\rho_{p}(w)$ is the average momentum of the system.
The effective velocity is given by the solution of the linear integral
equation of the form

$\displaystyle\centering
v^{\text{eff}}(r)=r+(T\rho_{p}v^{\text{eff}})(r)-(T\rho_{p})v^{\text{eff}}(r),\@add@centering$
(13)

with the operator $T$ being the integral of kernel ($K$) which takes into
account the two-particle scattering shift. It is defined as an integral
transform on a function $\psi(w)$ as,

$\displaystyle\centering T\psi(w)=\int
K(w-w^{\prime})\psi(w^{\prime})~{}dw^{\prime},\@add@centering$ (14)

where the kernel $K$ depends on the microscopic model and can be computed by
solving the large time asymptotic in two-particle scattering.

It is convenient to use a non-linear transformation to go to a coordinate
frame, where the hydrodynamic equation becomes simplier. This is done by
transformation to normal mode DOS $\rho_{n}$,

$\displaystyle\centering
1-T\rho_{n}(r)={{(1+T\rho_{p}(r)})^{-1}}.\@add@centering$ (15)

Key to this description is to introduce a dressing transformation that
transforms the real valued function $\psi$ as

$\displaystyle\centering\psi^{dr}=(1-T\rho_{n})^{-1}\psi.\@add@centering$ (16)

In this notation, the transformation between the normal mode and the effective
velocity is given as

$\displaystyle\centering\rho_{p}=\rho_{n}[1]^{dr},~{}~{}v^{eff}=\frac{[w]^{dr}}{[1]^{dr}},\@add@centering$
(17)

and the normal mode hydrodynamic equation is

$\displaystyle\centering
r\partial_{t}\rho_{n}+(v^{eff}-q_{1})\partial_{x}\rho_{n}=0,\@add@centering$
(18)

It is convenient to denote the average over DOS as $\left\langle
f\right\rangle_{p}=\int f(w)\rho_{p}(w)dw$. In this notation, the conserved
quantities are given as the integral identities $\left\langle
r\right\rangle_{p}=1$ and $\left\langle rw^{k}\right\rangle_{p}=I_{k}$. We are
interested in the correlation function of the conserved quantities,

$\displaystyle\centering
C_{\alpha\beta}=[\mathbf{C}]_{\alpha\beta}=\left\langle
r(w^{\alpha}-I^{\alpha}[1]^{dr})(w^{\beta}-I^{\beta}[1]^{dr})\right\rangle_{p}.\@add@centering$
(19)

Thermal states: For thermal states, the $\rho_{n}$ is the solution of the
thermodynamic Bethe-Ansatz (TBA) like equation characterized by the chemical
potential $\mu$ a function of the external pressure and temperature determined
at equilibrium by condition $\int n(w)dw=P$,

$\displaystyle\centering
V(w)-2(T\rho_{n})(w)+\log(\rho_{n}(w))-\mu=0,\@add@centering$ (20)

where $V(w)=\frac{1}{2}\beta w^{2}$ for the initial state in Gibbs
equilibrium. Note that for a classical system there is no TBA, but this
equation is obtained in Toda chain in a similar spirit using a nontrivial
mapping to random matrix theory with the kernel given by $K(x)=\log(x)$, which
is determined from the two-particle scattering. The solution to this highly
non-linear root equation is often supplemented with solving the corresponding
Fokker-Planck equation whose steady state solution satisfies the TBA equation
above,

$\displaystyle\centering\partial_{t}\rho(w,t)=\partial_{w}\left(V^{\prime}(w)\rho(w,t)-2P\int
K^{\prime}(w-v)\rho(v,t)\rho(w,t)dv\right)+\partial_{w}^{2}\rho(w,t),\@add@centering$
(21)

where $K^{\prime}$ denotes the derivative of the integral kernel $K$ which is
specific to the model. As a consistency check, the solution to Eq. (21) is
inserted into the TBA Eq. (20), gives the chemical potential $\mu$. At
equilibrium, $\mu$ turns out to be exactly the free energy when we set $w\to
0$ which acts as a cross-check for the results with the explicit formula
available kundu2016 ; takayama1986 . At this point, it is convenient to
introduce the quasi-energy parameterization of the normal mode obtained from
the solution of the Fokker-Planck solution of TBA by
$\epsilon(w)=-\log[\rho_{n}(w)]$. Numerically, this step is carried out by
first obtaining $\epsilon$ from the Fokker-Planck equation and then plugging
it into the TBA Eq. (20) and finding the roots of the equation. The average
momentum in terms of dressed variable is

$\displaystyle\centering\left\langle p\right\rangle_{\beta,P}=\left\langle
r\right\rangle_{\beta,P}\int\rho_{n}[w]^{dr}dw.\@add@centering$ (22)

The correlations in Eq. (19) in case of momentum correlation is conjectured as
[Eq. 3.17 in spohn2021 ],

$\displaystyle\centering C_{pp}\approx\left\langle
r\right\rangle_{\beta,P}\int dw\rho_{p}(w)\delta(x-t{v^{eff}(w)}{\left\langle
r\right\rangle_{\beta,P}^{-1}})([v]^{dr})^{2}(w).\@add@centering$ (23)

Other higher-order correlations can be computed appropriately by generalising
this formula.

##  III Toda chain

We now focus on applying the methods described on the Toda chain toda2012
which is a model of classical particles with exponentially decaying nearest-
neighbour interactions. This is an example of a classical interacting
integrable particle system spohn2018  which is solvable by Lax pairs as
described before. A large Toda lattice is used as a discrete model for
solitonic waves while a small chain of three particles can be mapped to the
famous Hennon Heiles problem. The history of understanding dynamics in Toda
lattice is decades old: The equilibrium dynamical properties in large Toda
chains were extensively studied using various methods. Numerically,
correlations have been explored in the literature in kundu2016  using
molecular dynamics simulations. In an integrable system of $N$ particles, an
equivalent number of conserved quantities is enough to describe the dynamics
completely. 222However, some recent results challenge this assumptiondhar2021
.

The Toda chain is defined as
$H({p_{i},r_{i}})=\sum_{i}e_{i}=\sum_{i}\frac{p_{i}^{2}}{2}+e^{-r_{i}}$ with
$r_{i}=q_{i+1}-q_{i}$. The equations of motion (e.o.m.) are:

$\displaystyle\dot{r_{i}}$ $\displaystyle=$ $\displaystyle p_{i+1}-p_{i},$
$\displaystyle\dot{p_{i}}$ $\displaystyle=$ $\displaystyle
e^{-r_{i-1}}-e^{-r_{i}}.$ (24)

These e.o.m. along with periodic boundary conditions $r_{N+1}=r_{1}$ and
$r_{0}=r_{N}$ can be cast to a Lax matrix as in Eq. (1), with $\lambda=0$ and
$L=L^{+}=L^{-}$. The symmetric matrix $L$ is a tri-diagonal matrix with
diagonal elements $b_{i}=\frac{p_{i}}{2}$ and off-diagonal elements
$a_{i}=\frac{1}{2}e^{-r_{i}/2}$. The conserved quantities are trace of the
powers of Lax matrix, denoted as $Q_{n}$ where $Q_{1}=\sum_{i=1}^{N}p_{i}$ and
$Q_{2}=\sum_{i=1}^{N}e_{i}$. In shastry2010 , explicit higher-order conserved
quantities are remunerated. Importantly, in real space, the conserved
quantities have increasing non-local support in space(see Fig. 1). However,
this is not the only way to view the conserved quantities of the system; for
example, the eigenvalues of the Lax matrix are also conserved under evolution.
It has been suggested that the spacing between the eigenvalues of the Lax
matrix provides a natural transition to normal mode frequencies in
goldfriend2018 ; goldfriend2019 . These conserved quantities are an adiabatic
invariant for small non-integrable perturbations grava2020 . Prior attempts
have been made to understand the dynamics of the Toda chain. The inverse
scattering method is useful in an isolated system and provides a way to
address not only the solitonic properties but also the large-scale
hydrodynamic properties, as discussed earlier. Thermalization properties and
the conundrum of the local observable are related to Gibbs or General Gibbs is
addressed in baldovin2021 . The majority of the works address single-point
functions and finite temperature dynamical spatio-temporal properties although
difficult to access has been attempted through non-interacting soliton gas
analogy diederich1981  and through quantum-classical correspondence with
quantum Toda chain cuccoli1993 . Earlier attempts at dynamic properties of the
Toda chain are reviewed in cuccoli1994 . After a gap in interest in Toda
dynamics, recent development in the integrable hydrodynamic theory of Toda
chain in spohn2016 ; doyon2019 , where the dynamics of single-point functions
starting from domain wall initial conditions was addressed mendl2022 .
Accurate numerical results for the dynamics of two-point functions were
discussed in kundu2016  and analytical results were only available in the
limiting cases when the model was solvable. As a by product, this work fills
the gap in theoretical description of the numerical results in spatio-temporal
correlation functions described in kundu2016  from the theory developed in
spohn2021 ; doyon2019 .

Results: Fig. (2) summaries the main numerical results of computing the
momentum spatio-temporal correlations at Gibbs thermal equilibrium. The system
described by Eq. (24) is solved together with the initial conditions at
equilibrium described by temperature ($\beta^{-1}=T$) and pressure ($P$) in a
Gibbs ensemble with $Z=\prod_{i}e^{-\beta(e_{i}+Pr_{i})}$. The noisy grey line
is a microscopic numerical simulation of the macroscopic number of $N=2048$
particles with Toda interaction with an average of $10^{7}$ initial conditions
sampled from the Gibbs distribution with $T=1,P=1$. The black dashed line is
obtained by numerically solving the GHD equation in Eq. (23) with the input
from the TBA equation for computing the effective velocity of the quasi-
particle and analytically computed thermal expectation values. Exact numerical
details of the relevant computation is given in mendl2022 , and a code is
available at mendlcode . The circle and the stars are data from numerically
solving the Euler scale hydrodynamic equation at Gibbs equilibrium as in Eq.
(LABEL:eq:_scale-free) for small system sizes with $N=8$ and $N=16$. As
anticipated earlier in the introduction, surprisingly, even for the small
system sizes, the spatiotemporal correlations have a good agreement with
correlations in a large system sizes. This is indeed a surprising observation,
and it make us wonder about the question of taking a thermodynamic limit,
which serves as the motivation of the problem. Note that the positions for
each of the dots corresponding to solving the Euler scale equation can be
interpreted as “generalized” velocities of the system in the ballistic scaling
limit. When the dimensions of the Lax matrix are of order $N$, there are $N+1$
velocities, which overlaps with the thermodynamic solution of the GHD. It is
interesting to note that these generalized velocities are neither symmetrical
nor are they equally spaced. The peak of the curves are not captured well with
this method, however, the possibility of the effects of taking of even-odd
number of particles cannot be ruled out. The relative density of state of the
distribution of the Lax eigenvalues are shown in Fig. (3), which shows that
the eigenvalues are closer to the edges and more spread out in the bulk.

Numerical method for direct Euler scale dynamics: The numerical method used to
simulate the hydrodynamics of the Euler scale is by Monte Carlo sampling of
the Lax matrix. The trick here is to do the Monte Carlo sampling from the Lax
matrix itself corresponding to the Gibbs ensemble of the system at a specific
temperature and pressure. For the Toda system, it can be shown that in Gibbs
distribution: i.e. when the system is prepared in
$P(\\{p_{i},r_{i}\\})=Z^{-1}e^{-\beta(H(\\{p_{i},r_{i}\\})+Pr_{i})}$ ($Z$ is
the normalization), the elements $a_{i}$ follow the chi-square distribution
(Fig. 4 for numerical sampling), $P_{a}(x)=\frac{2(4\beta)^{\beta
P}}{\Gamma(\beta P)}x^{2\beta P-1}e^{-4\beta x^{2}},$ and the momentum follow
the Gaussian distribution. The eigenspectrum of the Lax matrix is used to
compute the Euler-scale hydrodynamics along with the equilibrium correlations
of conserved quantities. In our simulations, we used a $10^{9}$ sample of the
Lax matrix and checked the convergence of the eigenvalues of the matrix
$\mathbf{A}$. The equilibrium correlations of the conserved quantities are
more difficult objects to deal with analytically. As we mentioned before, the
higher-order conserved quantities have increasing non-local support in real
space, and the convergence of equilibrium correlations becomes difficult.
However, here we are interested in the system at Gibbs equilibrium, where the
relevant conserved quantities that enter the distribution are energy and
stretch, which are both extremely local in space. This makes the connected
correlation of the higher-order conserved quantities negligible compared to
the correlations of energy and the stretch variables. Interestingly,
truncating the Lax matrix to include only the first three conserved quantities
which appear in the Gibbs ensemble: the energy, momentum, and the stretch
reproduce to $5\%$ error on the motion of the peaks of the dynamical
fluctuations kundu2016  with the non-zero eigenvalues $\lambda=\pm 0.8833$ and
we conclude that the correlations of the higher non-local conserved quantities
are important to capture the dynamical correlations of the more local
conserved quantities.

Comparison and discussion and limitations of results

We reiterate that these results are strictly for correlations in the ballistic
scaling limit. There are some intuitive results on the nature of diffusive
corrections on the ballistic scale, which we will not discuss here doyon2019 .
Although the method seems general, it is by no means more efficient
numerically than solving the GHD equations directly. The GHD framework
provides a general analytic framework for working with set generic initial
conditions. On the other hand, the Euler scale hydrodynamics is a more brute-
force approach to the problem using Monte Carlo sampling of the initial data
and is limited to small system sizes.

##  IV Conclusions

We have discussed the adequacy of using integrable hydrodynamics to describe
small classical (and quantum) systems. We show that although the hydrodynamic
theory is developed for a large number of particles in thermodynamic limit,
practically, using an independent direct method discussed here, a relatively
small system size in a classical interacting integrable system is already
quite accurate in describing the results in thermodynamic limit. We have
tested this for a short-ranged interacting system; however, questions related
to the finite-size hydrodynamics in a long-ranged system (for example in
Colegaro Sutherland model) remains to be explored.

Acknowledgements: I thank Herbert Spohn, Abhishek Dhar, Manas Kulkarni for
discussions. I thank Christian Mendl for his correspondence and GHD code. The
conference SPLDS2022 in Ponte-e-Mousson on the occasion of the birthday of
Malte Henkel in 2022 served as a transition point for writing this draft.
After completion of the article, I was informed that a similar study involving
MD of is being addressed in grava22 , due to a conflict of interest, we
decided to publish the articles separately. I thank Aurelia Chenu for reading
and improving the quality of the paper. This work was partially funded by the
Luxembourg National Research Fund (FNR, Attract grant 15382998).

## References

  * (1) Herbert Spohn.  Fluctuating hydrodynamics approach to equilibrium time correlations for anharmonic chains.  In Thermal Transport in Low Dimensions, pages 107–158. Springer International Publishing, 2016.
  * (2) Toshiya Kinoshita, Trevor Wenger, and David S. Weiss.  A quantum newton's cradle.  Nature, 440(7086):900–903, apr 2006.
  * (3) Vir B Bulchandani, Xiangyu Cao, and Joel E Moore.  Kinetic theory of quantum and classical toda lattices.  Journal of Physics A: Mathematical and Theoretical, 52(33):33LT01, jul 2019.
  * (4) Herbert Spohn.  Interacting and noninteracting integrable systems.  Journal of Mathematical Physics, 59(9):091402, sep 2018.
  * (5) 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(7):1115–1122, jul 1969.
  * (6) Nikos Theodorakopoulos.  Finite-temperature excitations of the classical toda chain.  Physical Review Letters, 53(9):871–874, 1984.
  * (7) Avijit Das, Kedar Damle, Abhishek Dhar, David A. Huse, Manas Kulkarni, Christian B. Mendl, and Herbert Spohn.  Nonlinear fluctuating hydrodynamics for the classical XXZ spin chain.  Journal of Statistical Physics, 180(1-6):238–262, oct 2019.
  * (8) Bruno Bertini.  Transport in out-of-equilibrium chains: Exact profiles of charges and currents.  Physical Review Letters, 117(20), 2016.
  * (9) G. A. El.  Kinetic equation for a dense soliton gas.  Physical Review Letters, 95(20), 2005.
  * (10) Vir B. Bulchandani.  Solvable hydrodynamics of quantum integrable systems.  Physical Review Letters, 119(22), 2017.
  * (11) Olalla A. Castro-Alvaredo.  Emergent hydrodynamics in integrable quantum systems out of equilibrium.  Physical Review X, 6(4), 2016.
  * (12) Paola Ruggiero.  Quantum generalized hydrodynamics.  Physical Review Letters, 124(14), 2020.
  * (13) S.R.S. Varadhan.  On the derivation of conservation laws for stochastic dynamics.  In Paul H. Rabinowitz and Eduard Zehnder, editors, Analysis, et Cetera, pages 677–694. Academic Press, 1990.
  * (14) Herbert Spohn.  Large scale dynamics of interacting particles.  Springer Science & Business Media, 2012.
  * (15) Herbert Spohn.  Hydrodynamic equations for the toda lattice.  01 2021.
  * (16) Benjamin Doyon.  Generalized hydrodynamics of the classical toda system.  Journal of Mathematical Physics, 60(7):073302, jul 2019.
  * (17) Jacopo De Nardis.  Hydrodynamic diffusion in integrable systems.  Physical Review Letters, 121(16), 2018.
  * (18) Jacopo De Nardis, Denis Bernard, and Benjamin Doyon.  Diffusion in generalized hydrodynamics and quasiparticle scattering.  SciPost Phys., 6:49, 2019.
  * (19) M. Schemmer, I. Bouchoule, B. Doyon, and J. Dubail.  Generalized hydrodynamics on an atom chip.  Physical Review Letters, 122(9), mar 2019.
  * (20) Aritra Kundu.  Equilibrium dynamical correlations in the toda chain and other integrable models.  Physical Review E, 94(6), 2016.
  * (21) B. Davies and V. E. Korepin.  Higher conservation laws for the quantum non-linear schroedinger equation, 2011.
  * (22) Abhishek Dhar, Aritra Kundu, and Keiji Saito.  Revisiting the mazur bound and the suzuki equality.  Chaos, Solitons and Fractals, 144:110618, 2021.
  * (23) Herbert Spohn.  Nonlinear fluctuating hydrodynamics for anharmonic chains.  Journal of Statistical Physics, 154(5):1191–1227, feb 2014.
  * (24) B. Sriram Shastry.  Dynamics of energy transport in a toda ring.  Physical Review B, 82(10), 2010.
  * (25) Keiji Saito, Masaru Hongo, Abhishek Dhar, and Shin ichi Sasa.  Microscopic theory of fluctuating hydrodynamics in nonlinear lattices.  Physical Review Letters, 127(1), jul 2021.
  * (26) Christian B. Mendl and Herbert Spohn.  High-low pressure domain wall for the classical Toda lattice.  SciPost Phys. Core, 5:2, 2022.
  * (27) Gleb Arutyunov.  Elements of Classical and Quantum Integrable Systems.  UNITEXT for Physics. Springer International Publishing.
  * (28) https://github.com/cmendl/toda-domainwall.
  * (29) Mendl. et.al.  To be published.
  * (30) H. Takayama and M. Ishikawa.  Classical thermodynamics of the toda lattice: –as a classical limit of the two-component bethe ansatz scheme–.  Progress of Theoretical Physics, 76(4):820–836, oct 1986.
  * (31) Morikazu Toda.  Theory of nonlinear lattices, volume 20.  Springer Science & Business Media, 2012.
  * (32) T. Goldfriend and J. Kurchan.  Fluctuation theorem for quasi-integrable systems.  EPL (Europhysics Letters), 124(1):10002, oct 2018.
  * (33) Tomer Goldfriend.  Equilibration of quasi-integrable systems.  Physical Review E, 99(2), 2019.
  * (34) T. Grava, A. Maspero, G. Mazzuca, and A. Ponno.  Adiabatic invariants for the FPUT and toda chain in the thermodynamic limit.  Communications in Mathematical Physics, 380(2):811–851, sep 2020\.
  * (35) Marco Baldovin, Angelo Vulpiani, and Giacomo Gradenigo.  Statistical mechanics of an integrable system.  Journal of Statistical Physics, 183(3):1–16, 2021.
  * (36) S. Diederich.  A conventional approach to dynamic correlations in the toda lattice.  Physics Letters A, 85(4):233–235, 1981.
  * (37) A. Cuccoli, M. Spicci, V. Tognetti, and R. Vaia.  Dynamic correlations of the classical and quantum toda lattices.  Physical Review B, 47(13):7859–7868, apr 1993.
  * (38) A. Cuccoli, R. Livi, M. Spicci, V. Tognetti, and R. Vaia.  Thermodynamics of the toda chain.  International Journal of Modern Physics B, 08(18):2391–2446, aug 1994.

##  V Appendix

Density of states: The relative density of states of the finite-dimensional
Lax matrix for Toda used to calculate the correlation function is shown in
Fig. 3.

Figure 3: The denisty of states $1/(\lambda_{i+1}-\lambda_{i})$ of the
distribution of sorted eigenvalues of the Lax matrix for the Toda lattice in
Gibbs equilibrium at temperature $T=1$ and pressure $P=1$. The sampling of the
Lax matrix at equilibrium is achieved through transformation of the
probability distribution. The diagonal elements of the Lax matrix are Gaussian
to sample while the off-diagonal elements is shown in Fig. 4

.

Monte-Carlo simulation of the Lax matrix:

Figure 4: Monte-carlo sampling of the distribution of the off-diagonal
elements of Lax matrix as given by equation Eq. III. We use this to compute
the Lax matrix’s eigenvalues and the correlation functions.

Numerical values in matrix $\mathbf{A}$:

The matrix $\mathbf{A}$ for $N=16$ truncated to $4\times 4$ matrix:

$\mathbf{A}=\begin{bmatrix}-0.0002&-0.9993&0.0017&-0.007&\dots\\\ \
-0.3357&0.0021&0.4819&-0.0114\\\ 0.0002&-0.2333&-0&1.2495\\\
-0.234&0.0038&-0.0556&-0.0313\\\ \vdots&&&&\ddots\end{bmatrix}$


# Plasma harmonic generation for highly efficient Breit-Wheeler pair creation

Suo Tang<EMAIL_ADDRESS>College of Physics and Optoelectronic
Engineering, Ocean University of China, Qingdao, Shandong, 266100, China

###### Abstract

Observation of efficient Breit-Wheeler electron-positron pair creation is one
of the main goals for modern-day strong laser-particle experiments. We propose
exploring this process by colliding plasma-generated harmonics, from currently
available laser pulses, with a beam of GeV photons. We show that the creation
yield is enhanced for more than one order of magnitude, compared with the
yield obtained without high-harmonic generation. The robustness of the yield
enhancement is demonstrated by considering multiple interaction parameters and
the potential photon source for the upcoming experiments. Moreover, we show
that due to the evident field asymmetry in the plasma-generated harmonics, the
spin-polarization of the created pairs could be as high as $60\%$.

Breit-Wheeler (BW) pair creation—the decay of light quanta into an electron-
positron pair, is one of the most intriguing strong-field quantum
electrodynamic (QED) phenomena Breit and Wheeler (1934); Ritus (1985); Di
Piazza et al. (2012); Fedotov et al. (2023); Gonoskov et al. (2022) and one of
the main goals for modern-day strong laser-particle experiments Joshi et al.
(2018); Jacobs (2022); Borysova (2021); Abramowicz et al. (2021); Naranjo et
al. (2021); Chen et al. (2022); Macleod (2022); Turner et al. (2022). The
current high-power laser facilities in operation or in development have
entered the multipetawatt regime and pave the way to measure the process in
the nonlinear and nonperturbative regime Sung et al. (2017); Kiriyama et al.
(2018); Shen et al. (2018); Danson et al. (2019); Yoon et al. (2021). This
process has been investigated theoretically in different type of laser fields
Nikishov and Ritus (1964); Jansen and Müller (2013); Meuren et al. (2016);
Jansen et al. (2016); Titov et al. (2018); Ilderton (2020, 2019a); Wistisen
(2020); Seipt and King (2020); Mercuri-Baron et al. (2021); Golub et al.
(2021); Tang (2021); Tang and King (2021); Tang (2022); Gao and Tang (2022);
Borysov et al. (2022) to improve the creation yield, tailor the spectra of the
created particles, and recently suggested to generate polarized positron beams
in ultraintense finely-tuned asymmetric fields Seipt et al. (2019); Chen et
al. (2019); Wan et al. (2020). A robust scheme to generate brilliant, highly
polarized positron beams is crucial for testing the physics beyond the
standard model Moortgat-Pick et al. (2008); Vauth and List (2016).

The typical scenario for observing the BW pair creation is to collide a beam
of high-energy probe photons with a strong laser field. In the center-of-mass
frame of the collision, there is an energy threshold $2m$ must be satisfied,
where $m$ is the electron’s rest mass and the natural units $\hbar=c=1$ is
used throughout. To overcome this threshold, one can improve either the energy
of the probe photon, or the intensity of the laser pulse to provide more laser
photons for the creation. The yield of the process thus depends mainly on two
factors: i) the parameter $\eta=k\cdot\ell/m^{2}$ characterizing the center-
of-mass energy of the collision and the number of laser photons $n\sim 2/\eta$
needed for the creation; ii) the classical intensity
$\xi_{l}=e\mathcal{E}_{l}/m\omega_{l}$ of the laser pulse characterizing the
nonlinearity of the creation, where $\ell^{\mu}$ is the photon momentum,
$k^{\mu}$, $\mathcal{E}_{l}$ and $\omega_{l}$ are respectively the laser’s
wavevector, intensity and carrier frequency of the laser, and $e$ is the
positron charge.

In the upcoming laser-particle experiments, _e.g._ LUXE Abramowicz et al.
(2021) and FACET II-E320 Chen et al. (2022), high-energy photons with
$E_{\gamma}\sim 10~{}\textrm{GeV}$ and laser pulses with the intensity
$\xi_{l}\sim O(1)$ and carrier frequency $\omega_{l}\sim 1~{}\textrm{eV}$
would be utilized. The BW yield of the potential measurement in these
campaigns would be limited as $\eta\ll 1$ and $\eta\xi\sim 1$. To improve the
particle yield, we propose to collide the high-energy photons with the plasma-
generated harmonics in the reflection from the currently available laser pulse
interacting with an overdense plasma Lichters et al. (1996); Baeva et al.
(2006); Teubner and Gibbon (2009); Tang and Kumar (2019). The energy parameter
is considerably improved by the rich high-order harmonics, and due to the
locking phase of the low-order harmonics Gordienko et al. (2004); Gonoskov et
al. (2011); Bulanov et al. (2013), the pulse amplitude is effectively
magnified. The creation yield could thus be significantly enhanced compared
with that in the initial pulse without the harmonic generation. Moreover, the
created particles could be highly polarized by the asymmetric harmonic field
when proper oblique incidence of the laser pulse is considered.

Figure 1: Sketch for efficient Breit-Wheeler pair creation using plasma
generated harmonics: a strong laser pulse irradiates an overdense plasma with
the incident angle $\theta$ and gets reflected to collide with a beam of high-
energy probe photons at the collision point (CP). The contour plots are the
electric field of the incident pulse with $\xi_{l}=3$ and reflected pulse
resulting from a 2D PIC simulation for the incidence $\theta=45^{\circ}$. In
the insets, the comparisons between the incident and reflected pulse at the
pulse center $\Delta_{c}=0$ are given: (a) the electric fields and (b) the
corresponding frequency spectra. The detail of the simulation is given in the
main text.

The sketch of our proposal is shown in Fig. 1: a linearly polarized laser
pulse, irradiating an overdense plasma with the incident angle $\theta$, gets
reflected at the plasma surface and then collides nearly head-on with a beam
of high-energy photons at the collision point. The periodically oscillating
current driven by the incident pulse at the plasma surface, modulates the
incident field and emits rich high-order harmonics into the reflected field
Teubner and Gibbon (2009). The electric fields of the incident and reflected
pulse are compared in Fig. 1 as an example for the incident angle
$\theta=45^{\circ}$, resulting from a 2D particle-in-cell (PIC) simulation
using the code EPOCH Arber et al. (2015): The incident laser pulse,
$\xi(\phi,\Delta_{c})=\xi_{l}\sin(\phi)\sin^{2}\left[\phi/(2N)\right]\exp(-\Delta^{2}_{c}/\sigma^{2})$,
is in the duration $0<\phi<2N\pi$, where $\xi_{l}=3$, $N=16$, $\phi=k\cdot x$
is the laser phase of the carrier frequency, $\sigma=6\lambda_{l}$ corresponds
to the full-width-at-half-maximum about $10\lambda_{l}$ and the laser
wavelength is $\lambda_{l}=0.8\mu\textrm{m}$; The plasma target has an
overdense bulk $n_{e}=15n_{c}$ and a pre-gradient $n_{c}\exp(d/L)/2$, where
$L=\lambda_{l}/11$, $d$ denotes the distance to the plasma surface,
$0<d<L\ln(2n_{e}/n_{c})$, $n_{c}=m\omega^{2}_{l}/(4\pi e^{2})$ is the plasma
critical density. As shown, the amplitude of the reflected pulse is
effectively magnified, see Fig. 1 inset (a) for the electric field at the
pulse center $\Delta_{c}=0$, and rich high-order harmonics of the incident
frequency are superposed in the reflected pulse, see Fig. 1 inset (b). As we
will show later, this effective amplitude magnification and rich harmonic
generation could result in orders of magnitude enhancement of the particle
yield. The other main feature of the reflected pulse as shown in Fig. 1 inset
(a) is the evident field asymmetry, which would then lead to high spin
polarization for the created particles.

The full QED calculation for the BW process has been done and outlined as
follows: the positron spectrum produced by a high-energy photon in the scaled
potential $a(\phi)=-\int^{\phi}_{0}\xi(\phi^{\prime})\mathrm{d}\phi^{\prime}$
can be written as

$\displaystyle\frac{\mathrm{d}\textrm{P}}{\mathrm{d}s}$
$\displaystyle=\frac{i\alpha}{2\pi\eta}\int\mathrm{d}\varphi\int\frac{\mathrm{d}\vartheta}{\vartheta}~{}e^{\frac{i\vartheta\Lambda}{2\eta
s(1-s)}}\left(1-\vartheta^{2}\langle\xi\rangle^{2}h_{s}\right)\,,$ (1)

where $\alpha=e^{2}\approx 1/137$ is the fine structure constant, $s=k\cdot
q/k\cdot\ell$ is the fraction of the lightfront momentum taken by the
positron, $h_{s}=(2s^{2}-2s+1)/[4s(1-s)]$, $\varphi$ ($\vartheta$) is the
average (interference) phase and $\Lambda=1-\langle a\rangle^{2}+\langle
a^{2}\rangle$ is the Kibble mass Brown and Kibble (1964); Dinu et al. (2016);
Ilderton et al. (2019); Ilderton (2019b). The window average $\langle
f\rangle$ is calculated as $\langle
f\rangle=\frac{1}{\vartheta}\int^{\varphi+\vartheta/2}_{\varphi-\vartheta/2}f(\phi)\mathrm{d}\phi$.
In the collision geometry, the plasma reflected pulse propagates in the
$z$-direction with the magnetic (electric) field in the $y$ ($x$)-direction.
The spin of the created positrons is polarized in the direction of the
magnetic field as

$\displaystyle S_{y}(s)$
$\displaystyle=\frac{1}{\mathrm{d}\textrm{P}/\mathrm{d}s}\frac{\alpha}{2\pi\eta}\int\mathrm{d}\varphi\int\mathrm{d}\vartheta~{}e^{\frac{i\vartheta\Lambda}{2\eta(1-s)s}}\frac{\langle\xi\rangle}{2s}\,,$
(2)

with the total polarization $S_{y}$ acquired by integrating Eq. (2) without
$\mathrm{d}\textrm{P}/\mathrm{d}s$ but normalized with the total yield P. From
Eq. (2), we can infer that the spin of the created positrons in the pulse with
symmetric field must be unpolarized as $S_{y}\propto\langle\xi\rangle\approx
0$ Seipt et al. (2018); Sorbo et al. (2018). The detailed derivation and
calculation of Eqs. (1) and (2) can be found in Ref. Tang (2022). Here, the
plane-wave approximation is used as the scale of the field transverse
variation is much larger than the electron Compton wavelength and the
particles’ transverse displacement in the near head-on collision Di Piazza
(2015, 2016, 2017, 2021).

Figure 2: Breit-Wheeler pair creation in the plasma reflected pulse. (a)
Positron spectra created by the probe photons
($E_{\gamma}=5,10~{}\textrm{GeV}$) colliding head-on with the incident and
reflected pulses at the center $\Delta_{c}=0$. (b) Positron yield _vs_
transverse misalignment $\Delta_{c}$ in the incident and reflected pulses. The
yield of the linear Breit-Wheeler (LBW) process in the reflected pulse is also
given. (c) Contributions of the different frequency components in the
reflected pulse to the total spectrum created by a $10~{}\textrm{GeV}$ photon
at $\Delta_{c}=-3.5\lambda_{l}$. The spectrum in the initial pulse is plotted
for comparison. (d) Double differential spectrum
$\mathrm{d}^{2}\textrm{P}_{l}/(\mathrm{d}s\mathrm{d}\omega)$ and energy
spectrum $\mathrm{d}\textrm{P}_{l}/\mathrm{d}s$ of the LBW pair creation. The
magenta dashed line is its frequency threshold, $\omega_{*}=\omega_{l}/[2\eta
s(1-s)]$, and the vertical black dashed lines denote the boundary of the
contribution from different order of harmonics $\omega/\omega_{l}=17,18,19,20$
for $E_{\gamma}=10~{}\textrm{GeV}$. The collision between the plasma harmonics
and probe photons is $30\lambda_{l}$ away from the plasma surface.

In Fig. 2 (a), we compare the positron spectra created by the high-energy
photons ($E_{\gamma}=5,10~{}\textrm{GeV}$, corresponding to
$\eta=0.06,~{}0.12$ respectively) in the head-on collision with the incident
and reflected pulses at $\Delta_{c}=0$: the spectrum created in the reflected
pulse is much higher and decays much slower than that in the incident pulse.
This means significant yield enhancement with the plasma reflected pulse as
shown in Fig. 2 (b), in which the positron yield with the change of the
transverse misalignment $\Delta_{c}$ is given. For the $10~{}\textrm{GeV}$
photon, the enhancement at $\Delta_{c}=0$ is about $30$ times, and for
$E_{\gamma}=5~{}\textrm{GeV}$, the yield is enhanced for more than three
orders of magnitude and even slightly larger than that created by the
$10~{}\textrm{GeV}$ photon in the incident pulse.

This significant yield enhancement can be attributed to the rich harmonics
contained in the reflected pulse. In Fig. 2 (c), we show the positron spectra
created by the different frequency components in the reflected pulse by
filtering out the other frequency components artificially. For the relatively
lower-order harmonics $\omega/\omega_{l}<2/\eta$, which come from the Doppler
frequency upshift during the backward acceleration of the plasma surface
Teubner and Gibbon (2009), the creation is triggered by absorbing more than
one harmonic photons, $n\sim 2\omega_{l}/(\eta\omega)>1$, _i.e._ nonlinear BW
process. The spectrum of this nonlinear process, shown as the blue dashed line
in Fig. 2 (c), is much higher than the total spectrum (green dotted-dashed
line) obtained in the incident pulse, due to two factors: i) the considerable
decrease of the photon number absorbed for the creation, and ii) the effective
increase of the reflected field amplitude superposed by these phase-locking
low-order harmonics Quéré et al. (2008). For the higher-order harmonics
$\omega/\omega_{l}>2/\eta$, which come from the synchrotron radiation of the
well-accelerated surface layer Mikhailova et al. (2012) and improve the
collision energy parameter considerably, $\eta\to\omega\eta/\omega_{l}>2$, the
creation can thus be triggered with only one harmonic photon, _i.e._ linear BW
process, which is forbidden in the incident pulse with the normal frequency.
The spectrum of this linear process [blue dotted line in Fig. 2 (c)] could
also be much higher than the total spectrum obtained in the incident pulse.

We point out that the total spectrum [blue solid line in Fig. 2 (c)] acquired
in the reflected pulse can be divided numerically into the contributions from
the linear and nonlinear processes, _i.e._
$\textrm{P}\approx\textrm{P}_{\omega<2\omega_{l}/\eta}+\textrm{P}_{\omega>2\omega_{l}/\eta}$
as the intensity of high-order harmonics ($\omega/\omega_{l}>2/\eta$) is too
weak to affect the nonlinear process and the low-order harmonics
($\omega/\omega_{l}<2/\eta$) are excluded from the linear process by the
frequency threshold [see Eq. (3)]. The nonlinear process is more important in
the energy regime $s\to 0.5$ and determines the height of the total spectrum,
while the linear process is more important in the regime $s\to 1$ (and $s\to
0$) resulting in the slow spectral decay at the high (low) energy tail. Around
the pulse center $\Delta_{c}\to 0$ where the laser intensity is high, the pair
creation is dominated by the nonlinear BW process as shown in Fig. 2 (b). The
linear process is dominant in the region away from the pulse center, where the
multiphoton absorption is strongly suppressed by the lower intensity, and thus
results in the effective pair creation in a much broader transverse region
than that in the incident pulse, also shown in Fig. 2 (b). For lower-energy
photons, the importance of the linear BW process is higher as the nonlinear
process is limited with the absorption of more harmonic photons.

In Fig. 2 (c), a clear step-structure is observed in the spectrum of the
linear BW process. This structure results directly from the harmonic structure
in the frequency spectrum of the reflected pulse, and can be clearly seen via
Eq. (1) after doing the perturbative expansion ($\xi\to 0$) and keeping only
$\xi^{2}$ terms as:

$\displaystyle\frac{\mathrm{d}\textrm{P}_{l}}{\mathrm{d}s}$
$\displaystyle=\frac{\alpha\omega_{l}}{\pi\eta}\int^{\infty}_{\omega_{*}}\frac{\mathrm{d}\omega}{\omega^{2}}|\tilde{\xi}(\omega)|^{2}\left(h_{s}+\omega_{*}\frac{\omega-\omega_{*}}{\omega^{2}}\right)\,,$
(3)

which gives the positron spectrum [red solid line in Fig. 2 (d)] matching
exactly the full QED calculation for $\omega/\omega_{l}>2/\eta$ in Fig. 2 (c),
where
$\tilde{\xi}(\omega)=\int\mathrm{d}\phi\xi(\phi)e^{i\phi\omega/\omega_{l}}$,
and $\omega_{\ast}=\omega_{l}/[2\eta s(1-s)]$ is the threshold frequency
required to create a pair [magenta dashed line in Fig. 2 (d)]. From Eq. (3),
we can see that the harmonic peaks could bring the horizontal fringes into the
double differential spectrum,
$\mathrm{d}^{2}\textrm{P}_{l}/(\mathrm{d}s\mathrm{d}\omega)\propto|\tilde{\xi}(\omega)|^{2}$,
shown in Fig. 2 (d). After integrating over the pulse frequency, a clear step-
structure is imprinted into the positron spectrum as lower-order harmonics
cannot contribute at larger $s$, see the boundaries (vertical black dashed
lines) of the contribution from different order of harmonics in Fig. 2 (d).
This harmonic structure may provide the idea to observe the linear BW process
by colliding the high-energy photon with the plasma reflected pulse after
filtering out its low-frequency components. (To filter out low-frequency
components, one can refer to Tsakiris et al. (2006); Hörlein et al. (2010).)

In Fig. 2, the horizontal axis of the spectra are restricted in the region
$0.5<s<1$ due to the symmetry, $\textrm{P}(s)=\textrm{P}(1-s)$ in Eq. (1). The
slight asymmetry at $\Delta_{c}=0$ in Fig. 2 (b) for the yield in the
reflected pulse is because of the mild deformation of the plasma surface. The
convergence of the harmonics due to plasma denting surface Vincenti et al.
(2014), is not critical as the steep plasma pre-gradient and intermediately
intense laser pulse with broad transverse profile are applied. The collision
point between the reflected pulse and photon beam is $30\lambda_{l}$ away from
the plasma reflection point, and no notable difference is observed when the
collision moves further.

Figure 3: Robustness of the yield enhancement with the plasma generated
harmonics by considering multiple interaction parameters: (a) laser incident
angle $\theta$; (b) plasma pre-gradient $L$; (c) plasma bulk density $n_{e}$;
(d) plasma material. Fully ionized materials: C, Ca, Sr, Ba, Au with
respectively the charge number $Z=6,20,38,56,79$, are used; (e) laser
intensity $\xi_{l}$; (f) photon energy $E_{\gamma}$ and (g) Compton photon
source scattered by a $16.5~{}\textrm{GeV}$ electron colliding head-on with
the incident pulse. 1D PIC simulations are performed for harmonic generation.
If not specified in each panels, the following parameters: $\xi_{l}=3$,
$\theta=45^{\circ}$, $L=\lambda_{l}/8$, $n_{e}=15n_{c}$ and $Z=79$ are used.
In (e), the plasma bulk density increases with the laser intensity as
$n_{e}=5\xi_{l}n_{c}$. The corresponding results in the incident pulse without
harmonics are given for comparison.

The robustness of the yield enhancement with the plasma reflected pulse is
tested in Fig. 3 for multiple interaction parameters. The plasma harmonic
generation is simulated with series of 1D PIC setups Bourdier (1983) for
weakly focused pulses $\sigma\to\infty$. The positron yield in the plasma
reflected pulse is significantly enhanced for a broad range of incident angle
$\theta$ and plasma pre-gradient $L$ as shown in Figs. 3 (a) and (b) with the
peak enhancement at around $\theta=45^{\circ}$ and $L=\lambda_{l}/10$, at
which the efficiency of the harmonic generation is maximized Thaury and Quéré
(2010); Bhadoria et al. (2022). For the highly overdense plasma $n_{e}\gg
n_{c}$ with a proper pre-gradient, its bulk density and target material may
not affect the efficiency of harmonic generation considerably, but only affect
the peak widening of high-order harmonics via the hole-boring effect Tang et
al. (2017). The yield enhancement is thus nearly stable for the high bulk
density in Fig. 3 (c) and different target materials in Fig. 3 (d). For the
relatively lower density $n_{e}\sim\xi_{l}n_{c}$ in Fig. 3 (c), the yield
enhancement becomes smaller as the harmonic generation is reduced with the
lower plasma current. In Figs. 3 (e) and (f), the level of the yield
enhancement with the plasma harmonics is declined for the larger laser
intensity and higher photon energy. This may be because i) the decreasing of
the photon number needed for the creation due to the harmonics becomes less
important, since more laser photons are provided with the larger laser
intensity and less laser photons are needed for the higher photon energy, and
ii) the importance of the linear BW contribution from the high-order harmonics
becomes weaker. However, in the intermediate intensity region $\xi_{l}\sim
O(1)$ with the moderate photon energy $E_{\gamma}\lesssim 10~{}\textrm{GeV}$,
the yield enhancement with the plasma harmonics is still significant as shown
in Fig. 3 (e): at $\xi_{l}=5$, the yield is enhanced by about $5$ times for
$E_{\gamma}=10~{}\textrm{GeV}$ and about $40$ times for
$E_{\gamma}=5~{}\textrm{GeV}$. This implies that the plasma-generated
harmonics have the great advantage for the current laser-particle experiments,
in which the dominant energy of the potential photon sources is much smaller
than $10~{}\textrm{GeV}$ Abramowicz et al. (2021); Naranjo et al. (2021); Chen
et al. (2022); MacLeod et al. (2023). In Fig. 3 (g), the spectrum,
$\mathrm{d}\textrm{P}_{\gamma}/\mathrm{d}E_{\gamma}$ of the Compton photon
source scattered by a $16.5~{}\textrm{GeV}$ electron colliding with the
incident pulse ($\xi_{l}=3$) is considered King and Tang (2020); Tang et al.
(2020); Tang and King (2023). As shown, with the plasma reflected pulse, the
energy spectrum,
$\textrm{P}(E_{\gamma})~{}\mathrm{d}\textrm{P}_{\gamma}/\mathrm{d}E_{\gamma}$
of the created positrons is significantly improved with the total yield
enhanced for more than $20$ times. The final yield could be further improved
if the laser pulse for the Compton scattering is also replaced with the
plasma-generated harmonics.

Figure 4: Spin polarization of the created positrons by the high-energy
photons with $E_{\gamma}=5,~{}10~{}\textrm{GeV}$ in the plasma reflected pulse
with the incident angle $\theta=45^{\circ}$. (a) Spectrum of the positron’s
spin polarization at $\Delta_{c}=0$. (b) Variation of the total spin
polarization with the change of the misalignment $\Delta_{c}$. The same
parameters are used as in Fig. 2.

The other advantage with the plasma-generated harmonics is to create the
positron beam with the high spin polarization. This is obtained due to the
evident field asymmetry as $S_{y}\propto\langle\xi\rangle$ in Eq. (2), when
the laser oblique incidence is applied for the harmonic generation, see Fig. 1
inset (a) for the incident angle $\theta=45^{\circ}$. The generation of the
harmonics in the oblique incidence is pronounced in one of the half laser
cycle in which the oscillation of the surface current is promoted by the laser
electric field, and suppressed in the other half cycle as the current
oscillation is weakened by the electric field Teubner and Gibbon (2009). As
shown in Fig. 4 for $\theta=45^{\circ}$, the spin of the created positrons is
highly polarized in the $+y$ direction, in which the magnetic field is
magnified [The electric field is magnified in the $+x$ direction, see Fig. 1
inset (a)]. In Fig. 4 (a), the spin polarization of the created positron with
the peak higher than $80\%$ around $s=0.2$, decreases with the increase of the
energy as $S_{y}(s)\propto s^{-1}$ in Eq. (2), and at the low (high) energy
limit $s\to 0$ ($s\to 1$), the spin polarization goes to zero because of the
dominance of the linear BW contribution. This can be seen via the perturbative
approximation, $S_{y}(s)\propto\xi^{3}\to 0$ at $\xi\to 0$. As the pairs are
more probable to be created around $s=0.5$, the total polarization of the
created positrons around the pulse center is about $62\%$ shown in Fig. 4 (b).
In the region away from the pulse center, the spin polarization decays to zero
quickly. As the asymmetry of the harmonic field depends on the laser incident
angle, the positrons’ spin polarization may be simply controlled by adjusting
the laser incidence.

In conclusion, we detailed the exploration of the Breit-Wheeler pair creation
using the plasma generated harmonics. The rich harmonics can not only improve
the yield of the nonlinear Breit-Wheeler process significantly, but also open
an effective channel for the linear Breit-Wheeler process. In the currently
available region of laser intensity and photon energy, the particle yield
could be improved for more than one order of magnitude. This proposal could be
realized at upcoming experiments and used as a robust scheme to generate
brilliant, highly polarized positron beams. We point out that the efforts, to
improve the harmonic-generation efficiency Edwards and Mikhailova (2016);
Yeung et al. (2017); Cantono et al. (2018); Li et al. (2022), can also be
devoted here to further enhance the particle yield.

###### Acknowledgements.

The author thanks Dr. A. Ilderton and Dr. B. King for useful discussion and
careful reading of the manuscript. The author acknowledges the support from
the Natural Science Foundation of China, Grant No.12104428. The work was
carried out at Marine Big Data Center of Institute for Advanced Ocean Study of
Ocean University of China.

## References

  * Breit and Wheeler (1934) G. Breit and J. A. Wheeler, Phys. Rev. 46, 1087 (1934), URL https://link.aps.org/doi/10.1103/PhysRev.46.1087.
  * Ritus (1985) V. I. Ritus, J. Russ. Laser Res. 6, 497 (1985), URL https://doi.org/10.1007/BF01120220.
  * Di Piazza et al. (2012) A. Di Piazza, C. Müller, K. Z. Hatsagortsyan, and C. H. Keitel, Rev. Mod. Phys. 84, 1177 (2012), URL https://link.aps.org/doi/10.1103/RevModPhys.84.1177.
  * Fedotov et al. (2023) A. Fedotov, A. Ilderton, F. Karbstein, B. King, D. Seipt, H. Taya, and G. Torgrimsson, Physics Reports 1010, 1 (2023), ISSN 0370-1573, advances in QED with intense background fields, URL https://www.sciencedirect.com/science/article/pii/S0370157323000352.
  * Gonoskov et al. (2022) A. Gonoskov, T. G. Blackburn, M. Marklund, and S. S. Bulanov, Rev. Mod. Phys. 94, 045001 (2022), URL https://link.aps.org/doi/10.1103/RevModPhys.94.045001.
  * Joshi et al. (2018) C. Joshi, E. Adli, W. An, C. E. Clayton, S. Corde, S. Gessner, M. J. Hogan, M. Litos, W. Lu, K. A. Marsh, et al., Plasma Physics and Controlled Fusion 60, 034001 (2018), URL https://dx.doi.org/10.1088/1361-6587/aaa2e3.
  * Jacobs (2022) R. Jacobs, arXiv e-prints (2022), eprint 2206.00403, URL https://doi.org/10.48550/arXiv.2206.00403.
  * Borysova (2021) M. Borysova, Journal of Instrumentation 16, C12030 (2021), URL https://doi.org/10.1088/1748-0221/16/12/c12030.
  * Abramowicz et al. (2021) H. Abramowicz et al., Eur. Phys. J. Spec. Top. 230, 2445 (2021), URL https://doi.org/10.1140/epjs/s11734-021-00249-z.
  * Naranjo et al. (2021) B. Naranjo, G. Andonian, N. Cavanagh, A. D. Piazza, A. Fukasawa, E. Gerstmayr, R. Holtzapple, C. Keitel, N. Majernik, S. Meuren, et al., THPAB270 (2021), ISSN 2673-5490, URL https://jacow.org/ipac2021/papers/thpab270.pdf.
  * Chen et al. (2022) Z. Chen, S. Meuren, E. Gerstmayr, V. Yakimenko, P. H. Bucksbaum, and D. A. R. and, in _Optica High-brightness Sources and Light-driven Interactions Congress 2022_ (Optica Publishing Group, 2022), p. HF4B.6, URL http://opg.optica.org/abstract.cfm?URI=HILAS-2022-HF4B.6.
  * Macleod (2022) A. J. Macleod, Journal of Physics: Conference Series 2249, 012022 (2022), URL https://dx.doi.org/10.1088/1742-6596/2249/1/012022.
  * Turner et al. (2022) M. Turner, S. S. Bulanov, C. Benedetti, A. J. Gonsalves, W. P. Leemans, K. Nakamura, J. van Tilborg, C. B. Schroeder, C. G. R. Geddes, and E. Esarey, Eur. Phys. J. D 76, 205 (2022), eprint 2210.09214, URL https://doi.org/10.1140/epjd/s10053-022-00535-y.
  * Sung et al. (2017) J. H. Sung, H. W. Lee, J. Y. Yoo, J. W. Yoon, C. W. Lee, J. M. Yang, Y. J. Son, Y. H. Jang, S. K. Lee, and C. H. Nam, Opt. Lett. 42, 2058 (2017), URL https://opg.optica.org/ol/abstract.cfm?URI=ol-42-11-2058.
  * Kiriyama et al. (2018) H. Kiriyama, A. S. Pirozhkov, M. Nishiuchi, Y. Fukuda, K. Ogura, A. Sagisaka, Y. Miyasaka, M. Mori, H. Sakaki, N. P. Dover, et al., Opt. Lett. 43, 2595 (2018), URL https://opg.optica.org/ol/abstract.cfm?URI=ol-43-11-2595.
  * Shen et al. (2018) B. Shen, Z. Bu, J. Xu, T. Xu, L. Ji, R. Li, and Z. Xu, Plasma Physics and Controlled Fusion 60, 044002 (2018), URL https://dx.doi.org/10.1088/1361-6587/aaa7fb.
  * Danson et al. (2019) C. N. Danson, C. L. Haefner, J. Bromage, T. J. Butcher, J.-C. Chanteloup, E. Chowdhury, A. D. Galvanauskas, L. A. Gizzi, J. Hein, D. I. Hillier, et al., High Power Laser Science and Engineering 7, e54 (2019), URL https://www.cambridge.org/core/journals/high-power-laser-science-and-engineering/article/petawatt-and-exawatt-class-lasers-worldwide/D85E3F002CBFC286C08076E127BB5F5C.
  * Yoon et al. (2021) J. W. Yoon, Y. G. Kim, I. W. Choi, J. H. Sung, H. W. Lee, S. K. Lee, and C. H. Nam, Optica 8, 630 (2021), URL https://opg.optica.org/optica/abstract.cfm?URI=optica-8-5-630.
  * Nikishov and Ritus (1964) A. I. Nikishov and V. I. Ritus, Sov. Phys. JETP 19, 529 (1964).
  * Jansen and Müller (2013) M. J. A. Jansen and C. Müller, Phys. Rev. A 88, 052125 (2013), URL https://link.aps.org/doi/10.1103/PhysRevA.88.052125.
  * Meuren et al. (2016) S. Meuren, C. H. Keitel, and A. Di Piazza, Phys. Rev. D 93, 085028 (2016), URL https://link.aps.org/doi/10.1103/PhysRevD.93.085028.
  * Jansen et al. (2016) M. J. A. Jansen, J. Z. Kamiński, K. Krajewska, and C. Müller, Phys. Rev. D 94, 013010 (2016), URL https://link.aps.org/doi/10.1103/PhysRevD.94.013010.
  * Titov et al. (2018) A. I. Titov, H. Takabe, and B. Kämpfer, Phys. Rev. D 98, 036022 (2018), URL https://link.aps.org/doi/10.1103/PhysRevD.98.036022.
  * Ilderton (2020) A. Ilderton, Phys. Rev. D 101, 016006 (2020), URL https://link.aps.org/doi/10.1103/PhysRevD.101.016006.
  * Ilderton (2019a) A. Ilderton, Phys. Rev. D 100, 125018 (2019a), URL https://link.aps.org/doi/10.1103/PhysRevD.100.125018.
  * Wistisen (2020) T. N. Wistisen, Phys. Rev. D 101, 076017 (2020), URL https://link.aps.org/doi/10.1103/PhysRevD.101.076017.
  * Seipt and King (2020) D. Seipt and B. King, Phys. Rev. A 102, 052805 (2020), URL https://link.aps.org/doi/10.1103/PhysRevA.102.052805.
  * Mercuri-Baron et al. (2021) A. Mercuri-Baron, M. Grech, F. Niel, A. Grassi, M. Lobet, A. D. Piazza, and C. Riconda, New Journal of Physics 23, 085006 (2021), URL https://dx.doi.org/10.1088/1367-2630/ac1975.
  * Golub et al. (2021) A. Golub, S. Villalba-Chávez, H. Ruhl, and C. Müller, Phys. Rev. D 103, 016009 (2021), URL https://link.aps.org/doi/10.1103/PhysRevD.103.016009.
  * Tang (2021) S. Tang, Phys. Rev. A 104, 022209 (2021), URL https://link.aps.org/doi/10.1103/PhysRevA.104.022209.
  * Tang and King (2021) S. Tang and B. King, Phys. Rev. D 104, 096019 (2021), URL https://link.aps.org/doi/10.1103/PhysRevD.104.096019.
  * Tang (2022) S. Tang, Phys. Rev. D 105, 056018 (2022), URL https://link.aps.org/doi/10.1103/PhysRevD.105.056018.
  * Gao and Tang (2022) Y. Gao and S. Tang, Phys. Rev. D 106, 056003 (2022), URL https://link.aps.org/doi/10.1103/PhysRevD.106.056003.
  * Borysov et al. (2022) O. Borysov, B. Heinemann, A. Ilderton, B. King, and A. Potylitsyn, Phys. Rev. D 106, 116015 (2022), URL https://link.aps.org/doi/10.1103/PhysRevD.106.116015.
  * Seipt et al. (2019) D. Seipt, D. Del Sorbo, C. P. Ridgers, and A. G. R. Thomas, Phys. Rev. A 100, 061402 (2019), URL https://link.aps.org/doi/10.1103/PhysRevA.100.061402.
  * Chen et al. (2019) Y.-Y. Chen, P.-L. He, R. Shaisultanov, K. Z. Hatsagortsyan, and C. H. Keitel, Phys. Rev. Lett. 123, 174801 (2019), URL https://link.aps.org/doi/10.1103/PhysRevLett.123.174801.
  * Wan et al. (2020) F. Wan, R. Shaisultanov, Y.-F. Li, K. Z. Hatsagortsyan, C. H. Keitel, and J.-X. Li, Physics Letters B 800, 135120 (2020), ISSN 0370-2693, URL https://www.sciencedirect.com/science/article/pii/S0370269319308421.
  * Moortgat-Pick et al. (2008) G. Moortgat-Pick, T. Abe, G. Alexander, B. Ananthanarayan, A. Babich, V. Bharadwaj, D. Barber, A. Bartl, A. Brachmann, S. Chen, et al., Physics Reports 460, 131 (2008), ISSN 0370-1573, URL https://www.sciencedirect.com/science/article/pii/S0370157308000136.
  * Vauth and List (2016) A. Vauth and J. List, International Journal of Modern Physics: Conference Series 40, 1660003 (2016), eprint https://doi.org/10.1142/S201019451660003X, URL https://doi.org/10.1142/S201019451660003X.
  * Lichters et al. (1996) R. Lichters, J. Meyer-ter Vehn, and A. Pukhov, Physics of Plasmas 3, 3425 (1996), eprint https://doi.org/10.1063/1.871619, URL https://doi.org/10.1063/1.871619.
  * Baeva et al. (2006) T. Baeva, S. Gordienko, and A. Pukhov, Phys. Rev. E 74, 046404 (2006), URL https://link.aps.org/doi/10.1103/PhysRevE.74.046404.
  * Teubner and Gibbon (2009) U. Teubner and P. Gibbon, Rev. Mod. Phys. 81, 445 (2009), URL https://link.aps.org/doi/10.1103/RevModPhys.81.445.
  * Tang and Kumar (2019) S. Tang and N. Kumar, Plasma Physics and Controlled Fusion 61, 025013 (2019), URL https://doi.org/10.1088%2F1361-6587%2Faaf378.
  * Gordienko et al. (2004) S. Gordienko, A. Pukhov, O. Shorokhov, and T. Baeva, Phys. Rev. Lett. 93, 115002 (2004), URL https://link.aps.org/doi/10.1103/PhysRevLett.93.115002.
  * Gonoskov et al. (2011) A. A. Gonoskov, A. V. Korzhimanov, A. V. Kim, M. Marklund, and A. M. Sergeev, Phys. Rev. E 84, 046403 (2011), URL https://link.aps.org/doi/10.1103/PhysRevE.84.046403.
  * Bulanov et al. (2013) S. V. Bulanov, T. Z. Esirkepov, M. Kando, A. S. Pirozhkov, and N. N. Rosanov, Physics-Uspekhi 56, 429 (2013), URL https://dx.doi.org/10.3367/UFNe.0183.201305a.0449.
  * Arber et al. (2015) T. D. Arber, K. Bennett, C. S. Brady, A. Lawrence-Douglas, M. G. Ramsay, N. J. Sircombe, P. Gillies, R. G. Evans, H. Schmitz, A. R. Bell, et al., Plasma Physics and Controlled Fusion 57, 113001 (2015), URL https://dx.doi.org/10.1088/0741-3335/57/11/113001.
  * Brown and Kibble (1964) L. S. Brown and T. W. B. Kibble, Phys. Rev. 133, A705 (1964), URL https://link.aps.org/doi/10.1103/PhysRev.133.A705.
  * Dinu et al. (2016) V. Dinu, C. Harvey, A. Ilderton, M. Marklund, and G. Torgrimsson, Phys. Rev. Lett. 116, 044801 (2016), URL https://link.aps.org/doi/10.1103/PhysRevLett.116.044801.
  * Ilderton et al. (2019) A. Ilderton, B. King, and D. Seipt, Phys. Rev. A 99, 042121 (2019), URL https://link.aps.org/doi/10.1103/PhysRevA.99.042121.
  * Ilderton (2019b) A. Ilderton, Phys. Rev. D 99, 085002 (2019b), URL https://link.aps.org/doi/10.1103/PhysRevD.99.085002.
  * Seipt et al. (2018) D. Seipt, D. Del Sorbo, C. P. Ridgers, and A. G. R. Thomas, Phys. Rev. A 98, 023417 (2018), URL https://link.aps.org/doi/10.1103/PhysRevA.98.023417.
  * Sorbo et al. (2018) D. D. Sorbo, D. Seipt, A. G. R. Thomas, and C. P. Ridgers, Plasma Physics and Controlled Fusion 60, 064003 (2018), URL https://dx.doi.org/10.1088/1361-6587/aab979.
  * Di Piazza (2015) A. Di Piazza, Phys. Rev. A 91, 042118 (2015), URL https://link.aps.org/doi/10.1103/PhysRevA.91.042118.
  * Di Piazza (2016) A. Di Piazza, Phys. Rev. Lett. 117, 213201 (2016), URL https://link.aps.org/doi/10.1103/PhysRevLett.117.213201.
  * Di Piazza (2017) A. Di Piazza, Phys. Rev. A 95, 032121 (2017), URL https://link.aps.org/doi/10.1103/PhysRevA.95.032121.
  * Di Piazza (2021) A. Di Piazza, Phys. Rev. D 103, 076011 (2021), URL https://link.aps.org/doi/10.1103/PhysRevD.103.076011.
  * Quéré et al. (2008) F. Quéré, C. Thaury, J.-P. Geindre, G. Bonnaud, P. Monot, and P. Martin, Phys. Rev. Lett. 100, 095004 (2008), URL https://link.aps.org/doi/10.1103/PhysRevLett.100.095004.
  * Mikhailova et al. (2012) J. M. Mikhailova, M. V. Fedorov, N. Karpowicz, P. Gibbon, V. T. Platonenko, A. M. Zheltikov, and F. Krausz, Phys. Rev. Lett. 109, 245005 (2012), URL https://link.aps.org/doi/10.1103/PhysRevLett.109.245005.
  * Tsakiris et al. (2006) G. D. Tsakiris, K. Eidmann, J. M. ter Vehn, and F. Krausz, New Journal of Physics 8, 19 (2006), URL https://dx.doi.org/10.1088/1367-2630/8/1/019.
  * Hörlein et al. (2010) R. Hörlein, Y. Nomura, P. Tzallas, S. G. Rykovanov, B. Dromey, J. Osterhoff, Z. Major, S. Karsch, L. Veisz, M. Zepf, et al., New Journal of Physics 12, 043020 (2010), URL https://dx.doi.org/10.1088/1367-2630/12/4/043020.
  * Vincenti et al. (2014) H. Vincenti, S. Monchocé, S. Kahaly, G. Bonnaud, P. Martin, and F. Quéré, Nature communications 5, 3403 (2014), URL https://doi.org/10.1038/ncomms4403.
  * Bourdier (1983) A. Bourdier, The Physics of Fluids 26, 1804 (1983), eprint https://aip.scitation.org/doi/pdf/10.1063/1.864355, URL https://aip.scitation.org/doi/abs/10.1063/1.864355.
  * Thaury and Quéré (2010) C. Thaury and F. Quéré, Journal of Physics B: Atomic, Molecular and Optical Physics 43, 213001 (2010), URL https://dx.doi.org/10.1088/0953-4075/43/21/213001.
  * Bhadoria et al. (2022) S. Bhadoria, T. Blackburn, A. Gonoskov, and M. Marklund, Physics of Plasmas 29 (2022), ISSN 1070-664X, 093109, eprint https://pubs.aip.org/aip/pop/article-pdf/doi/10.1063/5.0087854/16621798/093109_1_online.pdf, URL https://doi.org/10.1063/5.0087854.
  * Tang et al. (2017) S. Tang, N. Kumar, and C. H. Keitel, Phys. Rev. E 95, 051201 (2017), URL https://link.aps.org/doi/10.1103/PhysRevE.95.051201.
  * MacLeod et al. (2023) A. J. MacLeod, P. Hadjisolomou, T. M. Jeong, and S. V. Bulanov, Phys. Rev. A 107, 012215 (2023), URL https://link.aps.org/doi/10.1103/PhysRevA.107.012215.
  * King and Tang (2020) B. King and S. Tang, Phys. Rev. A 102, 022809 (2020), URL https://link.aps.org/doi/10.1103/PhysRevA.102.022809.
  * Tang et al. (2020) S. Tang, B. King, and H. Hu, Physics Letters B 809, 135701 (2020), ISSN 0370-2693, URL https://www.sciencedirect.com/science/article/pii/S0370269320305049.
  * Tang and King (2023) S. Tang and B. King, Phys. Rev. D 107, 096004 (2023), URL https://link.aps.org/doi/10.1103/PhysRevD.107.096004.
  * Edwards and Mikhailova (2016) M. R. Edwards and J. M. Mikhailova, Phys. Rev. Lett. 117, 125001 (2016), URL https://link.aps.org/doi/10.1103/PhysRevLett.117.125001.
  * Yeung et al. (2017) M. Yeung, S. Rykovanov, J. Bierbach, L. Li, E. Eckner, S. Kuschel, A. Woldegeorgis, C. Rödel, A. Sävert, G. Paulus, et al., Nature Photonics 11, 32 (2017), URL https://doi.org/10.1038/nphoton.2016.239.
  * Cantono et al. (2018) G. Cantono, L. Fedeli, A. Sgattoni, A. Denoeud, L. Chopineau, F. Réau, T. Ceccotti, and A. Macchi, Phys. Rev. Lett. 120, 264803 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.120.264803.
  * Li et al. (2022) B. Y. Li, F. Liu, M. Chen, F. Y. Wu, J. W. Wang, L. Lu, J. L. Li, X. L. Ge, X. H. Yuan, W. C. Yan, et al., Phys. Rev. Lett. 128, 244801 (2022), URL https://link.aps.org/doi/10.1103/PhysRevLett.128.244801.