# Moment-preserving Monte-Carlo Coulomb collision method for particle codes

Justin Ray Angus<EMAIL_ADDRESS>Yichen Fu  Vasily Geyko  Dave Grote  David
Larson  Lawrence Livermore National Laboratory, Livermore, CA 94550, USA

###### Abstract

Binary-pairing Monte-Carlo methods are widely used in particle-in-cell codes
to capture effects of small angle Coulomb collisions. These methods preserve
momentum and energy exactly when the simulation particles have equal weights.
However, when the interacting particles are of varying weight, these physical
conservation laws are only preserved on average. Here, we 1) extend these
methods to weighted particles such that the scattering physics is correct on
average, and 2) describe a new method for adjusting the particle velocities
post scatter to restore exact conservation of momentum and energy. The
efficacy of the model is illustrated with various test problems.

††journal: Manuscript for submission to the Journal of Computational Physics
(LLNL-JRNL-867233)

##  1 Introduction

Particle-in-cell (PIC) methods are a numerical approach for solving the phase-
space continuity law known as the Vlasov equation, which governs collisionless
kinetic processes in a plasma. The Vlasov equation is a subset of the more
general Boltzmann equation that additionally includes short-range collisional
processes. Binary-pairing Monte-Carlo collision (MCC) methods are a popular
numerical approach for capturing collisional effects in a PIC code.
Collisional processes in a fully ionized, ideal plasma are dominated by
cumulative small angle Coulomb collisions as governed by the Landau-Fokker-
Planck (LFP) equation [1].

The binary-pairing algorithms by Takizuka and Abe [2] (TA77) and Nanbu [3]
(N97) are two standard MCC methods for modeling the LFP equation [4]. These
algorithms are nearly identical, with the only difference being the formula
used to determine the center-of-mass (CM) scattering angle for a pair of
particles. For similar levels of accuracy, one can use a somewhat larger time
step with N97 compared to TA77 [4]. Both TA77 and N97 identically conserve
momentum and energy, consistent with the physical properties of the collision
operator. However, these methods, as originally formulated, are limited to
equally weighted particles. Extension of these methods to arbitrarily weighted
particles has been considered by Miller [5], Nanbu [6] and most recently by
Higginson et al [7]. In particular, it is shown in Ref. [7] how the weighted-
particle method by Nanbu in Ref. [6] produces erroneous results in certain
scenarios, and an alternative formulation that produces the correct scattering
physics on average for weighted particles is given. However, each of these
weighted-particle methods only preserve momentum and energy on average after
many scattering cycles.

Preservation of momentum and energy on average may work fine for homogeneous
problems. However, for non-homogeneous problems, where particles transport to
different regions on the grid with different physical properties, it may be
difficult to understand the error introduced by non-exact conservation of
momentum or energy during the application of the collision method.
Furthermore, there is increasing interest in fully implicit PIC methods that
are exactly energy conserving [8, 9]. Recently, these exactly energy
conserving methods have been coupled with models for collisions in Refs. [10,
11]. For these methods, it would not be ideal to throw away the exact
conservation of energy property when using weighted particles.

In this work, a binary-pairing MCC method for Coulomb collisions is presented
that 1) produces the correct scattering physics for weighted particles on
average and 2) maintains exact conservation of momentum and energy after each
scattering cycle. The method has two parts. First, a generalization of the
TA77 and N97 methods to arbitrarily weighted particles is derived starting
from the general $N{\times}N$ method. The second part of the algorithm is a
new method to adjust the velocity of the particles post-scatter such that
exact momentum and energy conservation are restored. The method to restore
energy conservation is based on inelastic scattering dynamics and works for
both non-relativistic and relativistic particles.

The remainder of the paper is outlined as follows. Starting with a general
discussion of the Boltzmann collision integral, the basics of small-angle
Coulomb collisions in a plasma relevant to binary-pairing MCC methods is
discussed in Sec. 2. This is followed by a derivation of the weighted-particle
extension of the TA77 and N97 methods in Sec. 3. The moment-correction method
is described in Sec. 4. The efficacy of the method is illustrated with various
test problems in Sec. 5. The paper concludes with some further discussion and
a summary in Sec. 6.

##  2 The Boltzmann collision integral

The Boltzmann equation governing the phase-space distribution function,
$f=f\left(\textbf{x},\textbf{v},t\right)$, for some arbitrary species is

$\frac{\partial f}{\partial
t}+\nabla_{\textbf{x}}\cdot\left(\textbf{v}f\right)+\nabla_{\textbf{v}}\cdot\left(\textbf{a}f\right)=\frac{\partial
f}{\partial t}\bigg{|}_{C},$ (1)

where x and v are the independent real-space and velocity-space coordinates,
respectively, a is the acceleration associated with long range, or collective,
forces, and the term on the right of the equal sign is the Boltzmann collision
integral. This term represents abrupt changes in velocity space owing to
short-range interactions between particles. The left-hand side of Eq. 1 is the
Vlasov equation, the phase-space continuity law for $f$. Particle-in-cell
(PIC) methods are a numerical technique for solving this equation. Monte Carlo
methods are commonly used in PIC codes to model the Boltzmann collision
integral on the right-hand side (RHS), which is the topic of the present
paper.

This Boltzmann collision integral can be written as a combination of sinks and
sources in velocity space (see Appendix B of Ref. [12]). For species $\alpha$
undergoing collisions with species $\beta$, this term can be written as

$\frac{\partial f_{\alpha}}{\partial t}\bigg{|}_{C}=\frac{\partial
f_{\alpha}}{\partial t}\bigg{|}_{\text{in}}-\frac{\partial
f_{\alpha}}{\partial t}\bigg{|}_{\text{out}},$ (2)

where

$\frac{\partial f_{\alpha}}{\partial t}\bigg{|}_{\text{out}}\equiv
f_{\alpha}\int\int
f_{\beta}|\textbf{v}_{\alpha}-\textbf{v}_{\beta}|Id\Omega_{\alpha}d^{3}v_{\beta},\
\text{and }\frac{\partial f_{\alpha}}{\partial
t}\bigg{|}_{\text{in}}\equiv\int\int
f^{\prime}_{\alpha}f^{\prime}_{\beta}|\textbf{v}^{\prime}_{\alpha}-\textbf{v}^{\prime}_{\beta}|I^{\prime}d\Omega^{\prime}_{\alpha}d^{3}v^{\prime}_{\beta}$
(3)

are the velocity-space sink and source of particles, respectively. The
quantity
$I=I\left(|\textbf{v}_{\alpha}-\textbf{v}_{\beta}|,\theta_{\alpha}\right)$
here is the differential cross section that gives the probability that a
particle from species $\alpha$ scatters into the solid angle
$d\Omega_{\alpha}=d\phi_{\alpha}\sin\theta_{\alpha}d\theta_{\alpha}$. The
total cross section is defined as
$\sigma(|\textbf{v}_{\alpha}-\textbf{v}_{\beta}|)\equiv\int
Id\Omega_{\alpha}$. The prime in the latter expression in Eq. 3 denotes post-
scattered values. The sink term describes particles of species $\alpha$
leaving a specific region of velocity space by scattering into all other
locations in velocity space. The source term describes particles entering a
specific region of velocity space by scattering out of all other regions of
velocity space.

In the analysis below, use is made of the Klimintovich representation of the
species distribution functions:

$f_{\alpha}=\sum_{i=1}^{N_{\alpha}}\frac{w_{i}}{\Delta
V}\delta\left(\textbf{v}-\textbf{v}_{i}\right),\ \text{and}\
f_{\beta}=\sum_{j=1}^{N_{\beta}}\frac{w_{j}}{\Delta
V}\delta\left(\textbf{v}-\textbf{v}_{j}\right),$ (4)

where $N_{\alpha}$ is the number of macro-particles (or just particles) used
to describe species $\alpha$ in a given cell of volume $\Delta V$, $w_{i}$ is
the weight of particle $i\in\alpha$, and similarly for species $\beta$. The
spatial dependence of the distribution functions is neglected for simplicity.
Furthermore, for notational simplicity, the subscript $i$ is used to refer to
a specific particle belonging to species $\alpha$ and the subscript $j$ is
used to refer to a specific particle belonging to species $\beta$.

###  2.1 Small angle Coulomb collisions

An expansion of the Boltzmann collision integral (Eqs. 2-3) for small angle
Coulomb collisions in an ideal plasma results in the LFP collision operator
(see Ch. 11 of [13] and Ch. 13 of Ref. [1]). This operator is of the drag-
diffusion type with nonlinear drag and diffusion coefficients. For each
particle pair, this process is a diffusive one for the CM scattering angle,
$\Theta_{ij}$. The mean square CM scattering angle for a single collision of
particle $i\in\alpha$ with particle $j\in\beta$ is computed by averaging
$\Theta$ over all possible impact parameters (see Sec. 3.3 of Ref. [12]) and
can be expressed as

$\left<\Theta_{ij}^{2}\right>\big{|}_{1}=\frac{2\pi b^{2}_{0,ij}}{\pi
b^{2}_{\max}}\ln\Lambda.$ (5)

Here, $b_{\max}=\lambda_{De}$ is the maximum impact parameter for screened
Coulomb collisions with $\lambda_{De}$ the plasma Debye length,
$\ln\Lambda=\ln\left(2\lambda_{De}/b_{0,ij}\right)$ is the Coulomb logarithm,
$b_{0,ij}/2=q_{\alpha}q_{\beta}/\left(4\pi\epsilon_{0}\mu_{\alpha\beta}u_{ij}^{2}\right)$
is the minimum impact parameter with
$\mu_{\alpha\beta}=m_{\alpha}m_{\beta}/\left(m_{\alpha}+m_{\beta}\right)$ the
reduced mass, and $u_{ij}\equiv|\textbf{v}_{i}-\textbf{v}_{j}|$. The total
number of collisions of particle $i$ with particle $j$ in time interval
$\Delta t$ is $\pi b^{2}_{\max}n_{j}u_{ij}\Delta t$, where $n_{j}=w_{j}/\Delta
V$ is the density associated with particle $j$. The total mean square CM angle
for particle $i\in\alpha$ interacting with particle $j\in\beta$ during time
interval $\Delta t$ is

$\left<\Theta_{ij}^{2}\right>_{\text{all
}b}=\left<\Theta_{ij}^{2}\right>\big{|}_{1}\pi b^{2}_{\max}\frac{w_{j}}{\Delta
V}u_{ij}\Delta t=2\pi b^{2}_{0,ij}\ln\Lambda\frac{w_{j}}{\Delta
V}|\textbf{v}_{i}-\textbf{v}_{j}|\Delta t.$ (6)

Being a diffusive process, the distribution of scattering angles $\Theta_{ij}$
after time $\Delta t$ is Gaussian with variance equal to that given in Eq. 6.
The polar scattering angle $\Theta_{ij}$ is thus chosen by sampling from such
a distribution. The total change of velocity of particle $i\in\alpha$ is
obtained by sequentially performing this process for all $j\in\beta$. To see
how this method connects with the LFP equation, the expected value for the
change in velocity of particle $i\in\alpha$ is derived below.

The change in $\textbf{u}_{ij}$ parallel to itself after scattering through a
small polar angle $\Theta_{ij}$ is $\Delta
u_{||}=-2u_{ij}\sin^{2}\left(\Theta_{ij}/2\right)\approx-
u_{ij}\Theta_{ij}^{2}/2$. After averaging $\Delta\textbf{u}_{ij}$ over the
uniformly distributed azimuthal angle $\phi$, the parallel component is the
only non-zero component [13] and we can write

$\left<\Delta\textbf{u}_{ij}\right>_{\text{all }b,\text{ all
}\phi}=-\textbf{u}_{ij}\frac{1}{2}\left<\Theta^{2}_{ij}\right>_{\text{all
}b}.$ (7)

The change in velocity of particle $i\in\alpha$ in the lab frame is related to
$\Delta\textbf{u}_{ij}$ by
$m_{\alpha}\Delta\textbf{v}_{i}=\mu_{\alpha\beta}\Delta\textbf{u}_{ij}$. The
final step in computing the expected value of $\Delta\textbf{v}_{i}$ in time
interval $\Delta t$ is to sum Eq. 7 for all $j\in\beta$:

$\left<\Delta\textbf{v}_{i}\right>=-\frac{\mu_{\alpha\beta}}{m_{\alpha}}\sum_{j=1}^{N_{\beta}}\textbf{u}_{ij}\frac{1}{2}\left<\Theta^{2}_{ij}\right>_{\text{all
}b}=-\frac{q^{2}_{\alpha}q^{2}_{\beta}\ln\Lambda}{4\pi\epsilon_{0}^{2}m_{\alpha}\mu_{\alpha\beta}}\sum_{j=1}^{N_{\beta}}\frac{w_{j}}{\Delta
V}\frac{\textbf{u}_{ij}}{u^{3}_{ij}}\Delta t.$ (8)

After dividing by $\Delta t$, Eq. 8 is the well-known drag coefficient for the
LFP equation (see Eq. 13.28 of Ref. [1]). The tensor diffusion coefficient
($\frac{1}{2}\left<\Delta v_{i}\Delta v_{k}\right>/\Delta t$), can
additionally be obtained following a similar process [13]. Note that each term
in the sum in Eq. 8 scales linearly with the normalized scattering length
$s_{ij}\equiv\frac{1}{2}\left<\Theta^{2}_{ij}\right>_{\text{all }b}$. The same
is true for the diffusion coefficient.

##  3 Monte Carlo method for weighted-particle Coulomb collisions

A method to extend binary-pairing algorithms for Coulomb collisions to
weighted particles is described in this section. This is done for the
$N{\times}N$ binary pairing method and the commonly used (and computationally
practical) order $N$ binary pairing methods by Takizuka and Abe [2] and Nanbu
[3]. For these latter methods, the scattering is done in such a way that the
expected value of the change in velocity for each particle from each species
is equal to that from the $N{\times}N$ method (i.e., Eq. 8) on average after
many pairings. This same constraint is used here to first extend the
$N{\times}N$ scattering method to include weighted particles. Then, the
reduction to an order $N$ binary-pairing method for weighted-particles follows
exactly as it was first done in TA77. For reference, examples of the order
$N{\times}N$ and order $N$ pairing methods considered in this work are
illustrated in Fig. 1. The order $N$ pairing method used here is the same as
that described in TA77.

Fig. 1: Example pairings for order $N^{2}$ and order $N$ pairing methods. The
left two panels are for inter-species collisions with $N_{\alpha}=5$ and
$N_{\beta}=2$. The right two panels are for intra-species collisions with
$N_{\alpha}=5$. The superscript in parentheses denote the number of times a
particle has previously being paired. Note that when $N$ is odd, the first
three particles for the order $N$ pairing method are done in an $N{\times}N$
way.

###  3.1 Order $N{\times}N$ binary pairing method

The discussion in Sec. 2.1 focused on the transport of a single particle
$i\in\alpha$ interacting with species $\beta$. For numerical implementation,
we seek to have the correct transport for each particle of each species. The
normalized scattering length for particle $i\in\alpha$ interacting with
particle $j\in\beta$ via cumulative small angle Coulomb collisions in an ideal
plasma, and vice versa, are define as

$s_{ij}\equiv\frac{1}{2}\left<\Theta_{ij}^{2}\right>_{\text{all
}b}=\frac{q_{\alpha}^{2}q^{2}_{\beta}\ln\Lambda}{4\pi\epsilon_{0}^{2}\mu^{2}_{\alpha\beta}u_{ij}^{3}}\frac{w_{j}}{\Delta
V}\Delta t,\ \text{and}\
s_{ji}\equiv\frac{1}{2}\left<\Theta_{ji}^{2}\right>_{\text{all
}b}=\frac{q_{\alpha}^{2}q^{2}_{\beta}\ln\Lambda}{4\pi\epsilon_{0}^{2}\mu^{2}_{\alpha\beta}u_{ij}^{3}}\frac{w_{i}}{\Delta
V}\Delta t.$ (9)

If the particles are equally weighted ($w_{i}=w_{j}$ for all $i$ and $j$),
then $s_{ij}=s_{ji}$ for all $i$ and $j$. One can then loop over all binary
pairings (see Fig. 1) and scatter each particle in the pair using the same
common scattering angles. This approach obtains the correct scattering physics
while also ensuring exact conservation of momentum and energy. However, when
the interacting particles have different weights, then $s_{ij}\neq s_{ji}$.
The reason is that macro-particles $i$ and $j$ are interacting with a
different number of physical particles. One could use independent scattering
angles for each particle in the pair. The scattering physics would be correct
for each particle, but this approach requires sampling from two different
Gaussian distribution functions for each binary pair.

We seek a weighted-particle implementation of the $N{\times}N$ method that
only requires sampling from one Gaussian for each binary pair. This is
achieved using the following normalized scattering length for the scattering
pairs:

$s\left(\alpha_{i},\beta_{j}\right)=\frac{q_{\alpha}^{2}q^{2}_{\beta}\ln\Lambda}{4\pi\epsilon_{0}^{2}\mu^{2}_{\alpha\beta}u_{ij}^{3}}\frac{w_{\max}}{\Delta
V}\Delta t,\ \ \ \left(\text{for order $N{\times}N$ pairings}\right)$ (10)

where $w_{\max}=\max(w_{i},w_{j})$ is used for the target density. Eq. 10 is
the correct value as given in Eq. 9 for the lower weight particle, but it is
too large by $w_{\max}/\text{min}(w_{i},w_{j})$ for the larger weight
particle. To compensate for this, the velocity vector for the larger weight
particle is only updated with a probability $w_{\text{min}}/w_{\max}$. For
simplicity, assume that $w_{j}>w_{i}$. Then, after many pairings between
particle $i\in\alpha$ and $j\in\beta$, the average value of the normalized
scattering length for the larger weight particle is

$\left<s\left(\alpha_{i},\beta_{j}\right)\frac{w_{i}}{w_{j}}\right>=\frac{q_{\alpha}^{2}q^{2}_{\beta}\ln\Lambda}{4\pi\epsilon_{0}^{2}\mu^{2}_{\alpha\beta}u_{ij}^{3}}\frac{w_{i}}{\Delta
V}\Delta t=s_{ji},$ (11)

which is the correct value as given in Eq. 9. This approach is referred to as
the rejection method and it is straight-forward to show that, in addition to
achieving the correct scattering physics, momentum and energy are also
preserved on average [6].

###  3.2 Order $N$ binary pairing method

$N{\times}N$ methods are accurate, but the problem is that their computational
cost scales with $N^{2}$. The first particle method for Coulomb collisions
with computational cost of order $N$ is the seminal work by Takizuka and Abe,
which uses $N_{\max}=\max\left(N_{\alpha},N_{\beta}\right)$ pairs per time
step. The pairs are formed randomly and are such that each particle from the
larger number population scatter exactly once per time step. The particles
from the smaller number population scatter $N_{\max}/N_{\text{min}}$ times per
time step on average. The normalized scattering length used in this method to
the determine the scattering angle for each binary pair is computed using the
value for the pair as used for the $N{\times}N$ method given by Eq. 10
multiplied by the number of particles for the population with the least number
of particles: $N_{\text{min}}=\text{min}\left(N_{\alpha},N_{\beta}\right)$.
This is expressed as

$s\left(\alpha_{i},\beta_{j}\right)=\frac{q_{\alpha}^{2}q^{2}_{\beta}\ln\Lambda}{4\pi\epsilon_{0}^{2}\mu^{2}_{\alpha\beta}u_{ij}^{3}}\frac{w_{\max}N_{\text{min}}}{\Delta
V}\Delta t,\ \ \ \left(\text{for order $N$ pairings}\right)$ (12)

To prove that this method gives the correct normalized scattering length on
average for each particle we assume, without loss of generality, that
$N_{\beta}>N_{\alpha}$. First, we consider equal weight particles as is
considered in TA77. On average after many scatterings, the expected value of
the change in velocity of particle $j\in\beta$ (the larger density population)
is

$\left<\Delta\textbf{v}_{j}\right>=-\frac{\mu_{\alpha\beta}}{m_{\beta}}\frac{1}{N_{\alpha}}\sum_{i=1}^{N_{\alpha}}\textbf{u}_{ij}s\left(\alpha_{i},\beta_{j}\right)=-\frac{\mu_{\alpha\beta}}{m_{\beta}}\sum_{i=1}^{N_{\alpha}}\textbf{u}_{ij}s_{ji},$
(13)

which is the correct value with $s_{ji}$ defined in the second expression in
Eq. 9. The factor of $1/N_{\alpha}$ in the second expression is the
probability that particle $j\in\beta$ is paired with particle $i\in\alpha$.
The expected value of the change in velocity of particle $i\in\alpha$ (the
smaller density population) is

$\left<\Delta\textbf{v}_{i}\right>=-\frac{\mu_{\alpha\beta}}{m_{\alpha}}\frac{1}{N_{\beta}}\sum_{j=1}^{N_{\beta}}\textbf{u}_{ij}s\left(\alpha_{i},\beta_{j}\right)\frac{N_{\beta}}{N_{\alpha}}=-\frac{\mu_{\alpha\beta}}{m_{\alpha}}\sum_{j=1}^{N_{\beta}}\textbf{u}_{ij}s_{ij},$
(14)

which is the correct value with $s_{ij}$ defined in the first expression of
Eq. 9. The $N_{\beta}/N_{\alpha}$ factor in the second expression is present
because each particle from species $\alpha$ is paired $N_{\beta}/N_{\alpha}$
times on average for each scattering event (see Fig. 1). Thus, the scattering
physics is correct for all particles on average when using $wN_{\min}/\Delta
V$ for the target density. Since the particles are equally weighted here,
momentum and energy are conserved exactly for each binary pair.

For weighted particles, Eq. 12 is used, but the velocity of the larger weight
particle is only updated with probability $w_{\text{min}}/w_{\max}$, same as
it is for the $N{\times}N$ method. The expected value of the change in
velocity of particle $j\in\beta$ from the larger number population is

$\left<\Delta\textbf{v}_{j}\right>=-\frac{\mu_{\alpha\beta}}{m_{\beta}}\frac{1}{N_{\alpha}}\sum_{i=1}^{N_{\alpha}}\textbf{u}_{ij}\frac{q_{\alpha}^{2}q^{2}_{\beta}\ln\Lambda}{4\pi\epsilon_{0}^{2}\mu^{2}_{\alpha\beta}u_{ij}^{3}}\frac{w_{\max}N_{\alpha}}{\Delta
V}\Delta
tP_{j}\left(\alpha_{i},\beta_{j}\right)=-\frac{\mu_{\alpha\beta}}{m_{\beta}}\sum_{i=1}^{N_{\alpha}}\textbf{u}_{ij}s_{ji}.$
(15)

For the last step, we made use of the fact that
$P_{j}\left(\alpha_{i},\beta_{j}\right)w_{\max}=w_{i}$ for all cases, where
$P_{j}\left(\alpha_{i},\beta_{j}\right)$ is the probability that particle
$j\in\beta$ scatters when paired with particle $i\in\alpha$. The expected
value of the change in velocity of particle $i\in\alpha$ from the smaller
number population is

$\left<\Delta\textbf{v}_{i}\right>=-\frac{\mu_{\alpha\beta}}{m_{\alpha}}\frac{1}{N_{\beta}}\sum_{j=1}^{N_{\beta}}\textbf{u}_{ij}\frac{q_{\alpha}^{2}q^{2}_{\beta}\ln\Lambda}{4\pi\epsilon_{0}^{2}\mu^{2}_{\alpha\beta}u_{ij}^{3}}\frac{w_{\max}N_{\alpha}}{\Delta
V}\Delta
tP_{j}\left(\alpha_{i},\beta_{j}\right)\frac{N_{\beta}}{N_{\alpha}}=-\frac{\mu_{\alpha\beta}}{m_{\alpha}}\sum_{j=1}^{N_{\beta}}\textbf{u}_{ij}s_{ij},$
(16)

where again for the last step use is made of the fact that
$P_{i}\left(\alpha_{i},\beta_{j}\right)w_{\max}=w_{j}$ for all cases, where
$P_{i}\left(\alpha_{i},\beta_{j}\right)$ is the probability that particle
$i\in\alpha$ scatters when paired with particle $j\in\beta$.

In summary, the weighted-particle version of TA77 presented here is the same
as the equal weight particle version with two changes. First, $w_{\max}$ is
used for the density in the normalized scattering length for each binary pair.
Second, the velocity of the larger weight particle is updated with probability
$w_{\min}/w_{\max}$. The pairing method used here (see Fig. 1), which is the
same as in TA77, can also be used for the N97 method. The only difference
between these two methods is the formula that takes the normalized scattering
length and returns the scattering angle. These methods asymptote to the same
variance as $\Delta t$ goes to zero [7].

Everything here is described for inter-species scattering. For intra-species
scattering of a species with $N_{\alpha}$ particles in a cell, the pairing is
done as described by TA77 and for weighted-particles one uses Eq. 12 with
$N_{\text{min}}=N_{\alpha}-1$. When $N_{\alpha}$ is odd, the first three
particles are paired in an $N{\times}N$ way and the normalized scattering
length is reduced by a factor of two for these pairings to maintain the
correct scattering length.

A slightly different method for extending the methods of TA77 and N97 to
include weighted particles while preserving the scattering physics on average
is described in Ref. [7]. In that work, the normalized scattering length is
computed using Eq. 12, but with $N_{\min}$ replaced by $N_{\max}$, and
$w_{\max}$ replaced by $w_{\max}/\mathcal{D}_{ij}$, where $\mathcal{D}_{ij}$
is duplicity factor for the particle pair defined by the number of times the
particle from the lower number population is paired during a scattering event.
Since the average duplicity factor per binary pair is
$\left<\mathcal{D}_{ij}\right>=N_{\max}/N_{\min}$, the method described in
this work is like that from Ref. [7] when using the average duplicity factor
for each particle pair.

##  4 Moment correction method

The Monte Carlo method for weighted-particle Coulomb collisions described in
the previous section only preserves momentum and energy on average after many
pairings. Several methods have been used before to restore exact momentum
and/or energy conservation after application of non-conservative scattering
algorithm, such as the models by Jones [14] and Lemons [15]. In the method by
Jones, the post-scattered velocity of each particle is linearly transformed by
shifting and rescaling the post-scattered velocity of each particle The method
by Lemons is similar, but the linear transformation is applied to the change
in velocity from scattering rather than to the velocity itself. In these
methods, each particle is shifted and rescaled by the same common values,
which are readily computed using the laws of conservation of momentum and
energy for non-relativistic particles. For relativistic particles, the shift
parameter is still readily computed, but a different formula for rescaling the
velocities needs to be used [16]. The formula used in Ref. [16] is for like-
species collisions and requires the energy-correction factor to be of the form
$1+\alpha$ with $|\alpha|\ll 1$. The parameter $\alpha$ is obtained by
expanding the kinetic energy to first order in $\alpha$ and is therefore not
exact. A weighted-particle scattering method is used by Sentoko and Kemp in
Ref. [17] that preserves energy exactly and works for relativistic particles,
but momentum is only preserved on average.

Here, a new method is proposed for restoring momentum and energy conservation.
Momentum is restored by applying a linear shift to the particles velocities
like that described above. The method to restore energy conservation is based
on kinematics for zero-angle inelastic scattering of binary pairs. Being based
on scattering kinematics, this method works for both non-relativistic and
relativistic particles. The method works as follows. First, the weighted-
particle scattering algorithm is applied to two interacting species. The
particle velocity before and post-scatter are denoted, respectively, as,
$\textbf{v}_{p}^{bs}$ and $\textbf{v}_{p}^{ps}$. Second, momentum conservation
is restored by shifting the particle velocities as

$\textbf{v}_{p}^{pm}=\textbf{v}_{p}^{ps}-\textbf{B}w_{p}\ \Rightarrow\
\textbf{B}=\frac{\sum_{p}m_{p}\left(\textbf{v}^{ps}_{p}-\textbf{v}_{p}^{bs}\right)}{\sum_{p}m_{p}w_{p}},$
(17)

where $m_{p}=m_{s}w_{p}$ is the mass of particle $p\in s$ and
$\textbf{v}_{p}^{pm}$ is the post-momentum corrected velocity of particle $p$.
For relativistic particles, one shifts the proper velocity
$\textbf{u}_{p}=\gamma_{p}\textbf{v}_{p}$ with
$\gamma_{p}=1/\sqrt{1-|\textbf{v}_{p}|^{2}/c^{2}}$. The difference between the
momentum correction used here and that used in Refs. [14, 15] is that here the
shift is biased to the larger weighted particles. After momentum conservation
has been restored, the total violation in energy conservation in a cell is
computed as

$\Delta E=\sum_{p}m_{p}\frac{1}{2}|\textbf{v}^{pm}_{p}|^{2}-KE^{bs},$ (18)

where $KE^{bs}=\sum_{p}m_{p}|\textbf{v}_{p}^{bs}|^{2}/2$ is the total kinetic
energy in a cell before scattering. The final step is to adjust the energy of
the particles such that $\Delta E$ is re-absorbed. As mentioned above, the
method we use is based on zero-angle inelastic scattering kinematics of binary
pairs. It is instructive to first go over this process. It is done here for
non-relativistic particles. An extension to relativistic particles is given in
Appendix A.

The post-collision velocity of two particles undergoing scattering in the CM
frame can be written as

$\displaystyle\textbf{v}^{\prime}_{1}$ $\displaystyle=$
$\displaystyle\textbf{v}_{1}+\frac{\mu_{R}}{m_{1}}\left(\textbf{v}^{\prime}_{R}-\textbf{v}_{R}\right),$
(19) $\displaystyle\textbf{v}^{\prime}_{2}$ $\displaystyle=$
$\displaystyle\textbf{v}_{2}-\frac{\mu_{R}}{m_{2}}\left(\textbf{v}^{\prime}_{R}-\textbf{v}_{R}\right),$
(20)

where $\mu_{R}=m_{1}m_{2}/(m_{1}+m_{2})$ is the reduced mass,
$\textbf{v}_{R}=\textbf{v}_{1}-\textbf{v}_{2}$ is the relative velocity, and
the prime superscript denotes post scatter. The energy conservation law for an
inelastic event with change in energy $U$ is

$\frac{1}{2}\mu_{R}|\textbf{v}^{\prime}_{R}|^{2}=\frac{1}{2}\mu_{R}|\textbf{v}_{R}|^{2}-U.$
(21)

For zero-angle scattering, the post-scatter relative velocity vector is
obtained directly from Eq. 21 and is given as

$\textbf{v}^{\prime}_{R}=\textbf{v}_{R}\sqrt{1-\frac{U}{0.5\mu_{R}|\textbf{v}_{R}|^{2}}}.$
(22)

The weighted-scattering methods described in Sec. 3 are such that the correct
scattering physics is preserved on average. When adjusting the distribution
function post-scatter to restore momentum and energy conservation, one should
minimize how much the distribution function is perturbed to maintain this
property. Here, the energy is corrected by forming binary pairs and adjusting
the CM energy using Eqs. 19-20. To try and minimize the perturbation, the
potential $U$ used for a given binary pair is chosen to be some small
percentage of the CM energy of the binary pair, that is $U=\text{sign}(\Delta
E)f_{E}0.5\mu_{R}|\textbf{v}_{p}|^{2}$ where $f_{E}$ is a user-specified small
parameter referred to as the energy-correction factor. If $\Delta E<0$, then
energy needs to be added back to the particles to restore energy conservation
and vice versa. The binary pairs are formed randomly. After each binary pair,
$\Delta E$ is adjust accordingly and the loop proceeds until $\Delta E=0$. To
restore energy conservation exactly, $|U|$ is chosen as:

$U=\text{sign}\left(\Delta E\right)\min\left(|\Delta
E|,f_{E}0.5\mu_{R}|\textbf{v}_{p}|^{2}\right).$ (23)

The free parameter $f_{E}$ in this model is the percentage that one allows the
CM energy of the binary pairs to change. If $f_{E}$ is too large, then the
distribution function may be perturbed too much, and the scattering physics
can be disrupted. However, more binary pairs are needed to restore energy
conservation for smaller values of $f_{E}$, which means the computational cost
of the algorithm increases. Typical values of $f_{E}$ used here are between
$0.02$ and $0.05$.

Another important things to note about the application of Eqs. 19, 20, and 22
used in this context is that the masses used are the weighted particle ones
($m_{1}=m_{s}w_{1}$ and $m_{2}=m_{s}w_{2}$). For this reason, it can be useful
to sort the particles by weight before forming the binary pairs such that the
particles are more likely to be paired with a similar weight particle.
Furthermore, this ensures that the energy-correction is biased to the larger
weight particles. This is preferable because the error in momentum and energy
conservation associated with each application of the weighted-scattering
algorithm reside in the larger weight particles, which only scatter
probabilistically. Another way to say this is that it is preferable to
minimize adjusting the velocity of the lower weight particles, which are
typically scattered correctly at each application of the weighted-scattering
algorithm.

The full weighted-scattering algorithm including steps to restore momentum and
energy conservation is outlined in Algorithm 1. This algorithm works for both
order $N$ TA-like pairings and for the $N{\times}N$ method. The only
difference is in the formula used for the normalized scattering length for the
binary pairs. Algorithm 1 is also written to be valid for intra-species and
inter-species collisions. For inter-species scattering, one needs to choose
how to partition the net violation in energy conservation between the two
species. In this work, $\Delta E$ is partitioned as follows

$\Delta
E_{s}=\frac{\bar{w}_{s}E_{s}^{pm}}{\bar{w}_{\alpha}E_{\alpha}^{pm}+\bar{w}_{\beta}E_{\beta}^{pm}}\Delta
E,\ \text{for}\ s=\alpha,\beta.$ (24)

In Eq. 24, $\bar{w}_{s}=\sum_{p\in s}w_{p}/N_{s}$ is the average particle
weight of species $s$ in the cell and $E_{s}^{pm}$ is the total energy in the
cell of species $s$ post momentum correction. The weighting of $\Delta E_{s}$
with $\bar{w}_{s}$ is for the purpose of biasing the energy correction to
larger weight particles.

Algorithm 1 Moment-preserving weighted Coulomb collision method.

1:Compute the initial energy and momentum of each species in the cell.

2:Formulate random binary pairs as described in Fig. 1.

3:Loop over pairs and for each pair:

4: Compute the scattering angle using $s(\alpha_{i},\beta_{j})$ given in Eq.
12 for TA-like or Eq. 10 for $N{\times}N$.

5: Use the formula from TA77 or N97 to compute the CM scattering angle from
$s(\alpha_{i},\beta_{j})$

6: Update $\textbf{v}_{p}$ for the smaller-weight particle of the pair.

7: Update $\textbf{v}_{p}$ for the larger-weight particle of the pair with
probability $w_{\min}/w_{\max}$.

8:Restore momentum conservation by shifting the velocity of each particle
using Eq. 17.

9:Compute the net violation in energy conservation $\Delta E$ using Eq. 18.

10:For inter-species scattering, partition $\Delta E$ using Eq. 24.

11:(Optional step) Sort particles of each species in descending order by
weight.

12:For each species $s$, loop over particles forming binary pairs and for each
pair:

13: Compute $\textbf{v}^{\prime}_{R}$ in Eq. 22 using $U$ from Eq. 23.

14: Modify $\textbf{v}_{p}$ of each particle using Eqs. 19-20.

15: Adjust net energy violation: $\Delta E_{s}\mathrel{-}=U$.

16: Break loop when $\Delta E_{s}=0$.

##  5 Numerical tests

Results from several numerical tests using the moment-preserving weighted-
particle Monte-Carlo scattering algorithm are presented in this section. The
first test considered is Test 1 from Ref. [7], which considers two populations
within the same species relaxing to a common velocity and temperature. Same as
in Ref. [7], simulations of this test are performed using a variety of
different weight ratios for the two populations to rigorously verify the
accuracy of the weighted-particle scattering method. The second test is Test 2
from Ref. [7], which is like Test 1, but includes three populations belonging
to two different species to test inter-species scattering. The third test is
an electron-ion thermalization test for a fully ionized carbon plasma. Unless
stated otherwise, the particles are sorted in descending order of particle
weights prior to adjusting the velocities to restore energy conservation.

###  5.1 Test 1

Consider two populations, call them $A$ and $B$, of fully ionized carbon 12
atoms ($A=12$, $Z=6$) with $m_{i}=Am_{u}-Zm_{e}$, where $m_{u}$ is the atomic
mass unit and $m_{e}$ is the electron mass. The initial density, temperature,
and velocity of population $A$ are $n_{A}=10^{19}/\text{cm}^{3}$, $T_{A}=500$
eV, and $U_{A}=+655$ km/s. Those for population $B$ are
$n_{B}=10^{20}/\text{cm}^{3}$, $T_{B}=500$ eV, and $U_{A}=0$ km/s. Momentum
conservation gives $U_{f}=U_{A}n_{A}/\left(n_{A}+n_{B}\right)=59.55$ for the
equilibrium value for the velocity of each species. Energy conservation gives
$T_{f}=1.97$ keV for each species. The simulations are performed in quasi 0D
using a 1D grid with 180 grid cells. The number of particles in Table 1 below
summarizing the numerical setup for this test refer to the number of particles
per cell. A fixed value of $\Lambda=10$ is used for the Coulomb logarithm. The
particle positions are not advances and there are no fields involved. The
simulations are velocity-space only and the results presented below are
averaged values on the grid.

Table 1: Test 1 simulation parameters. |  | $N_{A}$ | $N_{B}$ | $N_{A}{:}N_{B}$ | $w_{A}{:}w_{B}$ |  $\Delta t$ [fs] | $f_{E}$
---|---|---|---|---|---|---|---
T1a |  | 400 | 4,000 | 1:10 | 1:1 | 50 | 0.05
T1b |  | 400 | 400 | 1:1 | 1:10 | 50 | 0.05
T1c |  | 4,000 | 400 | 10:1 | 1:100 | 15 | 0.02
T1d |  | 200 | 8,000 | 1:40 | 4:1 | 50 | 0.05

The numerical parameters for tests T1a-T1d are summarized in Table. 1. One
thing to note is that the number of particles per cell used here is 10 times
less than that considered in Ref. [7]. Results for tests T1a (equal particle
weights) and tests T1b (equal number of particles) are shown in Fig. 2. The
results from test T1a are considered the baseline solution. Results from test
T1b shown in this figure, where the ratio of particles weights for the low-
density population to the high-density populations is 1:10, are
indistinguishable from test T1a. The relaxation of the mean velocity and
temperature profiles shown in Fig. 2 can additionally be compared with the
same results presented in Fig. 3 of Ref. [7].

Fig. 2: Mean velocity and temperature profiles for populations $A$ and $B$ for
tests T1a and T1b. The solid black and solid blue curves are for populations
$A$ and $B$, respectively, from test T1a (equally weighted particles). The red
$\times$ markers and the dashed yellow curves are from test T1b, where the
particle weights of the denser population B are 10$\times$ larger than that
for the lower density population A.

The violation in total x-momentum and energy for test T1b with and without the
moment-correction algorithm applied are shown in Fig. 3. Although not shown,
the velocity and temperature relaxation results shown in Fig. 2 for test T1b
are similar independent of using the moment-correction algorithm. The relative
error in momentum and energy conservation is between $0.5-1.5$% when not
correcting the moments and using unequal weight particles. Momentum and energy
conservation are preserved to machine precision with weighted particles when
applying the moment-correction method, as seen in the right panel of Fig. 3.

Fig. 3: Conservation of x-momentum (left panel) and energy (middle and right
panels) for test T1b with the moment correction method applied (dashed red
curves) and without (yellow curves). Results are compared with those from test
T1a, which uses equal particle weights and thus naturally conserves momentum
and energy. The right panel shows the relative change in total energy for
tests T1a and T1b using the moment correction method on a scale where it is
evident that the violation in energy conservation is at machine precision.

Results from all four test cases described in Table. 1 are shown in Fig. 4.
The relaxation of the velocity and temperature profiles for both species match
well for all cases. These results show 1) that the weighted scattering
algorithm achieves the correct scattering physics for a wide range of
weighted-particle scenarios, and 2) the moment-correction method can be
applied without significantly altering the scattering physics. It is worth
noting here that the original method for doing weighted-particle Coulomb
scattering by Nanbu in Ref. [6] does not achieve the correct relaxation for
any of the unequal weight particles tests for this problem (see Fig. 6 of Ref.
[7]).

Fig. 4: Mean velocity and temperature profiles of populations A and B for
tests T1a - T1d. See Table 1 for simulation parameters. The dashed-dotted
black curves are the equilibrium solutions.

Time profiles for the number of binary pairs required to restore energy
conservation for test T1d are shown in Fig. 5 for various values of the
energy-correction factor $f_{E}$. Results are shown in the middle panel of
this figure when using $10\times$ more particles. The number of binary pairs
needed to restore energy conservation scales approximately linear with $f_{E}$
and approximately as $\sqrt{N}$. The effects of $f_{E}$ and $N$ on $T_{A}$
from test T1d are shown in Fig. 6. The zoomed in inserts in this figure show
the solution converging to the equally weighted particle solution as
$f_{E}\rightarrow 0$ and with increasing $N$.

Fig. 5: Time profile of the number of binary pairs needed to restore energy
conservation for test T1d. The left and middle panels show results for various
values of the energy-correction factor $f_{E}$. The results in the middle
panel are obtained using $10\times$ more particles than used for the left
panel. The right panel shows how these results differ for this test if one
chooses not to sort the particles by weight prior to applying the energy
correction. Fig. 6: Long-time profiles for the temperature of the low-density
population A from T1d simulations with various values of the energy-correction
factor. The results in the right two panels are from simulations using
$10\times$ the number of particles compared to those for the left two panels.

###  5.2 Test 2

The second test considered is similar to Test 1, but populations A and B are
now separated into different species, species $\alpha$ and species $\beta$,
and a third population C is added to species $\beta$. Population C is a low-
density population with the following initial conditions:
$n_{C}=10^{18}/\text{cm}^{3}$, $T_{C}=500$ eV, and $U_{A}=-655$ km/s. The
simulation parameters for Test 2 are summarized in Table 2. A variety of
different weight combinations for the different populations are considered to
rigorously verify the weighted-particle scattering algorithm. Tests T2a-T2e
are the same as those considered in Ref. [7], but with $10\times$ less
particles, same as for Test 1. Sorting is used for $\alpha-\beta$ and
$\beta-\beta$ scattering but is not used for $\alpha-\alpha$ scattering since
each of the particles in species $\alpha$ have equal weights.

Table 2: Test 2 simulation parameters. |  | $N_{A}$ | $N_{B}$ | $N_{C}$ | $N_{A}{:}N_{B}{:}N_{C}$ | $N_{\alpha}{:}N_{\beta}$ | $w_{A}{:}w_{B}{:}w_{c}$ |  $\Delta t$ [fs] | $f_{E}$
---|---|---|---|---|---|---|---|---|---
T2a |  | 400 | 4,000 | 200 | 2:20:1 | 1:10.5 | 1:1:1 | 50 | 0.05
T2b |  | 400 | 4,000 | 800 | 1:10:2 | 1:12 | 4:4:1 | 50 | 0.05
T2c |  | 400 | 800 | 4,000 | 1:2:10 | 1:12 | 20:100:1 | 50 | 0.05
T2d |  | 4,000 | 500 | 100 | 40:5:1 | 40:6 | 1:80:20 | 50 | 0.05
T2e |  | 4,000 | 100 | 500 | 40:1:5 | 40:6 | 1:400:4 | 50 | 0.05
Fig. 7: Mean velocity and temperature profiles for populations A, B, and C for
tests T2a - T1e. See Table 2 for simulation parameters. The dashed-dotted
black curves are the equilibrium solutions.

Time-profiles of the mean velocity and temperature of populations A, B, and C
for tests T2a-T2e are shown in Fig. 7. Each of the quantities for each of the
populations shown in Fig. 7 are observed to relax to the steady state
solutions. Furthermore, the relaxation to the steady state values is observed
to agree well for all particle weighting combinations considered. There is a
small, but noticeable deviation of the $U_{x}$ profile for population $B$ for
test T2e, which is an extreme case where $w_{A}{:}w_{B}{:}w_{C}=1{:}400{:}4$.
Furthermore, there are only 100 particles per cell for population $B$ in this
test. This small discrepancy goes away with increasing the number of particles
per cell.

Fig. 8: Conservation of x-momentum (left panel) and energy (middle and right
panels) for test T2a-T2e.

The ability to maintain rigorous conservation of momentum and energy for Tests
T2a-T2e is shown in Fig. 8. The relative change in total x-momentum and energy
is shown in this figure to be on the order of machine-precision for each of
these tests.

###  5.3 Test 3

For the final test, electron-ion thermalization of a fully ionized carbon
plasma is considered ($A=12$, $Z=6$). The electron and ion distributions are
initialized as non-drifting Maxwellians with the following conditions:
$n_{i}=1{\times}10^{23}/$cm3, $T_{i}=50$ eV, $n_{e}=6{\times}10^{23}/$cm3, and
$T_{e}=150$ eV. The equilibrium solution for the temperature is $T_{f}=136$
eV. The simulations are performed using a 2D periodic domain on a
$40{\times}40$ grid with uniform spacing equal to $1.3285$ nm. For all
simulations, $\Delta t=3.54{\times}10^{-4}$ fs is used for the time step and
$\ln\Lambda=3$ is used as the Coulomb logarithm. This time step is roughly
$100$ times smaller than the characteristic electron collision time:
$\tau_{e}[\text{s}]=3.44{\times}10^{5}T_{e}^{3/2}[\text{eV}]/\left(Zn_{e}[1/\text{cm}^{3}]\Lambda\right)$.

Simulation results for this test problem are shown in Fig. 9. Here, time
profiles for the average electron and ion temperatures on the 2D grid are
shown for three cases. The number of ion particles per cell is set to
$N_{i}=512$ for each of these simulations. The results in the left panel are
obtained using equally weighted particles (T3a), in which case the number of
electron particles per cell is $N_{e}=6N_{i}=3072$. The simulations
corresponding to the middle and right panels use $N_{e}=N_{i}=512$, in which
chase $w_{e}=6w_{i}$. The results in the middle and right panels are obtained
without (T3b) and with (T3c) the moment-correction algorithm, respectively.
The electron and ion temperature profiles from each of these simulations are
observed to match well. Furthermore, the relaxation of $T_{e}$ and $T_{i}$
shown in Fig. 9 are also observed to agree well with results obtained from a
0D Spitzer model:

$\frac{dT_{e}}{dt}=-\nu_{T}\left(T_{e}-T_{i}\right);\
\frac{dT_{i}}{dt}=-Z\frac{dT_{e}}{dt},$ (25)

where $\nu_{T}=2m_{e}/M_{i}/\tau_{e}$ is the thermalization rate.

Fig. 9: Species temperature relaxation results for a fully ionized carbon
plasma. In each figure, the dashed black curve is the equilibrium solution,
and the dotted curves are obtained from a 0D Spitzer thermalization model as
described in the text. Equally weighted particles are used for the results
shown in the left panel. The results shown in the middle and right panels are
from simulations that use an equal number of simulation particles for both
species, without (middle panel) and with (right panel) the moment correction
method applied.

The ability to use weighted particles while conserving momentum and energy for
test T3c is shown in Fig. 10. The left two panels in this figure show the
change in total momentum in each direction. The amplitude of the change in
momentum for test T3c shown in the middle panel is about 2-3 times larger than
that for the simulation with equal weights shown in the left panel. In both
cases, the change in momentum is close to machine precision and behaving like
a random walk. Similarly, the relative change in energy for test T3c also
behaves like a random walk, but with a step size that is about $10{\times}$
larger than for test T1a where equal weight particles are used. The reason for
this is unclear. It may be due to the large electron-ion mass ratio. In either
case, energy is conserved nearly to machine-precision for weighted-particles
with the moment correction algorithm applied.

Fig. 10: Change in total momentum (left and middle panels) and relative change
in total energy (right panel) correspondig to test T3a and T3c shown in Fig.
9.

##  6 Discussion and summary

###  6.1 Analysis of perturbations to distribution functions

It is worthwhile to analyze how the energy-correction scheme described above
perturbs the species distribution functions. For simplicity, assume $\Delta E$
after momentum is restored is negative. Using Eq. 23 in Eq. 22, we can write

$\textbf{v}^{\prime}_{R}-\textbf{v}_{R}=\textbf{v}_{R}\sqrt{1+f_{E}}-\textbf{v}_{R}\approx\frac{f_{E}}{2}\textbf{v}_{R}.$
(26)

It is additionally assumed that the particles are sorted by weight such that
$m_{1}\approx m_{2}$. Then, using Eq. 26, the post-energy-correction particle
velocities in Eqs. 19-20 become

$\displaystyle\textbf{v}^{\prime}_{1}$ $\displaystyle=$
$\displaystyle\textbf{v}_{1}+\frac{f_{E}}{4}\left(\textbf{v}_{1}-\textbf{v}_{2}\right),$
(27) $\displaystyle\textbf{v}^{\prime}_{2}$ $\displaystyle=$
$\displaystyle\textbf{v}_{2}-\frac{f_{E}}{4}\left(\textbf{v}_{1}-\textbf{v}_{2}\right).$
(28)

We analyze Eqs. 27-28 in several different scenarios. First, assume the
velocity of the two particles have similar magnitudes, but are in opposite
directions. In this case, the velocity of both particles is scaled by
$1+f_{E}/2$. If the velocity of the two particles is similar in magnitude and
are in the same direction, then there is essential no modification to the
particle velocities because the CM energy is approximately zero. Another case
to consider is where the particle speeds are much different, e.g.,
$|\textbf{v}_{1}|\ll|\textbf{v}_{2}|$. In this scenario, the velocity of the
fast particle is modified by $1+f_{E}/4$, while the speed of the slower
particle gets enhanced to a value as large as $|\textbf{v}_{2}|f_{E}/4$. In
other words, for this scenario the fast particle is only perturbed slightly
while the slow particle can be perturbed significantly. We comment that more
strongly perturbing slower particles is preferred as perturbations to the
lower velocity part of the distribution function typically has less of an
effect on the physics.

###  6.2 Discussion of potential failures of the energy-correction algorithm

The energy-correction method works by adjusting the weighted-CM energy of
random binary pairs after the momentum correction has been applied. If the
violation of energy conservation, $\Delta E$ is large in magnitude compared to
the combined weighted-CM energy of all binary pairs, then the method can fail.
We find that this can occur when there are very few particles in a cell (e.g.,
$N$ order 10), and at least one of the particles has a weight that is much
larger than the others. When the large weight particle is determined to
scatter, it can cause a relatively large violation in energy conservation in
the cell that cannot be re-absorbed by adjusting the weighted-CM energy of the
very few binary pairs in the cell. We have only encountered this problem with
the number of particles in cell being order $10$. For such a low number of
particles in a cell, one cannot expect physically accurate results for
scattering with weighted particles. The simplest solution is to do nothing and
live with a momentum-conserving but not energy-conserving solution. The next
simplest solution is to save the momentum vectors of the particles prior to
scattering and, if the method fails, reset the particles velocities and skip
doing scattering in this cell during this time step.

###  6.3 Computational cost of momentum-correction algorithm

The moment-correction algorithm increases the computational cost of the
overall scattering algorithm. For test T1d, but with $10{\times}$ more
particles ($N=82,000$), the relative computational cost for no moment-
correction, moment-correction without sorting, and moment-correction with
sorting, are is 1.00, 1.23, and 1.72, respectively. The computational cost of
applying the moment-correction algorithm itself increases the cost by $23$%.
The cost increases by an additional $40$% when sorting the particles. We are
currently using the standard std::sort() function, which may not always have
the ideal computational cost of order $N\log\left(N\right)$. In the future, we
will seek out more efficient sorting algorithms to reduce the computational
cost. Furthermore, it is not always necessary to sort all of the particles
(see Fig. 5). The computational burden of sorting can also be greatly reduced
by only doing partial sorting.

###  6.4 Summary

A method for doing weighted-particle Coulomb collisions in a plasma that
preserves scattering physics while exactly conserving momentum and energy is
described in this manuscript. We first show how to extend $N{\times}N$
scattering methods for equally weighted particles to arbitrarily-weighted
particles such that the normalized scattering length is correct for each
particle on average. Then, it is straightforward to deduce an expression for
the scattering length of binary pairs for order $N$ methods that is a direct
generalization of that given by Takizuka and Abe [2] to include weighted
particles.

The weighted-particle scattering method alone only preserves momentum and
energy on average. A new method is introduced to rigorously restore exact
momentum and energy conservation post scatter. The method, which is based on
inelastic scattering dynamics, works by forming random binary pairs and making
small adjustments to the weighted CM energy of each pair. The method works for
both non-relativistic and relativistic particles. The efficacy of the entire
algorithm is illustrated with various test problems, including those in which
the original Nanbu method for weighted-particle Coulomb scattering [6] is
known to produce incorrect results [7].

While this work focused on small angle Coulomb collisions, the extension to
weighted particles and the method to restore momentum and energy conservation
are not specific to Coulomb collisions in any way. These methods work equally
well for other collision processes too, such as those governed directly by the
Boltzmann collision integral discussed at the beginning of Sec. 2 where the
normalized path length is used as a probability of a collision occurring
rather than as a variance for a Gaussian distribution function from which a
scattering angle is sampled from.

##  Appendix A Extension of moment-correction method to relativistic particles

Restoring momentum conservation for relativistic particles is done the same as
it is for non-relativistic particles, The only difference is that the proper
velocity ($\textbf{u}_{p}=\gamma_{p}\textbf{v}_{p}$) is shifted rather than
the actual velocity. The method for correcting for energy conservation, on the
other hand, must be done a bit differently. For relativistic particles, a
Lorentz transformation of the four-momentum vector must be done to go between
the lab and center-of-momentum frames (relativistic analog of the center-of-
mass frame). The CM velocity for a relativistic treatment is defined as

$\textbf{v}_{\text{CM}}=\frac{\textbf{p}_{1}+\textbf{p}_{2}}{\gamma_{1}m_{1}+\gamma_{2}m_{2}}=\frac{\textbf{p}_{1}+\textbf{p}_{2}}{E_{tot}/c^{2}}.$
(29)

One could attempt the same procedure as before; transform to the weighted-CM
frame, adjust the magnitude of the particle velocities to account for an
inelastic energy sink/source, then transform back to the lab frame. However,
there is a subtlety with respect to inelastic processes in relativistic
mechanics that needs to be mentioned. The potential energy of an inelastic
process in relativistic mechanics is realized via a change in particle mass.
This is well known for fusion processes, but it is the same for any inelastic
process. This means that the mass of the particles change in an inelastic
process such that the total energy,
$E_{tot}=\left(\gamma_{1}m_{1}+\gamma_{2}m_{2}\right)c^{2}$, does not change.
In other words, if we transform to the CM frame, alter the kinetic energy but
not the mass of the particles, and then transform back to the lab frame, we
will find that 1) momentum is not conserved and 2) the change in kinetic
energy does not match that in the CM frame. The reason is that the CM frame is
altered since $E_{tot}$ in Eq. 29 has changed.

To perform the energy-correction method for relativistic particles we write
the change in total energy and momentum as follows

$\displaystyle E^{\prime}_{tot}$ $\displaystyle=$ $\displaystyle E_{tot}-U,$
(30) $\displaystyle\textbf{p}^{\prime}_{1}$ $\displaystyle=$
$\displaystyle\textbf{p}_{1}+\alpha\Delta\textbf{u},$ (31)
$\displaystyle\textbf{p}^{\prime}_{2}$ $\displaystyle=$
$\displaystyle\textbf{p}_{2}-\alpha\Delta\textbf{u},$ (32)

where $\Delta\textbf{u}\equiv\textbf{u}_{1}-\textbf{u}_{2}$ and
$U\equiv\text{sign}\left(\Delta E\right)f_{E}KE_{cm}$ is defined the same as
before with $KE_{cm}$ the kinetic energy in the weighted-CM frame. From
Lorentz invariance of the product of 4-vectors, the total energy in the CM
frame is $E^{2}_{cm}=E^{2}_{tot}-|\textbf{p}_{tot}|^{2}c^{2}$ where
$\textbf{p}_{tot}=\textbf{p}_{1}+\textbf{p}_{2}$. The kinetic energy in the CM
frame is $KE_{cm}=E_{cm}-\left(m_{1}+m_{2}\right)c^{2}$. The scalar quantity
$\alpha$ in Eqs. 31-32 is computed using the following relativistic
relationship between energy and momentum:

$\displaystyle\frac{E_{tot}-U}{c^{2}}=m_{1}\gamma^{\prime}_{1}+m_{2}\gamma^{\prime}_{2}=\sqrt{m_{1}^{2}+|\textbf{p}^{\prime}_{1}|^{2}/c^{2}}+\sqrt{m_{2}^{2}+|\textbf{p}^{\prime}_{2}|^{2}/c^{2}}.$
(33)

After some algebraic manipulation, Eq. 33 can be written as the following
quadratic equation for $\alpha$:

$\frac{4}{c^{4}}\left[A^{2}|\Delta\textbf{u}|^{2}c^{2}-\left(\textbf{p}_{tot}\cdot\Delta\textbf{u}\right)^{2}\right]\alpha^{2}+\frac{4}{c^{2}}\Delta\textbf{u}\cdot\left[D\textbf{p}_{tot}-2A^{2}\textbf{p}_{2}\right]\alpha+\frac{4}{c^{4}}A^{2}E_{2}^{2}-D^{2}=0,$
(34)

where $A\equiv\left(E_{tot}-U\right)/c^{2}$ and $D\equiv
A^{2}+\left(E_{2}^{2}-E_{1}^{2}\right)/c^{4}$. The solution of Eq. 34 has two
roots. The positive root is the smaller of the two and is the solution that
converges to the same expression used for the non-relativistic treatment in
the limit where the particle velocities are small compared to the speed of
light.

## Acknowledgments

This work was performed under the auspices of the U.S. Department of Energy by
Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344 and
was supported by the LLNL-LDRD Program under Project No. 23-ERD-007.

## References

  * Bellan [2008] P. M. Bellan, Fundamentals of plasma physics, Cambridge university press, 2008\.
  * Takizuka and Abe [1977] T. Takizuka, H. Abe,  A binary collision model for plasma simulation with a particle code,  Journal of Computational Physics 25 (1977) 205 – 219.
  * Nanbu [1997] K. Nanbu,  Theory of cumulative small-angle collisions in plasmas,  Phys. Rev. E 55 (1997) 4642–4652.
  * Wang et al. [2008] C. Wang, T. Lin, R. Caflisch, B. I. Cohen, A. M. Dimits,  Particle simulation of coulomb collisions: Comparing the methods of takizuka and abe and nanbu,  Journal of Computational Physics 227 (2008) 4308 – 4329.
  * Miller and Combi [1994] R. H. Miller, M. R. Combi,  A coulomb collision algorithm for weighted particle simulations,  Geophysical Research Letters 21 (1994) 1735–1738.
  * Nanbu and Yonemura [1998] K. Nanbu, S. Yonemura,  Weighted particles in coulomb collision simulations based on the theory of a cumulative scattering angle,  Journal of Computational Physics 145 (1998) 639 – 654.
  * Higginson et al. [2020] D. P. Higginson, I. Holod, A. Link,  A corrected method for coulomb scattering in arbitrarily weighted particle-in-cell plasma simulations,  Journal of Computational Physics 413 (2020) 109450.
  * Chen et al. [2011] G. Chen, L. Chacón, D. Barnes,  An energy- and charge-conserving, implicit, electrostatic particle-in-cell algorithm,  Journal of Computational Physics 230 (2011) 7018–7036.
  * Markidis and Lapenta [2011] S. Markidis, G. Lapenta,  The energy conserving particle-in-cell method,  Journal of Computational Physics 230 (2011) 7037 – 7052.
  * Angus et al. [2022] J. R. Angus, A. Link, A. Friedman, D. Ghosh, J. D. Johnson,  On numerical energy conservation for an implicit particle-in-cell method coupled with a binary monte-carlo algorithm for coulomb collisions,  Journal of Computational Physics 456 (2022) 111030.
  * Angus et al. [2023] J. R. Angus, W. Farmer, A. Friedman, D. Ghosh, D. Grote, D. Larson, A. Link,  An implicit particle code with exact energy and charge conservation for electromagnetic studies of dense plasmas,  Journal of Computational Physics 491 (2023) 112383.
  * Lieberman and Lichtenberg [2005] M. A. Lieberman, A. J. Lichtenberg, Principles of Plasma Discharges and Material Processing, 245, 2nd ed., John Wiley and Sons, Inc., 2005.
  * Schmidt [1979] G. Schmidt, Physics of High Temperature Plasmas, 2nd ed., Academic Press, New York, San Francisco, and London, 1979.
  * Jones et al. [1996] M. E. Jones, D. S. Lemons, R. J. Mason, V. A. Thomas, D. Winske,  A grid-based coulomb collision model for pic codes,  Journal of Computational Physics 123 (1996) 169–181.
  * Lemons et al. [2009] D. S. Lemons, D. Winske, W. Daughton, B. Albright,  Small-angle coulomb collision model for particle-in-cell simulations,  Journal of Computational Physics 228 (2009) 1391–1403.
  * Cohen et al. [2013] B. I. Cohen, A. M. Dimits, D. J. Strozzi,  A grid-based binary model for coulomb collisions in plasmas,  Journal of Computational Physics 234 (2013) 33–43.
  * Sentoku and Kemp [2008] Y. Sentoku, A. Kemp,  Numerical methods for particle simulations at extreme densities and temperatures: Weighted particles, relativistic collisions and reduced currents,  Journal of Computational Physics 227 (2008) 6846 – 6861.


# Optical Two-dimensional Coherent Spectroscopy of Cold Atoms

Danfu Liang  Department of physics, Florida International University, Miami,
FL. 33199, USA  Lexter Savio Rodriguez  Department of physics, Florida
International University, Miami, FL. 33199, USA  Haitao Zhou  Department of
physics, Florida International University, Miami, FL. 33199, USA  Yifu Zhu
Department of physics, Florida International University, Miami, FL. 33199, USA
Hebin Li  Department of physics, Florida International University, Miami, FL.
33199, USA  Corresponding author<EMAIL_ADDRESS>
###### Abstract

We report an experimental demonstration of optical 2DCS in cold atoms. The
experiment integrates a collinear 2DCS setup with a magneto-optical trap
(MOT), in which cold rubidium (Rb) atoms are prepared at a temperature of
about 200 $\mu$K and a number density of $10^{10}$ cm-3. With a sequence of
femtosecond laser pulses, we first obtained one-dimensional second- and
fourth-order nonlinear signals and then acquired both one-quantum and zero-
quantum 2D spectra of cold Rb atoms. The capability of performing optical 2DCS
in cold atoms is an important step toward optical 2DCS study of many-body
physics in cold atoms and ultimately in atom arrays and trapped ions. Optical
2DCS in cold atoms/molecules can also be a new avenue to probe chemical
reaction dynamics in cold molecules.

††journal: ol

Ultrafast femtosecond lasers usually are not used for studies in cold atoms
since the two systems are apparently incompatible in their characteristic time
scales. A few experiments have utilized ultrafast lasers to demonstrate
ultrafast manipulation of quantum states and entanglement in trapped ions [1,
2] and neutral atoms [3, 4], enabling the possibility of atom-based quantum
gate operation at picosecond or shorter time scales. The high intensity of
femtosecond pulses also provides an opportunity to perform ultrafast nonlinear
spectroscopy in cold atoms which has been largely overlooked. Here we report
the implementation of an advanced ultrafast spectroscopic technique, namely
optical two-dimensional coherent spectroscopy (2DCS), in a cold atom cloud
prepared in a magneto-optical trap (MOT).

Originated in nuclear magnetic resonance [5], the concept of multi-dimensional
Fourier-transform spectroscopy has been implemented in the optical region
using femtosecond lasers and developed into a powerful tool to study energy
level structures, couplings, and dynamics in a variety of complex systems such
as proteins [6], photosynthetic systems [7, 8], semiconductor quantum wells
[9, 10, 11, 12, 13], quantum dots [14, 15, 16, 17], 2D materials [18, 19],
perovskites [20, 21], atomic vapors [22, 23, 24, 25, 26, 27, 28, 29, 30, 31,
32, 33, 34, 35], and weakly-bound molecules on helium nanodroplets in a
molecular beam [36]. Particularly, double-quantum 2DCS provids an extremely
sensitive background-free detection of dipole-dipole interactions in both
potassium (K) [25, 30] and rubidium (Rb) [27, 30] atomic vapors. The technique
can be extended to multi-quantum 2DCS which probe multi-atom correlated
coherent states (Dicke states) with a scalable and deterministic number of
atoms, up to eight atoms [30, 32] in a K atomic vapor. The collective states
of higher excited states ($D$ state), in addition to $P$ state, has also been
observed [34] by using double-quantum 2DCS. These studies have demonstrated
optical 2DCS as a powerful technique to study many-body interactions and
correlations in atomic systems.

However, atomic vapors at room or higher temperatures pose extra challenges in
the quantitative and deterministic study of many-body interactions and
correlations due to the presence of thermal motion. For example, although the
number density and the mean interatomic separation are known at a given cell
temperature, the exact separation and number of atoms have a statistical
distribution but are not deterministic for each shot of the experiment. The
thermal motion also speeds up the decoherence of many-body states, making it
difficult to form and observe coherent states with a large number of atoms. In
contrast, cold atoms provide a well-controlled and isolated environment with
virtually no thermal motion. Especially, recent advances in trapping neutral
atoms in an array of optical tweezers [37, 38, 39, 40, 41, 42] open the
possibility of studying many-body physics in an atom array with a
deterministic number of atoms and spatial distribution. Demonstrating optical
2DCS in a cold atom cloud is an essential first step toward this goal. Recent
progress in optical 2DCS made it feasible for 2DCS studies in cold atoms. The
sensitivity has been improved such that the dipole-dipole interaction can be
detected at a density ($10^{8}$ cm-3) lower than a typical cold-atom density
[30]. Optical frequency-comb-based 2DCS has achieved a frequency resolution
that is sufficient to resolve all hyperfine levels in Rb and K atoms [43, 44].

Figure 1: (a) Schematic of the experimental setup integrating a collinear 2DCS
setup with an Rb MOT. AOM: acousto-optic modulator. QWP: quarter wave plate.
PD: photodetector. (b) Relevant energy levels of Rb atoms at the $D_{2}$
transition and the frequencies of the lasers used in the experiment. The
hyperfine splittings are not to scale.

In this letter, we demonstrate the implementation of optical 2DCS in cold Rb
atoms by combining a collinear 2DCS setup with an MOT. Rubidium atoms are
laser cooled and trapped in the MOT at a temperature of $\sim 200$ $\mu$K. The
cold atom cloud includes approximately one million atoms at a number density
of $10^{10}$ cm-3. With the excitation pulses provided by the collinear 2DCS
setup, we first performed one-dimensional (1D) pump-probe and transient four-
wave mixing (FWM) spectroscopy. We then demonstrated optical 2DCS and acquired
one-quantum and zero-quantum 2D spectra of cold Rb atoms.

The experimental schematic is illustrated in Fig. 1(a). The system integrates
a collinear 2DCS setup (left) with an MOT (right). Four femtosecond pulses
derived from the collinear 2DCS setup are incident on a cold atom cloud
prepared in the MOT to perform 2DCS measurements. The MOT is constructed in an
octagonal glass cell in which the pressure is maintained $<10^{-8}$ torr. An
Rb dispenser provides both Rb isotopes at the natural abundance. Our
experiment cools and traps 85Rb atoms. The relevant energy levels and the
laser cooling scheme are shown in Fig.1(b). The cooling laser with an output
power of 200 mW is tuned to the transition
$|5^{2}S_{1/2},F=3\rangle\rightarrow|5^{2}P_{3/2},F^{\prime}=4\rangle$. The
repump laser with an output power of 30 mW is tuned to the transition
$|5^{2}S_{1/2},F=2\rangle\rightarrow|5^{2}P_{3/2},F^{\prime}=3\rangle$. The
cooling laser is divided into three beams and circularly polarized. The three
beams converge at the center of the glass cell in three orthogonal directions,
two horizontal and one vertical. Each beam is reflected back by a mirror so
there are two counter-propagating beams in each dimension for Doppler cooling.
The repump laser is also incident on the crossing of cooling laser beams to
pump atoms out of $|5^{2}S_{1/2},F=2\rangle$ state to maintain the cooling
cycle. The laser-cooled atoms are trapped in a magnetic trap generated by a
pair of anti-Helmholtz coils. In our experiment, the Rb atoms are cooled to a
temperature of $\sim$ 200 $\mu$K and the cold atom gas has a number density of
$10^{10}$ cm-3.

With cold atoms prepared in the MOT, optical 2DCS experiment is performed by
using a femtosecond laser tuned to excite the $D_{2}$ transition from
$|0\rangle=|5^{2}S_{1/2}\rangle$ to $|1\rangle=|5^{2}P_{3/2}\rangle$. The
femtosecond pulse is about $200$ fs in duration at a repetition rate of $78$
MHz. The laser spectrum has a central wavelength of $780$ nm and a bandwidth
of $2.6$ nm (the standard deviation). The collinear 2DCS setup consists of
three nested Mach-Zehnder interferometers [45], as shown in Fig. 1(a). The
femtosecond pulse is split into four pulses ($A,B,C,$ and $D$) which each go
through an acousto-optic modulator (AOM) and are subsequently combined into
one beam. The three time delays ($\tau$, $T$ and $t$) between the four pulses
are controlled by three translation delay stages. The four pulses are
modulated by AOMs at slightly different frequencies of $\Omega_{A}=80.107$
MHz, $\Omega_{B}=80.104$ MHz, $\Omega_{C}=80.0173$ MHz, and $\Omega_{D}=80$
MHz. At the same time, a continuous-wave (CW) reference laser copropagates
with the pulse laser and is modulated by the AOMs in each arm of the
interferometer. The combined reference laser beam is detected by a
photodetector (PD 1) to provide beat notes of the AOM frequencies. The
femtosecond pulses are incident on the cold atom cloud and generate nonlinear
fluorescence signals which are detected by another photodetector (PD 2) and a
lock-in amplifier. The nonlinear signal induced by a specific pulse sequence
can be measured by lock-in detection with a corresponding reference frequency
obtained by mixing beat notes of the CW reference laser beams in a digital
wave mixer.

Figure 2: Pulse sequences for (a) two-pulse pump-probe spectroscopy and (b)
four-pulse excitation for optical 2DCS. Double-side Feynman diagrams
describing contributing quantum pathways by (c) two-pulse excitation and (d)
four-pulse excitation.

Two different sequences of femtosecond pulses are used in the experiment. A
two-pulse sequence shown in Fig. 2(a) is used for pump-probe spectroscopy and
a four-pulse sequence shown in Fig. 2(b) is used for optical 2DCS. The pulses
with an asterisk represent that the pulse is considered conjugated,
corresponding to a negative AOM modulation frequency. The pump-probe
spectroscopy can be performed with any pair of pulses while the other two
pulses are blocked. This ensures a good pair-wise pulse overlap in both
spatially and temporally before starting the four-pulse experiment. Using
pulses $A^{*}$ and $B$ as an example, the lock-in amplifier detects the
fluorescence signal that is modulated at the reference frequency
$\Omega_{R0}=-\Omega_{A}+\Omega_{B}$. This signal is resulted from the
excitation pathway represented by the double-sided Feynman diagram in Fig.
2(c). The four-pulse sequence in Fig. 2(b) performs a fourth-order excitation
and the corresponding excitation pathways are represented by double-sided
Feynman diagrams in Fig. 2(d). In these processes, the first pulse $A^{*}$
creates a first-order coherence between states $|0\rangle$ and $|1\rangle$,
which evolves during time delay $\tau$. The second pulse $B$ converts the
first-order coherence into a population in either $|0\rangle$ or $|1\rangle$
depending on the relative phase. The third pulse $C$ generates a third-order
coherence between states $|1\rangle$ and $|0\rangle$, which evolves during
time delay $t$. Finally, the fourth pulse $D^{*}$ converts the third-order
coherence into a population in state $|1\rangle$ or $|0\rangle$. The
population in $|0\rangle$ does not emit fluorescence so the corresponding
pathways are not shown. The fourth-order population in state $|1\rangle$
decays and emits a fourth-order fluorescence signal, which can be detected by
a lock-in amplifier referenced to the frequency
$\Omega_{R1}=-\Omega_{A}+\Omega_{B}+\Omega_{C}-\Omega_{D}$. This fourth-order
nonlinear signal can be measured as a function of $\tau$, $T$, or/and $t$.
When the signal is measured with two varying time delays, a 2D spectrum can be
generated by Fourier transforming the time-domain signal with respect to two
time delays.

Figure 3: One-dimensional time-domain spectra for (a) two-pulse pump-probe
experiment and (b) the fourth-order signal due to four-pulse excitation. Real
part $X$, imaginary part $Y$, and amplitude $R=\sqrt{X^{2}+Y^{2}}$ are
plotted. (c) Frequency-domain spectra obtained by Fourier transforming the
time-domain pump-probe (blue) and fourth-order (red) signals.

We first performed one-dimensional (1D) scans with two and four pulses. In the
two-pulse experiment, the signal is measured by lock-in detection at the
reference frequency $\Omega_{R0}=-\Omega_{A}+\Omega_{B}$ as the time delay
between $A^{*}$ and $B$ is scanned. In the four-pulse experiment, the signal
measured at the reference frequency
$\Omega_{R1}=-\Omega_{A}+\Omega_{B}+\Omega_{C}-\Omega_{D}$ as the time delay
$t$ is scanned while the other two delays are fixed at $\tau=500$ fs and
$T=500$ fs. The obtained 1D time-domain spectra are shown in Fig. 3(a) for
two-pulse excitation and Fig. 3 for four-pulse excitation. The $X$ and $Y$
components are real and imaginary parts, respectively, of the lock-in
amplifier output and the magnitude of the signal is $R=\sqrt{X^{2}+Y^{2}}$.
Both 1D spectra reveal the dynamics of coherence $\rho_{01}$. The amplitude of
oscillations shows virtually no decay within the time window of 10 ps since
the coherence time is expected to be much longer at the ns scale. The signal
does not oscillate at the optical carrier frequency, but at a reduced
frequency $\nu^{*}=|\nu_{sig}-\nu_{ref}|$ which is the difference between the
signal frequency $\nu_{sig}$ and the CW reference laser frequency $\nu_{ref}$.
Fourier transforming these time domain signals gives the frequency-domain
spectra in Fig. 3(c), where the blue curve is the pump-probe spectrum and the
red curve is the fourth-order nonlinear spectrum. The frequency axis shows the
optical frequency calculated from the reduced frequency as
$\nu_{sig}=\nu^{*}+\nu_{ref}$. Both frequency-domain spectra have a resonance
at $384.2$ THz which matches the $D_{2}$ transition frequency. The linewidth
of the spectra is limited by the scanned delay time of 10 ps. The amplitude of
the pump-probe spectrum is greater than that of the fourth-order signal since
the pump-probe experiment measures a second-order nonlinear response.

Figure 4: (a) One-quantum 2D spectrum and (b) zero-quantum 2D spectrum of cold
Rb atoms. The amplitude of the spectra is plotted with the maximum normalized
to one.

Optical 2DCS measurements can be performed by scanning two time delays while
the third time delay is fixed. A one-quantum 2D spectrum, as shown in 3(a),
was obtained by scanning time delays $\tau$ and $t$ with $T=500$ fs. Fourier
transform of two time axes generates two frequency axes: absorption frequency
$\omega_{\tau}$ and emission frequency $\omega_{t}$ corresponding to $\tau$
and $t$, respectively. The resulting 2D peak is located on the diagonal line
at the $D_{2}$ transition frequency of $384.2$ THz. A zero-quantum 2D
spectrum, as shown in 3(b), was obtained by scanning time delays $T$ and $t$
with $\tau=500$ fs. The mixing frequency $\omega_{T}$ and the emission
frequency $\omega_{t}$ correspond to time delays $T$ and $t$, respectively.
The zero-quantum 2D peak has an emission frequency of $384.2$ THz while the
mixing frequency is zero since the population only has an exponential decay
with a lifetime of $26.2$ ns during time delay $T$. Both 2D spectra are
plotted with 20 contours with no visible background noise, demonstrating
excellent signal-to-noise ratio. The spectral resolution is $0.1$ THz limited
by the scanned time delay of $10$ ps so the hyperfine structures are not
resolved. The spectral resolution can be improved by scanning longer delays.
However, it is not practical to have a long enough delay to achieve a
sufficient frequency resolution ($29$ MHz) to fully resolve hyperfine
structures. Optical 2DCS with such frequency resolution is possible by
utilizing femtosecond frequency combs [43, 44]. The obtained 2D spectra of
cold atoms are consistent with the previous experiments [22, 35] performed in
hot Rb atomic vapor, demonstrating the capability of perform optical 2DCS in
cold atoms.

In summary, we have experimentally demonstrated the feasibility of optical
2DCS in cold atom clouds prepared in an MOT setup. We first probed 1D second-
order and fourth-order nonlinear signals to confirm that high-order nonlinear
signals can be generated in cold atoms by femtosecond laser pulses while
maintaining a stable cold atom cloud. We then extended the experiment to 2DCS
by scanning two time delays. Both one-quantum and zero-quantum 2D spectra were
acquired and the results are consistent with previous experiments in atomic
vapors. This measurement is extremely robust in probing a specific high-order
nonlinear process in cold atoms. This demonstration is an important first step
towards the goal of using optical 2DCS to study many-body physics in cold
atoms and ultimately in atom arrays trapped by optical tweezers. Similarly,
the technique also makes it possible to realize optical 2DCS experiments that
have been proposed to study quantum entanglements in trapped ions [46, 47].
Ultrafast spectroscopy has been traditionally used to study chemical
reactions, optical 2DCS in cold atoms/molecules will open a new avenue to
study reaction dynamics in ultracold molecular chemical reactions [48].

Funding. National Science Foundation (PHY-2216824). D.L. acknowledges the
support by a Dissertation Year Fellowship from FIU.

Disclosures. The authors declare no conflicts of interest.

## References

  * [1] J. Mizrahi, C. Senko, B. Neyenhuis, K. G. Johnson, W. C. Campbell, C. W. S. Conover, and C. Monroe, Phys. Rev. Lett. 110, 203001 (2013).
  * [2] J. D. Wong-Campos, S. A. Moses, K. G. Johnson, and C. Monroe, Phys. Rev. Lett. 119, 230501 (2017).
  * [3] Y. Song, H.-g. Lee, H. Kim, H. Jo, and J. Ahn, Phys. Rev. A 97, 052322 (2018).
  * [4] Y. Chew, T. Tomita, T. P. Mahesh, S. Sugawa, S. de Léséleuc, and K. Ohmori, arXiv:2111.12314 (2021).
  * [5] R. Ernst, G. Bodenhausen, and A. Wokaun, _Principles of Nuclear Magnetic Resonance in One and Two Dimensions_ (Oxford Science Publications, Oxford, U.K., 1987).
  * [6] P. Hamm, M. Lim, W. F. DeGrado, and R. M. Hochstrasser, Proc. Natl. Acad. Sci. U.S.A 96, 2036 (1999).
  * [7] J. Cao, R. J. Cogdell, D. F. Coker, H.-G. Duan, J. Hauer, U. Kleinekathöfer, T. L. C. Jansen, T. Mančal, R. J. D. Miller, J. P. Ogilvie, V. I. Prokhorenko, T. Renger, H.-S. Tan, R. Tempelaar, M. Thorwart, E. Thyrhaug, S. Westenhoff, and D. Zigmantas, Science Advances 6, eaaz4888 (2020).
  * [8] H.-G. Duan, V. I. Prokhorenko, R. J. Cogdell, K. Ashraf, A. L. Stevens, M. Thorwart, and R. J. D. Miller, Proceedings of the National Academy of Sciences 114, 8493 (2017).
  * [9] S. T. Cundiff, A. D. Bristow, M. Siemens, H. Li, G. Moody, D. Karaiskaj, X. Dai, and T. Zhang, IEEE J. Sel. Top. Quantum Electron. 18, 318 (2012).
  * [10] X. Li, T. Zhang, C. Borca, and S. Cundiff, Phys. Rev. Lett. 96, 057406 (2006).
  * [11] G. Nardin, G. Moody, R. Singh, T. M. Autry, H. Li, F. Morier-Genoud, and S. T. Cundiff, Phys. Rev. Lett. 112, 046402 (2014).
  * [12] R. Singh, T. M. Autry, G. Nardin, G. Moody, H. Li, K. Pierz, M. Bieler, and S. T. Cundiff, Phys. Rev. B 88, 45304 (2013).
  * [13] D. Turner, P. Wen, D. Arias, K. Nelson, H. Li, G. Moody, M. Siemens, and S. Cundiff, Phys. Rev. B 85, 201303 (2012).
  * [14] G. Moody, R. Singh, H. Li, I. A. Akimov, M. Bayer, D. Reuter, A. D. Wieck, A. S. Bracker, D. Gammon, and S. T. Cundiff, Phys. Status Solidi (b) 250, 1753 (2013).
  * [15] G. Moody, R. Singh, H. Li, I. Akimov, M. Bayer, D. Reuter, A. Wieck, and S. Cundiff, Solid State Commun 163, 65 (2013).
  * [16] G. Moody, R. Singh, H. Li, I. A. Akimov, M. Bayer, D. Reuter, A. D. Wieck, and S. T. Cundiff, Phys. Rev. B 87, 045313 (2013).
  * [17] G. Moody, R. Singh, H. Li, I. A. Akimov, M. Bayer, D. Reuter, A. D. Wieck, A. S. Bracker, D. Gammon, and S. T. Cundiff, Phys. Rev. B 87, 41304 (2013).
  * [18] G. Moody, C. Kavir Dass, K. Hao, C.-H. Chen, L.-J. Li, A. Singh, K. Tran, G. Clark, X. Xu, G. Berghäuser, E. Malic, A. Knorr, and X. Li, Nat. Commun. 6, 8315 (2015).
  * [19] M. Titze, B. Li, X. Zhang, P. M. Ajayan, and H. Li, Phys. Rev. Materials 2, 054001 (2018).
  * [20] F. Thouin, S. Neutzner, D. Cortecchia, V. A. Dragomir, C. Soci, T. Salim, Y. M. Lam, R. Leonelli, A. Petrozza, A. R. S. Kandada, and C. Silva, Phys. Rev. Materials 2 (2018).
  * [21] M. Titze, C. Fei, M. Munoz, X. Wang, H. Wang, and H. Li, J. Phys. Chem. Lett. 10, 4625 (2019).
  * [22] P. F. Tian, D. Keusters, Y. Suzaki, and W. S. Warren, Science 300, 1553 (2003).
  * [23] X. Dai, A. D. Bristow, D. Karaiskaj, and S. T. Cundiff, Phys. Rev. A 82, 052503 (2010).
  * [24] X. Dai, M. Richter, H. Li, A. D. Bristow, C. Falvo, S. Mukamel, and S. T. Cundiff, Phys. Rev. Lett. 108, 193201 (2012).
  * [25] X. Dai, M. Richter, H. Li, A. D. Bristow, C. Falvo, S. Mukamel, and S. T. Cundiff, Phys. Rev. Lett. 108, 193201 (2012).
  * [26] H. Li, A. D. Bristow, M. E. Siemens, G. Moody, and S. T. Cundiff, Nat. Commun. 4, 1390 (2013).
  * [27] F. Gao, S. T. Cundiff, and H. Li, Opt. Lett. 41, 2954 (2016).
  * [28] B. Lomsadze and S. T. Cundiff, Phys. Rev. Lett. 120, 233401 (2018).
  * [29] S. Yu, M. Titze, Y. Zhu, X. Liu, and H. Li, Opt. Express 27, 28891 (2019).
  * [30] S. Yu, M. Titze, Y. Zhu, X. Liu, and H. Li, Opt. Lett. 44, 2795 (2019).
  * [31] M. Binz, L. Bruder, L. Chen, M. F. Gelin, W. Domcke, and F. Stienkemeier, Opt. Express 28, 25806 (2020).
  * [32] D. Liang and H. Li, The Journal of Chemical Physics 154, 214301 (2021).
  * [33] S. Yu, Y. Geng, D. Liang, H. Li, and X. Liu, Opt. Lett. 47, 997 (2022).
  * [34] D. Liang, Y. Zhu, and H. Li, Phys. Rev. Lett. 128, 103601 (2022).
  * [35] J. Yan, S. Revesz, D. Liang, and H. Li, Phys. Rev. A 105, 052810 (2022).
  * [36] L. Bruder, U. Bangert, M. Binz, D. Uhl, R. Vexiau, N. Bouloufa-Maafa, O. Dulieu, and F. Stienkemeier, Nature Communications 9, 4823 (2018).
  * [37] A. M. Kaufman, B. J. Lester, C. M. Reynolds, M. L. Wall, M. Foss-Feig, K. R. A. Hazzard, A. M. Rey, and C. A. Regal, Science 345, 306 (2014).
  * [38] B. J. Lester, N. Luick, A. M. Kaufman, C. M. Reynolds, and C. A. Regal, Physical Review Letters 115, 073003 (2015).
  * [39] H. Kim, W. Lee, H.-g. Lee, H. Jo, Y. Song, and J. Ahn, Nature Communications 7, 13317 (2016).
  * [40] D. Barredo, S. de Léséleuc, V. Lienhard, T. Lahaye, and A. Browaeys, Science 354, 1021 (2016).
  * [41] M. Endres, H. Bernien, A. Keesling, H. Levine, E. R. Anschuetz, A. Krajenbrink, C. Senko, V. Vuletic, M. Greiner, and M. D. Lukin, Science 354, 1024 (2016).
  * [42] A. M. Kaufman and K.-K. Ni, Nature Physics 17, 1324 (2021).
  * [43] B. Lomsadze and S. T. Cundiff, Science 357, 1389 (2017).
  * [44] B. Lomsadze and S. T. Cundiff, Phys. Rev. Lett. 120, 233401 (2018).
  * [45] G. Nardin, T. M. Autry, K. L. Silverman, and S. T. Cundiff, Opt. Express 21, 28617 (2013).
  * [46] F. Schlawin, M. Gessner, S. Mukamel, and A. Buchleitner, Phys. Rev. A 90, 023603 (2014).
  * [47] M. Gessner, F. Schlawin, H. Häffner, S. Mukamel, and A. Buchleitner, New Journal of Physics 16, 092001 (2014).
  * [48] M.-G. Hu, Y. Liu, D. D. Grimes, Y.-W. Lin, A. H. Gheorghe, R. Vexiau, N. Bouloufa-Maafa, O. Dulieu, T. Rosenband, and K.-K. Ni, Science 366, 1111 (2019).