#  Resolving Multiphoton Processes with High-Order Anisotropy Ultrafast X-ray
Scattering

Adi Natan<EMAIL_ADDRESS>Stanford PULSE Institute, SLAC National
Accelerator Laboratory
2575 Sand Hill Road, Menlo Park, CA 94025  Aviad Schori  Stanford PULSE
Institute, SLAC National Accelerator Laboratory
2575 Sand Hill Road, Menlo Park, CA 94025  Grace Owolabi  Department
Electrical Engineering and Computer Science, Howard University, Washington DC
20059, United States of America  James P. Cryan  Stanford PULSE Institute,
SLAC National Accelerator Laboratory
2575 Sand Hill Road, Menlo Park, CA 94025  Linac Coherent Lightsource, SLAC
National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA 94025,
USA  James M. Glownia  Linac Coherent Lightsource, SLAC National Accelerator
Laboratory, 2575 Sand Hill Road, Menlo Park, CA 94025, USA  Philip H.
Bucksbaum  Stanford PULSE Institute, SLAC National Accelerator Laboratory
2575 Sand Hill Road, Menlo Park, CA 94025  Department of Physics, Stanford
University, Stanford, CA 94305  Department of Applied Physics, Stanford
University, Stanford, CA 94305

###### Abstract

We present first results on ultrafast X-ray scattering of strongly driven
molecular Iodine and analysis of high-order anisotropic components of the
scattering signal, up to four-photon absorption. We discuss the technical
details of retrieving high fidelity high-order anisotropy components, and
outline a method to analyze the scattering signal using Legendre
decomposition. We use simulated anisotropic scattering signals and Fourier
analysis to map how anisotropic dissociation motions can be extracted from the
various Legendre orders. We observe multitude dissociation and vibration
motions simultaneously arising from various multiphoton transitions. We use
the anisotropy information of the scattering signal to disentangle the
different processes and assign their dissociation velocities on the Angstrom
and femtosecond scales de-novo.

Capturing motions in space and time at the atomic scale is fundamental to the
understanding of chemical reactions and structural dynamics of molecules of
different complexities. Some of the emerging tools that allow such studies are
ultrafast scattering modalities, primarily by X-rays and relativistic
electrons, that were made feasible in recent years. In these studies, the
motions from excited molecules are captured in a pump-probe scheme, where the
excitation pump pulse is often an ultra-short linearly polarized optical laser
pulse with a duration shorter than the typical timescales of motion of
interest, as illustrated in Fig 1. The scattering signal is then usually
integrated over angle for improved fidelity and subtracted from the scattering
signal of the unexcited system, to allow tracing changes of signal positions
and infer structural dynamics Ihee  _et al._ (2010). This approach was used
successfully to demonstrate coherent motions and dynamics in molecules in the
gas phase Minitti  _et al._ (2015); Stankus  _et al._ (2019), as well as
structural changes in molecules in solution after the electronic excitation
Biasin  _et al._ (2016); Van Driel  _et al._ (2016); Chollet  _et al._ (2015).
In these experiments, the optical pulse parameters were carefully chosen to
ensure that only single-photon absorption is taking place and that the
molecular system is photoexcited to a specific electronic state. This is often
done by varying the pump intensity in the pump-probe setup and finding
conditions where linearity of the excitation, as manifested by the scattering
signal is achieved. The motion is inferred by measuring the scattering
difference as a function of delay and applying modeling and simulations. Using
angle integrated scattering signal is well justified as it captures all types
of motions that take place in the photoexcited system. However, angle-
dependent signals can be dramatically attenuated if only the isotropic
component is being analyzed. In many cases, there is an inherent anisotropy in
the scattering signal when a sample is excited by linearly polarized light due
to an optically induced dipole moment transition. This interaction creates
geometric alignment in the ensemble, and can be used to filter and enhance the
specific processes under study, such as in the case of a single-photon
absorption process Glownia  _et al._ (2016); Biasin  _et al._ (2018); Haldrup
_et al._ (2019), as well as perturbative two-photon excitation Baskin and
Zewail (2006); Yang  _et al._ (2018).

Higher orders of anisotropy play a significant role in understanding and
probing cases where the molecular system is in the presence of multi-photon
absorption and strong laser fields, such as dissociation due to bond softening
Bucksbaum  _et al._ (1990), above-threshold dissociation McKenna  _et al._
(2008), quantum coherent control Natan  _et al._ (2012), and light-induced
conical intersections Natan  _et al._ (2016). In addition, the interaction of
ultrashort pulses with molecules with anisotropic polarizability will generate
non-adiabatic (or impulsive) alignment Rosca-Pruna and Vrakking (2001). The
broad bandwidth of an ultrashort pulse creates rotational wavepackets that
evolve and rephase at periodic time delays, forming molecular alignment,
manifested by high order anisotropy in the sample under field-free conditions.
Molecular alignment is often used to probe diverse phenomena in the molecular
frame, such as, polyatomic vibrational dynamics and fragmentation Bisgaard
_et al._ (2009), ultrafast molecular frame electronic coherences Makhija  _et
al._ (2020), laser-induced rotational dynamics and control Karamatskos  _et
al._ (2019); Larsen, Wendt-Larsen, and Stapelfeldt (1999), Auger decay of
double core-hole statesCryan  _et al._ (2010), high-harmonic generation from
inner valence orbitals McFarland  _et al._ (2008) and diffractive imaging in
the molecular frame Yang  _et al._ (2015, 2016); Kierspel  _et al._ (2020);
Xiong, Wilkin, and Centurion (2020); Kierspel  _et al._ (2015).

Figure 1: Schematic description of a pump-probe experimental setup and
relevant coordinates presented in the text. A linearly polarized ultrashort
optical laser pulse pumps the molecular sample. The sample is then probed at
some time delay by an ultrashort X-ray pulse via scattering on a 2D detector
array with a finite range of scattering angles. The anisotropy induced in the
sample following photoexcitation is manifested by the anisotropy of the
scattering pattern.

For the case of an ensemble of diatomic molecular Iodine that will be
discussed in this work, multi-photon absorption will excite wavepackets from
its ground state to a multitude of states and pathways Brown  _et al._ (1999);
Pastirk, Lozovoy, and Dantus (2001); Fang and Gibson (2008), including three
dissociation limits, as well as higher bound Rydberg and ion-pair (IP) states,
with crossing occurring between them Bogomolov  _et al._ (2014); Lukashov,
Petrov, and Pravilov (2018). Such excitation will create time-dependent
angular distributions that will carry information regarding the number of
photons absorbed, the symmetry of the states involved, and the type of motion
that is taking place. Fig 2 describes some of the relevant potential curves
and states mentioned.

In a single-photon absorption ($\lambda=520nm$) process, the dynamics is
limited to transitions from the ground $X0^{+}_{g}$ state to the first
dissociation limit $I(^{3}P_{3/2})+I(^{3}P_{3/2})$ via states such as the
$A1_{u},C(B^{{}^{\prime\prime}})1_{u},B^{{}^{\prime}}(0_{u})$, or to vibration
motion via excitation to the bound $B0^{+}_{u}$ state. Upon absorbing two
photons, non-resonant and resonant Raman processes can take place to excite a
vibration wavepacket at the ground state via $X\leftarrow M\leftarrow X$ type
transition with $M$ the allowed symmetry intermediate state for the resonant
case. Furthermore, two-photon absorption can lead to dissociation at all the
dissociation limits, for example by exciting to the $C1_{g}$ and $0^{+}_{g}$
states. At three-photon absorption, excitation of Rydberg and IP states become
accessible, nonlinear Raman processes can take place to excite and mix various
lower energy states, and the ionization threshold is reached via four-photon
absorption

The excited rovibrational motion will undergo rotational dephasing that will
limit the delay time window such induced anisotropy can be detected, while
prompt dissociation along a single dissociative state will preserve the
anisotropy.

Figure 2: Potential energy curves of several valence, Rydberg and ion-pair
(IP) states of the Iodine molecule Lukashov, Petrov, and Pravilov (2018).
Arrows represent the energy of the laser excitation (520 nm) for the case of
single and multi-photon absorption.

Here, we present the first results for the case of a strongly driven molecular
Iodine vapor, where multiphoton processes take place, beyond the perturbative
single or two-photon absorption processes. We show that we can analyze and
retrieve time-dependent high-fidelity high-order anisotropy information of an
ultrafast X-ray scattering signal and assign motions of various excitation
processes that take place simultaneously de-novo.

##  I Theory

In this section we shall discuss how to link between the spatial-temporal
information of a photoexcited molecular system from a wavepacket perspective,
with its observed scattering pattern, using some results that have been
derived similarly before in other studies Ben-Nun, Cao, and Wilson (1997);
Lorenz, Møller, and Henriksen (2010a, b). We start with the assumption that we
have solved a time-dependent Schrodinger equation (TDSE) and have the exact
charge density of an excited system. We express the charge density by an
optically excited nuclear wave packet that is propagated along some electronic
state.

We then obtain expressions for the scattering pattern this wavepacket
generates. We will be interested in the general case of an n-photon absorption
that induces a non-trivial anisotropy and will result in up to $2n$ order
Legendre polynomials terms in the scattering signal. We also simplify the
treatment assuming both the X-ray and optical laser pulse are co-propagating
along $\vec{k}_{0}$, and that the angle between the incident X-ray beam and
the laser polarization is $\pi/2$, allowing the Legendre decomposition
approach in angle space. Scattering from time-evolving systems is generally
inelastic Lorenz, Møller, and Henriksen (2010a); Simmermacher  _et al._
(2019), however, assuming typical experimental conditions, we can use the so-
called static approximation Lorenz, Møller, and Henriksen (2010a) and replace
the general electronic scattering operator with its elastic expression for the
differential cross-section for scattering:

$\frac{d\sigma}{d\Omega}=\sigma_{T}\int
d\vec{R}\rho(\vec{R},t)|F(\vec{q},\vec{R})|^{2}$ (1)

where $\sigma_{T}$ is the Thomson cross-section, $\rho$ is the charge density,
and the squared molecular form factor $|F|^{2}$ is invariant under space
inversion $R\mapsto-R$. If we further simplify the treatment and consider only
a single atom pair, as in the case of diatomic Iodine, we can then write the
differential cross-section:

$\frac{d\sigma}{d\Omega}=\sigma_{T}f_{1}(q)f_{2}(q)\int
d\vec{R}\rho(\vec{R},t)(e^{\imath\vec{q}\vec{R}}+e^{-\imath\vec{q}\vec{R}})$
(2)

with $f_{i}(q)$ the ith atomic form factor. We can expand the scattering
exponential term using the plane wave expansion:

$e^{\imath\vec{q}\vec{R}}=\sum_{n=0,2,\ldots}(2n+1)(-1)^{n/2}P_{n}(cos\theta_{qR})j_{n}(qR)$
(3)

where $P_{n}$ are Legendre polynomials, $j_{n}$ are spherical Bessel
functions, and $\theta_{qR}$ is the angle between $\vec{q}$ and $\vec{R}$. We
only need to sum over the even orders as the components with odd polynomials
are anti-symmetric under space inversion and will cancel when we will
calculate the differential cross-section.

We use the spherical harmonics addition theorem to expand the Legendre
polynomials to express $\theta_{qR}$ in terms of the experimentally measured
scattering angles $(\theta_{q},\phi_{q}$), and $(\theta,\phi)$, the angle
between the laser polarization and molecular axis and its corresponding
azimuth:

$P_{n}(cos\theta_{qR})=\frac{4\pi}{2n+1}\sum\limits_{j=-k}^{k}Y_{jk}^{*}(\theta_{q},\phi_{q})Y_{jk}(\theta,\phi)$
(4)

Using the scattering symmetry, we can integrate over $\phi_{q}$ and use the
expressions above to arrive to:

$\frac{d\sigma}{d\Omega}=4\pi\sigma_{T}\sum_{n=0,2,\ldots}(-1)^{n/2}P_{n}(\cos\theta_{q})S_{n}(q)$
(5)

$\begin{split}S_{n}(q)=f_{1}(q)f_{2}(q)\int_{-1}^{1}d\cos(\theta)\int_{0}^{\infty}dR\int_{-\infty}^{\infty}dt\Xi(t-\tau)R^{2}\\\
|\psi(R,\cos(\theta),\tau)|^{2}P_{n}(\cos(\theta))j_{n}(qR)\end{split}$ (6)

where $\Xi(\tau)$ is the X-ray pulse intensity envelope, and
$\psi(R,\cos(\theta),\tau)$ is the molecular wavepacket. The nth order process
will be manifested both by the intensity distribution on the detector via the
scattering angle $\theta_{q}$, and via the anisotropy curves $S_{n}(q)$. A
schematic description of the relevant coordinates and experimental approach is
described in Fig 1. The subsequent analysis of the scattering pattern given by
Eq 5 will use a Legendre decomposition over the detector angle $\theta_{d}$ to
recover $S_{n}(q)$.

##  II Methods

The experimental procedure is described in detail in a previous studyGlownia
_et al._ (2016). In short, time-resolved X-ray scattering was performed at the
Linac Coherent Light Source (LCLS) free-electron laser (FEL) facility, SLAC
National Accelerator Laboratory, using the hard X-ray pump-probe (XPP)
instrument. Molecular Iodine vapor with a column density of $\sim
10^{18}cm^{-2}$ was excited by a 520 nm, 40 $\mu$J, 50 fs, optical pulses of
$\sim 5\times 10^{11}W/cm^{2}$, and probed by a 9 keV, 2 mJ, 40 fs X-ray
pulses at a variable pump-probe time delay provided by the LCLS. Approximately
$10^{7}$ X-ray photons per pulse were scattered onto a 2.3 megapixel array
detector Hart  _et al._ (2012), with $\sim 50$ photons per pulse per pixel.

Initial processing of the raw scattering data included detector corrections
that are mentioned in a previous study Glownia  _et al._ (2016). The analysis
includes a "dark" detector correction, where an average image of the detector
without incident x-rays is subtracted from the raw data, following single-
pixel detector corrections due to their non-linear response. The images are
then corrected for polarizationHura  _et al._ (2000) of the LCLS pulses, as
well as scattering geometry for the case of a plane detector Bösecke and Diat
(1997). Subsequently, the scattering patterns were corrected for absorption
artifacts due to the scattering cell geometry and upstream beam noise by
imposing an isotropy conditionVan Driel  _et al._ (2016) that is expected for
the unexcited scattering signal from Iodine thermal ground state. The images
were sorted according to a jitter correction timing tool and averaged to
obtain a pump-probe delay resolution of 20 fs. While the resolution of such
time delay binning exceeds the excitation pulse duration, it is appropriate
for the case of the multi-photon excitation observed which has a shorter
effective duration, as well as to the jitter correction resolution Harmand
_et al._ (2013); Glownia  _et al._ (2019). We improve the signal to noise
ratio by analyzing the data from several pump-probe scans that are combined
and sorted according to the pump-probe time delay.

Figure 3: (a) The LCLS 2.3 megapixel array detector (CSPAD Hart  _et al._
(2012)) showing the difference scattering signal $\Delta
I^{(\tau)}=I^{(\tau)}-\langle I^{\tau<0}\rangle$ for a pump-probe delay of
$\tau=240fs$. (b) The scattering signal is four-folded and averaged and (c)
transformed to polar detector coordinates. The usable $q$ range for analysis
was $1.2\text{\AA}^{-1}<q<4.3\text{\AA}^{-1}$ due to the scattering cell
geometry and mask used. (d) Applying Legendre decomposition to the polar
representation captures its anisotropy orders as well as filters for higher
angle dependent noise and missing signal. The reconstruction of the polar
signal here was done using Eq. 7 where we used Legendre orders up to $P_{12}$.
(e) We use the Legendre coefficients to obtain the anisotropy curves $S_{n}$
(Eq. 8) and show the estimated standard error per curve (shaded areas).

A representative time-binned difference signal is shown in Fig 3a. We use the
cylindrical symmetry of the scattering signal to increase its fidelity without
loss of angle resolution by four-folding and averaging the scattering image
quadrants. The averaging is weighted by the number of detector pixels that
contribute from each quadrant. Four-folding the detector image also helps
reduce the effective missing data that appears as gaps between the ASIC
elements of the array detector. The remaining missing data points in Fig 3b
are excluded from the later Legendre decomposition analysis. The signal is
then transformed and binned in the ($q$,$\theta_{d}$) polar coordinates.

We decompose the signal in each $q$ bin to even order Legendre basis up to the
relevant significant order (Fig 3d):

$I(q,\theta_{d})=\beta_{0}(q)\sum_{n=0,2,\ldots}\widetilde{\beta}_{n}(q)P_{n}(\cos\theta_{d})$
(7)

With the radial intensity $\beta_{0}(q)$ and the normalized detector
anisotropy terms $\widetilde{\beta}_{n}(q)=\beta_{n}(q)/\beta_{0}(q)$. The fit
is done only on $\theta_{d}$ values that contain signal. The relation between
the nth order anisotropy curve $S_{n}(q)$ in Eq. 6 and the corresponding
$\beta_{n}(q)$ term is given by

$S_{n}(q)=\frac{\beta_{0}(q)\widetilde{\beta}_{n}(q)}{\cos^{n}(\theta_{q})}=\frac{\beta_{n}(q)}{\left(1-\frac{q^{2}}{4|k_{0}|^{2}}\right)^{n/2}}\text{,
n=0,2,\ldots}$ (8)

where $|k_{0}|$ is the length of the wave vector of the incoming X-ray beam.
We note that while $\beta_{0}(q,\tau)$ has units of intensity or the average
number of scattered photons per $q$, the higher-order
$\widetilde{\beta}_{n}(q,\tau)$ are dimensionless and represent ratios between
the relevant angle components that dictate the degree of anisotropy.

Additional care is needed in centering the images before they are being four-
fold averaged and transformed to polar coordinates. The residual error to
$S_{n}$ due to centering inaccuracy was simulated for the detector binning
used in the analysis. For example, in the case of a pure isotropic signal
$S_{0}$, we find that a centering error at the single bin level $dq\sim
4\times 10^{-3}\text{\AA}^{-1}$ will create an artificial residual signal for
the higher-order $S_{n>0}$ that amounts to $\sim 10^{-4}S_{0}$ across the $q$
range used. Deviation from the correct image center beyond the single bin
resolution will create artificial anisotropy that will obscure and distort the
weak signals of the high anisotropy orders.

We find the center of the images by analyzing time delays where scattering can
only be isotropic and validate our approach by limiting its residual
anisotropy below the single-bin limit discussed. We analyze the contrast
signal:

$I_{c}=\frac{\langle I^{(\tau<0)}\rangle-\langle I^{(\tau\gg
0)}\rangle}{\langle I^{(\tau<0)}\rangle+\langle I^{(\tau\gg 0)}\rangle}$

Where $\langle I^{(\tau<0)}\rangle$ is obtained by averaging the scattering
signal before the laser excitation, containing only isotropic scattering, and
$\langle I^{(\tau\gg 0)}\rangle$ is the averaged signal at time delays after
rotational dephasing took place and most the dynamics have equilibrated. We
then apply to $I_{c}$ a series of intensity thresholds, each yields a
scattering signal $\tilde{I_{c}}$ that will be distributed symmetrically
around the center. We find the centers of these distributions using a Random
Sample Consensus (RANSAC) algorithm Torr and Zisserman (2000), an iterative
method adapted to robustly fit circles in the presence of noise.

Briefly, data points in $\tilde{I_{c}}$ are classified as outliers or inliers,
and the fitting procedure ignores the outliers. The classification is done by
randomly sampling a small subset of points to estimate the model parameters.
The condition for points to be considered inliers is given by
$|\sqrt{(x_{i}-x_{c})^{2}+(y_{i}-y_{c})^{2}}-r|<d$, where $x_{i},y_{i}$ are
data points of the subset , $x_{c},y_{c},r$ are the estimated model parameters
(center and radius), and $d$ is the distance fit tolerance. The random
sampling repeat until the fraction of inliers over the total number of points
that share the same model parameters exceeds a threshold. We then take the
trimmed mean of the estimated centers of $\tilde{I_{c}}$ as the center of the
scattering image, and obtain for the case of isotropic signal, a residual
$<10^{-4}S_{0}$ for the higher-order $S_{n>0}$, indicating centering accuracy
at the single bin level.

The values and estimated standard errors for the measured anisotropy curves
$S_{n}(q)$ are shown in Fig 3e. We calculate the standard error by first
measuring the experimental weighted sample variance of each $\Delta
I(q,\theta_{d})$ element at each time bin delay. The variance for the $k^{th}$
time delay difference signal $\Delta I^{(k)}(q,\theta_{d})$ is obtained by:

$Var(I^{(k)})=\frac{1}{N^{(k)}}\sum_{m=1}^{N^{(k)}}|w(q,\theta_{d})\big{(}\Delta
I^{(k)}(q,\theta_{d})-\Delta I^{(k)}_{m}(q,\theta_{d})\big{)}|^{2}$

where $N^{(k)}$ is the number of images recorded for that time delay bin,
$w(q,\theta_{d})$ is the statistical weight each element in $(q,\theta_{d})$
has due to the four-folding and detector pixel binning, $\Delta
I^{(k)}=I^{(k)}-\langle I^{\tau<0}\rangle$ is the weighted average intensity
difference of the $k^{th}$ time bin, and $\Delta I_{m}$ expresses the $m^{th}$
image that belongs to the $k^{th}$ time bin. The measured variance is then
translated to a weights vector in the weighted least squares Legendre fitting
process. The standard error obtained from the fit is calculated for each $q$
bin and is propagated to each anisotropy order according to Eq. 8.

Because we will be measuring a time-dependent signal $S_{n}(q,\tau)$ we would
like to use its anisotropy information as a function of the temporal pump-
probe delay to recover and disentangle the different multiphoton processes it
captures. An approach we have recently introduced to successfully characterize
one and two-photon interaction of isotropic scattering signals Bucksbaum  _et
al._ (2020); Ware  _et al._ (2019a) employs a temporal Fourier-transform to
obtain frequency-resolved X-ray-scattering signals. In this approach, time-
periodic vibration motions in $S(q,\tau)$ will appear as peaks in the
frequency domain $S(q,f)$ , while ballistic dissociation motions will appear
as lines where the dissociation velocity is linearly proportional to $v=2\pi
f/q$, where $f$ is the Fourier-transformed frequency coordinate.

The reason for such behavior can be understood if we model motion by a simple
outgoing charge density described by $\delta(R(t))=\delta(R-v\tau)$. The time
dependence in the scattering signal is then weighted by the spherical Bessel
function according to the $qR(t)$ product. The Fourier transform of such
function $\int d\tau e^{\imath\omega\tau}j_{n}(qR(\tau))$ will result in an
exponential integral function that has a maxima when $qR-\omega R/v=0$. A
detailed derivation of this result can be found in Ref Ware  _et al._ (2019b).
We shall discuss implementing such an approach in recovering the anisotropy
information for the case of multiple order contributions.

Figure 4: (top) Simulated difference signal of the anisotropy curves for the
case of dissociation with a $cos^{8}(\theta)$ angular distribution. The
highest order of anisotropy that is obtained is $S_{8}$, this term also
captures the strongest signal in the $(q,\tau)$ domain. However, the angular
distribution modeled has also projections on all lower $S_{k<n}$ orders.
(middle) Fourier transforming the temporal domain and plotting the magnitude
of each order reveals the dissociation signature as straight lines along $2\pi
f=qv$. (bottom) Analysis of how the dissociation signal of order $n$
propagates to lower orders $k<n$ as a function of $q$. In each plot, we
normalized the $k<n$ order signal by the leading order $n$ (marked dotted line
at unity). The most right-side plot shows the case for $n=8$ that is obtained
from the middle row of the figure. This type of signal propagation mapping
from higher to lower orders allows tracing the origins of each anisotropy for
the general case of multiple order contributions.

Consider an anisotropic dissociation motion modelled by the charge density:

$|\psi(R,\cos(\theta),\tau)|^{2}=\delta(R-(R_{0}+v\tau))\cos^{n}(\theta)$

where $R_{0}$ is the initial position at the instance of dissociation, $v$ is
the velocity of the dissociation of anisotropy order $n$. We use Eq. 5 and Eq.
6 to simulate the scattering pattern on a detector at each pump-probe delay
and truncate the scattering signal in the $q$-range similar to typical
experimental conditions, $1\text{\AA}^{-1}<q<4\text{\AA}^{-1}$. We subtract
the time delayed scattering signal from the stationary signal $\Delta
I=I(q,\tau)-I(q,\tau=0)$, and decompose the difference signal in Legendre
polynomials to obtain the $\widetilde{\beta}_{n}(q,\tau)$ coefficients as
shown in Eq 7. We then obtain the difference anisotropy curves $\Delta
S_{n}(q,\tau)$ using Eq. 8, and Fourier transform them to obtain the
dissociation signature in $(q,f)$ space.

For example, in Fig 4 we demonstrate the case for a pure anisotropy of order
$n=8$. We present the time-dependent anisotropy difference curves $\Delta
S_{n}(q,\tau)$ and their temporal Fourier transforms. We verify that we obtain
up to the 8th order Legendre coefficients without any higher orders
contributing. We also observe that the highest anisotropy order measured
($\Delta S_{8}$ in the example) has the strongest contribution both in the
both temporal and frequency domains. The analysis captures the expected
behavior where the dissociation velocity is obtained via $v=2\pi f/q$,
however, we also note that the same dissociation signature is present at lower
anisotropy orders.

The reason for the appearance of a similar dissociation line in the lower
orders is because of the nature of the Legendre polynomials used in the
decomposition. A scattering signal from a $cos^{n}$ distribution will not only
be captured by $P_{n}$ order polynomial but have projections to all $k<n$
order Legendre polynomials. In this example, a $cos^{8}$ distribution will be
decomposed to a Legendre series of even orders up to $n=8$. For a general case
where several orders contribute, one may choose to only analyze the $n=2$
order signal instead of all orders. This order includes besides the $cos^{2}$
contribution the projections of all higher orders, thus obtaining a
representation of the total anisotropy in the sample. Doing so, however,
hinders the assignment of specific processes with a particular degree of
anisotropy, as well as attenuates higher-order contributions by at least an
order of magnitude.

Instead, we show we can map the way the signal of each anisotropy order is
projected among lower orders to trace how different order contributions are
observed given typical experimental sampling in $q$ and $\tau$. In Fig 4 we
also show a calculation of the expected projection as function of $q$ for
orders up to $n=8$. We normalize the projections by the leading order and
obtain typical ranges of the magnitude of signal propagation, as summarized in
Table 1. We will use the information of that mapping to disentangle
dissociation processes among the orders and uncover the contribution of each
order.

Table 1:  The ratio ranges of $\mathcal{P}|S_{k}(q,f)|/|S_{n}(q,f)|$, where $\mathcal{P}|S_{k}|$ is the $k^{th}$ order projection of $|S_{n}(q,f)|$ with $k=0,2,..n-2$. The numbers in each cell present the value limits obtained for the range $1\text{\AA}^{-1}<q<4\text{\AA}^{-1}$. The values in the table are obtained from the curves seen at the bottom row of Fig 4 .  | $\mathcal{P}S_{0}$ | $\mathcal{P}S_{2}$ | $\mathcal{P}S_{4}$ | $\mathcal{P}S_{6}$
---|---|---|---|---
$S_{2}$ | 0.65 - 2.98 | 1 | - | -
$S_{4}$ | 0.68 - 2.62 | 0.15 - 0.46 | 1 | -
$S_{6}$ | 0.49 - 0.96 | 0.09 - 0.25 | 0.13 - 0.9 | 1
$S_{8}$ | 0.36 - 1.07 | 0.09 - 0.15 | 0.08 - 0.34 | 0.13 - 1.23

##  III Results and discussion

Figure 5: The experimental anisotropy curves $\Delta S_{n}$ obtained from
Legendre decomposition of the measured scattering signal. We observe
anisotropy up to $n=8$ order. The $n=0$ term is the angle averaged scattering
signal that captures all dynamics. Time-dependent oscillations of several
frequencies across the entire $q$ range can be seen (most noticeable at
$2<q<3$), caused by various single and multiphoton vibration excitations.
Higher $\Delta S_{n}$ terms capture various dissociation processes that are
manifested by $q$-dependent modulations curves, the longer the dissociation
distance the faster the modulation in $q$. The $q$-dependence nature of the
amplitude in the various orders is due to the spherical Bessel $j_{n}(qR)$
dependence as seen in Eq. 6. The drop of signal in the higher-order terms as a
function of delay is due to rotational dephasing for the case of the molecules
that are excited to bound states. However, excitations that lead to
dissociation carry the anisotropic nature of the signal regardless of
dephasing as the atoms produced in the dissociation process will move outward
along the straight line defined by the internuclear axis of the molecule.

We have applied the Legendre decomposition as described in Eq. 7 to the
measured scattering difference signal $\Delta
I^{(\tau)}=I^{(\tau)}(q,\theta_{d})-\langle I^{(\tau<0)}(q,\theta_{d})\rangle$
to derive the experimental anisotropy curves $\Delta S_{n}(q,\tau)$ using Eq.
8. The result is seen in Fig 5 for the different anisotropy orders. The
decomposition was done up to order $n=12$ and we have found that for the
experimental conditions used it was sufficient to only consider orders of up
to $n=8$, indicating that multiphoton processes up to 4-photon absorption are
taking place. The signal at $\Delta S_{10}(q,\tau)$ has a negligible
contribution for the analysis, as its observable signal at $q\sim
1.2\text{\AA}^{-1}$ does not allow information to be extracted, and may be
related to the limit of anisotropy detection at this range due to the
asymmetry of the inner hole mask of the detector.

The entire excitation dynamics are captured in the angle averaged $\Delta
S_{0}(q,\tau)$ term, where we primarily observe two types of signals: Time
periodic oscillations across the entire $q$ range related to various
vibrational excitations, and, $q$-dependent modulations that change their
radius of curvature as a function of time and indicate dissociative motion.
The dissociation signature is also evident in the higher $\Delta
S_{n}(q,\tau)$ terms, as this type of signal is not sensitive to rotational
dephasing where fragmentation preserves the molecules axis.

Figure 6: The (top) magnitude and (bottom) real part of the measured Fourier
transformed anisotropy curves $\Delta S_{n}(q,f)$. For visualization we use a
logarithmic scale for the magnitude. The $|\Delta S_{0}(q,f)|$ term captures
all dynamics, including many multiphoton vibrationaly excited states that
appear as peaks along $f$. These peaks broaden and attenuate in higher orders
due to rotational dephasing that limits their time periodic sampling. Several
dissociation pathways that appear as diagonal lines are captured in the
$S_{2}$ and $S_{4}$ orders, as well as weaker dissociation in the $S_{6}$ and
$S_{8}$ orders. To analyze the dissociation signatures we also use the
information of the real or imaginary part of $S_{n}(q,f)$ that allows better
contrast vs the noise floor, as seen for the terms $\Re{S_{n}(q,f)}$ scaled by
the Fourier vector that appear below the magnitude information. (bottom left)
Analyzing the positions of the dissociation signatures in different anisotropy
orders allows to identify the pathways and the number of photons that
participated (see text). The line colors code for the different anisotropy
orders, the slopes correspond to the average dissociation velocities: (a) 7.9
(b) 16.1 (c) 17.5 (d) 19.1 (e) 21.6 (f) 15.7 (g) 19.8 (h) 21.6 (i) 15.5 (j)
20.4 (k) 19.6 $\text{\AA}/ps\pm 0.5\text{\AA}/ps$.

We then Fourier transform the anisotropy curves to obtain the frequency-
resolved $\Delta S_{n}(q,f)$ terms as seen in Fig 6. We analyzed the magnitude
and real part of each term and obtained the estimated dissociation velocities
for the different anisotropy orders via a linear model estimation using a
RANSAC algorithm. For vibration frequencies detected, the estimate was based
on the center of mass of the peak found in frequency. We observe that the
dissociation velocities take place at a range of $8-22\text{\AA}/ps$. We
analyze the signal strengths of each dissociation line starting from the
highest order measured and consider the way the signal can propagate to lower
orders to resolve its origin and the pathway that takes place. We also observe
a multitude of vibration excitation that is manifested as peaks in frequency
at a wide $q$ range, most noticeably in the $|\Delta S_{0}(q,f)|$ term.
Analysis of the vibrational excitation is limited in the high order anisotropy
terms due to rotational dephasing that restricts the effective time-window to
sample periodic motions. As a result, for Iodine, slower vibrations will only
have two to three periods before the signal will be lost to dephasing. This
broadens the detected peak widths and restricts resolving details in
frequency.

To consider the way the signal of higher-order anisotropy terms is projected
to lower orders, we start by analyzing the highest observable order. The
observed dissociation signal of $|\Delta S_{8}(q,f)|$ indicates that a
4-photon transition to a repulsive _gerade_ excited state occurred. From the
estimated dissociation velocity of $19.6\pm 0.5\text{\AA}/ps$ we calculate a
kinetic energy release (KER) of $2.5\pm 0.1eV$ that corresponds to a
transition to the $C1_{g}$ state that dissociates to the second dissociation
limit: $I(^{3}P_{3/2})+I(^{3}P_{1/2})$. This transition is accomplished via a
3-photon absorption to Rydberg states following a 1-photon transition back to
the $C1_{g}$ state, as well as a Raman excitation to the ground state
following 2-photon absorption. Furthermore, 4-photon absorption from the
ground state can populate the 6d-Rydberg series which lies above the
dissociation limit of the first-tier IP states. While these states are known
to give rise to several fragmentation channels, these processes take place at
a much later time compared to the time delays probed here. Previous studies
that measured the fragmentation of these Rydberg states Donovan, Lawley, and
Ridley (2015); Bogomolov  _et al._ (2014) show that there was no anisotropy
measured for the fragments. This result indicates that dissociation happened
at a timescale longer than the rotation period of Iodine, at least an order of
magnitude larger than the range of time delays probed here. Around the range
$1.5\text{\AA}^{-1}<q<2.5\text{\AA}^{-1}$ where the dissociation signal was
observed, we deduce that $\sim 30\%$ of the signal in $|\Delta S_{8}(q,f)|$ is
projected to the $|\Delta S_{6}(q,f)|$. Similarly, the projections to $|\Delta
S_{4}(q,f)|$ and $|\Delta S_{2}(q,f)|$ are about $20\%$ and $10\%$
respectively.

The observed dissociation signal of the $|\Delta S_{6}(q,f)|$ term is related
to a 3-photon transition to an _ungerade_ repulsive excited state, as well as
to a possible projection from $n=8$ order. We measure two dissociation
signatures for this term, with estimated velocities of $15.5\pm
0.5\text{\AA}/ps$ and $19.1\pm 0.5\text{\AA}/ps$. The slower dissociation
velocity agrees with the expected KER of $\sim 1.6eV$ for nonlinear-Raman
transition to the $C(B^{{}^{\prime\prime}})1_{u}$ state that decays to the
first dissociation limit $I(^{3}P_{3/2})+I(^{3}P_{3/2})$. The second
dissociation signal with $19.1\text{\AA}/ps$ is interpreted as a projection of
the $n=8$ order of similar velocity. To arrive at that conclusion we sampled a
region of interest around $q\sim 2\text{\AA}^{-1}$ where the dissociation
signal is present and calculated the ratio of the median intensity
$\bar{|S_{6}|}/\bar{|S_{8}|}$ for this range. We obtain a ratio of $0.42$,
which is similar to the expected range of values for a projection from $n=8$
in that $q$ range. Considering the noise floor observed for these terms in Fig
6, we deduce that most, if not all the signal observed at that dissociation
channel is actually due to a projection of the $n=8$ order. We also observe
peaks around $1.4\pm 0.3$ and $2.1\pm 0.3$ THz corresponding to periods of 710
and 470 fs respectively. The first peak can be interpreted as a 3 photon
absorption process to high lying vibration levels $v=278-290$ of the
$D0^{+}_{u}$ IP state, as well as for the $v=136$ range of the $F0^{+}_{u}$ IP
state Hoy and Lipson (1990); Hiraya  _et al._ (1988). However, the limited
resolution of the approach used cannot resolve finer details. The second
detected peak is interpreted as a vibration wavepacket that was excited via
nonlinear Raman scattering to the bound $B(0^{+}_{u})$ state with an estimated
period of 450 fs.

For the $|\Delta S_{4}(q,f)|$ term we observe three dissociation channels with
velocities of $15.7\pm 0.5\text{\AA}/ps$, $19.8\pm 0.5\text{\AA}/ps$, and
$21.6\pm 0.5\text{\AA}/ps$. This order captures 2-photon transitions to a
_gerade_ excited state along with projections of the higher orders. Analyzing
the relative signal strengths we observe that high order projections are minor
here. The measured dissociation of $15.7\text{\AA}/ps$ correspond to a
2-photon transition to the $O^{+}_{g}$ state that dissociate via the 3rd limit
$I(^{3}P_{1/2})+I(^{3}P_{1/2})$ with a KER of $1.6eV$. The next dissociation
is more prominent and relates to excitation of the $C1_{g}$ state via the 2nd
limit $I(^{3}P_{3/2})+I(^{3}P_{1/2})$ with a KER of $2.6eV$. Further, we
measure another dissociation signal with a faster velocity of
$21.6\text{\AA}/ps$ that corresponds to a KER of $3.1eV$ indicating that we
excite the $a1_{g}$ state that decay to the 1st limit
$I(^{3}P_{3/2})+I(^{3}P_{3/2})$. We observe peaks in frequency at $1.75\pm
0.2$ THz, $2.3\pm 0.2$ THz corresponding to periods of 570 and 430 fs
respectively. We also observe a weaker peak around $6\pm 0.2$ THz which
indicates Raman excitation of the ground state at $\sim 167$ fs period. The
lower peaks around $2$ THz can only be explained by higher-order projections
that are more visible in lower orders, as will be elaborated for the $n=2$
order.

The lowest anisotropy order $n=2$ contains multiple dissociation channels that
coincide with the higher orders, as also shown in Fig 6 summary plot. In the
perturbative limit this order will capture 1-photon transitions to _ungerade_
excited states such as the bound $B(0^{+}_{u})$ state, the repulsive
$C(B^{\prime\prime})1_{u}$ and $A1_{u}$ states. Analyzing the signal strength
of the different channels concerning possible projections of the higher orders
we learn that the $n=4$ order contributes most of the signal measured along
the same dissociation channels with two exceptions. We measure a dissociation
line at the $16.1\text{\AA}/ps$ that correspond to 1-photon transition to the
$C(B^{{}^{\prime\prime}})1_{u}$ state. Additionally, the slowest dissociation
velocity we measure $7.9\text{\AA}/ps$ is related to a transition to the
lowest repulsive $A1_{u}$ state with a KER of 0.41 eV. Such slow dissociation
is at the limit of detection using the frequency-resolved approach as the line
that corresponds to it has a greater overlap with the different vibration
peaks, as seen around 1.5-2 THz.

The $n=2$ order also shows multiple strong peaks in frequency up to 6 THz,
that relate to vibration excitation, besides the strong and broad peak at 1.9
THz that correspond to the vibration motion of the $B(0^{+}_{u})$. We see this
effect also for the $n=4$ order lower frequency peaks. The ability to resolve
vibration modes in frequency depends on the way rotation motion is projected
to the different anisotropy orders. As a result, projections of higher to
lower orders for periodic motions depend on the rotation dynamics in addition
to the mathematical projection. This contrasts with to the way dissociation
signal projection propagates, as it is insensitive to rotational dephasing.
Consequently, the $S_{2}$ term effectively samples more time delay bins
compared to higher orders when considering periodic motions. The attenuation
and broadening of periodic signals due to rotational dephasing in the higher
orders restrict the projection analysis between the orders to the point it is
impractical to assign excitation pathways based on anisotropy.

The isotropic signal $|\Delta S_{0}(q,f)|$ combines all types of signals, and
in particular, can be used to identify all peaks in frequency that corresponds
to several Raman and nonlinear Raman excitations. While it can be challenging
to resolve some of these peaks because of the overlap they have with all the
diagonal dissociation lines, we conclude that the spectroscopy of vibrational
motion is mostly limited to the isotropic signal in the context of the
frequency-resolved X-ray scattering analysis presented here.

##  IV Conclusions

In conclusion, we are presenting the first results of the observation of high-
order anisotropy, up to 8th order, in ultrafast x-ray scattering. We present
an analysis method that uses a Legendre decomposition and Fourier analysis on
measured and simulated scattering signal. We map the way different anisotropy
orders project signal to lower orders and use it to resolve and interpret the
experimentally observed signals. We resolve many dissociation channels and
vibration modes that are excited by multiphoton transitions and assign the
processes based on the order of the anisotropy they appeared at. Leveraging
high-order anisotropy information of the scattering signal is limited by the
rotational timescale. In particular, it is limiting the ability to use it for
Fourier-resolved x-ray scattering of vibration motions, implying that direct
real-space approaches should be considered in this context. However, during
that rotational timescale, the anisotropy information provides a sensitive
filter to differentiate and trace different excitation pathways that take
place simultaneously. The information provided using this approach can be also
applied in other analysis schemes, as well as open the possibility to serve as
a prior when analyzing the isotropic time-dependent signal past the rotation
timescale.

## Conflicts of interest

There are no conflicts to declare.

## Acknowledgments

This work was supported by the U.S. Department of Energy, Office of Science,
Basic Energy Sciences, Chemical Sciences, Geosciences, and Biosciences
Division. The experiment described was carried out at the Linac Coherent Light
Source (LCLS) at the SLAC National Accelerator Laboratory. LCLS is an Office
of Science User Facility operated for the U.S. Department of Energy Office of
Science by Stanford University.

## References

  * Ihee  _et al._ (2010) H. Ihee, M. Wulff, J. Kim, and S.-i. Adachi, International Reviews in Physical Chemistry 29, 453 (2010).
  * Minitti  _et al._ (2015) M. Minitti, J. Budarz, A. Kirrander, J. Robinson, D. Ratner, T. Lane, D. Zhu, J. Glownia, M. Kozina, H. Lemke, _et al._ , Physical review letters 114, 255501 (2015).
  * Stankus  _et al._ (2019) B. Stankus, H. Yong, N. Zotev, J. M. Ruddock, D. Bellshaw, T. J. Lane, M. Liang, S. Boutet, S. Carbajo, J. S. Robinson, _et al._ , Nature chemistry 11, 716 (2019).
  * Biasin  _et al._ (2016) E. Biasin, T. B. Van Driel, K. S. Kjær, A. O. Dohn, M. Christensen, T. Harlang, P. Vester, P. Chabera, Y. Liu, J. Uhlig, _et al._ , Physical review letters 117, 013002 (2016).
  * Van Driel  _et al._ (2016) T. B. Van Driel, K. S. Kjær, R. W. Hartsock, A. O. Dohn, T. Harlang, M. Chollet, M. Christensen, W. Gawelda, N. E. Henriksen, J. G. Kim, _et al._ , Nature communications 7, 1 (2016).
  * Chollet  _et al._ (2015) M. Chollet, R. Alonso-Mori, M. Cammarata, D. Damiani, J. Defever, J. T. Delor, Y. Feng, J. M. Glownia, J. B. Langton, S. Nelson, and et al., Journal of Synchrotron Radiation 22, 503–507 (2015).
  * Glownia  _et al._ (2016) J. Glownia, A. Natan, J. Cryan, R. Hartsock, M. Kozina, M. Minitti, S. Nelson, J. Robinson, T. Sato, T. van Driel, _et al._ , Physical review letters 117, 153003 (2016).
  * Biasin  _et al._ (2018) E. Biasin, T. B. Van Driel, G. Levi, M. G. Laursen, A. O. Dohn, A. Moltke, P. Vester, F. B. Hansen, K. S. Kjaer, T. Harlang, _et al._ , Journal of synchrotron radiation 25, 306 (2018).
  * Haldrup  _et al._ (2019) K. Haldrup, G. Levi, E. Biasin, P. Vester, M. G. Laursen, F. Beyer, K. S. Kjær, T. B. Van Driel, T. Harlang, A. O. Dohn, _et al._ , Physical review letters 122, 063001 (2019).
  * Baskin and Zewail (2006) J. S. Baskin and A. H. Zewail, ChemPhysChem 7, 1562 (2006).
  * Yang  _et al._ (2018) J. Yang, X. Zhu, T. J. Wolf, Z. Li, J. P. F. Nunes, R. Coffee, J. P. Cryan, M. Gühr, K. Hegazy, T. F. Heinz, _et al._ , Science 361, 64 (2018).
  * Bucksbaum  _et al._ (1990) P. H. Bucksbaum, A. Zavriyev, H. G. Muller, and D. W. Schumacher, Physical review letters 64, 1883 (1990).
  * McKenna  _et al._ (2008) J. McKenna, A. Sayler, F. Anis, B. Gaire, N. G. Johnson, E. Parke, J. Hua, H. Mashiko, C. Nakamura, E. Moon, _et al._ , Physical review letters 100, 133001 (2008).
  * Natan  _et al._ (2012) A. Natan, U. Lev, V. S. Prabhudesai, B. D. Bruner, D. Strasser, D. Schwalm, I. Ben-Itzhak, O. Heber, D. Zajfman, and Y. Silberberg, Physical Review A 86, 043418 (2012).
  * Natan  _et al._ (2016) A. Natan, M. R. Ware, V. S. Prabhudesai, U. Lev, B. D. Bruner, O. Heber, and P. H. Bucksbaum, Physical Review Letters 116, 143004 (2016).
  * Rosca-Pruna and Vrakking (2001) F. Rosca-Pruna and M. Vrakking, Physical Review Letters 87, 153902 (2001).
  * Bisgaard  _et al._ (2009) C. Z. Bisgaard, O. J. Clarkin, G. Wu, A. M. Lee, O. Geßner, C. C. Hayden, and A. Stolow, Science 323, 1464 (2009).
  * Makhija  _et al._ (2020) V. Makhija, K. Veyrinas, A. E. Boguslavskiy, R. Forbes, I. Wilkinson, R. Lausten, S. P. Neville, S. T. Pratt, M. S. Schuurman, and A. Stolow, Journal of Physics B: Atomic, Molecular and Optical Physics 53, 114001 (2020).
  * Karamatskos  _et al._ (2019) E. T. Karamatskos, S. Raabe, T. Mullins, A. Trabattoni, P. Stammer, G. Goldsztejn, R. R. Johansen, K. Długołecki, H. Stapelfeldt, M. J. Vrakking, _et al._ , Nature communications 10, 1 (2019).
  * Larsen, Wendt-Larsen, and Stapelfeldt (1999) J. J. Larsen, I. Wendt-Larsen, and H. Stapelfeldt, Physical review letters 83, 1123 (1999).
  * Cryan  _et al._ (2010) J. P. Cryan, J. M. Glownia, J. Andreasson, A. Belkacem, N. Berrah, C. I. Blaga, C. Bostedt, J. Bozek, C. Buth, L. F. DiMauro, _et al._ , Physical review letters 105, 083004 (2010).
  * McFarland  _et al._ (2008) B. K. McFarland, J. P. Farrell, P. H. Bucksbaum, and M. Gühr, Science 322, 1232 (2008).
  * Yang  _et al._ (2015) J. Yang, J. Beck, C. J. Uiterwaal, and M. Centurion, Nature communications 6, 1 (2015).
  * Yang  _et al._ (2016) J. Yang, M. Guehr, T. Vecchione, M. S. Robinson, R. Li, N. Hartmann, X. Shen, R. Coffee, J. Corbett, A. Fry, _et al._ , Nature communications 7, 11232 (2016).
  * Kierspel  _et al._ (2020) T. Kierspel, A. Morgan, J. Wiese, T. Mullins, A. Aquila, A. Barty, R. Bean, R. Boll, S. Boutet, P. Bucksbaum, _et al._ , The Journal of Chemical Physics 152, 084307 (2020).
  * Xiong, Wilkin, and Centurion (2020) Y. Xiong, K. J. Wilkin, and M. Centurion, Physical Review Research 2, 043064 (2020).
  * Kierspel  _et al._ (2015) T. Kierspel, J. Wiese, T. Mullins, J. Robinson, A. Aquila, A. Barty, R. Bean, R. Boll, S. Boutet, P. Bucksbaum, _et al._ , Journal of Physics B: Atomic, Molecular and Optical Physics 48, 204002 (2015).
  * Brown  _et al._ (1999) E. J. Brown, I. Pastirk, B. I. Grimberg, V. V. Lozovoy, and M. Dantus, The Journal of chemical physics 111, 3779 (1999).
  * Pastirk, Lozovoy, and Dantus (2001) I. Pastirk, V. V. Lozovoy, and M. Dantus, Chemical Physics Letters 333, 76 (2001).
  * Fang and Gibson (2008) L. Fang and G. Gibson, Physical review letters 100, 103003 (2008).
  * Bogomolov  _et al._ (2014) A. S. Bogomolov, B. Grüner, S. A. Kochubei, M. Mudrich, and A. V. Baklanov, The Journal of chemical physics 140, 124311 (2014).
  * Lukashov, Petrov, and Pravilov (2018) S. Lukashov, A. Petrov, and A. Pravilov, Cham: Springer International Publishing (2018).
  * Ben-Nun, Cao, and Wilson (1997) M. Ben-Nun, J. Cao, and K. R. Wilson, The Journal of Physical Chemistry A 101, 8743–8761 (1997).
  * Lorenz, Møller, and Henriksen (2010a) U. Lorenz, K. B. Møller, and N. E. Henriksen, Phys. Rev. A 81, 023422 (2010a).
  * Lorenz, Møller, and Henriksen (2010b) U. Lorenz, K. B. Møller, and N. E. Henriksen, New Journal of Physics 12, 113022 (2010b).
  * Simmermacher  _et al._ (2019) M. Simmermacher, N. E. Henriksen, K. B. Møller, A. M. Carrascosa, and A. Kirrander, Physical review letters 122, 073003 (2019).
  * Hart  _et al._ (2012) P. Hart, S. Boutet, G. Carini, M. Dubrovin, B. Duda, D. Fritz, G. Haller, R. Herbst, S. Herrmann, C. Kenney, _et al._ , in  _X-Ray free-electron lasers: beam diagnostics, beamline instrumentation, and applications_ , Vol. 8504 (International Society for Optics and Photonics, 2012) p. 85040C.
  * Hura  _et al._ (2000) G. Hura, J. M. Sorenson, R. M. Glaeser, and T. Head-Gordon, The Journal of Chemical Physics 113, 9140 (2000).
  * Bösecke and Diat (1997) P. Bösecke and O. Diat, Journal of Applied Crystallography 30, 867 (1997).
  * Harmand  _et al._ (2013) M. Harmand, R. Coffee, M. R. Bionta, M. Chollet, D. French, D. Zhu, D. Fritz, H. Lemke, N. Medvedev, B. Ziaja, _et al._ , Nature Photonics 7, 215 (2013).
  * Glownia  _et al._ (2019) J. M. Glownia, K. Gumerlock, H. T. Lemke, T. Sato, D. Zhu, and M. Chollet, Journal of synchrotron radiation 26, 685 (2019).
  * Torr and Zisserman (2000) P. H. Torr and A. Zisserman, Computer vision and image understanding 78, 138 (2000).
  * Bucksbaum  _et al._ (2020) P. H. Bucksbaum, M. R. Ware, A. Natan, J. P. Cryan, and J. M. Glownia, Physical Review X 10, 011065 (2020).
  * Ware  _et al._ (2019a) M. R. Ware, J. M. Glownia, A. Natan, J. P. Cryan, and P. H. Bucksbaum, Philosophical Transactions of the Royal Society A 377, 20170477 (2019a).
  * Ware  _et al._ (2019b) M. R. Ware, J. M. Glownia, N. Al-Sayyad, J. T. O’Neal, and P. H. Bucksbaum, Physical Review A 100, 033413 (2019b).
  * Donovan, Lawley, and Ridley (2015) R. J. Donovan, K. P. Lawley, and T. Ridley, The Journal of Chemical Physics 142, 204306 (2015).
  * Hoy and Lipson (1990) A. Hoy and R. Lipson, Chemical physics 140, 187 (1990).
  * Hiraya  _et al._ (1988) A. Hiraya, K. Shobatake, R. Donovan, and A. Hopkirk, The Journal of chemical physics 88, 52 (1988).


# Dynamical Stability of Slip-stacking Particles

Jeffrey Eldred  Department of Physics, Indiana University, Bloomington,
Indiana 47405, USA  Fermi National Accelerator Laboratory, Batavia, Illinois
60510, USA  Robert Zwaska  Fermi National Accelerator Laboratory, Batavia,
Illinois 60510, USA

###### Abstract

We study the stability of particles in slip-stacking configuration, used to
nearly double proton beam intensity at Fermilab. We introduce universal area
factors to calculate the available phase space area for any set of beam
parameters without individual simulation. We find perturbative solutions for
stable particle trajectories. We establish Booster beam quality requirements
to achieve 97% slip-stacking efficiency. We show that slip-stacking dynamics
directly correspond to the driven pendulum and to the system of two standing-
wave traps moving with respect to each other.

###### pacs:

29.20.dk, 02.60.Cb, 37.10.Jk, 43.25.Uv

## Introduction

Slip-stacking is integral to high-intensity operation at Fermilab and will
likely play a central role in upgrades to the accelerator complex PIP Mariani
Galambos . Particle loss in the slip-stacking process is a limiting factor on
ultimate performance PIP  Brown . Single-particle dynamics associated with
slip-stacking contribute directly to the particle losses. This paper analyzes
these dynamics at depth, both analytically and numerically. Our numerical
results completely characterize the stable phase-space boundary. We use these
results to recommend an upgrade to the Fermilab Booster that would
substantially reduce slip-stacking losses.

Our analytical results provide insight into slip-stacking by presenting a
perturbative general solution and new parameteric resonances. These results
should also be of interest to the greater field of dynamical mathematics
because, as we demonstrate, the dynamics of slip-stacking are isomorphic to
the well-studied dynamics of the driven pendulum. The analysis in this paper
is also intended to facilitate application of slip-stacking to other
accelerators and non-accelerator systems with analogous dynamics.

## Background

Slip-stacking is a particle accelerator configuration that permits two high-
energy particle beams of different momenta to use the same transverse space in
a cyclic accelerator. The two beams are longitudinally focused by two sets of
rf cavities with a small frequency difference between them. Each frequency is
tuned to the momentum of one of the beams.

The two azimuthal beam distributions are manipulated as a consequence of their
difference in rf frequency. The two beams injected on separated portions of
azimuth with a small frequency difference will overlap gradually, allowing
injection Brown . When the cyclic accelerator is filled, the azimuthal
distribution of the beams will coincide at a certain tune and can then be
accelerated simultaneously. The accelerating rf cavities operate at the
average frequency, capturing both beams as one. The potential beam intensity
of a synchrotron is doubled through the application of this technique.

A preliminary study explored the beam dynamics in a 2-rf system Mills . The
slipping of bunched beams was first demonstrated at the CERN SPS Boussard  but
the emittance growth led to unacceptable particle losses. Fermilab has
subsequently implemented slip-stacking operationally since 2004 MacLachlan
SeiyaBC Brown . Initially, the higher beam intensity was used to increase
antiproton production for proton-antiproton collider experiments Cd .
Subsequently, slip-stacking was applied to neutrino production for Neutrinos
at Main Injector (NuMI) experiments Minos Minerva1 Minerva2 Nova .

Beam-loading effects can impact the effectveness of slip-stacking and were
addressed in the Main Injector by the development of a beam-loading
compensation system with -14dB feedback and -20dB feedforward SeiyaB DeyKK Dey
. The beam-loading effects on slip-stacking in the Recycler will be an order
of magnitude weaker than in the Main Injector and can be compensated if
necessary. The typical beam-loading voltage is $\sim$2kV DeyKK  compared to a
typical rf cavity voltage of 90kVDey . In the Main Injector the $R_{sh}/Q$ of
the rf cavities is $100\Omega$ DeyKK , while in the Recycler the $R_{sh}/Q$ is
$13\Omega$ Madrak . This paper focuses on the constraints on the stable phase-
space area from the single-particle dynamics of the two-rf system; direct
space-charge effects are an order of magnitude weaker.

For the particular case of slip-stacking at Fermilab, the difference between
the two rf frequencies must be equal to the product of the harmonic number of
the Booster rf and the cycle rate of the Fermilab Booster. The cycle-rate of
the Booster is 15-Hz and therefore $\Delta f=h_{B}f_{B}=1260$ Hz. A possible
upgrade to a 20-Hz Booster is also analyzed, for which $\Delta f=1680$ Hz. A
20-Hz Booster would enable slip-stacking buckets with substantially greater
phase-space area. However, the slip-stacking cyclic accelerator (either the
Main Injector or the Recycler) must be able to simultaneously accommodate
beams in a range of momentum corresponding to their frequency difference.

## Single-rf Longitudinal Stability

The motion of a single particle under the influence of a single stationary rf
cavity is described in terms of its phase space coordinates $\phi$ and
$\delta$. The phase $\phi$ is the phase of the particle relative to the
resonating electromagnetic field in the rf cavity. $\delta$ is the fractional
deviation from the reference momentum; $\delta=0$ corresponds to a particle
whose revolution frequency $f_{rev}$ is a subharmonic of the frequency of the
rf cavity $f_{rf}=hf_{rev}$. The phase-slip factor $\eta$ is used to describe
how the revolution period $T$ changes with $\delta$ and is given by
$\displaystyle\eta\delta=\Delta T/T$ (see SYLee ).

The equations of motion associated with the trajectory of a single particle
under the influence of a single stationary rf cavity SYLee  are given by:

$\dot{\phi}=2\pi
f_{rev}h\eta\delta,~{}\dot{\delta}=f_{rev}\frac{eV}{\beta^{2}E}\sin(\phi).$
(1)

$V$ is the effective voltage of the rf cavity, $e$ is the charge of the
particle, $\beta=v/c$ is the velocity fraction of the speed of light, $E$ is
the total energy of the particle. Let $V_{\delta}$ be equal to
$\frac{eV}{\beta^{2}E}$, the maximum change in $\delta$ during a single
revolution.

For small $\phi$, Eq. 1 has a stable solution known as a synchrotron
oscillation

$\phi=\rho\sin(\omega_{s}t+\psi)$ (2)

where the synchrotron frequency $\omega_{s}$ is given by
$\displaystyle\omega_{s}=2\pi f_{rev}\sqrt{\frac{V_{\delta}h\eta}{2\pi}}$. The
amplitude $\rho$ and the initial phase $\psi$ are set by initial conditions.
More generally, stable oscillatory motion is bound in phase and momentum by
the separatrix

$\delta=\pm\frac{2}{h|\eta|}\frac{\omega_{s}}{\omega_{rev}}\cos\left(\frac{\phi}{2}\right).$
(3)

The equations of motion given in Eq. 1 are isomorphic to that of a simple
pendulum. The stable region of phase space within the separatrix is referred
to as the rf bucket and the region outside of the separatrix is referred to as
slipping with respect to the bucket.

In contrast, the dynamics of slip-stacking are explicitly time-dependent and
there is no simple separatrix delineating the bucket boundary. The lack of a
clearly defined bucket confounds beam operation as the incoming particles
cannot be conventionally inserted into the bucket. We broaden the term bucket
to include cases without a separatrix: a particle trajectory is in a
particular rf bucket if the particle phase with respect to the rf cavity is
bounded and averages to zero.

## Slip-stacking and the Driven Pendulum

The equations of motion for a single particle under the influence of two rf
cavities with identical voltage and different frequencies are:

$\displaystyle\dot{\phi_{A}}$ $\displaystyle=2\pi f_{rev}h\eta\delta_{A}$
$\displaystyle\dot{\delta_{A}}$
$\displaystyle=2f_{rev}V_{\delta}\sin(\phi_{A})\cos\left(\frac{\omega_{\phi}t+\phi_{D}}{2}\right).$
(4)

$\omega_{\phi}$ is what we refer to as the phase-slipping frequency, the
angular frequency separation between the two rf cavities
$\omega_{\phi}=2\pi\Delta f$. $\phi_{A}$ is the average of the phases and
$\phi_{D}$ is the difference between the phases for the two rf cavities at
$t=0$. Without loss of generality, we eliminate $\phi_{D}$ with $t\rightarrow
t-\phi_{D}/\omega_{\phi}$. Applying a substitution to Eq. 4 yields the
corresponding second-order equation of motion:

$\ddot{\phi_{A}}=-2\omega_{s}^{2}\sin(\phi_{A})\cos\left(\frac{\omega_{\phi}t}{2}\right).$
(5)

The corresponding Hamiltonian is then given by:

$H=\pi
f_{rev}h\eta\delta_{A}^{2}+f_{rev}V_{\delta}\cos(\phi_{A})\cos\left(\frac{\omega_{\phi}t}{2}\right).$
(6)

This Hamiltonian leads to nonlinear, resonant, and chaotic phase-space
trajectories. We find that this Hamiltonian is isomorphic to that of a
pendulum under a sinusoidal driving force (in the absence of gravity). The
driven pendulum is a type of nonlinear Mathieu equation that is a subject of
ongoing research in computational mathematics Broer Xu . A canonical form for
the driven pendulum may be parameterized as Broer :

$\frac{d^{2}}{d\tau^{2}}x=(a+b\cos(\tau))\sin(x),~{}a\geq 0,~{}b\geq 0.$ (7)

Eq. 7 is obtained from Eq. 5 under the substitution

$t\rightarrow(2/\omega_{\phi})\tau,\quad\phi_{A}\rightarrow x+\pi,\quad
b=8\left(\frac{\omega_{s}}{\omega_{\phi}}\right)^{2}.$ (8)

The parameter $a$ corresponds to the force of gravity and $a=0$ in this
analogue.

The accelerator literature Mills Boussard MacLachlan  has identified the
importance of the slip-stacking parameter

$\alpha_{s}=\omega_{\phi}/\omega_{s}$ (9)

as the criterion for effective slip-stacking. The parameter $b$ in Eq. 8 is a
function of $\alpha_{s}$, indicating that all nontrivial dynamics of slip-
stacking depend only on $\alpha_{s}$. For example, if one slip-stacking
configuration has phase-slipping frequency $\omega_{\phi}$ and another
configuration with the same $\alpha_{s}$ has phase-slipping frequency
$\omega_{\phi}^{\prime}$ then the second phase space diagram is isomorphic to
the first where the $\delta$ axis must be scaled by
$\omega_{\phi}^{\prime}/\omega_{\phi}$.

## General Perturbative Solution

We analyze the dynamics within a slip-stacking bucket by shifting the origin
from one synchronized to the average frequency and phase of the two rf
cavities to one synchronized to the rf cavity with higher frequency:

$\dot{\phi}=2\pi
f_{rev}h\eta\delta,~{}\dot{\delta}=f_{rev}V_{\delta}[\sin(\phi)+\sin(\phi-\omega_{\phi}t)].$
(10)

Substituting one equation into the other we have:

$\ddot{\phi}=-\omega_{s}^{2}[\sin(\phi)+\sin(\phi)\cos(\omega_{\phi}t)-\cos(\phi)\sin(\omega_{\phi}t)].$
(11)

Eq. 11 can be expanded into powers of $\phi$ to consider the small $\phi$
perturbation:

$\displaystyle\ddot{\phi}=-\omega_{s}^{2}\Bigg{\\{}$
$\displaystyle\sum_{k=0}^{\infty}\frac{(-1)^{k}}{(2k+1)!}\phi^{2k+1}[1+\cos(\omega_{\phi}t)]$
$\displaystyle-\sum_{k=0}^{\infty}\frac{(-1)^{k}}{(2k)!}\phi^{2k}\sin(\omega_{\phi}t)\Bigg{\\}}.$
(12)

We use the Poincare-Lindstedt method (see Ch. 2 of Lichtenberg ) to find the
perturbative solution to Eq. 12 as a linear combination of oscillatory terms.
We substitute the case $\phi=0$ into Eq. 12 and we solve to generate
$\phi=-\alpha_{s}^{2}\sin(\omega_{\phi}t)$. Next we use
$\phi=-\alpha_{s}^{2}\sin(\omega_{\phi}t)$ to generate
$A_{n}\sin(n\omega_{\phi}t)$ terms. The coefficients $A_{n}$ are of the order
$\alpha_{s}^{-2n}$ and do not depend on the initial coordinates of the
particle. These terms form the particular solution:

$\phi_{p}=\sum_{n=1}^{\infty}A_{n}\sin(n\omega_{\phi}t).$ (13)

There is no stable equilibrium point inside of the bucket. The particular
solution is analogous to a moving bucket center; we term the trajectory of a
particle at the moving bucket center to be quasi-synchronous because the
frequency spectrum will depend only on harmonics of $\omega_{\phi}$.

We continue the perturbation to the general case by using
$\phi=\phi_{p}{+\rho\sin(\omega_{s}t+\psi)}$, the sum of the particular
solution (Eq. 13) and the small-oscillation single-rf solution (Eq. 2). Using
$\phi=\phi_{p}+\rho\sin(\omega_{s}t+\psi)$ generates terms of the form
$B_{m,n}{\sin[m(1+\sigma)\omega_{s}t+n\omega_{\phi}t+m\psi]}$. The shift in
the synchrotron oscillation frequency $\sigma$ is a necessary contribution to
the coefficient of ${\sin[(1+\sigma)\omega_{s}t+\psi]}$ in Eq. 12 to
counterbalance the contribution made by the cross-multiplication of higher
order $B_{m,n}{\sin[m(1+\sigma)\omega_{s}t+n\omega_{\phi}t+m\psi]}$ terms. For
any integer $m>0$ and any integer $n$, the coefficients $B_{m,n}$ are of the
order $\rho^{m}\alpha_{s}^{-2|n|}$; except when $m$ is even and $n=0$, in
which case the coefficients $B_{m,0}$ are of the order
$\rho^{m}\alpha_{s}^{-2}$. Writing out the full perturbative solution, we
have:

$\displaystyle\phi=$
$\displaystyle\sum_{n=1}^{\infty}A_{n}\sin(n\omega_{\phi}t)$
$\displaystyle+\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}B_{m,n}\sin[m(1+\sigma)\omega_{s}t+n\omega_{\phi}t+m\psi]$
$\displaystyle+\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}B_{m,-n}\sin[m(1+\sigma)\omega_{s}t-n\omega_{\phi}t+m\psi].$
(14)

The trajectory of a particle in a slip-stacking rf bucket is referred to as a
rotating solution in the driven-pendulum literature. The particular solution
was previously obtained by Zhang and Ma Zhang . An alternate perturbative
approach for the general solution in implicit form is given in Lenci . We are
first to find a general and explicit solution.

The perturbative solution for the small oscillations around the moving bucket
center can be expressed in coefficients up to order $\alpha_{s}^{-4}$ and
$\rho\alpha_{s}^{-2}$. The derivation shown in the appendix leads to the
equations of motion:

$\displaystyle\phi=$ $\displaystyle
A_{1}\sin(\omega_{\phi}t)+A_{2}\sin(2\omega_{\phi}t)$
$\displaystyle+\rho\sin[(1+\sigma)\omega_{s}t+\psi]$
$\displaystyle+B_{1,1}\sin[(1+\sigma)\omega_{s}t+\omega_{\phi}t+\psi]$
$\displaystyle+B_{1,-1}\sin[(1+\sigma)\omega_{s}t-\omega_{\phi}t+\psi].$ (15)
$\displaystyle\delta=$ $\displaystyle\frac{1}{2\pi f_{rev}h\eta}\dot{\phi}.$
(16) $\displaystyle A_{1}=$ $\displaystyle-\frac{1}{\alpha_{s}^{2}-1}.$ (17)
$\displaystyle A_{2}=$
$\displaystyle\frac{1}{(2\alpha_{s})^{2}-1}\left(\frac{A_{1}}{2}\right).$ (18)
$\displaystyle B_{1,\pm 1}=$
$\displaystyle\frac{\alpha_{s}^{-1}}{\alpha_{s}\pm
2}\left(\frac{\rho}{2}\right).$ (19) $\displaystyle\sigma=$
$\displaystyle\frac{3}{4}\alpha_{s}^{-4}.$ (20)

The parameters $\rho$ and $\psi$ are determined by initial conditions, shown
explicitly in the appendix. We are the first to discover and calculate
$\sigma$, the synchrotron frequency shift from a slip-stacking perturbation.
Generally, particles within a slip-stacking bucket will undergo synchrotron
oscillations at a higher frequency than the corresponding single-rf bucket.

Substituting the perturbative terms from Eq. 14 into Eq. 11 indicates that a
new parametric resonance will occur wherever
$m\omega_{s}(1+\sigma)=n\omega_{\phi}$. For example, the
$\rho^{m}\alpha_{s}^{-2(n-1)}\sin[m(1+\sigma)\omega
t-(n-1)\omega_{\phi}t+m\psi]$ term will be multiplied by
$\cos(\omega_{\phi}t)$ in Eq. 11 and lead to a growth term proportional to
$\rho^{m}\alpha_{s}^{-2(n-1)}\sin(m\psi)$. The case where
$m\omega_{s}=\omega_{\phi}$ was previously investigated by Mills Mills . An
analytical description of the stable phase-space boundary would require a
complete determination of the cases in which parametric resonances lead to
particle loss; this may be the subject of future work.

## Stability Maps & Area Factors

The size and shape of slip-stacking buckets determine which portion of an
injected beam distribution is lost. In application, lost particles migrate to
an incorrect azimuthal location and consequently collide into the beampipe
during injection, extraction, or acceleration Brown . We map the stability of
initial particle positions by integrating the equations of motion for each
position. The integration is iterated for a sufficiently large number of
revolutions (at least 30 synchrotron periods). The stability of the particle
is tested after every phase-slipping period. A particle is considered lost if
its phase with respect to each of the first rf cavity, the second rf cavity,
and the average of the two rf cavities, is larger than a certain cut-off (we
used $3\pi/2$). The remaining particles therefore belong to one of four stable
regions shown in Fig. 1: one for the higher frequency, one for the lower
frequency, one for the average frequency and average phase, and one with
average frequency but $\pi$ offset from the average phase. These two stable
regions at the average frequency are the original examples of dynamic
stabilization Butikov . Fig. 2 shows the stability of initial coordinates in
the higher bucket for $\alpha_{s}=3.6$ and $\alpha_{s}=4.1$, in which the
effects from slip-stacking resonances are evident. The supplemental material
SM  shows the stability maps of the higher slip-stacking bucket for values of
$\alpha_{s}$ from 2 to 8 in descending 0.1 increments.

We find some trajectories are “metastable” because they lead to particle loss
only after thousands of revolutions. The stable phase-space area as a function
of time is shown in Fig. 3 for several values of $\alpha_{s}$.

Figure 1: Stability of initial coordinates for $\alpha_{s}=4.4$. The color
corresponds to the number of synchrotron periods a particle with the
corresponding initial coordinates survives before it is lost. The two large
stable regions correspond to the higher and lower rf buckets where beam is
injected and maintained. The two stable regions along the $\delta_{A}=0$ axis
are created by the interaction between the two rf cavities.

Figure 2: Stability of initial coordinates for selected values of
$\alpha_{s}$. The color corresponds to the number of synchrotron periods a
particle with the corresponding initial coordinates survives before it is
lost. On the left, $\alpha_{s}=3.6$ and four resonance islands can be seen due
to the $\omega_{s}(1+\sigma)=4\omega_{\phi}$ resonance. On the right,
$\alpha_{s}=4.1$ and five resonance islands can be seen due to the
$\omega_{s}(1+\sigma)=5\omega_{\phi}$ resonance. Figure 3: The stable area of
the slip-stacking bucket relative to a single rf bucket, is plotted on a log
scale and plotted over time. Each curve corresponds to a simulation with a
different value of $\alpha_{s}$ with $\alpha_{s}=$ 2.0, 2.5, 3.0, 3.5, 4.0,
4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0 (going from the bottom line to the top
line). The rapid losses at the beginning corresponds to regions of phase space
in which particles rapidly slip by both slip-stacking buckets. In the next
phase, the metastable particle loss occurs asymptotically.

The bucket area is computed as the product of the total number of ultimately
surviving points and the cell area. We define the slip-stacking area factor
$F(\alpha_{s})=\mathcal{A}_{s}/\mathcal{A}_{0}$ as the ratio of the slip-
stacking bucket area to that of a single-rf bucket with the same rf voltage
and frequency. The area factor follows the notation of Lee (Ch. 3.II of SYLee
) for accelerating beams, in which the ratio of running bucket area to
stationary bucket area is used. Particles in the bucket are described by Eq.
14 with finite coefficients, therefore the bucket area is conserved.
Consequently $F(\alpha_{s})$ does not depend on the initial rf phase
difference used to generate the stability map. We write the phase space area
($\phi\cdot\delta$ units) using $F(\alpha_{s})$:

$\mathcal{A}_{s}=\mathcal{A}_{0}F(\alpha_{s})=\frac{16}{h|\eta|}\frac{\omega_{s}}{\omega_{rev}}F(\alpha_{s}).$
(21)

Fig. 4(a) plots the numerically derived slip-stacking area factor
$F(\alpha_{s})$. Using Fig. 4(a) with Eq. 21 provides the first method for
calculating the slip-stacking stable phase-space area without requiring each
case to be simulated individually. $F(\alpha_{s})$ increases rapidly above
$\alpha_{s}\approx 3$ and asymptotically approaches 1. $F(\alpha_{s})$ has
several local minimum where resonances are crossed; this loss of area occurs
when large amplitude trajectories have a parametric resonance and therefore
does not occur at precise integer values of $\alpha_{s}$.

Figure 4: (a) The slip-stacking area factor as a function of $\alpha_{s}$. As
$\alpha_{s}$ increases the distance between the rf buckets becomes greater,
the buckets become more independent, and the slip-stacking bucket area
approaches the single-rf bucket area.
(b) The modified slip-stacking area factor as a function of $\alpha_{s}$. The
modified slip-stacking area factor is maximized near $\alpha_{s}=6.2$.

In application, slip-stacking is tuned to maximize stable phase-space area
while holding $\omega_{\phi}$ constant. The value of $\omega_{\phi}$ is
generally constrained by gross features of the accelerators, for example the
harmonic number and cycle time. The slip-stacking parameter $\alpha_{s}$ is
tuned through changing $\omega_{s}$ which is proportional to the square root
of the applied rf voltage. Furthermore $\omega_{s}$ changes the bucket area by
both the slip-stacking area factor $F(\alpha_{s})$ and the single-rf bucket
area, so there is an optimal voltage in which phase space area is maximized.
We rewrite Eq. 21 to separate the parameters that are held constant from those
dependent on $\alpha_{s}$:

$\mathcal{A}_{s}=\frac{16}{h|\eta|}\frac{\omega_{\phi}}{\omega_{rev}}\left(\frac{F(\alpha_{s})}{\alpha_{s}}\right)=\frac{16}{h|\eta|}\frac{\omega_{\phi}}{\omega_{rev}}Z(\alpha_{s}).$
(22)

This modified area factor $Z(\alpha_{s})$ is graphed in Fig. 4(b).
$Z(\alpha_{s})$ is maximal near $\alpha_{s}=6.2$ and when considering other
optimization criteria 5.5 to 7 is is a practical tuning range for
$\alpha_{s}$.

## Injection Efficiency, Emittance, and Aspect Ratio

The stability maps can also be used to analyze injection scenarios, by
weighting the (scaled) stability maps according to a distribution that
represents the number of incoming particles injected into that region of
phase-space. We used this technique to identify the greatest longitudinal
emittance an incoming Gaussian-distributed beam could have and still achieve
97% injection efficiency at its optimal value of $\alpha_{s}$. The
longitudinal beam emittance is given in Eq. 23 below:

$\epsilon=\pi\sigma_{p}\sigma_{T},~{}\epsilon_{97\%}=2.17^{2}\pi\sigma_{p}\sigma_{T}$
(23)

The current accelerator upgrade proposal, Proton Improvement Plan II (PIP-II)
PIP , defines a minimum 97% slip-stacking efficiency required to maintain
current loss levels while increasing intensity. Fig. 5 shows the 97%
longitudinal emittance as a function of aspect ratio and demonstrates the
consequences of a mismatched injection into a slip-stacking bucket. Fig. 6
shows the optimal value of $\alpha_{s}$ as a function of aspect ratio. The
optimal value of $\alpha_{s}$ determines the optimal rf cavity voltage, shown
in Fig. 7. These results were obtaining using parameter values specific to
slip-stacking in the Fermilab Recycler (see Table. 1).

Figure 5: The maximum 97% emittance at 97% efficiency (at an optimal value of $\alpha_{s}$) is shown as a function of aspect ratio. The bottom line (black) is for the 15-Hz Booster cycle-rate (status quo) and top line (red) is for 20-Hz Booster cycle-rate (proposed upgrade). The vertical dashed lines represent the nominal aspect ratios given in Table 1. Figure 6: The optimal slip-stacking parameter $\alpha_{s}$ for maximum emittance (at 97% efficiency) is shown as a function of aspect ratio. Values of $\alpha_{s}$ greater than 8 are not evaluated. The bottom line (black) is for the 15-Hz Booster cycle-rate (status quo) and top line (red) is for 20-Hz Booster cycle-rate (proposed upgrade). The vertical dashed lines represent the nominal aspect ratios given in Table 1. Figure 7: The optimal Recycler rf voltage for maximum emittance (at 97% efficiency) is shown as a function of aspect ratio. The bottom line (black) is for the 15-Hz Booster cycle-rate (status quo) and top line (red) is for 20-Hz Booster cycle-rate (proposed upgrade). The vertical dashed lines represent the nominal aspect ratios given in Table 1. Recycler Kinetic Energy ($E$) | 8 GeV
---|---
Recycler Reference RF freq. ($f$) | 52.8 MHz
Recycler Harmonic number ($h$) | 588
Recycler Phase-slip factor ($\eta$) | -8.6*$10^{-3}$
Maximum Recycler RF Voltage ($V$) | 2 $\times$ 150 kV
Booster harmonic number ($h_{B}$) | 84
Booster cycle rate ($f_{B}$) | 15/20 Hz
Difference in Recycler RF freq. ($\Delta f$) | 1260/1680 Hz
Nominal Booster emittance ($\epsilon_{97\%}$) | 0.12 eV$\cdot$s
Nominal Booster Aspect Ratio | 3.00 MeV/ns
Nominal Recycler Aspect Ratio (100 kV) | 1.06 MeV/ns
Nominal Recycler Aspect Ratio (57 kV) | 0.80 MeV/ns
Table 1: Recycler and Booster parameters used in analysis.

A nominal value for the Booster emittance is 0.12 eV$\cdot$s SeiyaCD . The
Fermilab Booster uses bunch rotation via quadrupole excitation Yang Ahrens ,
with parameters that are actively tuned to minimize losses. With bunch
rotation, the aspect ratio of at least 2.6 MeV/ns is achievable at extraction
from the Booster SeiyaCD . At Recycler rf cavity voltage $V_{0}=$ 100kV, the
slip-stacking parameter for the Recycler is $\alpha_{s}(V_{0})\approx 4.39$
for a 15-Hz Booster cycle-rate and $\alpha_{s}(V_{0})\approx 5.86$ for a 20-Hz
Booster cycle-rate. For other voltages, the Recycler slip-stacking parameter
is given by $\alpha_{s}(V)=\alpha_{s}(V_{0})\sqrt{V/V_{0}}$.

We examine the 97% efficiency benchmark not only for a 15-Hz Booster cycle-
rate but also for a proposed 20-Hz Booster cycle-rate. A 20-Hz Booster cycle-
rate would necessarily increase the phase-slipping frequency (rf frequency
separation) by a factor of 4/3 and therefore the bucket height would increase
by a factor of 4/3 (and have the same $\alpha_{s}$). Slip-stacking losses
could be reduced by a factor 4-10 for the same emittance Eldred .
Alternatively, the 97% efficiency benchmark is achievable at an emittance up
to a factor of $\sim$1.7 greater. Other implications of a 20-Hz Booster cycle-
rate are discussed in a Fermilab technical memo Eldred . A 20-Hz Booster
cycle-rate is clearly superior for high-intensity operation.

The scaling symmetry used to analyze the 20-Hz Booster cycle-rate can
generalized. An optimization at phase-slipping frequency $\omega_{\phi}$ and
aspect ratio $r$ is equivalent to an optimization at phase-slipping frequency
$\omega_{\phi}^{\prime}$ and aspect ratio
$(\omega_{\phi}^{\prime}/\omega_{\phi})r$. The same optimal slip-stacking
parameter would be obtained at a higher synchrotron frequency
$(\omega_{\phi}^{\prime}/\omega_{\phi})\omega_{s}$, increasing the rf voltage
at $(\omega_{\phi}^{\prime}/\omega_{\phi})r$ to
$(\omega_{\phi}^{\prime}/\omega_{\phi})^{2}V$.

## Application to Other Physical Systems

In general, the dynamics discussed in this paper apply to any system governed
by (nearly) identical sinusoidal potentials moving with respect to each other
(or equivalently, a sinusoidal potential oscillating in amplitude). This is
relevant to standing wave traps which are used in optical and acoustic physics
and are instances of a controllable sinusoidal potential. Optical lattices are
a type of standing wave trap used in ultracold atomic physics. Two optical
lattices moving with respect to each other could occupy the same transverse
space yet focus two groups of atoms with independent momenta. We know of no
such experiment that utilizes this optical slip-stacking laser configuration,
but we believe it may be relevant to at least two optical applications: trap
accumulation Davies  and atomic collisions Jaksch . If $M$ is the mass of the
atom, $V_{H}$ is the potential barrier height, and $v$ is the relative
velocity between the two standing wave traps, then the optical slip-stacking
parameter is given by:

$\alpha_{s}=v\sqrt{\frac{M}{2V_{H}}}.$ (24)

The relative velocity $v$ can be calculated from the frequency difference in
the standing waves $v=\Delta f$ Clairon . Our approach suggests a value of
$\alpha_{s}$ equal to at least 5 for efficient stacking.

Acoustic standing waves can be used to trap small spheres or droplets in a
sinusoidal potential originally described by Gor’kov Gorkov . This technique
has grown in sophistication and application Foresti Shields . Eq. 24 would
also determine the stability of objects (with mass $M$) in a possible acoustic
slip-stacking configuration.

## Conclusion

In summary, we have provided a framework for addressing both the trajectory
and stability of particles in a slip-stacking potential. We introduce the
slip-stacking area factor $F(\alpha_{s})$ and the modified area factor
$Z(\alpha_{s})$ as tools to calculate the stable slip-stacking bucket area for
any combination of accelerator parameters. We introduce the quasi-synchronous
particle trajectory and provide a perturbative solution near it. We describe a
series of new parametric resonances in slip-stacking. We provide a general
method for analyzing slip-stacking injection scenarios and describe the
implications for the operation and upgrades of the Fermilab Booster. We
identify for the first time how the dynamics of slip-stacking correspond to
the driven pendulum and moving standing wave traps.

## Acknowledgments

This work is supported in part by grants from the US Department of Energy
under contract DE-FG02-12ER41800 and the National Science Foundation NSF
PHY-1205431. Special thanks to SY Lee for providing a crucial mentoring role
immediately prior to the beginning of this research.

##  Appendix A Derivation of Perturbative Solution

In this paper we show that the second-order equation of motion for a single
particle in a slip-stacking bucket is given by Eq. 11 and the perturbative
solution is given by Eq. 14. Recall also, that the coefficients $A_{n}$ are of
order $\alpha_{s}^{-2n}$ and $B_{m,n}$ are of order
$\rho^{m}\alpha_{s}^{-2|n|}$ except the $B_{2k,0}$ coefficients which are of
the order $\rho^{2k}\alpha_{s}^{-2}$. The coefficient $B_{1,0}$ is defined to
be equal to $\rho$. The parameters $\rho$ and $\psi$ are set by initial
conditions.

For clarity, we adopt a short-hand notation for the oscillatory terms as
follows:

$\displaystyle\sin[m(1+\sigma)\omega_{s}t+n\omega_{\phi}t+m\psi]\equiv
s_{m,n},$
$\displaystyle\cos[m(1+\sigma)\omega_{s}t+n\omega_{\phi}t+m\psi]\equiv
c_{m,n}.$ (25)

In this appendix, we explicitly obtain the perturbative solution up to order
$\alpha_{s}^{-4}$ and $\rho\alpha_{s}^{-2}$ (or equivalently, all coefficients
up to order $\alpha_{s}^{-5}$ with $\rho\sim\alpha_{s}^{-3}$). Therefore, we
use the coefficients $A_{1}$, $A_{2}$, $\rho$, $B_{1,1}$ and $B_{1,-1}$; all
other coefficients are neglected at this precision. We start by assuming a
solution form and will demonstrate it to be self-consistent:

$\phi=A_{1}s_{0,1}+A_{2}s_{0,2}+\rho s_{1,0}+B_{1,1}s_{1,1}+B_{1,-1}s_{1,-1}.$
(26)

It will be sufficient to substitute this expression into the form of Eq. 11
expanded up to second order in $\phi$:

$\omega_{\phi}^{-2}\ddot{\phi}=-\alpha_{s}^{-2}\left[\phi(1+c_{0,1})-s_{0,1}+\frac{1}{2}\phi^{2}s_{0,1}\right].$
(27)

Here we have divided both sides by $\omega_{\phi}^{2}$ to make the order of
the perturbation terms more explicit.

Next we can substitute the solution given in Eq. 26 into Eq. 27 to calculate
the coefficients. We write out the left-hand side (LHS):

$\displaystyle-A_{1}s_{0,1}-4A_{1}s_{0,2}$
$\displaystyle-\alpha_{s}^{-2}(1+\sigma)^{2}\rho s_{1,0}$
$\displaystyle-[\alpha_{s}^{-1}(1+\sigma)+1]^{2}B_{1,1}s_{1,1}$
$\displaystyle-[\alpha_{s}^{-1}(1+\sigma)-1]^{2}B_{1,-1}s_{1,-1}=$ (28)

We write out the right-hand side (RHS):

$\displaystyle=-\alpha_{s}^{-2}\bigg{(}$ $\displaystyle
A_{1}s_{0,1}+A_{1}s_{0,1}c_{0,1}+\frac{1}{2}A_{1}^{2}s_{0,1}s_{0,1}s_{0,1}$
$\displaystyle+A_{2}s_{0,2}+A_{2}s_{0,2}c_{0,1}-s_{0,1}$
$\displaystyle+\frac{1}{2}\cancel{A_{2}^{2}}s_{0,2}s_{0,2}s_{0,1}+\cancel{A_{1}A_{2}}s_{0,1}s_{0,2}s_{0,1}$
$\displaystyle+\rho s_{1,0}+\rho s_{1,0}c_{0,1}$
$\displaystyle+B_{1,1}s_{1,1}+B_{1,1}s_{1,1}c_{0,1}$
$\displaystyle+B_{1,-1}s_{1,-1}+B_{1,-1}s_{1,-1}c_{0,1}$
$\displaystyle+\frac{1}{2}\cancel{(\rho
s_{1,0}+B_{1,1}s_{1,1}+B_{1,1}s_{1,-1})^{2}}s_{0,1}$ $\displaystyle+A_{1}\rho
s_{0,1}s_{1,0}s_{0,1}+\cancel{A_{2}\rho}s_{0,2}s_{1,0}s_{0,1}$
$\displaystyle+\cancel{A_{1}B_{1,1}}s_{0,1}s_{1,1}s_{0,1}+\cancel{A_{1}B_{1,-1}}s_{0,1}s_{1,-1}s_{0,1}$
$\displaystyle+\cancel{A_{2}B_{1,1}}s_{0,2}s_{1,1}s_{0,1}+\cancel{A_{2}B_{1,-1}}s_{0,2}s_{1,-1}s_{0,1}\bigg{)}.$
(29)

where the crossed-out terms are higher order than the precision of this
analysis. We rewrite the RHS without these immediately negligible terms, using
the trigonometric product-to-sum rules, and grouping by the oscillatory term:

$\displaystyle=-\alpha_{s}^{-2}\bigg{[}$
$\displaystyle\left(A_{1}-1+\frac{1}{2}A_{2}+\frac{3}{8}A_{1}^{2}\right)s_{0,1}$
$\displaystyle+\left(\frac{1}{2}A_{1}+A_{2}\right)s_{0,2}$
$\displaystyle+\left(\rho+\frac{1}{2}B_{1,1}+\frac{1}{2}B_{1,-1}+\frac{1}{2}A_{1}\rho\right)s_{1,0}$
$\displaystyle+\left(\frac{1}{2}\rho+B_{1,1}\right)s_{1,1}+\left(\frac{1}{2}\rho+B_{1,-1}\right)s_{1,-1}$
$\displaystyle+\left(\frac{1}{2}A_{2}-\frac{1}{8}A_{1}^{2}\right)s_{0,3}$
$\displaystyle+\left(\frac{1}{2}B_{1,1}-\frac{1}{4}A_{1}\rho\right)s_{1,2}$
$\displaystyle+\left(\frac{1}{2}B_{1,-1}-\frac{1}{4}A_{1}\rho\right)s_{1,-2}\bigg{]}.$
(30)

We then equate Eq. 28 with Eq. 30 for all time. Each oscillatory term
corresponds to its own equation:

$\displaystyle-A_{1}$
$\displaystyle=-\alpha_{s}^{-2}\left(A_{1}-1+\frac{1}{2}A_{2}+\frac{3}{8}A_{1}^{2}\right).$
(31) $\displaystyle-4A_{2}$
$\displaystyle=-\alpha_{s}^{-2}\left(\frac{1}{2}A_{1}+A_{2}\right).$ (32)
$\displaystyle-\alpha_{s}^{-2}(1+\sigma)^{2}\rho$
$\displaystyle=-\alpha_{s}^{-2}\left(\rho+\frac{1}{2}B_{1,1}+\frac{1}{2}B_{1,-1}+\frac{1}{2}A_{1}\rho\right).$
(33) $\displaystyle-[\alpha_{s}^{-1}(1+\sigma)+1]^{2}B_{1,1}$
$\displaystyle=-\alpha_{s}^{-2}\left(\frac{1}{2}\rho+B_{1,1}\right).$ (34)
$\displaystyle-[\alpha_{s}^{-1}(1+\sigma)-1]^{2}B_{1,-1}$
$\displaystyle=-\alpha_{s}^{-2}\left(\frac{1}{2}\rho+B_{1,-1}\right).$ (35)
$\displaystyle 0$
$\displaystyle=-\alpha_{s}^{-2}\left(\frac{1}{2}A_{2}-\frac{1}{8}A_{1}^{2}\right).$
(36) $\displaystyle 0$
$\displaystyle=-\frac{1}{2}\alpha_{s}^{-2}B_{1,1}+\frac{1}{4}\alpha_{s}^{-2}A_{1}\rho.$
(37) $\displaystyle 0$
$\displaystyle=-\frac{1}{2}\alpha_{s}^{-2}B_{1,-1}+\frac{1}{4}\alpha_{s}^{-2}A_{1}\rho.$
(38)

Solving Eq. 32 for $A_{2}$ we obtain:

$A_{2}=\frac{1}{4-\alpha_{s}^{-2}}\left(\frac{A_{1}}{2}\right)=\frac{1}{(2\alpha_{s})^{2}-1}\left(\frac{A_{1}}{2}\right).$
(39)

Solving Eq. 31 for the linear $A_{1}$ terms we obtain:

$A_{1}=\frac{\alpha_{s}^{-2}}{1-\alpha_{s}^{-2}}\left(-1+\frac{1}{2}A_{2}+\frac{3}{8}A_{1}^{2}\right).$
(40)

Since $A_{1}$ of the order $\alpha_{s}^{-2}$ (as expected) then the
$\alpha_{s}^{-2}A_{2}$ term and the $\alpha_{s}^{-2}A_{1}^{2}$ are of order
$\alpha_{s}^{-6}$ and are neglected. Rewriting Eq. 31 to reflect this, we
have:

$A_{1}=-\frac{\alpha_{s}^{-2}}{1-\alpha_{s}^{-2}}=-\frac{1}{\alpha_{s}^{2}-1}.$
(41)

All expressions on the RHS side of Eq. 36 are of order $\alpha_{s}^{-6}$,
therefore we are self-consistent to exclude the $A_{3}$ term.

We solve Eq. 34 and Eq. 35 for $B_{1,\pm 1}$ and obtain:

$\displaystyle B_{1,\pm 1}$
$\displaystyle=\frac{\alpha_{s}^{-2}}{[\alpha_{s}^{-1}\pm
1]^{2}-\alpha_{s}^{-2}}\left(\frac{\rho}{2}\right)$
$\displaystyle=\frac{\alpha_{s}^{-2}}{1\pm
2\alpha_{s}^{-1}}\left(\frac{\rho}{2}\right)$
$\displaystyle=\frac{\alpha_{s}^{-1}}{\alpha_{s}\pm
2}\left(\frac{\rho}{2}\right).$ (42)

where $\sigma$ makes a negligible contribution to $B_{1,\pm 1}$.

All expressions on the RHS of Eq. 37 and Eq. 38 are of order
$\alpha_{s}^{-4}\rho$, therefore we are self-consistent to exclude the
$B_{1,2}$ and $B_{1,-2}$ terms.

For Eq. 33 we call the $\sigma^{2}$ terms negligible, subtract the bare $\rho$
term from each side and solve for $\sigma$ to obtain the shift in synchrotron
frequency:

$\sigma=\frac{1}{4}\left(\frac{B_{1,1}+B_{1,-1}}{\rho}+A_{1}\right).$ (43)

To calculate Eq. 43 first we must calculate:

$\displaystyle\frac{B_{1,1}+B_{1,-1}}{\rho}$
$\displaystyle=\frac{\alpha_{s}^{-1}}{2}\left(\frac{1}{\alpha_{s}+2}+\frac{1}{\alpha_{s}-2}\right)$
$\displaystyle=\frac{\alpha_{s}^{-1}}{2}\left(\frac{\alpha_{s}-2+\alpha_{s}+2}{\alpha_{s}^{2}-4}\right)$
$\displaystyle=\frac{1}{\alpha_{s}^{2}-4}.$ (44)

We substitute Eq. 44 into Eq. 43 to obtain:

$\displaystyle\sigma$
$\displaystyle=\frac{1}{4}\left(\frac{1}{\alpha_{s}^{2}-4}-\frac{1}{\alpha_{s}^{2}-1}\right)$
$\displaystyle=\frac{1}{4}\left(\frac{\alpha_{s}^{2}-1-\alpha_{s}^{2}+4}{(\alpha_{s}^{2}-4)(\alpha_{s}^{2}-1)}\right)$
$\displaystyle=\frac{3}{4}\alpha_{s}^{-4}.$ (45)

To summarize, we write these coefficients together:

$\displaystyle A_{1}$
$\displaystyle=-\frac{1}{\alpha_{s}^{2}-1}=-\alpha_{s}^{-2}(1+\alpha_{s}^{-2}).$
(46) $\displaystyle A_{2}$
$\displaystyle=\frac{1}{(2\alpha_{s})^{2}-1}\left(\frac{A_{1}}{2}\right)=-\frac{1}{8}\alpha_{s}^{-4}.$
(47) $\displaystyle B_{1,\pm 1}$
$\displaystyle=\frac{\alpha_{s}^{-1}}{\alpha_{s}\pm
2}\left(\frac{\rho}{2}\right).$ (48) $\displaystyle\sigma$
$\displaystyle=\frac{1}{4}\left(\frac{B_{1,1}+B_{1,-1}}{\rho}+A_{1}\right)=\frac{3}{4}\alpha_{s}^{-4}.$
(49)

## Derivation of $\rho$ and $\psi$ from Initial Conditions

At this the synchrotron amplitude $\rho$ and initial synchrotron phase $\psi$
are still undetermined, but we can express them in terms of the initial
coordinates. At time $t=0$ (which is fixed to a time in which the relative
phase between the rf cavities is zero), let $\delta=\delta_{0}$ and
$\phi=\phi_{0}$.

We have shown that $\phi$ takes the form give in Eq. 26. We evaluate this
expression for $\phi$ at $t=0$, leave the short-hand notation (Eq. 25), and
find $\phi_{0}$ in terms of $\rho$ and $\psi$:

$\displaystyle\phi_{0}=(\rho+B_{1,1}+B_{1,-1})\sin(\psi).$ (50)

We can calculate $\delta$ from our solution for $\phi$ by taking the
derivative:

$\displaystyle\delta$ $\displaystyle=\frac{1}{2\pi f_{rev}h\eta}$
$\displaystyle\\!\\!\\!\\!\dot{\phi}.$ (51) $\displaystyle\delta$
$\displaystyle=\frac{\omega_{\phi}}{2\pi
f_{rev}h\eta}\bigg{\\{}\\!\\!\\!\\!\\!\\!$ $\displaystyle
A_{1}c_{0,1}+2A_{2}c_{0,2}+\alpha_{s}^{-1}(1+\sigma)\rho c_{1,0}$
$\displaystyle+[\alpha_{s}^{-1}(1+\sigma)+1]B_{1,1}c_{1,1}$
$\displaystyle+[\alpha_{s}^{-1}(1+\sigma)-1]B_{1,-1}c_{1,-1}\bigg{\\}}.$ (52)

We evaluate this expression for $\delta$ at $t=0$, leave the short-hand
notation (Eq. 25), and find $\delta_{0}$ in terms of $\rho$ and $\psi$:

$\displaystyle\delta_{0}=\frac{\omega_{\phi}}{2\pi
f_{rev}h\eta}\bigg{\\{}A_{1}$ $\displaystyle+2A_{2}$ $\displaystyle+\bigg{[}$
$\displaystyle\left(\frac{1+\sigma}{\alpha_{s}}\right)(\rho+B_{1,1}+B_{1,-1})$
$\displaystyle+\alpha_{s}(B_{1,1}-B_{1,-1})\bigg{]}\cos(\psi)\bigg{\\}}.$ (53)

Next we solve Eq. 50 and Eq. 53 for $\Phi_{0}=\rho\sin(\psi)$ and
$\Delta_{0}=\rho\cos(\psi)$:

$\displaystyle\Phi_{0}=$
$\displaystyle\left(1+\frac{B_{1,1}+B_{1,-1}}{\rho}\right)^{-1}\phi_{0}.$ (54)
$\displaystyle\Delta_{0}=\alpha_{s}\bigg{[}$ $\displaystyle
1+\frac{B_{1,1}+B_{1,-1}}{\rho}+\alpha_{s}\frac{B_{1,1}-B_{1,-1}}{\rho}$
$\displaystyle+\sigma\left(1+\frac{B_{1,1}+B_{1,-1}}{\rho}\right)\bigg{]}^{-1}$
$\displaystyle\times\bigg{(}$ $\displaystyle\frac{2\pi
f_{rev}h\eta}{\omega_{\phi}}\delta_{0}-A_{1}-2A_{2}\bigg{)}.$ (55)

$\phi_{0}$ and $\delta_{0}$ have been translated by the initial position of
the bucket center and rescaled to obtain the expressions for $\Phi_{0}$ and
$\Delta_{0}$. Using $\Phi_{0}=\rho\sin(\psi)$ and $\Delta_{0}=\rho\cos(\psi)$,
the solution for $\rho$ and $\psi$ be written as:

$\displaystyle\rho$ $\displaystyle=\sqrt{\Phi_{0}^{2}+\Delta_{0}^{2}}.$ (56)
$\displaystyle\psi$
$\displaystyle=\operatorname{sgn}(\Phi_{0})\arccos\left(\frac{\Delta_{0}}{\rho}\right).$
(57)

Eq. 54 and Eq. 55 can be further simplified by writing the $B_{1,1}$ and
$B_{1,-1}$ terms explicitly in terms of $\alpha_{s}$. We calculate:

$\displaystyle\alpha_{s}\frac{B_{1,1}-B_{1,-1}}{\rho}$
$\displaystyle=\frac{1}{2}\left(\frac{1}{\alpha_{s}+2}-\frac{1}{\alpha_{s}-2}\right)$
$\displaystyle=\frac{1}{2}\left(\frac{\alpha_{s}-2-\alpha_{s}-2}{\alpha_{s}^{2}-4}\right)$
$\displaystyle=\frac{-2}{\alpha_{s}^{2}-4}.$ (58)

We apply Eq. 44 and Eq. 58 to Eq. 54 and Eq. 55 to obtain:

$\displaystyle\Phi_{0}=$
$\displaystyle\left(1+\frac{1}{\alpha_{s}^{2}-4}\right)^{-1}\\!\\!\\!\\!\phi_{0}=\frac{\alpha_{s}^{2}-4}{\alpha_{s}^{2}-3}\phi_{0}.$
(59) $\displaystyle\Delta_{0}=$
$\displaystyle\alpha_{s}\left[1+\frac{1}{\alpha_{s}^{2}-4}-\frac{2}{\alpha_{s}^{2}-4}+\sigma\left(1+\frac{1}{\alpha_{s}^{2}-4}\right)\right]^{-1}$
$\displaystyle\times\left(\frac{2\pi
f_{rev}h\eta}{\omega_{\phi}}\delta_{0}-A_{1}-2A_{2}\right).$ (60)
$\displaystyle\Delta_{0}=$
$\displaystyle\frac{\alpha_{s}(\alpha_{s}^{2}-4)}{(\alpha_{s}^{2}-5)+\sigma\alpha_{s}^{2}}\left(\frac{2\pi
f_{rev}h\eta}{\omega_{\phi}}\delta_{0}-A_{1}-2A_{2}\right).$ (61)

$\Phi_{0}$ and $\Delta_{0}$ are fully expanded as follows:

$\displaystyle\Phi_{0}=$
$\displaystyle\phi_{0}(1-\alpha_{s}^{-2}-3\alpha_{s}^{-4}).$ (62)
$\displaystyle\Delta_{0}=$ $\displaystyle\alpha_{s}\left(\frac{2\pi
f_{rev}h\eta}{\omega_{\phi}}\delta_{0}+\alpha_{s}^{-2}+\frac{5}{4}\alpha_{s}^{-4}\right)$
$\displaystyle\times\left(1+\alpha_{s}^{-2}+\frac{17}{4}\alpha_{s}^{-4}\right).$
(63)

Eq. 62 and Eq. 63 can then be substituted into Eq. 56 and Eq. 57 to obtain
$\rho$ and $\psi$ respectively. Then $\rho$ can be substituted into Eq. 48 to
obtain $B_{1,1}$ and $B_{1,-1}$.

## References

  * (1) P. Derwent et al. Proton Improvement Plan-II December 2013, 2013.
  * (2) C. Mariani (LBNE/DUSEL Collaborations), in Proceedings of the Neutrino Oscillation Workshop, Otranto, Italy, 2011, edited by P. Bernardini, G. Fogli, and E. Lisi.
  * (3) J. Galambos, M. Bai, and S. Nagaitsev, Snowmass Workshop on Frontier Capability Summary Report, 2013.
  * (4) B. C. Brown, P. Adamson, D. Capista, W. Chou, I. Kourbanis, D. K Morris, K. Seiya, G. H. Wu, and M. J. Yang, Phys. Rev. ST Accel. Beams 16, 071001 (2013).
  * (5) F. E. Mills, Brookhaven National Laboratory Report No. 15936, 1971.
  * (6) D. Boussard and Y. Mizumachi, IEEE Trans. Nucl. Sci. 26, 3623 (1979).
  * (7) J. A. MacLachlan, Fermi National Accelerator Laboratory Report No. 0711, 2001.
  * (8) K. Seiya, T. Berenc, B. Chase, W. Chou, J. Dey, P. Joireman, I. Kourbanis, J. Reid, and D. Wildman, in Proceedings of HB2006, Tsukuba, Japan, 2006, edited by Y. H. Chin, H. Yoshikawa, and M. Ikegami.
  * (9) T. Aaltonen et al. (CDF and D0 Collaborations), Phys. Rev. D 86, 092003 (2012).
  * (10) P. Adamson et al. (MINOS Collaboration), Phys. Rev. Lett. 110, 251801 (2013).
  * (11) L. Fields et al. (MINERvA Collaboration), Phys. Rev. Lett. 111, 022501 (2013).
  * (12) G. A. Fiorentini et al. (MINERvA Collaboration), Phys. Rev. Lett. 111, 022502 (2013).
  * (13) M. Muether, Nucl. in Proceedings of Neutrino Oscillation Workshop, Lecce, Italy, 2002, edited by P. Bernardini, G. Fogli, E. Lisi.
  * (14) K. Seiya et al., in Proceedings of Particle Accelerator Conference, 2005, edited by C. Horak.
  * (15) J. Dey, K. Koba, I. Kourbanis, and J. Reid, in Proceedings of Particle Accelerator Conference, 2007, edited by C. Petit-Jean-Genaz.
  * (16) J. Dey and I. Kourbanis, in Proceedings of Particle Accelerator Conference, 2005, edited by C. Horak.
  * (17) R. Madrak and D. Wildman, in Proceedings of North American Particle Accelerator Conference, 2013, edited by T. Satogata, C. Petit-Jean-Genaz, and V. Schaa.
  * (18) S. Y. Lee, Accelerator Physics, 3rd Ed. (World Scientific, Singapore, 2012).
  * (19) H. W. Broer, I. Hoveijn, M. van Noort, C. Simó, and G. Vegter, J. Dyn. Diff. Eq. 16, 897 (2004).
  * (20) X. Xu, M. Wiercigroch, and M. P. Cartmell, Chaos, Solitons Fractals 23, 1537 (2005).
  * (21) A. J. Lichtenberg and M. A. Lieberman, Regular and Chaotic Dynamics, 2nd Ed. (Springer, California, 1992).
  * (22) H. Zhang and T. W. Ma, Nonlinear Dyn. 70, 2433 (2012).
  * (23) S. Lenci, E. Pavlovskaia, G. Rega, and M. Wiercigroch, J. Sound Vib. 310, 243 (2008).
  * (24) E. I. Butikov, J. Phys. A 44, 295202 (2013).
  * (25) See supplemental material for complete series of stability maps.
  * (26) K. Seiya, B. Chase, J. Dey, P. Joireman, I. Kourbanis, and J. Reid, in Proceedings of CARE-HHH-APD Workshop BEAM’07, 2007, edited by W. Scandale and F. Zimmermann.
  * (27) X. Yang, A. I. Drozhdin, and W. Pellico, in Proceedings of Particle Accelerator Conference, 2007, edited by C. Petit-Jean-Genaz.
  * (28) L. A. Ahrens et al., in Proceedings of Particle Accelerator Conference, 1999, edited by A. Luccio and W. MacKay.
  * (29) J. Eldred and R. Zwaska, Fermi National Accelerator Laboratory Report TM-2587-APC, 2014.
  * (30) H. J. Davies and C. S. Adams, J. Phys. B 33, 4079 (2000).
  * (31) D. Jaksch, H. J. Briegel, J. I. Cirac, C. W. Gardiner, and P. Zoller, Phys. Rev. Lett. 82, 1975 (1999).
  * (32) A. Clairon, C. Salomo, S. Guellati, and W.D. Phillip, Europhys. Lett. 16, 165 (1991).
  * (33) L. P. Gor’kov, Sov. Phys. Dokl. 6, 773 (1962).
  * (34) D. Foresti and D. Poulikakos, Phys. Rev. Lett. 112, 024301 (2014).
  * (35) C. W. Shields, L. M. Johhson, L. Gao, and G. P. López, Langmuir 30, 3923 (2014).