# Stochastic laser cooling enabled by many body effects

Roie Dann and Ronnie Kosloff  The Institute of Chemistry, The Hebrew
University of Jerusalem, Jerusalem 91904, Israel<EMAIL_ADDRESS>
###### Abstract

A novel laser cooling mechanism based on many body effects is presented. The
method can be applicable for cooling a large class of atoms and molecules in
higher density than commonly excepted by existing methods. The cooling
mechanism relies on the collective encounters of particle and light.
Stochastic events between the particles and photons as well as a collective
effect give rise to energy transfer between these media. Such mechanism relies
on multiple light-matter encounters, therefore requiring a sufficient particle
density, $\rho\sim 10^{14}\text{cm}^{-3}$. This is an advantage for
experiments where high phase space density is required. A second tuning laser
can be added increasing the applicability to many types of atoms and
molecules. This tuning laser changes the inter-particle potential by inducing
an AC stark effect. As a result the required trapping density can be reduced
down to $\rho\sim 10^{6}\text{cm}^{-3}$. Simulations of phase space
distributions were performed comparing different particle densities, trap
potentials and light field intensity profiles. The modelling shows efficient
cooling rates up to $~{}10^{2}\text{K/s}$ for a dense ensemble of $~{}^{87}$Rb
atoms, and cooling rates up to $~{}6\cdot 10^{2}\text{K/s}$ when adding an
additional tuning source.

##  I Introduction

Atom-photon interactions has been a major research topic in physics. Its
origins can be traced to Kepler who as early as 1619, suggested that light may
have a mechanical effect, when observing that a comet’s tail is always
pointing away from the sun Kepler (1619). Later on Maxwell suggested a
phenomena known as "light pressure", pressure exerted on a surface when
exposed to electro-magnetic radiation Maxwell (1862). The topic was
revolutionized by two papers (1909 and 1916), when Einstein, following
Planck’s law of black body radiation, showed that light energy quanta must
carry a momentum set by $p=\frac{h}{\lambda}$ Einstein (1909, 1917a, 1917b).

Many experimental realizations exploiting the photon momentum have been
performed Metcalf and Van der Straten (2012). In the present study we exploit
this momentum transfer for cooling an ensemble of colliding atoms. The
frequency shift during the collision is the source of momentum exchange
between the particles and light.

The first suggestion for cooling atoms via photon-atom interactions was
proposed by H.E.D. Scovil in 1959 Scovil and Schulz-DuBois (1959). Scovil
pioneered a quantum thermodynamic approach to laser cooling thus introducing
the first quantum refrigerator. Further advancement in laser cooling was not
recorded until more than a decade later. In 1975, simultaneously and
independently of Scovil’s work, two groups of Wineland and Dehmelt, Wineland
and Dehmelt (1975) as well as Hansch and Schawlow Hänsch and Schawlow (1975)
introduced new theories for laser cooling. Wineland’s and Dehmelt’s work treat
the cooling of ions in an ion trap, and Hansch’s and Schawlow’s theory
concentrates on neutral atoms. The initial theory proposed by Wineland and
Dehmelt, known as the Doppler Cooling, involves energy transfer from the
atomic media to photons depending on the relative velocity of the atoms to the
light source. The Doppler Cooling theory predicts a minimum temperature known
as the ’Doppler limit’ which for Sodium and Rubidium atoms amounts to
$240\mu\text{K}$ and $146\mu\text{K}$, respectively Steck (2001, 2000). When
experimental studies based on the theory took place the cooling was
unexpectedly efficient and led to temperatures below the Doppler limit,
raising theoretical questions concerning the underlying mechanism.

A significant effort was aimed at extending or replacing the Doppler cooling
theory. Diverse theories were developed such as Raman cooling, cavity mediated
cooling, and Sisyphus cooling Kasevich and Chu (1992); Ketterle et al. (1992);
Diedrich et al. (1989); Kerman et al. (2000); Dalibard and Cohen-Tannoudji
(1985, 1989); Neuhauser et al. (1978); Leibfried et al. (2003); Boozer et al.
(2006); Metcalf and Straten (2007); Breeden and Metcalf (1981). All proposed
theories describe different mechanisms of energy transfer between atomic and
photon media, and take advantage of the electro-magnetic waves’ global
character, enabling transport of energy away from the atomic medium.

The Sisyphus cooling theory was proposed in 1989 by Cohen Tannoudji and Jean
DalibardDalibard and Cohen-Tannoudji (1985). The theory involves two
interfering laser beams creating a standing wave with a polarization gradient
oscillating spatially between three polarizations $\sigma^{+}$, $\pi$ and
$\sigma^{-}$. The periodic potential imposed on the atoms affects the ground
and excited states differently, resulting in a varying energy gap alternating
spatially. Such spatial dependence of the energy gap allows an average energy
transfer from the moving atoms via repetitive excitations. In conjunction with
the theoretical work, experimental research achieved nano Kelvin temperatures,
setting the stage for the materialization of Bose-Einstein condensates by
additional evaporative cooling Anderson et al. (1995); Davis et al. (1995);
Zwierlein et al. (2003).

In the present cooling theory we change the focus from a single atom picture
to a collective many body approach in a dense particle medium. A new cooling
mechanism is proposed converting kinetic to optical energy, based on a
stochastic modulation of the emission frequency due to relative motion of the
particles.

The present cooling theory is based on high particle density when collisions
occur. It differs from the well established laser cooling mechanism
appropriate for a sparse medium. The high density is advantageous for
experiments where a large phase space density is desired. It has been claimed
Gallagher and Pritchard (1989) that inter-atomic collisions in high densities
will lead to trap loss or heating. Nevertheless, collective effects of light
trapped in the particle medium have been reported in ultracold medium Walker
et al. (1990); Ketterle et al. (1993); Townsend et al. (1995), experimentally
demonstrating that dense cold samples are feasible. The trapped light
generates an internal pressure which leads to the expansion of the atomic
cloud Pruvost et al. (2000); Bachelard et al. (2016); Guerin et al. (2017).
More recent experiments Camara et al. (2014) have observed these collective
effects in very large magnetic-optical-traps. In cooling, the high density
effects should become our ally and not our foe.

Theoretical approaches to model these phenomena are based on continuum
hydrodynamical theories Pruvost et al. (2000); Rodrigues et al. (2016). We
will adopt such a hydrodynamic description for our cooling theory in addition
to absorption and emission properties which depend on the atomic properties
and density.

The theory is demonstrated on cooling Rubidium 87 atoms. Rb has been the
workhorse of ultra-cold atomic physics, due to favourable properties such as
convenient vapour pressure at room temperature and a large absorption cross
section and a large scattering length for S-wave scattering. The cooling
scheme can also be applied to a variety of atomic and molecular systems where
an approximate closed cycle transition can be found and the particle density
is dominated by two-body collisions. To expand the possible cooling candidates
we suggest adding an additional tuning laser which modifies the collision
parameters. This laser also enables to reduce the inter-particle density.
Sympathetic cooling of a mixture of species can also become applicable.

The Stochastic Cooling theory is introduced at the beginning of section II
followed by an explanation of the modelling technique, and derivation of the
different model variables in Sections III. Section IV presents the results of
the Stochastic laser cooling theory. An additional generalized scheme,
Enhanced Stochastic Cooling theory, which allows efficient cooling of non
alkali atoms and molecules is presented in Section V. Following the
theoretical method is a discussion and conclusions.

##  II Stochastic Cooling theory

###  II.1 Laser cooling scheme

The prerequisite of the scheme is a trapping potential able to confine and
isolate a dense ensemble of gas phase atoms or molecules. The cooling is based
on applying beams of light to the interior of the trap detuned below the
resonance frequency. This light diffuses out, trapped by the absorption and
emission process. Such events also shift the light’s frequency to the blue on
average. Eventually the light is emitted from the dilute exterior regions of
the atomic cloud, Fig. 1. If such a scenario can be maintained it is obvious
that on average the energy consumed to generate the blue shift in the
radiation frequency will be extracted from the kinetic energy of the atoms
leading to cooling.

Energy flow from the particle medium to the electromagnetic field will occur
only if a single photon will go through multiple excitation/dexcitation
cycles. This imposes a restriction on the absorption cross section and the
particle density. The particles should be stratified in the trap such that the
high density is in the centre. Simultaneously, the particles undergo
diffusional dynamics as a result of repetitive inter-particle collisions and
interaction with light. Overall, under these conditions, two different
mediums, particle and photons, are captured in the trap interacting with each
other by energy transfer and influencing their respective motion. We will
refer to the photonic medium as the ’light medium’ in the following. The
energy shift of the light is a stochastic many body effect incorporating
multiple absorption emission cycles. It combines the asymmetry in absorption
and emission line-shapes with the spectral dependence of the light trapping in
the medium, the mechanism is explained in detail in the following.

Figure 1: $~{}^{87}\mbox{Rb}$ atoms are confined by a MOT (green), a laser
source, detuned slightly from resonance (red), is applied to the particle
ensemble (represented in grey). The photons are confined to the trap due to
repeated absorption. They diffuse through the atomic media until they escape
the trap blue detuned with respect to the incident laser (shown in blue and
purple arrows). Inter-particle collisions modify the absorption cross section
and locally equilibrate the kinetic energy.

The photons captured in the trap undergo repetitive excitation cycles, for
each cycle there is a probability of an energy shift to the photon on account
of the particle’s kinetic energy. Due to the asymmetry of the absorption
probability function (incident light is detuned below the atomic line) a
’blue’ shifted photon will have a higher probability of being reabsorbed, Fig.
4, 2. While a ’red’ detuned photon probability to be reabsorbed is decreased,
i.e, photons which transfer energy to the particle medium will diffuse faster
through the particle medium until reaching low particle densities on the edge
of the trap and escaping. This causes an effective cut off for energy transfer
from the light medium to the particles. Alternatively, photons detuned to the
’blue’ which reduce the particles’ energy, undergo more excitation cycles,
allowing further energy transfer from the particle medium to the light, Fig.
2. The collective effect of dual dependence induces a net energy transfer
between the two media and efficient cooling.

Figure 2:  A photon of the laser frequency, $\omega_{L}$, is absorbed by a
$~{}^{87}\mbox{Rb}$ atom. In a random process the frequency of the emitted
photon, $\omega$, is shifted with respect to $\omega_{L}$. For a positive
energy shift (top part) the probability to absorb increases and conversely for
a negative energy shift (lower part).

Another characteristic of the dual dependence involves pressure broadening
Pichler et al. (1983): When the density of an atomic ensemble increases, a red
shift to the atom transition frequency occurs. This red shift is universal and
is caused by the larger polarizability of the excited state which enhances the
long range attractive van der Waals force. The particle density profile can
therefore induce a spatially varying optimal absorption frequency. For a
particle density which decreases radially, the resonance frequency increases
accordingly. As a result, a photon starting in the centre when shifted to the
blue and propagating outward will be reabsorbed due to new resonance
conditions. This process will be repeated until a blue photon will escape the
trap at the low density outer region.

An important issue is how to stratify the light in the opaque particle medium.
A few options can be envisioned. When the trap forces are large causing a
large gradient in density a far red detuned incident light will only be
absorbed at the centre. Once the light frequency is shifted to the blue it
becomes trapped by the particles. The only escape route is from the dilute
outer boundary. Another option is to use electromagnetic induced transparency
EIT to inject the light to the interior.

###  II.2 Energy transfer mechanism:

For a homogeneous atomic gas with the typical S$\rightarrow$ P transition, the
excited inter-atomic van der Waals potential scales as $C_{3}/r^{3}$ compared
to $C_{6}/r^{6}$ for the ground state potential, Cf. Fig. 3, Stwalley et al.
(1978). As a result, the absorption frequency varies with the relative
distance. Once a photon is absorbed, the atom spends an average lifetime
($27.7\cdot 10^{-8}\text{s}$) in the excited state. At this stage, the two
neighbouring atoms will undergo random relative motion, until their decay by
photo-emission. This random motion, accompanied by Doppler phenomena,
collisions, random electromagnetic fields and the natural linewidth of the
excited state, cause an energy shift of the emitted photon. Summarizing the
phenomena: For sufficient density, the relative motion causes a change in the
van der Waals potential energy at the expense of the emitted photon energy.
Following the description, a random energy shift requires random relative
motion between neighbouring atoms, this is indeed the case in the
semiclassical limit, where the relative motion is isotropic.

Figure 3:  $~{}^{87}\mbox{Rb}$ energy levels, lowest singlet ground state
$X^{1}\Sigma_{g}$ and $a^{3}\Sigma_{g}$ and one of the excited states
$0_{u}^{-}$. The excitations occur in the long range part of the potential
$r\approx 10^{3}$Å. (Energy in units of wavenumbers).

It is important to note that in the present modeling we neglected elastic
Rayleigh scattering which occurs in addition to the repeated
absorption/emission events. In typical cases, such as propagation of light
through biological tissue or planetary atmosphere, Chandrasekhar (2013); Dudko
and Weiss (2005), elastic scattering constitutes the main contribution, thus
influencing photon propagation. However, for photons near atomic absorption
line propagating in a particle medium with a large absorption cross section,
the absorption phenomena is the major contribution to the light propagation,
arising from a large difference in the typical lifetime of the two processes.

In the exposition the inter-atomic interactions are modeled by the four
electronic energy states of Rb2 molecules, the two ground states; the singlet
and triplet, $X\,^{1}\Sigma_{g}^{+}$ and $a\,^{3}\Sigma_{u}^{+}$
correspondingly, and $0_{g}^{-}$ and $1_{g}^{-}$ excited states Steck (2001).
The two ground states differ at close range distances, but for large distances
($r>100$ Å), the ground states’ singlet and triplet coalesce, scaling as van
der Waals interactions $\propto-1/r^{6}$. The excited state’s long range
potential scale as $\propto-1/{r^{3}}$ due to a degeneracy of the P state
Stwalley et al. (1978).

##  III Modelling methods for the combined particle and light media

###  III.1 Probabilistic analysis over phase space

The Stochastic Cooling theory was modelled by a probabilistic simulation where
the physics is embedded in terms of the dynamics of continuous probability
distribution functions (PDF) over phase space. This is the suitable
description for diffusional behaviour and dominant collective effects. Both
the light and particles are confined in the trap, and are described by the
position and momentum in the trap. It is important to note that for the
particles the momentum is proportional to the velocity while for the photons
it is linearly dependent on the energy.

A full stochastic model involves a 12 dimensional probability function,
including all the particle and photon degrees of freedom (DOF). Such a system
is computationally very demanding. However, if an isotropic environment is
assumed, all the axes are degenerate, and only 4 DOF (position and momentum
DOF for particles and light) are required. A 4D model is still computationally
challenging with respect to the desired accuracy. A solution to this problem
is achieved by comparing typical time scales characterizing both media. The
particle diffusion and thermalization rate is much faster than the resonant
photon diffusion rate. Such a separation of time scales effectively decouples
the two ensembles in a short time regime. This assumption allows us to break
down the general model to two separate phase space distributions, viz. to
particle and light media. This separation follows the mean field
approximation.

####  III.1.1 The Fokker Planck Eq. particle ensemble

##### Initialization:

In the initial stage, the particles, described by the PDF, $P$, are confined
in a trap with an initial temperature $T_{init}$. $P$ is propagated in time by
the Smoluchowski equation Einstein (1905); Sutherland (1905); Von Smoluchowski
(1906); Kolmogoroff (1931) until steady state is reached.

$\displaystyle\frac{\partial P({x_{par},p_{par}})}{\partial
t}=-\frac{\partial}{\partial
x}_{par}\left({\frac{p_{par}}{m}P}\right)+\frac{\partial}{\partial
p}_{par}\left({(V_{h.o}^{\prime}(x_{par})+\mu_{par}p_{par})P}\right)$ (1)
$\displaystyle+D_{par}\left({\mu_{par},T_{init}}\right)\frac{\partial^{2}P}{\partial^{2}p_{par}}$

where: $x_{par}$ and $p_{par}$ are the particle position and velocity,
correspondingly; $m$ is the particle mass; $\mu_{par}$ is the drag constant,
calculated from the experimental relaxation time Bloch et al. (2001).
$D_{par}(\mu_{par},T_{init})$ is the particle momentum diffusion function,
dependent on the drag constant and temperature.

The first term on the R.H.S. describes the coupling between the velocity and
location of the atoms. The collisions between the atoms transfer momentum
between the two particles, creating an overall diffusion in momentum which is
described by the last term. Balancing the diffusion is the trap’s potential,
associated to the term $\frac{\partial}{\partial
p_{par}}\left({V_{h.o}^{\prime}(x_{par})P}\right)$, a mixed term coupling the
confining force and the momentum and $\frac{\partial}{\partial
p_{par}}\left({\mu_{par}v_{par}P}\right)$ is the drag term originating from
particle collisions. An additional term can be added in extremely low
temperatures, where the particle de Broglie wavelength is in the order of the
mean inter-atomic distance and the particles’ scattering length should be
considered. Large scattering length should lead to an additional spatial
diffusion term. For the studied temperature regime this effect is negligible.
For particles in an harmonic trap potential the steady state distribution has
been shown to be a Gaussian with a variance dependent on the ratio between the
diffusion constant and the drag force, $D_{par}=\mu_{par}k_{B}T$, where
$k_{B}$ is the Boltzmann constant Einstein (1905); Smoluchowski (1952)

##### Coupling of the light to the particle medium:

The interacting particle light ensemble is modeled by the diffusion function,
$D_{par\,coupled}$. The change in the diffusion variable arises from an
average energy flow from the particles’ ensemble to the radiation field and
momentum alterations by photon absorption/emission processes. The variable
$D_{par\,coupled}$ is described in detail in Section III.2.

The radiation characteristics are described in detail in table IX.3.

####  III.1.2 The Fokker Planck Eq. for light

To describe the light medium we construct a second 2D probability distribution
function over phase space. The function is propagated in time with a F-P Eq.
derived from the ’Radiative Transfer Equation’ (RTE) Chandrasekhar (2013),
similar to the Photon Diffusion Eq. Furutsu and Yamada (1994); Durduran et al.
(1997), further details are given in the Appendix. These equations usually
describe light propagation in a scattering medium. While, in our study, we
treat excitation cycles as scattering events characterized by long interaction
times, resulting from the atomic decay time.

The dynamics of the light phase space are described by the following equation;

$\displaystyle\frac{\partial\phi(x,p_{l},t)}{\partial
t}=\frac{\partial}{\partial
x}\left({D_{x_{l}}(\rho_{par}(x),E_{photon}(p_{l}))\frac{\partial}{\partial
x}\phi(x_{l},p_{l},t)}\right)$ (2) $\displaystyle+\frac{\partial}{\partial
p_{l}}\left({D_{p_{l}}(\rho_{par}(x),E_{photon},T_{par})\frac{\partial}{\partial
p_{l}}\phi(x,p_{l},t)}\right)$

$x$ and $p_{l}$ are the position and momentum of the photon ensemble,
$E_{photon}$ is the photon energy, $T_{par}$ and $\rho_{par}$ are the
instantaneous particle temperature and density, respectively.

The equation has two diffusion terms, in space and momentum, describing the
diffusion in the particle medium and energy transfer. This is a similar
equation to the general Photon Diffusion equation Furutsu and Yamada (1994);
Durduran et al. (1997) but lacks any source or sink term, due to the fact that
for atoms there are no clear non-radiative processes. Loss mechanisms, such as
photo-association are also negligible for this case.

##### Diffusion functions explanation:

$D_{x_{l}}(\rho_{par}(x),E_{photon})$ is the light position diffusion
function; the photons’ propagation in the particle medium is described by a
diffusional movement, caused by repeated absorption/emission cycles and the
isotropic nature of spontaneous emission. The particle density sets the mean
distance between consecutive absorptions, the photons’ energy compared to the
transition line determines the probability of absorption, both affecting the
diffusion rate directly.

$D_{p_{l}}(\rho_{par}(x),E_{photon},T_{par})$ is the light momentum diffusion
function: The diffusion rate is determined by atom-atom interactions
influenced by the atomic density and velocity at temperature $T_{par}$.

The physics and interaction between both media is embedded in the properties
of these variables as a function of the different parameters. A complete
analysis follows.

###  III.2 Derivation of the diffusion variables

The light-matter momentum diffusion function, $D_{par\,coupled}$;

$D_{par\,coupled}=m\cdot\mu(k_{B}T_{inter}+E_{recoil}\cdot R)$ (3)

$D_{par\,coupled}$ is determined by accounting for all the different effects
influencing the particles’ energy or momentum, considering energy and momentum
conservation. Each term represents a different mechanism for energy transfer.
The $\mu\cdot E_{recoil}\cdot R$ term arises from the condition of a pressure
balance between radiation pressure and the particles’ momentum, when the
recoil temperature is achieved Walker et al. (1990). The light medium exerts a
constant radiation pressure on the particles by continuous absorption. The
effect is insignificant when the magnetic force of the trap is bigger than the
radiation pressure force, but should be considered when the particles are
cooled to a temperature where both forces are on the same scale. $E_{recoil}$
is the recoil energy and $R$ is a constant dependent on the ratio of photons
to particles or the intensity of the light.

The energy conservation between the ensembles is described by the drag
constant times the typical kinetic energy of a single particle, $\mu
k_{b}T_{int}$. The total energy of both media is kept constant by adjusting
the temperature variable, $T_{inter}$. The local energy transfer between the
particles to the radiation field is calculated from the net energy change due
to the interactions, as well as the total energy change arising from the
photon flow in and out of the trap.

A further term $(-D_{int})$ can be added to the diffusion coefficient. The
added term relates the instantaneous momentum transfer between the two media.
However, the additional term is negligible in the long range due to momentum
transfer accounted for in the energy transfer term, $m\cdot\mu
k_{B}T_{inter}$.

The present description does not account for quantum effects. At lower
temperatures the theory should be modified by adjusting the modeling
parameters, taking into account the particle wave characteristics. This can be
done by adding a position diffusional term describing the weak localization
due to collision diosi1993calderia .

####  III.2.1  Spatial diffusion amplitude,
$D_{x_{l}}(\rho_{par}(x),E_{photon}(p_{l}))$

The RTE derivation for photon propagation in a highly scattering medium
predicts the value of : $D_{x_{l}}=\frac{v}{3\mu_{s}^{\prime}}$ where
$\mu_{s}^{\prime}=\mu_{s}(1-g)$ , and $\mu_{s}^{-1}$ is the mean distance
between consecutive scattering events in the original derivation. In the case
above, where scattering events are neglected, $\mu_{s}^{-1}$ is the mean
distance between consecutive absorption events, $g$ is the scattering
anisotropy constant $\langle{cos(\theta)}\rangle$, which vanishes for
isotropic scattering Graaff and Ten Bosch (2000). Following the assumptions
mentioned, III.1.2, we derive a similar expression for $D_{x_{l}}$.

$D_{x_{l}}=\frac{l^{2}}{3\delta t}$ (4)

For the range of densities common for a MOT, the emission decay time, $\delta
t$, is the relevant time scale between adjacent excitations, and $l$ is the
mean distance between consecutive absorptions. For near resonance light,
$\delta t$ is the lifetime of the excited state, independent of the detuning
Guerin et al. (2017); Lagendijk and Van Tiggelen (1996); van Rossum and
Nieuwenhuizen (1999)

We assume a homogeneous medium for a small element in space. In such a medium
the probability distribution for a photon to cover a distance $y$ without
being absorbed by a particle is:

$P(y)=\sigma_{abs}ne^{-\rho_{par}(x)\sigma_{abs}y}$ (5)

$\sigma_{abs}(\nu)$ is the absorption cross section, $\nu$ is the photon
frequency and $x$ the position in the MOT. The mean free path is given by
$l=\langle{y}\rangle=\int_{0}^{\infty}yP(y)dy=\frac{1}{\rho_{par}\sigma_{abs}}$
Valeau et al. (2006). From Eq. 4 and 5 we obtain:

$D_{x_{l}}=\left[{3(\rho_{par}(x)\sigma_{abs}(\nu))^{2}\delta t}\right]^{-1}$
(6)

####  III.2.2 Momentum diffusion function,
$D_{p_{l}}(\rho_{par}(x),E_{photon},T_{par})$

The diffusion coefficient of light can be decomposed into a product of two
contributions: 1\. The probability function of a photon being absorbed by a
pair of interacting atoms, $G_{abs}(\rho_{par},\nu)$. 2\. The diffusion
function describing the diffusion rate in momentum ($\propto E$), caused by
the random energy shift of the absorbed photon,
$\mathcal{D}(\rho_{par},T_{par})$.

Assuming the absorption and energy transfer mechanisms are independent and the
diffusion function can be written as;

$D_{p_{l}}\left({\rho_{par}(x),E_{photon},T_{par}}\right)=G_{abs}(\rho_{par},\nu)\cdot\mathcal{D}(\rho_{par},T)$
(7)

##### Absorption probability function:

The details of the derivation of the absorption probability function,
$G_{abs}(\rho_{par},\nu)$, is shown in the appendix. Here, we present an
overall description and the results, Fig. 4.

The absorption probability is calculated employing a quantum description of
the absorption and emission process. For a low particle density, Cf Table
IX.1, the analysis can be restricted to a two particle interaction. The
calculation is then reduced to a three-body interaction, two neutral
$~{}^{87}\mbox{Rb}$ atoms and a light field characterizing a single photon.
The absorption cross section is solved assuming a weak field by the time
dependent perturbation theory. The free propagation is obtained by solving the
time dependent (TD) Schrödinger equation. The wave propagation is calculated
using the Chebychev polynomial expansion method with a Fourier grid, together
with a Gaussian Random Phase approach Tal-Ezer and Kosloff (1984); Koch et al.
(2009). The details of the propagation are described in Appendix IX.4.

Figure 4: The normalized absorption probability of $~{}^{87}\text{Rb}_{2}$ as
a function of the shift from the atomic transition line,
$G_{abs}(\rho_{par},\nu)$. The details of the calculation can be found in
Appendix IX.6. The function is used as an input for the momentum diffusion
function, $D_{p_{l}}$, using the experimental cross-section value to rescale
the probability function.

##### Energy transfer between the atom and radiation field:

There are a number of processes which cause a photonic energy shift: The
natural line broadening due to spontaneous emission Tannoudji et al. (1973);
Weisskopf and Wigner (1930), and Doppler phenomena Siegman (1986). However, in
high density the most dominant phenomena is pressure broadening and pressure
shift Rotondaro and Perram (1997); Peach (1981), which arise from atom-atom
interactions. The process is stochastic and can be viewed as a 1D random walk
on an energy axis. In such a case the momentum diffusion amplitude is the
squared mean momentum related to the energy shift per unit time. The function
can be calculated as, $\mathcal{D}(\Delta p)\propto\frac{var(F(\Delta
E))}{\delta t}$, where $var(F(\Delta E))$ is the variance of a specific
probability function, $F(\Delta E)$. $F(\Delta E)$ describes the probability
for a certain energy shift, between the emitted and absorbed photons. We show
here the main point of the derivation of $F(\Delta E)$.

We begin with change in energy due in a typical excitation:

$\Delta
E=-\frac{C_{3}}{r_{f}^{3}}+\frac{C_{6}}{r_{f}^{6}}+\frac{C_{3}}{r_{i}^{3}}-\frac{C_{6}}{r_{i}^{6}}=C_{3}\bigg{(}{\frac{1}{r_{i}^{3}}-\frac{1}{r_{f}^{3}}}\bigg{)}+C_{6}\bigg{(}\frac{1}{r_{f}^{6}}-\frac{1}{r_{i}^{6}}\bigg{)}$
(8)

where $r_{i}$ and $r_{f}$ are the relative distance between a pair of
functions at absorption and excitation times respectively, and $C_{3}$,
$C_{6}$ are the van der Waals potential’s constants.

Transforming to the centre of mass and relative velocity coordinates, the
velocity distribution is a Maxwell Boltzmann distribution of particles with a
reduced mass $\mu=m/2$ and kinetic energy of $E_{k}=\frac{p_{r}^{2}}{2\mu}$:

$f(v)=\sqrt{\frac{\mu}{2\pi k_{B}T}}e^{-\frac{\mu v^{2}}{2k_{B}T}}$ (9)

The initial relative density determines the average distance,
$r_{i}=(\rho(x))^{-1/3}$,($x$ is the position in the trap), and the final
distance is written in terms of the decay time and the initial distance,
$r_{f}=r_{i}+v\cdot\delta t$, where $r_{i}\gg v\cdot\delta t$ in the relevant
density and temperature range.

By expanding up to the first term in $v\cdot\delta t$, we obtain a relation
between the energy transferred and the particle’s relative velocity. From this
relation the energy transfer distribution function is obtained by a random
variable transformation for the Maxwell Boltzmann distribution.

$\Delta E=3(\rho(x))^{4/3}(C_{3}-2\rho(x)C_{6})\cdot v\cdot\delta t=C(x)v$
(10)

$C(x)=3(\rho(x))^{4/3}\delta t(C_{3}-2\rho(x)C_{6})$

The distribution function of the change in energy for a single excitation;

$F(\Delta E)=N_{norm}e^{-\frac{\mu(\Delta E)^{2}}{2C^{2}k_{B}T}}$

Making an ansatz in equation 7, the radiation field momentum diffusion
amplitude can be written as;

$\displaystyle
D_{p_{l}}(\rho_{par}(x),E_{photon},T_{par})=G_{abs}(\rho_{par},\nu)$ (11)
$\displaystyle\cdot\frac{[3(\rho(x))^{4/3}(C_{3}-2\rho(x)C_{6})]^{2}\delta
t\cdot k_{B}T_{par}}{\mu_{par}\cdot c^{2}}$

where $c$ is the speed of light.

The photon diffusion function is highly dependent on the density of the
particles and the spatial distribution of photons in the trap. The linear
temperature dependence demonstrates the fact that when the particles cool it
becomes harder to extract entropy.

Figure 5: A a schematic flow chart of the modelling method, as described in
section III.3

###  III.3 Final modelling summary

The Fokker-Plank equation is propagated by a Chebychev polynomial expansion
method for the evolution operator $U(t+\Delta t)=e^{-\hat{G}\Delta t}\rho(t)$,
where $\rho(t)$ is the modelled distribution function at time $t$ and the
propagator, $\hat{G}=\frac{\partial\rho(t)}{\partial t}$, is the corresponding
Fokker Planck operator Kosloff (1994); Johnson (2011). A Fourier method is
used to calculate the derivative terms in $\hat{G}$ operation. This scheme is
highly accurate and efficient. The two-phase space distribution of light and
particle ensembles are propagated simultaneously for small time laps,
transferring information about energy, momentum and density between the
models, Fig. 5. Absorbing boundary conditions are applied to the light density
function to account for the photons escaping the trap. In addition, new
photons are added with a frequency distribution corresponding to the laser
source. Such a scenario models a constant laser incident intensity.

##  IV Results A: Probabilistic analysis of the Stochastic Cooling

Following the evolution of the initial particle Gaussian distribution after a
transient time its phase space distribution is compressed. This is a signature
of cooling. Figure 6, (A,B plots) shows an increase of phase space density
after $6\mu\text{s}$. On the other hand, the light medium experiences a fast
broadening of the momentum distribution (time scale of $0.1\mu\text{s}$). For
low momentum, the distribution reaches a threshold due to a rapid decrease of
the absorption probability and a fast spatial diffusion of the low (red
detuned) momentum photons escaping the trap. While high momentum photons are
confined for longer periods of time in the particle medium. This occurs until
photons reach off-resonant frequencies, leading to a fast diffusion for
extremely high frequency blue detuned photons. The phenomena can be seen in
figure 6 D as off-resonance photons (large gaps between the resonant momentum,
$8.334901e^{-28}\frac{\text{Kgm}}{\text{s}}$; bright horizontal strips in the
figure) diffuse rapidly. At this stage of the process the cooling continues
and approaches a constant rate, (see Fig. 7), resulting from continuous
replacement of the diffused photons by laser light absorbed by the particle
medium.

The cooling rate is linear for the initial coupling with light but eventually
saturates because the energy transfer depends on the particle velocity,
$D_{p_{l}}\propto T_{par}$. The cooling rate slows down for low temperatures
until reaching the quantum regime where additional processes should be
incorporated in the model. An example of such an effect is a further
contribution to the spatial diffusion resulting from localization of the
particle wave packet due to collision with a neighbouring particle.

###  IV.1 Comparison of different trap potentials

The trap’s potential shape determines the particle density, which in turn
affects the probability of a photon to be absorbed by the particle medium. To
test the cooling sensitivity a set of models were studied with different
potentials: harmonic, linear, and quartic potentials. Comparison between
different trap shapes was made while keeping the potential energy at the
positions $x=\frac{L}{4}$, $x=\frac{3L}{4}$, almost identical. In addition, an
equal amount of particles was used for both simulations. The results are
presented in Figure 7. The most efficient cooling rate is predicted by the
harmonic trap,$1.45\cdot 10^{2}\text{K/s}$, which is almost by $50\%$ larger
than the particle cooling in a linear trap, $1.02\cdot 10^{2}\text{K/s}$. The
quartic $(\propto x^{4})$ potential shows a cooling rate of $37.1\text{K/s}$.
Due to the red shift of absorption with density the cooling is optimal when
there is a significant gradient in particle density such as in the harmonic
trap.

###  IV.2 Comparison of different densities

A direct connection between the average density in the trap and the cooling
rate was found. At low density a linear increase in the cooling rate is
observed Cf. Fig. 8. At higher densities the the cooling rate reduces. This is
in a density range which is considerably lower in comparison with the quantum
regime. In such a regime, the basic many body cooling phenomena should still
be valid but complimented by quantum corrections to the model.

Figure 6:  The particle phase space on the left (A,B) before coupling to the
radiation field (A) and after at time $t\approx 6*10^{-6}\text{s}$ at
$T=10^{-4}\text{K}$ (B). The right hand side represents the light phase space
before the coupling (C) and at time $t$ (D). The vertical axis of all the
figures describes the momentum, and the relevant scaling of the units is
presented on the left of the axis. Figure 7: Particle temperature as a
function of time, for different potentials. For a density of
$\rho=10^{14}\text{cm}^{-3}$; The trap potential:
$V_{harmonic}=\frac{1}{2}kx^{2}$ ; $V_{linear}=k|x|$;
$V_{quadratic}=\frac{1}{2}kx^{4}$. Figure 8:  The cooling rate in absolute
value as a function of the initial average particle density.

An estimation of the maximum cooling rate is obtained by considering a sphere
filled with a uniform gas of atoms with density $\rho$, and a typical
spontaneous emission lifetime of $\tau$. The bound for the cooling rate can be
calculated by noticing that only the atoms at the outer boundary of the sphere
effectively emit energy. The number of atoms per a unit area occupying the
outer shell is $N_{surf}\sim\rho^{2/3}$. On average the excited atoms out of
$N_{surf}$ will emit a blue shifted photon giving an energy difference of
$\hbar\Delta\omega$, with a rate of $\dot{Q}=\frac{\hbar\Delta\omega}{\tau}$.
Furthermore, if $1\%$ of the atoms on the surface are excited at a certain
instant, the upper bound to the cooling rate can be estimated as:

$R\sim 0.01\cdot\dot{Q}=0.01\frac{\hbar\Delta\omega}{k_{B}\tau}$ (12)

Inserting in equation 12 the data for $~{}^{87}\mbox{Rb}$:
$\rho=10^{13}\text{cm}^{-}3$, $\Delta\omega\leq 10^{7}\text{rad/s}$,
$\tau_{d}=27.7\text{ns}$ gives $R\sim 10^{6}\text{K/s}$. Where the
$\Delta\omega$ was estimated from Figure 8 considering the stochastic nature
of the process. Note that this is an upper bound, not taking into account the
stochastic nature of the cooling and possible heating sources. Comparing to
the current modeling, it predicts a smaller cooling rate of the order of
$R\sim 10^{3}\text{K/s}$.

The geometrical arguments can explain the dependence of the cooling rate on
the density. At low density the whole volume emits therefore a linear scaling
is expected as seen in Fig. 8. For high density asymptotic cooling rate should
scale as of $\rho^{2/3}$. For the data of the asymptotic cooling rate (Fig.
8), we obtain the scaling of $R\sim\rho^{0.67}$, in accordance with the
geometrical analysis.

##  V Enhanced Stochastic Cooling - An extension of the Stochastic Cooling
method

The Stochastic Cooling of $~{}^{87}\mbox{Rb}$ atoms, described in Section II,
utilizes the energy gap dependence on the inter-atomic distance. To generalize
this mechanism we propose a method applicable to different types of
constituents. The main idea incorporates additional control of the energy gap
between the ground and first excited potential energy surface.

The Stochastic Cooling method, (Sec. II and III) requires a spatial dependence
on the energy gap between the ground and excited states. $~{}^{87}\mbox{Rb}$
is a unique case, the ground and excited energy states scale differently with
the relative distance between the atoms, Cf. Fig. 3, resulting in an energy
gap with a sufficient spatial gradient. In the general case, both energy
states scale similarly, and the spatial gradient may not be sufficient to
achieve efficient cooling.

An enhancement of the spatial gradient between the energy states can be
induced by employing a second CW field with a frequency , $\omega_{1}$, in
resonance with the transition line between the excited state, $E_{e}$, to a
higher excited state, denoted by $E_{f}$. It is crucial that the frequency
$\omega_{1}$ should be different from the atomic transition, not affecting the
excitation process from the ground state.

##### Derivation:

The atomic energy gap between subsequent levels roughly scales as
$\propto\left({\frac{1}{n^{2}}-\frac{1}{(n+1)^{2}}}\right)$, demonstrating
that the first energy gap is much bigger than the other gaps. The big
difference between the energy gaps allows an explicit treatment of a two-
level-system coupled to an oscillating classical field Cohen-Tannoudji et al.
(1977). The solution is given in terms of the Rabi frequency , $\Omega$,
linearly dependent on the vector electric field amplitude of the laser,
$\vec{E_{\mathcal{L}}}$.

We will focus on two energy levels of the excited state, $\left|e_{1}\right>$
and the level of a higher excited state, $\left|e_{2}\right>$. The classical
radiation field induces an energy shift to the bare Hamiltonian levels, and
the new shifted states are given by

$E_{\pm}=\pm\frac{\hbar\sqrt{|\Omega|^{2}+\Delta^{2}}}{2}$ (13)

where $\Delta=\omega_{L}-(\omega_{f}-\omega_{h})$, can be neglected for a
resonant radiation, $\omega_{1}\equiv\omega_{f}-\omega_{e}\gg\Delta$.

$E_{\pm}=\pm\frac{\hbar\Omega}{2}$ (14)

For a classical radiation field, of frequency $\omega_{1}$, with a spatial
dependence, the intensity varies in the trap. For example, a high intensity
light focused at the centre of the trap will have a gradient toward lower
intensities at the edge of the trap. The intensity gradient results in a
spatial dependent Rabi frequency, $\Omega(x)$, where $x$ is the trap’s radial
coordinate. Concentrating on a pair of atoms in the trap, the Rabi frequency
can be written as a function of the inter-atomic distance, $r$. This leads to
an excited state $E_{e}$ which varies spatially as well, while the ground
state stays unperturbed by the classical EM field. This phenomena induces a
spatial energy dependent gap between the ground and excited states, with a
gradient depending on the EM intensity.

Once an energy gap is controlled by a tuning laser of frequency $\omega_{1}$,
a second laser of frequency $\omega_{0}=\omega_{e}-\omega_{g}$ is applied to
the particle medium. As in the Stochastic Cooling method, the laser of a
frequency $\omega_{0}$ induces excitations between the ground and excited
states. The relative random motion of the atoms in the excited state will
induce an average energy transfer and cooling. This process is controlled by
the electric field amplitude $\vec{E}_{\mathcal{L}}$ with a frequency
$\omega_{1}$.

Figure 9:  The excited state, $\left|e_{1}\right>$ , and the next energy
state, $\left|e_{2}\right>$, are coupled by a tuning laser of frequency
$\omega_{1}$ generating a Stark shift which is intensity dependent. The Stark
shift is enhanced at resonance conditions occurring at a specific inter-atomic
distances. As a result, the potential $\left|e_{1}\right>$ is modified and
with it the resonance conditions of the cooling laser
$\omega_{0}=\frac{E_{e}-E_{g}}{\hbar}$. By varying the intensity of the tuning
laser in the trap we obtain a gradient in the absorption probability of the
cooling laser.

The energy level configuration and the addition of a second tuning laser is
similar to the scenario utilized for Electromagnetic Induced Transparency,
(EIT). In EIT a combined AC-Stark splitting and quantum interference results
in a transparency at a frequency of the probe (cooling) laser. Similar
applications have been achieved for Rb and Pb. Such phenomena, under
favourable circumstances, will allow an easy penetration of the photons to a
partly opaque particle medium Xiao et al. (1995); Field et al. (1991).

The control of the energy gap gradient enables the Enhanced Stochastic Cooling
to operate at lower densities relative to the densities required for
Stochastic Cooling of Rubidium. For Stochastic Cooling of 87Rb, the high
densities are required because of the small dependence of the energy gap on
the inter-atomic distance. Once the gradient is engineered with an external
field, the energy gap can be shifted to longer inter-atomic distances. For
this method, the required density is bounded only by densities where the
photons are characterized by diffusional motion.

We present in table IX.4 the sufficient densities for applying Enhanced
Stochastic Cooling of different Alkali atoms appear .

Constituent |  Required density $[\text{cm}^{-3}]$
---|---
Rubidium |  $1.47\cdot 10^{6}$
Sodium |  $2.85\cdot 10^{6}$
Caesium |  $1.04\cdot 10^{6}$

###  V.1 Optimization of the tuning radiation field

In the following section we discuss how optimization of tuning laser intensity
and frequency affects the cooling.

####  V.1.1 Intensity variations of the tuning radiation field

The question arises what is the optimal field profile? The laser frequency is
determined by the gap $\omega_{1}=\omega_{f}-\omega_{e}$, but different
intensity profiles can be realized. Modern-day optics allow creating many
intensity profiles, utilizing optical holographic lenses, and state-of-the-art
optical devices. An optimization including all possible scenarios can be
complicated. However, we notice that by a similar derivation as in section
III.2.2 the diffusion constant $Dp_{l}$ is proportionate to the square of the
tuning laser’s electric field gradient. Proportionality suggests that a
largely varying field in the trap region will result in increased cooling to
the particle media, Cf. Fig. 11.

####  V.1.2 Frequency variation of the classical radiation field

The frequency of the second laser source can be tuned to induce cooling for a
specific atomic density. In the last section we consider a classical radiation
field, applied to the trap, with resonance frequency $\omega_{1}$ matching the
asymptotic transition $\left|e_{1}\right>\rightarrow\left|e_{2}\right>$.
However, the energy gap in resonance to the transition changes along the
inter-atomic distance, $r$. Modern experimental techniques allow accurate
control of the laser frequency and spatial intensity. This enables control of
the exact region of the inter-atomic distances which are coupled to the
cooling field. A laser with a frequency of $\omega_{1}(r_{i})$ will couple
between the ground and excited states in a region near $r_{i}$, inducing a
gradient in the energy gap between the states. The gradient will induce
cooling, originating from a pair of atoms with a certain inter-atomic
distance. Alternately, when averaging over the inter-atomic distances, the
laser frequency, $\omega_{1}(r_{i})$, will match an atomic medium of a density
$\rho(r_{i},x)$, (where $x$ is the radial component of the trap).

##  VI Results B: Enhanced Stochastic Cooling

Enhanced Stochastic Cooling was modelled on $~{}^{87}\mbox{Rb}$ atoms, in a
similar method as described in section III. A second tuning laser in
conjunction with the cooling laser is employed with a wavelength of $1475.6$
nm ($~{}^{2}P_{1/2}\rightarrow~{}^{2}D_{3/2}$). Once the light and particle
media were coupled, the two phase spaces were propagated in time, while
synchronizing the parameters after each time step. Energy transfer is assumed
to originate solely from coupling of the external field. The diffusion
functions are recalculated based on the derivation presented in Section V.

A number of different intensity profiles were studied, $E_{1}\propto
x^{2}$,$E_{2}\propto x^{4}$, and also a highly oscillating profile $E_{2}$. Cf
Fig. 10, 11. The sinusoidal profile shows the fastest cooling rate, $6.85\cdot
10^{2}\text{K/s}$ as a result of a large gradient. The other profiles $E_{2}$
and $E_{1}$ have an inferior cooling rate of, $1.33\cdot 10^{2}\text{K/s}$ and
$40\text{K/s}$, correspondingly.

Figure 10: Profiles of tuning light intensity: Black dashed $E_{1}=2\cdot
10^{14}x^{2}/2$, Red dash-dot $E_{2}=5\cdot 10^{12}x^{4}$, Blue solid
$E_{3}=4.96\cdot 10^{13}sin(200x)/200$ (MKS). The field profiles where
adjusted so that the maximum intensity in the trap will be equal. Figure 11:
Particle temperature as a function of time for different external field
profiles corresponding to figure 10.

##  VII Discussion

Cooling of neutral atoms via collective many-body interactions is an efficient
universal cooling scheme, applicable as a complementary method to prior
cooling methods, Kozyryev et al. (2016); Dieckmann et al. (1998), or
independently.

At sufficient particle density and large absorption cross sections photons are
trapped in the particle media. In such density regimes, the photon propagation
in the trap is characterized by diffusion. A single photon exhibits a large
number of excitation cycles, allowing energy and entropy transfer between
particles and photons. The mechanism proposed depends on the collective
behaviour of the particles and light media, giving rise to coordinated
dynamics. The cooling rate depends on particle density, on the density
gradient and the asymmetry in the spectrum between absorption and emission. On
average, an absorbed red photon will be emitted as a blue photon. Density of
similar values used in the demonstration have been achieved experimentally,
Camara et al. (2014); ferlaino12 et al. (2014); chomaz et al. (2017),
$\left({10^{13}-10^{15}\text{cm}^{-3}}\right)$. A simple experimental setup
for stochastic laser cooling requires the ability to change the trap potential
and with it the density. This will also allow to check the validity of the
theory.

The main thermodynamic principle unifying all cooling methods is an increase
of the total entropy of the joint particle and light ensembles. The energy
transfer from the particle medium to the light medium decreases the entropy of
the particle ensemble. This comes at the expense of the light entropy where
the constant radiation loss from the trap is the entropy generating mechanism.

Our modeling for Rubidium demonstrates that cooling can be achieved by the
’Stochastic Cooling method’. A single ’cooling laser’, coupling the ground and
lowest excited state ($\left|g\right>=X\,^{1}\Sigma_{g}^{+}$,
$a\,^{3}\Sigma_{u}^{+}$ ; $\left|e\right>=0_{g}^{-}$, $1_{g}^{-}$) is
sufficient. Modeling predicts an efficient cooling for a density range of
$(1-7)\cdot 10^{14}\text{cm}^{-3}$, cooling rates ranging between
$100-800\text{K/s}$. Furthermore, the model predicts that the asymptotic
cooling rate for a fixed trap volume will scale with the density as
$\rho^{2/3}$ and be linear at low density. This is in accordance with a simple
geometrical model.

This scheme can be extended to cool other constituents by adding a second
’tuning laser’ coupling the first state and a higher excited state. The
generalized cooling method, ’Enhanced Stochastic Cooling’ , allows universal
extension of cooling for different types of neutral atoms as well as
molecules. The experimental requirement is an additional CW laser. The rate of
cooling can be controlled by determining the values of the gradient of the
intensity of the tuning laser, $|\frac{\partial E_{tune}}{\partial x}|$, Merli
et al. (2009). The generalized scheme predicts efficient cooling rates that
can be maximized by choosing an intensity profile with a large spatial
gradient in the trap, Cf. Fig. 11.

The phenomena, enabling cooling of Rubidium (Stochastic Cooling), is related
to pressure line shift and pressure broadening, which arises at sufficiently
high densities Vdović et al. (2006). The pressure shift and inter-atomic
interactions allow energy transfer between particle translational degrees of
freedom to internal degrees of freedom and to the photonic medium. The
magnitude of the pressure shift influences the cooling rate directly. Similar
effects, arising for increased densities, can be seen in other condensed
matter phenomena, such as charge transfer to solvent and modifications to
absorption/emission spectra for liquid phase relative to a gas phase spectrum.
For cooling, higher particle densities increase the pressure shift and, in
turn, the gradient of the energy gap between the ground and excited states. As
a result of the stochastic nature of the process, a larger gradient leads to
faster cooling.

The mechanisms described, responsible for energy transfer from the particle to
the light medium, is valid for the semi-classical regime. At low temperatures,
the present theory should be modified by a quantum theory. The crossover
temperature is when the de Broglie thermal wavelength is in the range of the
mean particle distance ($T~{}>~{}\frac{h^{2}\rho^{2/3}}{2\pi mk_{B}}$). For
$~{}^{87}\text{Rb}$ this temperature is $\sim 10^{-6}\text{K}$ for a density
of $\rho=10^{14}\text{cm}^{-3}$. At low temperatures the asymmetry between red
and blue shift emission is larger. The reason is that the ground state density
is peaked at the attractive region of the van der Waals potential and the
emission is biased toward the outer turning point of the vibration of the
excited potential which is larger. Near the BEC limit, additional corrections
can be made based on particle wave characteristics. These are out of the scope
of this paper. The Enhanced Stochastic Cooling method is applicable to low
densities and lower temperature regimes. The details are presented in Table
IX.4.

Manipulation of cold molecules and cooling molecules to extremely low
temperatures has been one of the main focal points of the Atomic Molecular
Optical research field, Stuhl et al. (2012); Bochinski et al. (2003);
Narevicius et al. (2009); Zeppenfeld et al. (2012); Bethlem et al. (2000);
Bethlem and Meijer (2003); Jochim et al. (2003). The method proposed is
applicable to molecular cooling experiments. For sufficient densities, light
can be trapped for long time periods in the molecular medium, and an energy
transfer is predicted.

In contrast to the general simplicity of laser cooling atoms, the higher
number of degrees of freedom in molecules induces a complex internal energy
structure. These features complicate the cooling process due to additional
relaxation channels which leads to induced heating. For efficient cooling the
molecules need to have large diagonal Franck-Condon factors. This will allow
repeated electronic excitations while minimizing the excitation of vibrational
states. Low inelastic rates are favourable which result in heating and a fast
molecular loss rate. In addition, the energy transitions should match the
available laser cooling frequencies. A number of different constituents
qualify for efficient Stochastic cooling or Enhanced Stochastic cooling, ${\rm
OH}$,${\rm CaF}$ and ${\rm YO}$.

${\rm CaF}$ has been cooled to velocities of $\approx 10\pm 4\text{m/s}$,
which is below the capture velocity of a molecular MOT Hemmerling et al.
(2016), and has suitable electronic transitions. The cooling laser can be
applied, detuned slightly below the $X^{2}\Sigma_{1/2}^{+}$, $v=0$, $N=-1$,
$A^{2}\Pi_{1/2}$, $v=0,~{}J=+3/2$ transition of $606\text{nm}$, while the
coupling laser couples between $A^{2}\Pi_{1/2}$, $v=0~{},J=+3/2$ and
$C^{2}\Pi_{1/2}$, $v=0,~{}J=-1/2$ states of $729.5\text{nm}$. Additional
lasers and a magnetic field may be required in order to bring back dark
magnetic sub-states to the optical cycle and reduce population loss to excited
vibrational states. Such a scheme is envisaged to induce Enhanced Stochastic
cooling. A spatial gradient in the energy gap between the ground state and
excited state allows efficient cooling. The gradient is created by the
polarizability difference between the two states. This means that a dense
ensemble of ${\rm CaF}$ trapped in a MOT could be cooled further, towards sub-
millikelvin temperatures.

Other molecular candidates are ${\rm OH}$ radical and ${\rm YO}$ which have
been confined in a trap Sawyer et al. (2007); Stuhl et al. (2012); Hummon et
al. (2013) and studied in context with optical cooling Yeo et al. (2015);
Collopy et al. (2015); Bochinski et al. (2003). They have a suitable internal
energy structure with convenient optical transitions in the visible light.
Both show similarities to $\rm CaF$, allowing application of Stochastic
Cooling methods. Similarly to the case of $\rm CaF$, additional pumping lasers
and a magnetic field could be needed to ensure a closed cooling cycle. These
additional lasers prevent population trapping in dark states by re-pumping
back to the cooling cycle.

Almost all molecules are more polarizable in the excited state. As a result, a
gradient in the energy gap will arise and the Stochastic Cooling method is
therefore applicable since the Stochastic Cooling method is based on van der
Waals forces. Molecules with higher polarizability could be cooled more
efficiently. The molecular density should be in the range where two-body
elastic collisions are dominant over three-body inelastic ones. To conclude,
’Enhanced Stochastic Cooling’ can greatly increase the class of atoms and
molecules that can be cooled to sub-Kelvin temperatures.

##  VIII Acknowledgements

We thank Amikam Levi and Shimshon Kallush for fruitful discussions and Maya
Dann and Moshe Armon for their help. This work was partially supported by the
ISF - Israeli Science Foundation

##  IX Appendix

###  IX.1 Rubidium data table

D1$\,\mbox{transition }(5^{2}S_{1/2}\rightarrow 5^{2}P_{1/2})$ Steck (2001)

Wavelength (Vacuum) |  794.979 nm
---|---
Lifetime |  27.7 ns
Recoil Energy |  22.8236 kHz
Effective Far-Detuned Saturation Intensity |  $4.484\text{mW/cm}^{2}$
Effective Far-Detuned Resonant Cross section |  $1.082\cdot 10^{-9}\text{cm}^{2}$

###  IX.2 Model parameters

Density Range |  $10^{13}-10^{14}$ $\text{cm}^{-3}$
---|---
Trap Length |  1 mm
Cooling Laser Frequency |  $3.7711\cdot 10^{14}\text{s}^{-}1$
Particle Photon Ratio in The Trap |  1
Initial Temperature |  $0.01\text{K}$

###  IX.3 Pulse Parameters, Cf Sec. III.2.2

Variance |  $10^{-16}\text{s}^{2}$
---|---
Amplitude |  Normalized so the cross section will fit the experimental value

###  IX.4 Numerical methods

The dynamical equations to be solved for the particles and light Eq. 1 and 2
have the structure:

$\frac{\partial}{\partial t}F(x,p,t)={\bf O}F(x,p,t)$ (15)

where $F(x,p,t)$ is the probability function in phase space and ${\bf O}$ is a
differential operator defined in Eq. 1 and 2 . We can write a formal solution
for a short time step $\Delta t$:

$F(x,p,t+\Delta t)\approx e^{{\bf O}\Delta t}F(x,p,t)$ (16)

The exponent in Eq. 16 is expanded by a Chebychev polynomial of order $N$
Kosloff (1994):

$e^{{\bf O}\Delta t}\approx\sum_{k=0}^{N}C_{k}(\Delta t)T_{k}({\bf O})$ (17)

where $C_{k}$ are expansion coefficients (Bessel functions) and $T_{k}(x)$ is
the Chebychev polynomial of order $k$.

The following table presents the details of the Fourier-Chebychev numerical
scheme for solving coupled Fokker-Planck equations.

Grid size |  $10^{-3}$ m
---|---
Number of spatial grid points in position | 250
Number of spatial grid points in momentum | 500
Grid spacing |  $\Delta r=4\cdot 10^{-6}$ m
Grid spacing |  $\Delta p=4.8491\cdot 10^{-27}$ kg m/s
Order of Chebychev Polynomial (light) | 45
Grid size |  $2\cdot 10^{-3}$ m
Grid points of light medium position | 500
Grid points of light medium momentum | 100
Grid spacing (light) |  $\Delta r=4\cdot 10^{-6}$ m
Grid spacing (light) |  $\Delta p=4.8395\cdot 10^{-36}$ kg m/s
Typical time step (light) |  $10^{-9}$ s
Number of time steps | 6000
Order of Chebychev Polynomial (light) | 271

###  IX.5 Estimate the number of excitations for a single photon

The photon propagation through the atomic medium can be modeled as a 3D random
walk, resulting from repeated absorption/emission cycles. The square of the
distance that a photon reaches after N steps, or variance is;

$var_{3d}=N\varepsilon^{2}$ (18)

where $N$ is the number of absorption/emission cycles and $\varepsilon$ is the
length of each step between consecutive absorption events. Assuming a
spherical trap with a uniform density;

$\varepsilon=P_{abs}^{-1}\cdot\rho^{-\frac{1}{3}}$ (19)

where $P_{abs}=\rho^{\frac{2}{3}}\sigma$ is the absorption probability. (For
$\rho^{\frac{2}{3}}\sigma\leq 1$ the equality holds.)

To escape the trap, the photon has to reach a distance of $R$ (trap radius)
from the center of the trap.

$N=\frac{R^{2}}{\left({P_{abs}^{-1}\cdot\rho^{-\frac{1}{3}}}\right)^{2}}=R^{2}\rho^{-2}\sigma^{-2}$
(20)

$\rho(N)=\left({NR^{-2}\sigma^{2}}\right)^{-\frac{1}{2}}$ (21)

###  IX.6 Absorption probability function

We solve the transition probability between the ground and excited state of a
quantum system of two $~{}^{87}\mbox{Rb}$ atoms and a light field
characterizing a single photon.

The original Hamiltonian, with no coupling to a radiation field is:

$\hat{H}_{g/e}=\hat{T}+\hat{V}_{g/e}=\frac{\hat{P}}{2m}+V_{g/e}(\mathbf{r})$
(22)

In the presence of an electromagnetic field the two surfaces of the ground and
excited states are coupled by the interaction of the field and the dipole
momentum operator. The new Hamiltonian is written as:

$\hat{H}=\hat{H_{g}}\otimes\hat{P_{-}}+\hat{H_{e}}\otimes\hat{P_{+}}+\varepsilon(t)\hat{\mu}\otimes\hat{S}_{+}+\varepsilon(t)\hat{\mu}\otimes\hat{S}_{-}=\bigg{[}\begin{array}[]{cc}\hat{H_{e}}&\varepsilon(t)\hat{\mu}\\\
\varepsilon*(t)\hat{\mu}&\hat{H_{g}}\end{array}\bigg{]}$ (23)

The electromagnetic field is given by;

$\varepsilon(t)=\bar{\varepsilon}(t)e^{-i\omega_{L}t}+\bar{\varepsilon}*(t)e^{i\omega_{L}t}$
(24)

where $\omega_{L}$ is the laser carrier frequency and $\bar{\varepsilon}(t)$
the envelope of the pulse. After the rotating wave approximation, the
Hamiltonian reduces to;

$\hat{H_{S}}=\bigg{[}\begin{array}[]{cc}\hat{H_{e}}-\hbar\omega_{L}/2&\bar{\varepsilon}(t)\hat{\mu}\\\
\bar{\varepsilon}(t)\hat{\mu}&\hat{H_{g}}+\hbar\omega_{L}/2\end{array}\bigg{]}$
(25)

The amplitude of absorption of a photon is calculated considering a system
following dynamics governed by $\hat{H}_{S}$. For the basis states
$(\\{\psi_{k}\\})$ of $\hat{H}_{S}$ , the amplitude transfer from the
eigenstate $\left|\psi_{i}\right>$ to eigenstate $\left|\psi_{n}\right>$,
after time $t$, is given by the time dependent perturbation theory. Assuming a
weak field the amplitude, to good approximation, is given by the first order
term;

$b_{n}^{(1)}(t)=-\frac{i}{\hbar}\int_{0}^{t}{e^{i\omega_{ni}t}\hat{W_{ni}}\left({t^{\prime}}\right)dt}$
(26)

where $\hat{W}_{ni}(t)$ is the time perturbation term;
$\omega_{ni}=\frac{E_{n}-E_{i}}{\hbar}$. Defining

$c_{n}(t)=b_{n}(t)e^{-iE_{n}t/\hbar}$ (27)

With the help of Eq. 26

$c_{n}(t)=-\frac{i}{\hbar}\int_{0}^{t}{d\tau
e^{-iE_{n}(t-\tau)/\hbar}\hat{W_{ni}}(t^{\prime})e^{-iE_{i}\tau/\hbar}}$ (28)

For a perturbation $\hat{W_{ni}}(t^{\prime})=\bar{\varepsilon}(t)\hat{\mu}$,
and the following identity;
$e^{-iE_{n}(t-\tau)/\hbar}\left|\psi_{n}\right>=e^{-i\hat{H_{n}}(t-\tau)/\hbar}\left|\psi_{n}\right>$

$c_{n}(t)=-\frac{i}{\hbar}\int_{0}^{t}{d\tau
e^{-iE_{n}(t-\tau)/\hbar}\bar{\varepsilon}(t)\hat{\mu}(t^{\prime})e^{-iE_{i}\tau/\hbar}}$
(29)

Similarly, for the excited state, using Eq. 25.

$\left|\psi_{e}(t)\right>=-\frac{i}{\hbar}e^{-\frac{i}{\hbar}(\hat{H}_{e}-\hbar\omega_{L}/2)t}\int_{0}^{t}{d\tau
e^{\frac{i}{\hbar}\hat{H_{e}}\tau}e^{-i\hbar\omega_{L}\tau}\hat{\mu}\bar{\varepsilon}(\tau)e^{-\frac{i}{\hbar}\hat{H_{g}}\tau}\left|\psi_{g}(0)\right>}$
(30)

Assuming a narrow Gaussian pulse which centred at $t=0$ far away from the
source, the integral boundaries can be taken to infinity.

$\bar{\varepsilon}(\tau)=\frac{B}{\sqrt{2\pi\sigma_{t}^{2}}}e^{-\frac{\tau^{2}}{2\sigma_{t}^{2}}}$
(31)

where $B$ is an amplitude constant. Using the identity Ashkenazi et al. (1997)
(Sec. III)

$\sigma_{A}(\omega_{L})\propto\left<\psi_{i}\right|{\hat{A}}{\left|\psi_{i}\right>}\rightarrow\hat{A}=\hat{\mu}\int_{-\infty}^{\infty}d\tau
e^{\frac{i}{\hbar}(\hat{H_{e}}-\hat{H_{g}}-\hbar\omega_{L})\tau}\hat{\mu}$
(32)

Assuming the system is in the ground state at the initial time, the propagator
is given by;

$\hat{A}=\hat{\mu}\int_{-\infty}^{\infty}d\tau\bar{\varepsilon}(\tau)e^{\frac{i}{\hbar}(\hat{H_{e}}-\hbar\omega_{L})\tau}e^{-\frac{i}{\hbar}E_{g}\tau}\hat{\mu}$
(33)

When the transition dipole moment is constant in postion and momentum, the
solution of the integral gives;

$\hat{A}=\hat{\mu}^{2}Be^{-\frac{\sigma_{t}^{2}\Delta^{2}}{2}}$ (34)

where $\Delta=\frac{1}{\hbar}\hat{H_{e}}-\omega_{g}-\omega_{L}$ Defining,
$\alpha=\frac{\sigma_{t}^{2}}{2}$, as well as emitting the global phase from
the expression and defining $c=(\omega_{g}+\omega_{L})$ and decomposing the
Hamiltonian to the kinetic and potential terms, $H_{e}=\hat{T}+\hat{V}_{e}$.
The expression for the propagator is given by;

$\sigma_{A}\propto{\langle\psi_{g}|}{e^{-\frac{\sigma_{t}^{2}}{2}(\hat{G}+\hat{F}+\hat{K})}|\psi_{g}\rangle}$
(35)

where:

$\hat{G}=\hat{T}^{2}-2\hbar c\hat{T}$ (36)

$\hat{F}=\hat{V_{e}}^{2}-2\hbar c\hat{V_{e}}$ (37)

$\hat{K}=\hat{T}\hat{V_{e}}+\hat{V_{e}}\hat{T}$ (38)

Using the Zassenhaus formula to expand the exponent, we find that the high
order commutators can be neglected Magnus (1954).

Taking the first term in the expansion, the cross section can be summarized by
the expression;

$\sigma_{A}\propto{\langle\psi_{g}|}{e^{-\frac{\sigma_{t}^{2}}{2}\hat{G}}e^{-\frac{\sigma_{t}^{2}}{2}\hat{K}}e^{-\frac{\sigma_{t}^{2}}{2}\hat{F}}|\psi_{g}\rangle}$
(39)

Where $\hat{G}$ and $\hat{F}$ are the kinetic and potential energy terms
correspondingly, and $\hat{K}$ is a correlation term.

The solution is a Gaussian function with a variance of $2.89$
$(\text{nHz})^{2}$, Cf. Table IX.3, centred around $(V_{e}-V_{g})/\hbar$, for
large $r$’s the contribution of the van Der Waals interactions to the
probability to be absorbed is negligible. However, for a short range the
interaction will shift the resonance frequency towards lower frequencies in
comparison with the atomic transition line, influencing the optimized detuning
from resonance, $\Delta\omega_{optimise}$, used for optimal cooling.

##### Direct cross section calculation:

We can decompose the initial thermal state to random phase Gaussian wave
functions, when assuming a contribution of the kinetic term only. The
approximation is valid for the density and temperature regime in our
experiment, Cf. Table IX.1, where the ground state potential has a minor
effect on the wave function of the ground state.

$e^{-\frac{\hat{H}}{k_{B}T}}\approx e^{-\frac{p^{2}}{2mk_{B}T}}$ (40)

Each thermal Gaussian wave function, in the momentum representation, has a
temperature dependent standard deviation, $\sigma=\sqrt{mk_{B}T}$ and an added
random phase $G(p)=e^{-\frac{p^{2}}{2mk_{B}T}+ipR_{0}}$.

In the position representation this amounts to an ensemble of Gaussians
centred at different locations $\\{R_{0}\\}$. In the final stage of the
calculation all Gaussians are summed and averaged, the random phases cancel
one another constructing the asymptotic thermal state propagated in time.

The overall effect of the described calculation is equivalent to the following
process: Each Gaussian, centred at a different location, is coupled to an
electric field at time $\tau$, the EM field couples the ground and excited
states resulting in a population transfer to the excited state. The excited
state is then propagated until time $t$ to achieve a single realization. The
overall excited state is then achieved by integrating on all possible
transition times, $\tau$. The calculation converges to the first order time
perturbation term assuming a weak pulse.

This process is repeated for different laser frequency shifts, $\Delta\omega$,
and an absorption probability distribution function dependent on the laser
frequency shift is achieved.

###  IX.7 Relation between the cross section and the matrix element

Deriving the proportionality
$\sigma_{A}(\omega_{L})\propto{\langle\psi_{i}|}{\hat{A}|\psi_{i}\rangle}$.
The power can be written as

$P=\frac{dE}{dt}=\langle{\frac{dH}{dt}}\rangle$ (41)

Making an ansatz of Eq. 23

$P=\langle{\frac{d\varepsilon(t)}{dt}\hat{\mu}\otimes\hat{S}_{+}+\frac{d\varepsilon^{*}(t)}{dt}\hat{\mu}\otimes\hat{S}_{-}}\rangle=\langle{\frac{d\varepsilon(t)}{dt}\hat{\mu}\otimes\left|\psi_{e}\right>\left<\psi_{g}\right|+\frac{d\varepsilon^{*}(t)}{dt}\hat{\mu}\otimes\left|\psi_{g}\right>\left<\psi_{e}\right|}\rangle$
(42)

Inserting the density matrix expression;
$\rho=\frac{1}{2}(\left|\psi_{g}\right>\left<\psi_{g}\right|+\left|\psi_{e}\right>\left<\psi_{e}\right|)$
: The state is written as:

$\displaystyle\frac{d\varepsilon(t)}{dt}\left<\psi_{e}\right|\hat{\mu}\left|\psi_{g}\right>+\frac{d\varepsilon^{*}(t)}{dt}\left<\psi_{g}\right|\hat{\mu}\otimes\left|\psi_{e}\right>$
$\displaystyle=-2\text{Real}\left({\frac{d\varepsilon(t)}{dt}\left<\psi_{e}\right|\hat{\mu}\left|\psi_{g}\right>}\right)=-2\text{Real}\left({\frac{d\varepsilon(t)}{dt}\langle{\hat{\mu}\otimes\hat{S}_{+}}\rangle}\right)$

The power at time $t$ ;

$P(t)=-2\text{Real}\left({\frac{d\varepsilon(t)}{dt}\left<\psi_{e}(t)\right|\hat{\mu}\left|\psi_{g}(t)\right>}\right)\propto\left<\psi_{g}(t)\right|\hat{\mu}\int_{-\infty}^{\infty}d\tau
e^{\frac{i}{\hbar}(\hat{H_{e}}-\hat{H_{g}}-\hbar\omega_{L})\tau}\hat{\mu}\left|\psi_{g}(t)\right>$
(43)

The power is proportionate to the population change which has a linear
dependency on the cross section

$P=\hbar\omega_{0}\frac{dN_{e}}{dt}\propto\sigma(\omega_{L})$ (44)

combining Eq. 43 and 44 we get the desired relation;

$\sigma_{A}(\omega_{L})\propto\left<\psi_{g}(t)\right|\hat{\mu}\int_{-\infty}^{\infty}d\tau
e^{\frac{i}{\hbar}(\hat{H_{e}}-\hat{H_{g}}-\hbar\omega_{L})\tau}\hat{\mu}\left|\psi_{g}(t)\right>=\left<\psi_{g}(t)\right|\hat{A}\left|\psi_{g}(t)\right>$
(45)

###  IX.8 Energy transfer between the atom and radiation field and calculation
of $\mathcal{D}(\rho_{par},T_{par})$

The energy change due to a typical excitation is:

$\Delta
E=-\frac{C_{3}}{r_{f}^{3}}+\frac{C_{6}}{r_{f}^{6}}+\frac{C_{3}}{r_{i}^{3}}-\frac{C_{6}}{r_{i}^{6}}=C_{3}\left({\frac{1}{r_{i}^{3}}-\frac{1}{r_{f}^{3}}}\right)+C_{6}\left({\frac{1}{r_{f}^{6}}-\frac{1}{r_{i}^{6}}}\right)$
(46)

where $r_{i}$ is the inter-atomic distance for time $t$ when the photon is
absorbed and $r_{f}$ is the relative distance at time $t+\delta t$, when the
photon is emitted. Where $\delta t$ is the typical decay time for the Rubidium
87 D1 transition.

Transforming to the center of mass and relative velocity coordinates the
velocity distribution is a Maxwell Boltzmann distribution of particles with a
reduced mass $\mu=m/2$, a kinetic energy of
$E_{k}=\frac{\langle{{\vec{{p_{r}}}}^{2}}\rangle}{2\mu}$ , momentum
$\vec{p_{r}}=\mu\cdot\vec{v}$ , and relative velocity $\vec{v}$. The velocity
distribution for the relative particle,

$f(v)=\sqrt{\frac{\mu}{2\pi k_{B}T}}e^{-\frac{\mu v^{2}}{2k_{B}T}}$ (47)

The initial relative distance is assumed to be the mean distance for a density
$\rho(x)$ , where $x$ is the spatial position in the trap.

$r_{i}=(\rho(x))^{-1/3}$ (48)

The final relative atomic distance, $r_{f}$ , can be written in terms of the
relative velocity $v$ ; $r_{f}=r_{i}+v\cdot\delta t$. Making an ansatz of Eq.
48.

$\Delta E=C_{3}\left({\rho(x)-\frac{1}{(r_{i}+v\cdot\delta
t)^{3}}}\right)+C_{6}\left({\frac{1}{(r_{i}+v\cdot\delta
t)^{6}}-\rho^{2}(x)}\right)$ (49)

Since $r_{i}\gg v\cdot\delta t$ (in the density range discussed) we can expand
in a Taylor series up to the first term.

$\frac{1}{(r_{i}+v\cdot\delta
t)^{n}}\approx\frac{1}{r_{i}^{n}}\cdot\left({1-n\frac{v\cdot\delta
t}{r_{i}}}\right)$ (50)

The energy gap is reduced to;

$\Delta E=3(\rho(x))^{4/3}v\cdot\delta t(C_{3}-2\rho(x)C_{6})=C\cdot v$ (51)

$C=3(\rho(x))^{4/3}\delta t(C_{3}-2\rho(x)C_{6})$ (52)

In the first order approximation the energy change and the relative velocity
are linearly dependent. The distribution function in velocity translates to an
energy distribution function, for $E=\hbar(\omega_{f}-\omega_{i})$.

$f(E=\hbar\omega_{f})=N_{norm}e^{-\frac{\mu\hbar^{2}(\omega_{f}-\omega_{i})^{2}}{2C^{2}k_{B}T}}$
(53)

$N_{norm}=\sqrt{\frac{1}{\pi}\frac{\mu}{2C^{2}k_{B}T}}$ (54)

The variance of the function $f(\hbar\omega_{f})$ can be used to calculate the
light phase space diffusion variable of energy transfer,
$\mathcal{D}_{E}(\rho_{par},T_{par})=\frac{var(f(E))}{\delta t}$ , arising
from the interaction of the particle and photons, including only photons which
are absorbed.

$Var(f(E))=\frac{C^{2}k_{B}T}{\mu}$ (55)

$C=3(\rho(x))^{4/3}\delta t(C_{3}-2\rho(x)C_{6})$ (56)

$\mathcal{D}_{E}(\rho,T)=\frac{[3(\rho(x))^{4/3}(C_{3}-2\rho(x)C_{6})]^{2}\delta
t\cdot k_{B}T}{\mu}$ (57)

For the diffusion in momentum using the photon energy relation, $E=p\cdot c$ ,
the diffusion variable in momentum is given by:
$\mathcal{D}=\frac{1}{c^{2}}\mathcal{D}_{E}$

###  IX.9 Random phase approach for calculating the absorption probability
function

We use the Random phase approach as an efficient scheme for propagating a
thermal state $\hat{\rho}$. A thermal state is an incoherent state which
undergoes coherent time evolution. In this case a direct approach is a full
solution of the Liouville von Neumann equation in the Schrödinger picture,

$i\hbar\frac{\partial\rho}{\partial t}=[H,\rho]$ (58)

For a time evolution operator $\hat{U}(t)=e^{-\frac{i}{\hbar}\hat{H}t}$ the
dynamics can be captured by the equation
$\rho(t)=\hat{U}(t,0)\rho(0)\hat{U}^{\dagger}(t,0)$. When a wide range of
energy states are populated the direct solution of the initial state can be
difficult and time consuming. An alternative approach decomposes the initial
thermal state to random phase Gaussian wave functions. Time evolution can be
calculated on each realization and averaged to assemble the thermal state at
time $t$. A detailed description follows. For a high number of realizations
the random phases cancel each other leaving no effect on the the desired
calculation. This is the underlying principle of the method. For a general
random phase $e^{i\theta_{\alpha}}$ where $N\gg 1$ we can write the Cronicer
delta function as:

$\frac{1}{N}\sum_{k=1}^{N}e^{i(\theta_{\alpha}^{k}-\theta_{\beta}^{k})}=\delta_{\alpha\beta}$
(59)

$k$ labels a set random angle, each angle given for each basis state ,
$\alpha$ and $\beta$. If $\alpha=\beta$, $k_{\alpha}=k_{\beta}$ for all $k$ ,
we get the unity, for any other case the equality converges to zero as
$\frac{1}{\sqrt{N}}$. This characteristic allows a composition of the operator
with an arbitrary complete orthonormal basis $\\{\left|\alpha\right>\\}$ and
the random phases $\\{e^{i\theta_{\alpha}^{k}}\\}$ . We define a thermal
random wave function
$\left|\psi_{\alpha}^{k}\right>=e^{i\theta^{k}_{\alpha}}\left|\alpha\right>$
and an accumulated wave function
$\left|\Psi^{k}\right>=\sum_{\alpha}\left|\psi_{\alpha}^{k}\right>=\sum_{\alpha}e^{i\theta^{k}_{\alpha}}\left|\alpha\right>$

$\hat{1}=\frac{1}{N}\sum_{k=1}^{N}\left|\Psi^{k}\right>\left<\Psi^{k}\right|=\sum_{\alpha,\beta}\left|\alpha\right>\left<\beta\right|\frac{1}{N}\sum_{k=1}^{N}e^{i(\theta_{\alpha}^{k}-\theta_{\beta}^{k})}=\frac{1}{N}\sum_{k=1}^{N}\sum_{\alpha,\beta}\left|\psi_{\alpha}^{k}\right>\left<\psi_{\beta}^{k}\right|$
(60)

Therefore the thermal state at time $t=0$ is;

$\displaystyle\hat{\rho}=\hat{\rho}\cdot\hat{1}=\frac{1}{Z}e^{-\frac{\hat{H}\beta}{2}}e^{-\frac{\hat{H}\beta}{2}}\sum_{\alpha,\beta}\left|\alpha\right>\left<\beta\right|\frac{1}{N}\sum_{k=1}^{N}e^{i(\theta_{\alpha}-\theta_{\beta})}$
$\displaystyle=\frac{1}{Z}\frac{1}{N}\sum_{k=1}^{N}\sum_{\alpha,\beta}e^{i\theta_{\alpha}^{k}}e^{-\frac{E_{\alpha}\beta}{2}}\left|\alpha\right>\left<\beta\right|e^{-\frac{E_{\beta}\beta}{2}}e^{-i\theta_{\beta}^{k}}=\frac{1}{Z}\frac{1}{N}\sum_{k=1}^{N}\left|\varphi^{k}\right>\left<\varphi^{k}\right|$

while the thermal random wave functions are
$\left|\varphi^{k}\right>=\sum_{\alpha}e^{-\frac{E_{\alpha}\beta}{2}+i\theta_{\alpha}^{k}}\left|\alpha\right>$
, and the temperature dependence given by $\beta=\frac{1}{k_{B}T}$ .

The thermally averaged time dependent states , $\hat{\rho}(t)$ , can be
calculated by the same process, decomposed to time dependent thermal random
wave functions $\left|\varphi^{k}(t)\right>$. We obtain the thermal state
$\hat{\rho}(t)$ by propagating $N$ accumulated thermal random function,
$\left|\varphi^{k}(t)\right>=\hat{U}(t,0)\left|\varphi^{k}\right>$ , and
taking an average defined by equation IX.9. Taking a closer look at a single
thermal random wave function
$\left|\varphi_{\alpha}^{k}\right>=e^{-\frac{E_{\alpha}\beta}{2}+i\theta_{\alpha}^{k}}\left|\alpha\right>$,
$\\{\left|\alpha\right>\\}$ are chosen to be the momentum state basis. For
small potential energy the state can be written as;

$e^{-\frac{\hat{H}\beta}{2}+i\theta_{p}^{k}}\left|p\right>\approx
e^{-\frac{p^{2}\beta}{2m}+i\theta_{p}^{k}}\left|p\right>=e^{-\frac{p^{2}}{2mk_{B}T}+ipR_{0}}\left|p\right>$
(61)

having defined $\theta_{p}^{k}=pR_{0}$ in the second equalization.

The thermal random state, for the range of high kinetic energy or weak
interactions, is a thermal Gaussian with a variance $mk_{B}T$ and an
additional random phase. In the position representation the wave function has
a form of a Gaussian displaced by $R_{0}$,
$\left|\varphi_{r}^{k}\right>=e^{-\frac{1}{2}mk_{B}T(r-R_{0})^{2}}$ In the
position representation the random phase approach leads to a decomposition of
the initial thermal state to many thermal Gaussian wave functions centered
randomly in space. The validity of such an approximation for a two-body
interaction holds only for large $r$ where the potential is weak.

## References

  * Kepler (1619) J. Kepler (1619).
  * Maxwell (1862) J. C. Maxwell, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 23, 12 (1862).
  * Einstein (1909) A. Einstein, On a Heuristic Point of View about the Creation and Conversion of Light 1 On the Electrodynamics of Moving Bodies 10 The Development of Our Views on the Composition and Essence of Radiation 11 The Field Equations of Gravitation 19 The Foundation of the Generalised Theory of Relativity 22 p. 11 (1909).
  * Einstein (1917a) A. Einstein, also Physikalische Zeitschrift 18, 121 (1917a).
  * Einstein (1917b) A. Einstein, Physikalische Zeitschrift 18 (1917b).
  * Metcalf and Van der Straten (2012) H. J. Metcalf and P. Van der Straten, _Laser cooling and trapping_ (Springer Science & Business Media, 2012).
  * Scovil and Schulz-DuBois (1959) H. Scovil and E. Schulz-DuBois, Physical Review Letters 2, 262 (1959).
  * Wineland and Dehmelt (1975) D. Wineland and H. Dehmelt, in _Bulletin of the American Physical Society_ (AMER INST PHYSICS CIRCULATION FULFILLMENT DIV, 500 SUNNYSIDE BLVD, WOODBURY, NY 11797-2999, 1975), vol. 20, pp. 637–637.
  * Hänsch and Schawlow (1975) T. W. Hänsch and A. L. Schawlow, Optics Communications 13, 68 (1975).
  * Steck (2001) D. A. Steck, _Rubidium 87 d line data_ (2001).
  * Steck (2000) D. A. Steck, Report, Los Alamos National Laboratory, Los Alamos 124 (2000).
  * Kasevich and Chu (1992) M. Kasevich and S. Chu, Physical review letters 69, 1741 (1992).
  * Ketterle et al. (1992) W. Ketterle, A. Martin, M. A. Joffe, and D. E. Pritchard, Physical review letters 69, 2483 (1992).
  * Diedrich et al. (1989) F. Diedrich, J. Bergquist, W. M. Itano, and D. Wineland, Physical Review Letters 62, 403 (1989).
  * Kerman et al. (2000) A. J. Kerman, V. Vuletić, C. Chin, and S. Chu, Physical review letters 84, 439 (2000).
  * Dalibard and Cohen-Tannoudji (1985) J. Dalibard and C. Cohen-Tannoudji, JOSA B 2, 1707 (1985).
  * Dalibard and Cohen-Tannoudji (1989) J. Dalibard and C. Cohen-Tannoudji, JOSA B 6, 2023 (1989).
  * Neuhauser et al. (1978) W. Neuhauser, M. Hohenstatt, P. Toschek, and H. Dehmelt, Physical Review Letters 41, 233 (1978).
  * Leibfried et al. (2003) D. Leibfried, R. Blatt, C. Monroe, and D. Wineland, Reviews of Modern Physics 75, 281 (2003).
  * Boozer et al. (2006) A. Boozer, A. Boca, R. Miller, T. Northup, and H. J. Kimble, Physical review letters 97, 083602 (2006).
  * Metcalf and Straten (2007) H. J. Metcalf and P. Straten, _Laser cooling and trapping of neutral atoms_ (Wiley Online Library, 2007).
  * Breeden and Metcalf (1981) T. Breeden and H. Metcalf, Physical Review Letters 47, 1726 (1981).
  * Anderson et al. (1995) M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. Cornell, Science 269, 198 (1995).
  * Davis et al. (1995) K. B. Davis, M.-O. Mewes, M. R. Andrews, N. Van Druten, D. Durfee, D. Kurn, and W. Ketterle, Physical review letters 75, 3969 (1995).
  * Zwierlein et al. (2003) M. W. Zwierlein, C. A. Stan, C. H. Schunck, S. M. Raupach, S. Gupta, Z. Hadzibabic, and W. Ketterle, Physical review letters 91, 250401 (2003).
  * Gallagher and Pritchard (1989) A. Gallagher and D. E. Pritchard, Physical review letters 63, 957 (1989).
  * Walker et al. (1990) T. Walker, D. Sesko, and C. Wieman, Physical Review Letters 64, 408 (1990).
  * Ketterle et al. (1993) W. Ketterle, K. B. Davis, M. A. Joffe, A. Martin, and D. E. Pritchard, Physical review letters 70, 2253 (1993).
  * Townsend et al. (1995) C. Townsend, N. Edwards, C. Cooper, K. Zetie, C. Foot, A. Steane, P. Szriftgiser, H. Perrin, and J. Dalibard, Physical Review A 52, 1423 (1995).
  * Pruvost et al. (2000) L. Pruvost, I. Serre, H. T. Duong, and J. Jortner, Physical Review A 61, 053408 (2000).
  * Bachelard et al. (2016) R. Bachelard, N. Piovella, W. Guerin, and R. Kaiser, Physical Review A 94, 033836 (2016).
  * Guerin et al. (2017) W. Guerin, M. Rouabah, and R. Kaiser, Journal of Modern Optics 64, 895 (2017).
  * Camara et al. (2014) A. Camara, R. Kaiser, and G. Labeyrie, Physical Review A 90, 063404 (2014).
  * ferlaino12 et al. (2014) K. Aikawa, A. Frisch, M. Mark, S. Baier, A. Rietzler, R. Grimm, and F. Ferlaino, Physical Review Lett. 108, 210401 (2012).
  * chomaz et al. (2017) L. Chomaz, R. van Bijnen, D. Petter, G. Faraoni, S. Baier, J. Becher, M. Manfred, F. Waechtler, L. Santos, and F. Ferlaino, arXiv preprint arXiv:1705.06914, (2017).
  * Rodrigues et al. (2016) J. Rodrigues, J. Rodrigues, A. Ferreira, and J. Mendonça, Optical and Quantum Electronics 48, 1 (2016).
  * Pichler et al. (1983) G. Pichler, S. MiloSeviC, D. Viza, and R. Beuc, Journal of Physics B: Atomic and Molecular Physics 16, 4619 (1983).
  * Stwalley et al. (1978) W. C. Stwalley, Y.-H. Uang, and G. Pichler, Physical Review Letters 41, 1164 (1978).
  * Chandrasekhar (2013) S. Chandrasekhar, _Radiative transfer_ (Courier Corporation, 2013).
  * Dudko and Weiss (2005) O. K. Dudko and G. H. Weiss, Diffusion Fundamentals 2, 114 (2005).
  * Einstein (1905) A. Einstein, Annalen der physik 322, 549 (1905).
  * Sutherland (1905) W. Sutherland, Philosophical Magazine Series 6 9, 781 (1905), eprint http://dx.doi.org/10.1080/14786440509463331.
  * Von Smoluchowski (1906) M. Von Smoluchowski, Annalen der physik 326, 756 (1906).
  * Kolmogoroff (1931) A. Kolmogoroff, Mathematische Annalen 104, 415 (1931).
  * Bloch et al. (2001) I. Bloch, M. Greiner, O. Mandel, T. W. Hänsch, and T. Esslinger, Physical Review A 64, 021402 (2001).
  * Smoluchowski (1952) R. Smoluchowski, Physical Review 87, 482 (1952).
  * Furutsu and Yamada (1994) K. Furutsu and Y. Yamada, Physical Review E 50, 3634 (1994).
  * Durduran et al. (1997) T. Durduran, A. Yodh, B. Chance, and D. Boas, JOSA A 14, 3358 (1997).
  * (49) L. Diósi, Physica A 199, 517-526 (1993).
  * Graaff and Ten Bosch (2000) R. Graaff and J. Ten Bosch, Optics letters 25, 43 (2000).
  * Lagendijk and Van Tiggelen (1996) A. Lagendijk and B. A. Van Tiggelen, Physics Reports 270, 143 (1996).
  * van Rossum and Nieuwenhuizen (1999) M. v. van Rossum and T. M. Nieuwenhuizen, Reviews of Modern Physics 71, 313 (1999).
  * Valeau et al. (2006) V. Valeau, J. Picaut, and M. Hodgson, The Journal of the Acoustical Society of America 119, 1504 (2006).
  * Tal-Ezer and Kosloff (1984) H. Tal-Ezer and R. Kosloff, The Journal of chemical physics 81, 3967 (1984).
  * Koch et al. (2009) C. P. Koch, M. Ndong, and R. Kosloff, Faraday discussions 142, 389 (2009).
  * Tannoudji et al. (1973) C. C. Tannoudji, D. Bernard, and L. Franck (1973).
  * Weisskopf and Wigner (1930) V. F. Weisskopf and E. P. Wigner, Z. Phys. 63, 54 (1930).
  * Siegman (1986) A. E. Siegman, Mill Valley, CA 37, 208 (1986).
  * Rotondaro and Perram (1997) M. D. Rotondaro and G. P. Perram, Journal of Quantitative Spectroscopy and Radiative Transfer 57, 497 (1997).
  * Peach (1981) G. Peach, Advances in Physics 30, 367 (1981).
  * Kosloff (1994) R. Kosloff, Annual review of physical chemistry 45, 145 (1994).
  * Johnson (2011) S. G. Johnson, _Notes on fft-based differentiation_ (2011).
  * Cohen-Tannoudji et al. (1977) C. Cohen-Tannoudji, B. Diu, and F. Laloë, New York (1977).
  * Xiao et al. (1995) M. Xiao, Y.-q. Li, S.-z. Jin, and J. Gea-Banacloche, Physical Review Letters 74, 666 (1995).
  * Field et al. (1991) J. Field, K. Hahn, and S. Harris, Physical review letters 67, 3062 (1991).
  * Kozyryev et al. (2016) I. Kozyryev, L. Baum, K. Matsuda, B. Hemmerling, and J. M. Doyle, Journal of Physics B: Atomic, Molecular and Optical Physics 49, 134002 (2016).
  * Dieckmann et al. (1998) K. Dieckmann, R. Spreeuw, M. Weidemüller, and J. Walraven, Physical Review A 58, 3891 (1998).
  * Merli et al. (2009) A. Merli, F. Eimer, F. Weise, A. Lindinger, W. Salzmann, T. Mullins, S. Götz, R. Wester, M. Weidemüller, R. Ağanoğlu, et al., Physical Review A 80, 063417 (2009).
  * Vdović et al. (2006) S. Vdović, D. Sarkisyan, and G. Pichler, Optics communications 268, 58 (2006).
  * Stuhl et al. (2012) B. K. Stuhl, M. T. Hummon, M. Yeo, G. Quéméner, J. L. Bohn, and J. Ye, Nature 492, 396 (2012).
  * Bochinski et al. (2003) J. R. Bochinski, E. R. Hudson, H. J. Lewandowski, G. Meijer, and J. Ye, Physical review letters 91, 243001 (2003).
  * Narevicius et al. (2009) E. Narevicius, S. T. Bannerman, and M. G. Raizen, New Journal of Physics 11, 055046 (2009).
  * Zeppenfeld et al. (2012) M. Zeppenfeld, B. G. Englert, R. Glöckner, A. Prehn, M. Mielenz, C. Sommer, L. D. van Buuren, M. Motsch, and G. Rempe, Nature 491, 570 (2012).
  * Bethlem et al. (2000) H. L. Bethlem, G. Berden, F. M. Crompvoets, R. T. Jongma, A. J. Van Roij, and G. Meijer, Nature 406, 491 (2000).
  * Bethlem and Meijer (2003) H. L. Bethlem and G. Meijer, International reviews in physical chemistry 22, 73 (2003).
  * Jochim et al. (2003) S. Jochim, M. Bartenstein, A. Altmeyer, G. Hendl, S. Riedl, C. Chin, J. H. Denschlag, and R. Grimm, Science 302, 2101 (2003).
  * Hemmerling et al. (2016) B. Hemmerling, E. Chae, A. Ravi, L. Anderegg, G. K. Drayna, N. R. Hutzler, A. L. Collopy, J. Ye, W. Ketterle, and J. M. Doyle, Journal of Physics B: Atomic, Molecular and Optical Physics 49, 174001 (2016).
  * Sawyer et al. (2007) B. C. Sawyer, B. L. Lev, E. R. Hudson, B. K. Stuhl, M. Lara, J. L. Bohn, and J. Ye, Physical review letters 98, 253002 (2007).
  * Hummon et al. (2013) M. T. Hummon, M. Yeo, B. K. Stuhl, A. L. Collopy, Y. Xia, and J. Ye, Physical review letters 110, 143001 (2013).
  * Yeo et al. (2015) M. Yeo, M. T. Hummon, A. L. Collopy, B. Yan, B. Hemmerling, E. Chae, J. M. Doyle, and J. Ye, Physical review letters 114, 223003 (2015).
  * Collopy et al. (2015) A. L. Collopy, M. T. Hummon, M. Yeo, B. Yan, and J. Ye, New Journal of Physics 17, 055008 (2015).
  * Ashkenazi et al. (1997) G. Ashkenazi, U. Banin, A. Bartana, R. Kosloff, and S. Ruhman, Advances in Chemical Physics 100, 229 (1997).
  * Magnus (1954) W. Magnus, Communications on pure and applied mathematics 7, 649 (1954).


# A Unique Perfect Power Decagonal Number

Philippe Michaud-Rodgers

###### Abstract

Let $\mathcal{P}_{s}(n)$ denote the $n$th $s$-gonal number. We consider the
equation

$\mathcal{P}_{s}(n)=y^{m},$

for integers $n,s,y,$ and $m$. All solutions to this equation are known for
$m>2$ and $s\in\\{3,5,6,8,20\\}$. We consider the case $s=10$, that of
decagonal numbers. Using a descent argument and the modular method, we prove
that the only decagonal number $>1$ expressible as a perfect $m$th power with
$m>1$ is $\mathcal{P}_{10}(3)=3^{3}$.

##  1 Introduction

The $n$th $s$-gonal number, with $s\geq 3$, which we denote by
$\mathcal{P}_{s}(n)$, is given by the formula

$\mathcal{P}_{s}(n)=\frac{(s-2)n^{2}-(s-4)n}{2}.$

Polygonal numbers have been studied since antiquity [6, pp.1–39] and relations
between different polygonal numbers and perfect powers have received much
attention in the literature (see [7], for example, and references therein).
Kim, Park, and Pintér [7, Theorem 1.2] find all solutions to the equation
$\mathcal{P}_{s}(n)=y^{m}$ when $m>2$ and $s\in\\{3,5,6,8,20\\}$ for integers
$n$ and $y$. We extend this result (for $m>1$) to the case $s=10$, that of
decagonal numbers.

######  Theorem 1.

All solutions to the equation

$\mathcal{P}_{10}(n)=y^{m},\quad n,y,m\in\mathbb{Z},\leavevmode\nobreak\ m>1,$
(1)

satisfy $n=y=0,n=\mathinner{\\!\left\lvert y\right\rvert}=1,$ or $n=y=m=3$.

In particular, the only decagonal number $>1$ expressible as a perfect $m$th
power with $m>1$ is $\mathcal{P}_{10}(3)=3^{3}$.

We will prove Theorem 1 by carrying out a descent argument to obtain various
ternary Diophantine equations, to which one may associate Frey elliptic
curves. The difficulty in solving the equation $\mathcal{P}_{s}(n)=y^{m}$ for
a fixed value of $s$ is due to the existence of the trivial solution $n=y=1$
(for any value of $m$). We note that adapting our method of proof also works
for the cases $s\in\\{3,5,6,8,20\\}$ mentioned above, but will not extend to
any other values of $s$ (see Remark 3.1).

##  2 Descent and Small Values of $m$

We note that it will be enough to prove Theorem 1 in the case $m=p$, prime. We
write (1) as

$n(4n-3)=y^{p},\quad n,y\in\mathbb{Z},\leavevmode\nobreak\ p\text{ prime},$
(2)

and suppose that $n,y\in\mathbb{Z}$ satisfy this equation with $n\neq 0$.

### Case 1: $3\nmid n$.

If $3\nmid n$, then $n$ and $4n-3$ are coprime, so there exist coprime
integers $a$ and $b$ such that

$n=a^{p}\quad\text{ and }\quad 4n-3=b^{p}.$

It follows that

$4a^{p}-b^{p}=3.$ (3)

If $p=2$ we see that $(2a-b)(2a+b)=3$, so that $a=b=\pm 1$ and so
$n=\mathinner{\\!\left\lvert y\right\rvert}=1$. If $p=3$ or $p=5$, then using
Magma’s [5] Thue equation solver, we find that $a=b=1$ also.

### Case 2: $3\parallel n$.

Suppose $3\parallel n$ (i.e. $\mathrm{ord}_{3}(n)=1$). Then after dividing (2)
by $3^{\mathrm{ord}_{3}(y)p}$, we see that there exist coprime integers $t$
and $u$, with $3\nmid t$ such that

$n=3t^{p}\quad\text{ and }\quad 4n-3=3^{p-1}u^{p}.$

Then

$4t^{p}-3^{p-2}u^{p}=1.$ (4)

If $p=2$ we have $(2t-u)(2t+u)=1$ which has no solutions. If $p=3$, then
$4t^{3}-3u^{3}=1$, and using Magma’s Thue equation solver, we verify that
$u=t=1$ is the only solution to this equation, and this gives $n=y=3$. If
$p=5$, we use Magma’s Thue equation solver to see that there are no solutions.

### Case 3: $3^{2}\mid n$.

If $3^{2}\mid n$ then $3\parallel 4n-3$, and arguing as in Case 2, there exist
coprime integers $v$ and $w$, with $3\nmid w$ such that

$n=3^{p-1}v^{p}\quad\text{ and }\quad 4n-3=3w^{p}.$

So

$4\cdot 3^{p-2}v^{p}-w^{p}=1.$ (5)

If $p=2$, then as in Case 2 we obtain no solutions. If $p=3$ or $p=5$ then we
use Magma’s Thue equation solver to verify that there are no solutions with
$v\neq 0$.

##  3 Frey Curves and the Modular Method

To prove Theorem 1, we will associate Frey curves to equations (3), (4), and
(5), and level-lower to obtain a contradiction. We have considered the cases
$p=2,3,$ and $5$ in Section 2, and so we will assume $m=p$ is prime, with
$p\geq 7$.

We note that at this point we could directly apply [3, Theorem 1.2] to
conclude that the only solutions to (6) are $a=b=1$, giving $n=1$, and apply
[2, Theorem 1.2] that (7) and (8) have no solutions. The computations for (6)
are not explicitly carried out in [3], so for the convenience of the reader
and to highlight why the case $s=10$ is somewhat special, we provide some
details of the arguments.

### Case 1: $3\nmid n$.

We write (3) as

$-b^{p}+4a^{p}=3\cdot 1^{2},$ (6)

which we view as a generalised Fermat equation of signature $(p,p,2)$. We note
that the three terms are integral and coprime.

We suppose $ab\neq\pm 1$. Following the recipes of [3, pp. 26–31], we
associate Frey curves to (6). We first note that $b$ is odd, since
$b^{p}=4n-3$. If $a\equiv 1\pmod{4}$, we set

$E_{1}\mathrel{\mathop{\mathchar 58\relax}}Y^{2}=X^{3}-3X^{2}+3a^{p}X.$

If $a\equiv 3\pmod{4}$, we set

$E_{2}\mathrel{\mathop{\mathchar 58\relax}}Y^{2}=X^{3}+3X^{2}+3a^{p}X.$

If $a$ is even, we set

$E_{3}\mathrel{\mathop{\mathchar
58\relax}}Y^{2}+XY=X^{3}-X^{2}+\frac{3a^{p}}{16}X.$

We level-lower each Frey curve and find that for $i=1,2,3,$ we have
$E_{i}\sim_{p}f_{i}$, for $f_{i}$ a newform at level $N_{p_{i}}$, where
$N_{p_{1}}=36,N_{p_{2}}=72,$ and $N_{p_{3}}=18$. The notation $E\sim_{p}f$
means that the mod-$p$ Galois representation of $E$ arises from $f$. There are
no newforms at level $18$ and so we focus on the curves $E_{1}$ and $E_{2}$.
There is a unique newform, $f_{1}$, at level $36$, and a unique newform,
$f_{2}$, at level $72$.

The newform $f_{1}$ has complex multiplication by the imaginary quadratic
field $\mathbb{Q}(\sqrt{-3})$. This allows us to apply [3, Proposition 4.6].
Since $2\nmid ab$ and $3\nmid ab$, we conclude that $p=7$ or $13$, and that
all elliptic curves of conductor $2p$ have positive rank over
$\mathbb{Q}(\sqrt{-3})$. However, it is straightforward to check that this is
not the case for $p=7$ and $13$. We conclude that $E_{1}\not\sim_{p}f_{1}$.

Let $F_{2}$ denote the elliptic curve with Cremona label 72a2 whose isogeny
class corresponds to $f_{2}$. This elliptic curve has full two-torsion over
the rationals and has $j$-invariant $2^{4}\cdot 3^{-2}\cdot 13^{3}$. We apply
[3, Proposition 4.4], which uses an image of inertia argument, to obtain a
contradiction in this case too.

######  Remark 3.1.

The trivial solution $a=b=1$ (or $n=y=1$) corresponds to the case $i=1$ above.
The only reason we are able to discard the isomorphism $E_{1}\sim_{p}f_{1}$ is
because the newform $f_{1}$ has complex multiplication. The modular method
would fail to eliminate the newform $f_{1}$ otherwise. For each value of $s$,
we can associate to (1) generalised Fermat equations of signature $(p,p,2)$,
$(p,p,3)$, and $(p,p,p)$. We found we could only obtain newforms with complex
multiplication (when considering the case corresponding to the trivial
solution) when $s=3,6,8,10,$ or $20$. A similar strategy of proof also works
for $s=5$ using the work of Bennett [1, p. 3] on equations of the form
$(a+1)x^{n}-ay^{n}=1$ to deal with the trivial solution.

### Case 2: $3\parallel n$.

We rewrite (4) as

$4t^{p}-3^{p-2}u^{p}=1\cdot 1^{3},$ (7)

which we view as a generalised Fermat equation of signature $(p,p,3)$. The
three terms are integral and coprime. We suppose $tu\neq\pm 1$. Using the
recipes of [4, pp. 1401–1406], we associate to (7) the Frey curve

$E_{4}\mathrel{\mathop{\mathchar 58\relax}}Y^{2}+3XY-3^{p-2}u^{p}\,Y=X^{3}.$

We level-lower $E_{4}$ and find that $E_{4}\sim_{p}f$, where $f$ is a newform
at level $6$, an immediate contradiction, as there are no newforms at level
$6$.

### Case 3: $3^{2}\mid n$.

We rewrite (5) as

$-w^{p}+4\cdot 3^{p-2}v^{p}=1\cdot 1^{3},$ (8)

which we view as a generalised Fermat equation of signature $(p,p,3)$. The
three terms are integral and coprime. We suppose $vw\neq\pm 1$. The Frey curve
we attach to (8) is

$E_{5}\mathrel{\mathop{\mathchar 58\relax}}Y^{2}+3XY+4\cdot
3^{p-2}v^{p}\,Y=X^{3}.$

We level-lower, and find that $E_{5}\sim_{p}f$, where $f$ is a newform at
level $6$, a contradiction as in Case 2.

This completes the proof of Theorem 1.

## References

  * [1] M. Bennett.  Rational approximation to algebraic numbers of small height: the Diophantine equation $\,\mathinner{\\!\left\lvert ax^{n}+by^{n}\right\rvert}=1$.  _J. Reine Angew. Math._ , 535:1–49, 2001.
  * [2] M. Bennett.  Products of consecutive integers.  _Bull. Lond. Math. Soc._ , 36(5):683–694, 2004.
  * [3] M. Bennett and C. Skinner.  Ternary Diophantine equations via Galois representations and modular forms.  _Canad. J. Math._ , 56(1):23–54, 2004.
  * [4] M. Bennett, V. Vatsal, and S. Yazdani.  Ternary Diophantine equations of signature $(p,p,3)$.  _Compos. Math._ , 140(6):1399–1416, 2004.
  * Bosma et al. [1997] W. Bosma, J. Cannon, and C. Playoust.  The Magma algebra system. I. The user language.  _J. Symbolic Comput._ , 24(3-4):235–265, 1997.
  * [6] L. Dickson  _History of the Theory of Numbers. Volume II: Diophantine Analysis_.  Dover Publications, New York, 2005.
  * [7] D. Kim, K. Park, and A. Pintér.  A Diophantine problem concerning polygonal numbers.  _Bull. Aust. Math. Soc._ , 88(2):345–350, 2013.

Mathematics Institute, University of Warwick, CV4 7AL, United Kingdom

E-mail address<EMAIL_ADDRESS>