# Cesium atoms in cryogenic argon matrix

T. Battard(1) S. Lahs(1) C. Crépin(2) D. Comparat(1)
<EMAIL_ADDRESS>(1) Université Paris-Saclay, CNRS,
Laboratoire Aimé Cotton, 91405, Orsay, France.  (2) Université Paris-Saclay,
CNRS, Institut des Sciences Moléculaires d’Orsay, 91405, Orsay, France

###### Abstract

This paper presents both experimental and theoretical investigations into the
spectroscopy of dilute cesium (Cs) atoms within a solid argon (Ar) matrix at
cryogenic temperatures. This system is relevant for matrix isolation
spectroscopy and in particular for recently proposed methods for investigating
phenomena that extend beyond the standard model of particle physics. We record
absorption spectra at various deposition temperatures and examine the
evolution of these spectra post-deposition with respect to temperature
changes.

Taking advantage of Cs-Ar and Ar-Ar pairwise interaction potentials, we
conduct a stability study of trapping sites, which indicates a preference for
Td (tetrahedral, 4 vacancies) and Oh (cubic, 6 vacancies) symmetries. By
implementing a mean-field analysis of the long-range Cs(6s,6p)-Ar-Ar triple
dipole interaction, combined with a temperature-dependent shift in zero point
energy, we propose effective Cs(6s,6p)-Ar pairwise potentials. Upon
integrating these pairwise potentials with spin-orbit coupling, we achieve a
satisfactory agreement between the observed and simulated absorption line
positions. The observed line broadening is reasonably well reproduced by a
semi-classical thermal Monte Carlo approach based on Mulliken-type differences
between excited and ground potential curves. Additionally, we develop a
simple, first-order crystal field theory featuring only 6 interaction mode
coordinates. It uses the reflection approximation and incorporates quantized
(phonon) normal modes. This produces a narrow triplet structure but not the
observed amount of splitting.

##  I Introduction

Despite matrix isolation spectroscopy having been studied since the 1950s
barnes81matrix ; almond89 ; bondybey1996new , our understanding of the shape
and behavior of trapping environments within the matrix remains limited. This
is particularly true for even the simplest systems comprising single valence
electron atoms (alkali atoms) in rare gas matrix environments, where no
consensus has yet been reached regarding the trapping site with the host atom,
especially for heavier atoms crepin1999photophysics ; davis2018investigation ;
ozerov2021generic .

Gaining a more accurate understanding of the trapping site and the matrix’s
effect on the dopant would prove invaluable for precision spectroscopy
experiments, such as magnetometry, spin-dependent interactions, and
investigations into physics beyond the Standard Model. These investigations
may include searches for violations of fundamental symmetries (parity, time-
reversal) like electric dipole moment (EDM) particles or axion-like dark
matter candidates. Precision spectroscopy experiments utilize spin-induced
transitions, which are not directly affected by the matrix’s pure
electrostatic interaction at first order. However, a detailed study is
necessary to assess the precision that these experiments can attain
pryor1987artificial ; arndt1993can ; kozlov2006proposal ; atoms1018Hessel ;
upadhyay2019spin ; braggio2022spectroscopy ; budker2022quantum .

In this article, our focus is on studying cesium atoms trapped in an argon
matrix. Understanding this simple system will help to further characterize
more complex experiments. Among the (stable) alkali atoms, Cs is the heaviest,
making it highly sensitive to the effects of the electron EDM
safronova2018search . We have chosen Argon because it is the only rare gas
whose naturally occurring isotopes do not possess a nuclear spin that would
interact with the cesium spin.

Absorption (or transmission) spectra form the foundation for all optical
manipulations necessary for the above-proposed experiments. Intriguingly, only
a few experiments have studied Cs embedded in Argon ien1962electron ;
goldsborough1964electron , and only two optical transmission spectra have been
published weyhmann1965optical ; Kanda1971  (balling1983laser  gives only
absorption and emission frequencies). For Cs in Ar, two structures were
observed. They could arise from $6s\rightarrow 6p$ transition in two different
trapping site environments: the $6p$ triplet degeneracy being lifted either by
a low symmetry site weyhmann1965optical  or by a dynamical Jahn-Teller effect
in a cubic symmetry site Kanda1971 . Ref. balling1983laser  observed, in
addition, $6s\rightarrow 5d$ and $6s\rightarrow 7p$ transitions, along with
relaxation and fluorescence from $6p$ and $5d$ levels.

In this article, we conduct new experiments with the aim of reconciling the
discrepancies between previous measurements. We begin by presenting the
experimental setup, which enables precise Cs density measurements in highly
polycrystalline Ar samples. Then, using Cs-Ar and Ar-Ar pairwise interaction
potentials, we first conduct a trapping site stability study. By correcting
the pairwise interaction with effective third-order effects, we predict
absorption line positions for the different trapping sites found. We then
validate, using a simple crystal field model, whether the found linewidth
agrees with the proposed trapping sites. Finally, we provide conclusions and
potential improvements for further exploration.

##  II Experimental setup

Figure 1 illustrates the experimental setup. The experiment is performed under
a residual gas pressure of $\sim 10^{-7}\,$mbar in a two-stage cryostat
manufactured by Mycryofirm using a pulsed tube cryocooler (SHI RP-082E2S). The
cryogenic matrix is grown by depositing Ar gas onto a C-cut sapphire (Al2O3)
plate, which is $20\text{\,}\mathrm{m}\mathrm{m}$ in diameter and
$1\text{\,}\mathrm{m}\mathrm{m}$ in thickness, and embedded in a copper frame.
The cryogenic temperature of the sample holder, ranging from
$3\text{\,}\mathrm{K}$ to $30\text{\,}\mathrm{K}$, is regulated via resistive
heating. In order to minimize thermal disturbances caused by the kinetic
energy transfer of the Ar beam, the sample holder’s temperature is maintained
at a constant level via a PID loop.

###  II.1 Sample growth

Figure 1: (Color available online) Schematic representation of the cryostat
and the optical systems utilized to monitor the matrix growth rate and to
measure the absorption and fluorescence spectra. The arrows indicate the
propagation directions of the beams. The laser illumination of the matrix is
imaged on CCD camera 2 to display the interference pattern.

The thickness $L$ of the argon layer is continuously monitored via 2D
reflection Fizeau interferometry. A pattern formed by the interference between
the light reflected from the front and back sides of the matrix is captured by
a lens onto a CCD camera. The 2D imaging offers numerous advantages compared
to previous 1D transmission or reflection laser interferometry studies
hollenberg1962absolute ; groner1973measurement ; jiang1975absolute : it
features a background-free environment ensuring excellent contrast, it is less
restricted by beam size and spatial interference narrowing, and it enables the
monitoring of the whole sample thereby ensuring the quality of deposition.

For this Ar growth measurement, a DFB diode laser operating at
$\lambda=$852\text{\,}\mathrm{n}\mathrm{m}$$ was employed, providing the
advantage of also allowing the monitoring of the Cs flux through saturated
absorption.

The interference fringes at a specific point evolve as per
$\frac{2n_{\text{Ar}}}{\lambda}\frac{dL}{dt}$, where $n_{\text{Ar}}=1.3$ is
the refractive index of an argon bulk crystal at 852 nm sinnock1969refractive
.

A growth of approximately 3 fringes per minute is observed, corresponding to a
growth rate for an ideal argon bulk crystal of
$\frac{dL}{dt}=1\mu\mathrm{m/min}$. After about one hour of argon deposition,
a crystal is formed with a thickness of
$\sim$$50\text{\,}\mathrm{\SIUnitSymbolMicro m}$.

###  II.2 Cs density

To embed cesium into the matrix, we utilize Alfasource 3S dispensers from
Alfavakuo e.U., loaded with a CsBi25 alloy. The alkali purity of the cesium is
99.980%. The dispenser is situated outside the $50\text{\,}\mathrm{K}$ and
$4\text{\,}\mathrm{K}$ thermal shields and is activated by resistive current
heating. The Cs vapor is funneled to the substrate through a steel tube with
an inner diameter $d$ of $4\text{\,}\mathrm{m}\mathrm{m}$, which pierces both
thermal shields. Independently, the tube is heated to approximately
$40\text{\,}\mathrm{\SIUnitSymbolCelsius}$ via resistors to minimize blockage
from condensing Cs. We have extensively studied in Ref. hahn2022comparative
the behavior of cesium atom effusion through a collimating tube and
established that, irrespective of the tube material, Cs adheres to the wall.
This results in a well-collimated and homogeneous Cs atomic beam.

We monitor the Cs flux by recording a laser absorption spectrum while
adjusting the frequency of the low-intensity diode laser intersecting the
cesium beam. We assume that the Cs beam diameter aligns with the inner tube
diameter $d$ and its density $n$ is uniform. Confirmation of this is obtained
by scanning the laser beam through the atomic beam or by a simple 2D
absorption image of the Cs in the Ar matrix, where the cesium-doped zone
creates a visible spot that is about $5\text{\,}\mathrm{m}\mathrm{m}$ wide.
The velocity distribution is effusive hahn2022comparative :
$f_{z}(v_{z})=\left(\frac{m}{\sqrt{2}k_{B}T_{z}}\right)^{2}v_{z}^{3}e^{-\frac{mv_{z}^{2}}{2k_{B}T_{z}}}$
along $z$. The Cs beam is intersected by the laser at an angle of
$\alpha=$22\text{\,}\mathrm{\SIUnitSymbolDegree}$$ relative to the beam axis,
as shown in Fig. 1. Hence, the Beer-Lambert-Bouguer absorption is given by
$nd\sigma_{0}\int_{0}^{\infty}\frac{1}{1+\left(\frac{\omega_{L}-\omega_{\rm
Cs}-kv_{z}\sin\alpha}{\Gamma/2}\right)^{2}}f_{z}(v_{z})dv_{z}$ where
$\sigma_{0}=$1.4\text{\times}{10}^{-9}\text{\,}\mathrm{c}\mathrm{m}^{2}$$ is
the absorption cross section for an isotropic light polarization at the
resonant $6s_{1/2}(F=4)\rightarrow 6p_{3/2}(F=5)$ frequency $\hbar\omega_{\rm
Cs}$, which has a natural spontaneous decay rate of
$\Gamma=1/($30.4\text{\,}\mathrm{n}\mathrm{s}$)$ steck2019cesium . We fit the
experimental absorption profile to this formula while scanning the laser
angular frequency $\omega_{L}$. The results align with and are not highly
sensitive to $T_{z}\approx$400\text{\,}\mathrm{K}$$. We then determine the
beam density $n$ and consequently the Cs flux $n\pi(d/2)^{2}\bar{v}$, where
$\bar{v}=\sqrt{8k_{B}T_{z}/\pi m}$ is the average velocity. We only account
for the population in the hyperfine level $F=4$ (which represents only 9/16 of
the total Cs(6s) population in the beam, assuming Boltzmann-equidistribution
between the 9 times degenerate $F=4$ and 7 times $F=3$ levels). Ultimately, we
obtain a Cs(6s) flux of a few $10^{12}\ $at/s reaching the Ar matrix. We thus
anticipate a $\sim 10^{-4}$ atomic ratio of Cs in the argon matrix (density on
the order of $10^{18}$ cm-3) and a typical average internuclear distance
between Cs atoms on the order of $10\text{\,}\mathrm{n}\mathrm{m}$.

To prevent Cs contamination of the window and to generate well-defined
trapping environments for all the Cs atoms, we always initiate Cs deposition
only after a few micrometers of Ar have already been deposited.

###  II.3 Quality of the samples

Undoubtedly, the epitaxy of an fcc Ar crystal is unfeasible on the hexagonal
lattice of sapphire due to the incompatibility of their structural lattices
nepijko2005morphology , leading to the formation of several crystal defects,
including holes or cavities, as observed in Ref nepijko2005morphology .

In theory, Cs can be trapped in various defect types such as point defects
surrounded by vacancies or line-dislocation defects (including screw, edge,
etc.), surface defects (like stacking fault, tilt and twist grain boundaries),
or volume defects (e.g., pores, cracks, or alternate phases like hcp
inclusions) patterson2018solid .

It is known from literature that the grain size is inversely proportional to
the growth rate smith1970inert , and the deposition temperature significantly
influences the crystal quality. The crystals contain twins, stacking faults,
and dislocations venables1977rare . At $20\text{\,}\mathrm{K}$ and
$10\text{\,}\mathrm{K}$, the inter-dislocation distances are roughly
$8\text{\,}\mathrm{n}\mathrm{m}$ and $3\text{\,}\mathrm{n}\mathrm{m}$,
respectively kovalenko1970structure ; shakhsemampour1984magnetic . This
suggests the likely presence of small crystallites comprised of thousands of
atoms, separated by dislocations. Additionally, if a film is grown at a
condensation temperature greater than $2/3$ of the sublimation temperature
(that is $30\text{\,}\mathrm{K}$ for Ar), the fcc grains are about
$100\text{\,}\mathrm{n}\mathrm{m}$ in size. Conversely, at 1/3 of the
sublimation temperature (that is $10\text{\,}\mathrm{K}$ for Ar), the grain
size reduces to around $10\text{\,}\mathrm{n}\mathrm{m}$, with indications of
a minority hcp phase smith1970inert ; rudman1978rare ; song2013self .
Furthermore, below a critical temperature (18 K for Ar), atomic-scale cavities
appear in the lattice schulze1973density .

These studies align well with our experimental observations regarding the
quality of the Fizeau interference image during crystal growth, and the
reduction in light transmission intensity due to scattering on grain
boundaries. Specifically, as described in pollack1964solid , a sample turns
almost opaque when the Ar crystal grain diameter reaches around $0.1\lambda$.
We attribute this to exactly this increased scattering.

Although we attempted annealing up to $30\text{\,}\mathrm{K}$, the Ar crystal
quality showed no significant improvement. This outcome aligns with
expectations, as higher temperatures (closer to the triple point) would be
necessary tovbin2015temperature ; tovbin2017issues  to promote vacancy
formation or atom-vacancy exchange through diffusion tishchenko1982self ;
tishchenko1984kinetics ; leibin2021modeling . However, achieving such
temperatures in a free-standing sample in vacuum is not possible as
sublimation would occur prematurely due to the high vapor pressure of solid
argon.

##  III Absorption spectra

###  III.1 Experimental Spectra

In our absorption (or transmission) spectroscopy, we use an Avantes AvaLight-
HAL tungsten-halogen light source and an Ocean Optics QE65000 spectrometer.
The spectrometer was calibrated by measuring the spectra of Hg and Kr lamps.

To measure the absorption in the doped Ar crystal, we shine light from the
lamp through the sample. The position and shape of the light on the plate can
be adjusted with irises and mirrors and monitored via a camera.

We then compute the spectral absorbance of the sample due to Cs as
$A=\log_{10}\frac{I_{0}}{I_{\rm t}}$, where $I_{t}$ denotes the spectral
irradiance density recorded in the presence of the Cs-Ar sample, and $I_{0}$
is obtained prior to any Cs deposition.

The captured spectra from varying deposition temperatures are presented in
Fig. 2 a). Above 14K, the spectra become significantly broader, and we
encountered difficulties producing any spectra for deposition temperatures
above this.

Conversely, if Cs and Ar are deposited at a low temperature (8K in the case of
Fig. 2 b), the spectra remain visible even when the sample is subsequently
heated up to approximately 35K, a temperature at which argon begins to
sublimate under vacuum.

Figure 2: a) Absorbance spectra of Cs embedded in Ar for deposits conducted at
different temperatures. b) Evolution with temperature for a spectrum initially
created at 8K. The spectra are offset in ordinate for better visual clarity.
The positions of the 6s${}_{1/2},\rightarrow$6p1/2 (left) and
6s${}_{1/2},\rightarrow$6p3/2 (right) Cs gas-phase transitions are indicated
for reference.

The spectra show two triplet structures similar to those observed by
weyhmann1965optical ; Kanda1971 , attributed to two types of trapping sites.
For comparison, spectra of previous experiments were obtained at
$1.8\text{\,}\mathrm{K}$ in Ref. Kanda1971 , likely performed at
$4.2\text{\,}\mathrm{K}$ in weyhmann1965optical  and performed at
$10\text{\,}\mathrm{K}$ in balling1983laser . Typical positions and widths of
the two triplets structures at 8K are presented in Table 1. The differing
relative intensities of the two triplets in the various samples (Fig. 2 a)
imply that depending on the deposition conditions (deposit rate, sample
density, temperature, argon purity, etc.), one or the other trapping site is
favored. It is worth noting that we used slightly different Ar and Cs fluxes
for different samples. The flux estimates given in the previous section are
measured for the sample deposited at $8\text{\,}\mathrm{K}$, which will be the
one we focus on in the following. For other deposits the rates used can easily
vary by a factor of 3.

The temperature evolution shown in Fig. 2 b) shows a more pronounced spectral
deformation in the red site (the one between 11000 and 12500cm^-1) than in the
blue site (the one between 12500 and 14000cm^-1). The temperature changes
below 20K are mainly reversible. In the 6-14K temperature range corresponding
to a well-defined red triplet, there is a slight blue shift ($\sim 50$ cm-1)
of the two lower energy components from 6K to 14K, while the position of the
highest energy component remains constant. The bandwidths of this triplet
structure do not depend on the temperature. Conversely, a slight broadening of
the components attributed to the blue site can be inferred from the analysis
in Gaussian bands (less than 50 cm-1 from 6 to 20K), while their positions do
not shift with temperature. Heating beyond 20K typically leads to a drop in
the baseline (likely due to stronger light scattering in the Ar), along with
irreversible changes in the red triplet. The blue site, however, remains well-
defined until 28K. Due to the lack of reproducibility, it is not possible to
draw more precise conclusions about the evolution of temperature. These
informations should be handled with care.

line center (cm-1) | line center (nm) | FWHM (cm-1)
---|---|---
11 220 | 891 | 310
11 780 | 849 | 360
12 170 | 822 | 350
12 770 | 783 | 400
13 400 | 746 | 380
13 850 | 722 | 400
Table 1: Line center and width extracted from Fig. 2 b) at
$8\text{\,}\mathrm{K}$. Both triplets were fitted by a sum of 3 Gaussians.

More broadly, these types of triplet structures are well-documented in alkali
atoms trapped in rare gas matrices. They are not due to resolved phonon lines
but are created by homogeneous broadening and result from the lifting of the
degeneracy of the $p$ level by the crystal environment crepin1999photophysics
. However, due to the finite temperature, it is not straightforward to
differentiate between the splitting created by a low symmetry crystalline
field, as in dislocations or surface defects, and one created by a
(temperature-dependent) dynamical Jahn-Teller effect in a more symmetric
(cubic, tetrahedral, etc.) crystalline field. From the shape of the lines,
however, it already seems evident that the main splitting is on the order of
the fine structure splitting between $6p_{1/2}$ and $6p_{3/2}$ (of nearly
$554\text{\,}\mathrm{c}\mathrm{m}^{-1}$ in the gas phase). Then the static or
dynamical electric crystal field will lift the degeneracy between the
$|m|=3/2$ and $|m|=1/2$ of the $6p_{3/2}$ level.

###  III.2 Trapping sites

####  III.2.1 Potential curves

To explore potential trapping sites within an fcc lattice, we should conduct a
stability analysis.

Utilizing a Lennard-Jones potential with an equivalent depth of
$\sim$45\text{\,}\mathrm{c}\mathrm{m}^{-1}$$ and an equilibrium distance of
$\sim$0.55\text{\,}\mathrm{n}\mathrm{m}$$ as the X${}^{2}\Sigma^{+}_{1/2}$,
Cs-Ar potential blank2012m+ ; kobayashi2016ab ; medvedev2018ab , generic
calculations for the accommodation of an atom in fcc and hcp rare-gas solids
ozerov2019computational ; ozerov2021generic  suggest that the likely trapping
sites for Cs in an Ar fcc matrix could have 6, 8, or potentially 10 vacancies
with respective symmetries of Oh, C2v, C4v. Nonetheless, a 7-vacancy fcc
lattice (C3v symmetry), a 5-vacancy trapping site located in an hcp
environment, or a grain boundary akin to a stacking fault accommodation (D3h
symmetry) could also emerge davis2018investigation ; kleshchina2019stable .

Figure 3: a) Spin-orbit constant for Cs(6p)-Ar interaction. b) Cs(6p)-Ar
potential curves. c) Cs(6s)-Ar potential curves. d) Ar-Ar ground state
potential curves. The potential curves $V$ used in this study are compared to
the zero-point energy $V_{FH}$, and a third-order correction $V^{(3)}$. The
original potentials are taken from kobayashi2016ab . For the excited state, we
chose to plot only a third-order case with no cut-off (in analogy with the
ground state case, see text); i.e., the two-body potential $V$ (offset by
$E(6p)$) multiplied by
$\frac{1}{3}\frac{8\pi\rho}{9}\frac{C_{9}^{*}}{C_{6}^{*}}$.

However, as depicted in Fig. 3 c), the Lennard-Jones potential fails to
accurately reproduce the Cs-Ar ground state interaction. Consequently, we
conducted our study using more precise potentials.

For the $V_{\rm Ar-Ar}$ potentials, we use the simple analytical formula
provided by deiters2019two . For the Cs-Ar ground state $V_{\rm Ar-Cs}=V_{\rm
X\,^{2}\Sigma_{1/2}^{+}}$, we employed the potential curve from Ref.
kobayashi2016ab  with a third-order interpolation between points. Due to a
lack of data points in the long-range, we incorporated parts of the potential
from blank2012m+  for inter-nuclear distances exceeding
$2\text{\,}\mathrm{n}\mathrm{m}$. These potentials are presented in Fig. 3.
Although the Lennard-Jones approximation accurately represents the Ar-Ar
potential, it is considerably inaccurate for the Cs-Ar potential.

For comprehensive Ar-Ar and Cs-Ar interactions, we use a simple two-body
(pairwise potential) approximation. Here, the full interaction energy of a Cs-
doped crystal with $(N-n)$ Ar atoms and $n$ vacancies is given by the
following equation, using obvious notations:

$V=\sum_{i=1}^{N-n}V_{\rm Ar-Cs}(\bm{R}_{Cs-Ar_{i}})+\sum_{1\leq i<j\leq
N-n}V_{\rm Ar-Ar}(\bm{R}_{ij})$ (1)

We began by using a simple cubic grid with $N=4\times 7^{3}=1372$ atoms (the
factor 4 arises from the 4 atoms in the fcc primitive cell). Convergence is
achieved within one percent of the energy for this configuration. This number
of atoms enables us to simulate a small homogeneous chrystallite that would
sit between dislocations.

For a central Ar atom located at $\bm{R}_{0}$, the sum over all other atoms
$\sum_{i}V_{\rm Ar-Ar}(\bm{R}_{0}-\bm{R}_{i})$ reaches
$\approx$773\text{\,}\mathrm{c}\mathrm{m}^{-1}$$ (here, a larger grid is
necessary for convergence). This significantly overestimates the experimental
value for cohesive energy of $645\text{\,}\mathrm{c}\mathrm{m}^{-1}$
lotrich1997three ; schwerdtfeger2016towards . Similarly, the equilibrium
distance is established for a lattice constant of
$0.521\text{\,}\mathrm{n}\mathrm{m}$, compared to the correct value of
$0.531\text{\,}\mathrm{n}\mathrm{m}$ schwerdtfeger2016towards . To reproduce
the correct cohesive (also called atomization $-E_{\rm at}$) energy and the
equilibrium lattice constant, we scaled the Ar-Ar potential by $V=\alpha
V_{\rm Ar-Ar}(\beta R)$ with $\alpha=645/773$ and $\beta=0.521/0.531$. Other
variants of empirically modified potentials, such as $\alpha V_{\rm Ar-
Ar}(\beta+R)$ used in Ref. bezrukov2019empirically , yielded similar outcomes.
This method is a straightforward way to retain the simple two-body pairwise
model while ensuring the appropriate cohesive energy and equilibrium lattice
constant. This is achieved by incorporating many-body corrections and zero-
point energy effectively schwerdtfeger2016towards .

####  III.2.2 Site Stability

To evaluate the thermodynamic stability of a structure with $n$ vacancies, we
calculate a variant of the Gibbs (or enthalpy) free accommodation energy for
the structure lara2019effects ; kleshchina2017modeling ;
ozerov2019computational : $\Delta E_{N}(n)=E_{\rm Cs}(n,N)-E_{\rm
Ar}(N)(N-n)/N$. Here, $E_{\rm Cs}(n,N)$ signifies the total energy of the
system when a Cs atom is introduced into the crystal. Conversely, $E_{\rm
Ar}(N)$ represents the total energy of a pure Ar crystal comprising N atoms,
excluding the vacancies.

For the system to exhibit stability, it is crucial that the dissociation of
$M$ systems, each containing $n$ vacancies, into $M^{\prime}$ subsystems with
$n^{\prime}$ vacancies and $M^{\prime\prime}$ subsystems with
$n^{\prime\prime}$ vacancies be energetically unfavorable. This requirement
can be represented as: $M\Delta E_{N}(n)<M^{\prime}\Delta
E_{N}(n^{\prime})+M^{\prime\prime}\Delta E_{N}(n^{\prime\prime})$, with the
condition $n^{\prime\prime}M^{\prime\prime}+n^{\prime}M^{\prime}=nM$ ensuring
the conservation of atom numbers. For instance, two systems each with 5
vacancies can evolve into a more stable configuration of one system with 4
vacancies and another with 6 vacancies. Thus, for $n$ vacancies to form a
stable trapping site, $\Delta E_{N}(n)$ must fall within the convex hull of
the free energy curve as it varies with the number of vacancies tao2015heat ;
kleshchina2017modeling .

Despite its significance, $\Delta E_{N}(n)$ is just one of several vital
factors, such as migration and activation energies, needed for understanding
the trapping sites. The $\Delta E_{N}(n)$ formula is primarily beneficial for
small crystals, or in our case, a small single crystallite in a poly-crystal.
Alternatively, for $N\rightarrow\infty$, the thermodynamic limit is often
applied with $E_{\rm Ar}(N)/N\approx E_{\rm at}$, resulting in $\Delta
E^{\infty}_{N}(n)=E_{\rm Cs}(n,N)-E_{\rm Ar}(N)+nE_{\rm at}$ lam1985multiple ;
lam1986calculations ; sabochick1988atomistic ; zenia2016stability ;
lounis2016stability .

Figure 4: For a given number of vacancies $n$, the optimized configuration
found to lower $\Delta E(n)$ is displayed (calculation performed with $N=1372$
atoms). The figure illustrates both $\Delta E(n)$ and $\Delta E^{\infty}(n)$,
alongside a depiction of the final structure realized for the Cs atom and a
few of the nearest-neighbor Ar atoms.

In order to find the most stable (minimum energy) configuration at a fixed $n$
value, we do not explore the entire position space for the Cs atom and the
vacancies through varying random initial configurations. Instead, we
initialize with pre-established defect structures in fcc crystals, using
positions sourced from existing literature cotterill1966morse ;
damask1968defects ; johnson1973empirical ; lam1985multiple ;
lam1986calculations ; sabochick1988atomistic ; nemirovich2007vacancies ;
wang2013defect ; zenia2016stability ; lounis2016stability ; lara2019effects ;
peng2020solute , specifically footnoteTrapping .

In anticipation of the necessity for further studies, such as the Jahn-Teller
dynamical effect, we permit the Cs atoms and the initial few (typically $\sim
60$) layers of nearest neighboring Ar atoms to move and optimize their
position to minimize $\Delta E_{N}(n)$ as much as possible. We utilize the
Adaptive Moment Estimation (Adam) kingma2014adam  and its modified version,
Nesterov-accelerated Adaptive Moment Estimation (Nadam) dozat2016incorporating
, collectively considered as ”the best overall choice” among gradient descent
optimization algorithms ruder2016overview ; haji2021comparison . Upon nearing
the optimum, employing the second-order derivative offers a faster method to
adjust the equilibrium position of the normal coordinates
simons2003introduction . This procedure is iteratively applied, optimizing
positions until a minimum is attained. Figure 4 showcases results from such
optimized geometry. However, it should be noted that using initial fixed
(pure) fcc Ar and vacancy positions suffices to reproduce the results with
satisfactory accuracy for stable site predictions, eliminating the need for
atomic position optimization.

This study’s findings contrast from the previous one carried out using a
Lennard-Jones potential ozerov2019computational ; ozerov2021generic . Our
study indicates that the most stable trapping sites probably are Td with 4
vacancies and Oh with 6 vacancies. This seems consistent with our observation
of two triplet lines in absorption.

###  III.3 Simulation of the position of the absorption lines

####  III.3.1 Line Positions in the Pairwise Approach

Contrary to the observations reported by Ref. balling1983laser , our data, as
shown in Fig. 2, does not show evidence of the Cs$(6s\rightarrow 5d)$
transition. Only Cs$(6s\rightarrow 6p)$ triplet structures appear visible,
with minimal evidence of any other Cs levels involved. Also Cs dimers seem to
be absent at the low density and temperature we are operating at. We are
therefore confident that a simple pairwise approximation of the diatomic-in-
molecule (DIM) type of approach, involving only the Ar ground state and
Cs$(6s)$, Cs$(6p)$ levels, should yield accurate results to simulate the
Cs$(6s\rightarrow 6p)$ absorption spectra.

The detailed calculation is outlined in Appendix B. Briefly, the energy shift
of the $Cs(6s)$ is calculated in a similar manner as was done for the
stability study. Now, we include the 6p manifold, where the energies are the
eigenvalues of the Hamiltonian. These are calculated by summing up the
interaction $\langle
L^{\prime}M^{\prime}|\hat{V}_{\text{Cs,Ar}}(\bm{R}_{\text{Cs,Ar}}=\\{X,Y,Z\\})|LM\rangle$
between a Cs atom (with $|6pM\rangle$ levels quantized along an arbitrary but
fixed axis $z$) and each Ar atom. This interaction is defined by the matrix
obtained using Wigner rotations:

$\frac{1}{3}(2V_{\Uppi}(R)+V_{\Upsigma}(R))\bm{I}_{3}+\frac{1}{6}\frac{V_{\Upsigma}(R)-V_{\Uppi}(R)}{R^{2}}\bordermatrix{&|M=-1\rangle&|M=0\rangle&|M=1\rangle\cr&X^{2}+Y^{2}-2Z^{2}&3\sqrt{2}(X+iY)Z&-3(X+iY)^{2}\cr&3\sqrt{2}(X-iY)Z&-2(X^{2}+Y^{2}-2Z^{2})&-3\sqrt{2}(X+iY)Z\cr&-3(X-iY)^{2}&-3\sqrt{2}(X-iY)Z&X^{2}+Y^{2}-2Z^{2}}$
(2)

Here, we also include the spin-orbit $A^{\rm
SO}(R_{\text{Cs,Ar}})\bm{L}.\bm{S}$ term, where $R=\sqrt{X^{2}+Y^{2}+Z^{2}}$.

The Hund’s case (a) potential curves $V_{\Uppi}(R)$, $V_{\Upsigma}(R)$, and
the spin-orbit coefficient $A^{\rm SO}(R)$, where $R$ denotes the Cs-Ar
separation, are constructed as shown in Fig. 3. They are based on the most
recent $V_{\rm X\,^{2}\Upsigma_{1/2}^{+}}$, $V_{\rm A\,\Uppi_{1/2}}$, $V_{\rm
A\,\Uppi_{3/2}}$, $V_{\rm B\,^{2}\Upsigma_{1/2}^{+}}$ potentials
kobayashi2016ab . The long-range part is taken from the potentials of Ref.
blank2012m+ , which provide more points, but are calculated from a smaller
basis sets and hence, are less accurate at short range.

Given that the significant variation of the spin-orbit coefficient $A^{\rm
SO}(R)$ is limited to $R$ values smaller than the distance between the closest
neighbours heaven2010potential , as seen in Fig. 3 a), the use of a constant
spin-orbit coefficient $A^{\rm SO}$ yields minimal impact (a maximum shift
difference of $10\text{\,}\mathrm{c}\mathrm{m}^{-1}$ in the final
Cs(6s)-Cs(6p) transition lines). Hence, in the following analysis, we simply
use the $V_{\Uppi},V_{\Upsigma}$ potential curves and a constant spin-orbit
coefficient $A^{\rm SO}$ based on the experimental value.

Figure 5: This figure presents a comparison between an experimental
transmittance spectrum (represented by the solid pink lines), for the two
trapping sites, Oh with 6 vacancies, and Td with 4 vacancies. The lower
section indicates the position of the $0\text{\,}\mathrm{K}$ lines employing
pairwise potential, where dots correspond to the two-body cases with
potentials sourced from kobayashi2016ab . The influence of zero point energy
on these lines is marginal (at most $150\text{\,}\mathrm{c}\mathrm{m}^{-1}$).
However, the incorporation of the third-order effect (refer to the main text)
results in an additional shift that we predict to fall within the elliptical
regions.

The big dots on Figure 5 display the results for the Cs(6s)-Cs(6p) transition
line positions. In this $0\text{\,}\mathrm{K}$ theoretical computation, only
doublets emerge, not triplets as would be anticipated in highly symmetric
trapping sites, due to the lack of thermal or dynamical broadening.
Nevertheless, the experimental and calculated positions show a reasonable
degree of agreement.

It may seem surprising initially that the most shifted lines, relative to the
gas phase case, are those with 6 vacancies. This happens while the more
constrained structure with 4 vacancies appears to be hardly affected. The
observation that the red triplet is Td and the blue one is Oh is contrary to
what has been noted (or simulated) for other systems such as Na in Ar
ryan2010investigations ; jacquet2011spectroscopic . This discrepancy can be
qualitatively explained by a simplified model where Cs is surrounded by Ar
atoms in a symmetric spherical environment.

In this spherical context, the average interaction term becomes
$2V_{\Uppi}+V_{\Upsigma}$ (refer to matrix (2)). Consequently, a valuable
approach to study the shift of the $6s-6p$ line position involves examining
$2V_{\Uppi}+V_{\Upsigma}-V_{\rm X}$. This is demonstrated in Fig. 6, along
with the positions of the nearest neighbors in both the Td and the Oh
structures as found in our stability study.

Fig. 6 highlights that the blue shift, resulting from positive $\Delta
V=2V_{\Uppi}+V_{\Upsigma}-V_{\rm X}$ values, is more pronounced in the Oh case
than in the Td scenario. This difference mainly arises because, in the Oh
situation, the first 8 and second 24 nearest neighbors (the solvation-type
shell) each contribute to a blue shift of approximately
$100\text{\,}\mathrm{c}\mathrm{m}^{-1}$. Meanwhile, in the Td case, the first
shell almost causes no shift due to the cancellation of
$2V_{\Uppi}+V_{\Upsigma}$ and $V_{\rm X}$ shifts, and the second shell only
contains 12 atoms. The cumulative effects of these shifts, coupled with the
red shifts produced by shells at greater distances, offer a qualitative
explanation for the $\sim$2000\text{\,}\mathrm{c}\mathrm{m}^{-1}$$ blue shift
in the Oh scenario. Also, the cancellation effect results in Td lines
appearing near the gas phase location.

Figure 6: This figure illustrates $\Delta V=2V_{\Uppi}+V_{\Upsigma}-V_{\rm X}$
Cs-Ar potentials with an offset of the asymptotic values. The first $R_{1}$,
second $R_{2}$, and third $R_{3}$ nearest neighbor shells in the Td and Oh
scenarios are provided (for the Oh case, an example of the Cs and Ar positions
is also given).

In line with our approach for Ar-Ar potentials, to enhance the comparison of
the line positions, we should include the zero-point energy and Third-Order
Many-Body interactions. However, as these factors merely provide corrections,
we didn’t use them to re-optimize the atomic positions.

####  III.3.2 Third order (dipole-dipole-dipole type) correction

To the best of our knowledge, the third-order correction has not been
previously utilized for simulating absorption spectra crepin1999photophysics ;
jacquet2011spectroscopic . However, this correction is non-negligible,
particularly for ground states and even more so for excited states.

We currently lack accurate three-body Cs-Ar-Ar potential curves, leading us to
approximate their contribution diaz1980combination ; dridi2022development . As
elaborated in Appendix E, we employ two and three-body expressions for the
pure long dipolar range part to extrapolate an appropriate three-body
potential also applicable for shorter distances. Following that, we aggregate
it over all atoms. Unfortunately, the nearest neighbors play a significant
role jacquet2011spectroscopic , and it is in the area where our formula proves
to be the least accurate. Consequently, we opt for a simpler approximation,
akin to the approach used by Mutô-Stenschke-Marcelli-Wang-Sadus
stenschke1994effective ; dridi2022development ; akhouri2022thermodynamic . We
average the formula, employing a mean-field approach for the position of the
third Ar atom, resulting in an effective correction to the two-body Cs-Ar
potential.

These are straightforward density-dependent expressions that connect the
averaged three-body interactions to the two-body ones. This connection
facilitates merging both types of interactions into effective Cs-Ar two-body
potential curves. With this methodology, the final Cs-Ar potential curve is a
corrected version of the original one, offset by the effective three-body
contribution $\bar{E}_{12}^{(3)}$. Ultimately, we obtain an effective two-body
potential as described by Eqs. (35) and (37).

The influence of this adjustment can be qualitatively estimated from Eq. (34).
By utilizing $2V_{\Uppi}+V_{\Upsigma}-V_{\rm X}$ as a measure of the shift, we
observe that the third-order effect $\bar{E}_{12}^{(3)}$ displays a
$\frac{8\pi\rho}{3R^{6}}(C_{9}^{*}-C_{9})>0$ dependence. This implies a
probable net blue shift, where $\rho$ is the Ar density and $C_{9},C_{9}^{*}$
represent the triple dipole coefficient for the ground and excited states,
respectively. This effect can also be interpreted by examining the potential
curves provided in Fig. 3. In this figure, the depicted third-order excited
state potential curve is computed similarly to the ground state case. Hence,
we have
$V_{\Upsigma}^{(3)}=(V_{\Upsigma}-E(6p))\frac{1}{3}\frac{8\pi\rho}{9}\frac{C_{9}^{*}}{C_{6}^{*}}$
and
$V_{\Uppi}^{(3)}=(V_{\Uppi}-E(6p))\frac{1}{3}\frac{8\pi\rho}{9}\frac{C_{9}^{*}}{C_{6}^{*}}$.

####  III.3.3 Zero Point Energy Correction

In order to incorporate the quantum effects attributable to the zero-point
energy, we utilize the Feynman-Hibbs approach, equivalent to the Wigner-
Kirkwood method. This approach introduces a temperature-dependent effective
correction to a pair potential between ground state atoms, given as: $\Delta
V_{FH}(R)=\sigma_{0}^{2}(V^{\prime\prime}(R)+2V^{\prime}(R)/R)$, where $\mu$
is the reduced mass, and $\sigma_{0}=\sqrt{\frac{\hbar^{2}}{24\mu k_{B}T}}$
represents the Gaussian width of quantum particles, as derived from the
uncertainty principle berendsen2007simulating ; stroker2022thermodynamic .

####  III.3.4 Absorption Line Positions

Upon integrating the zero-point energy and third-order correction, we find
novel line positions for the absorption spectra depicted in Fig. 5. In this
computation, for the excited states, we applied the effective two-body
potential (two-body corrected with three-body), and for the ground state, we
additionally incorporated the Feynman-Hibbs zero-point energy correction,
$\Delta V_{FH}(R)$, assuming an experimental temperature of
$8\text{\,}\mathrm{K}$. This procedure was not used for the Ar-Ar interaction,
given that we already have a well-scaled potential. However, it is intriguing
to validate this method against the known Ar-Ar case. As Fig. 3 suggests, our
method of amending a two-body potential with $\Delta V_{FH}(R)+\bar{E}^{(3)}$
accurately replicates the scaled potential we used.

A key insight is that the zero-point energy is quite negligible (on the order
of $100\text{\,}\mathrm{c}\mathrm{m}^{-1}$), whereas the third-order
correction significantly impacts the outcome and thus should not be
disregarded. Indeed, the shift generated by the third-order effect is
extremely sensitive to the selected cut-off parameters, and it can be either
attractive or repulsive, contingent on their values. As such, the induced
shift should be viewed more as an indicator of the effect rather than a
precise quantitative result, which is why we chose to represent its effect as
an uncertainty in our findings. The ellipses in Fig. 5, acting as effective
error bars, have been roughly assessed using different cut-off values and
exponents in the power of the cut-off function (see Appendix E), along with
the potential curves $V_{\Upsigma}^{(3)}$ and $V_{\Uppi}^{(3)}$ shown in Fig.
3.

Despite the substantial uncertainty of the third order, we can observe a
qualitative agreement between theoretical line positions and experimental
ones, such as the Oh lines being blue-shifted compared to the Td ones.
However, the correct line separation between Oh and Td remains elusive. Using
theoretical potential curves from Ref. blank2012m+  yielded similar results
but they were blue-shifted by approximately
$1000\text{\,}\mathrm{c}\mathrm{m}^{-1}$ compared to results obtained using
the more precise and recent potentials from Ref. kobayashi2016ab ). Still,
this does not resolve the issue of the line separation between Oh and Td,
which consistently remains around $1000\text{\,}\mathrm{c}\mathrm{m}^{-1}$,
half of the value observed in the experiment. This discrepancy indicates the
need for further investigations to reconcile the differences between
experimental and simulated positions. High-accuracy Cs-Ar potential curves are
required, and better Cs-Ar-Ar third-order calculations, particularly for the
nearest neighbor distances, should be developed.

The line position alone is a good indicator of the trapping site but is not
always sufficient to confirm it. Even though the calculated positions for the
Td 4-vacancies and Oh 6-vacancies cases align with the observation, it is
prudent to also consider other stable trapping sites on the convex hull such
as the single vacancy case, where one Ar atom is replaced by a Cs atom, and
the 10-vacancies case. The simulation of the absorption lines for the single-
vacancy substitution site results in very red-shifted theoretical lines below
$10\,000\text{\,}\mathrm{c}\mathrm{m}^{-1}$, where no experimental lines are
observed, ruling out this case. However, the 10-vacancies case predicts line
positions at $13\,760\text{\,}\mathrm{c}\mathrm{m}^{-1}$ and
$13\,210\text{\,}\mathrm{c}\mathrm{m}^{-1}$.

To study the lineshape and not just the central position, we now need to
incorporate thermal broadening of the lines due to electron-lattice
interaction.

We have already provided the results for the 1 and 10 vacancies cases in Fig.
7, confirming that they can likely be dismissed due to the large shift for 1
vacancy and the overly large and inverted triplet splitting for 10 vacancies.
Henceforth, we will focus primarily on the 4 and 6 vacancies cases.

Figure 7: Comparison between the experimental absorption spectrum and
theoretical semi-classical Mulliken formula (4) for the 1 and 10 vacancy sites
at $10\text{\,}\mathrm{K}$.

###  III.4 Lattice interaction

We must also account for the fact that the nuclei are not stationary, due to
zero-point energy and temperature effects.

####  III.4.1 Electron-lattice Interaction for Non-equilibrium Positions

We employ pair-wise potentials to compute the electronic ground state energy
$V(\bm{R})$ for the positions of the Cs and $(N-n)$ Ar nuclei
$\bm{R}=\\{\bm{R}_{1},\ldots,\bm{R}_{N-n+1}\\}$ around the equilibrium point
$(\bm{R}_{0})$, as determined by the stability study.

Subsequently, as detailed in Appendix C, we diagonalize the mass-scaled
Hessian matrix, computed via finite difference, to identify the normal mode
coordinates $\bm{Q}^{\rm n}$. The Hessian matrix provides all the oscillation
frequencies; this method is a more accurate way of including the zero-point
effect than the previously used Feynman-Hibbs formula.

####  III.4.2 Semi-classical Approximation for the Transitions

To study the laser excitation of the $(6s)$ cesium atom toward the $6p$
manifold under the dipolar operator $\hat{d}$, we start by noting (using the
Beer-Lambert-Bouguer’s law) that the spectral density optical absorption
coefficient $A(E)$ for a photon of energy $E$ is given by the sum over all
initial vibronic levels $|i\rangle$, populated with probability $P_{i}$, with
state $|\Psi_{i}\rangle$ of energy $E_{i}$, towards all possible final ones
$|\Psi_{f}\rangle$ of energy $E_{f}$:

$A(E)\propto\sum_{f,i}P_{i}\left|\langle\Psi_{f}|\hat{d}|\Psi_{i}\rangle\right|^{2}\delta(E-(E_{f}-E_{i}))$
(3)

$A(E)dE$ is the absorption coefficient for a photon in the energy band
$(E,E+dE)$. Several approximations are derived in Appendix D. For instance, we
neglect the variation of the dipole $d_{if}(\bm{Q})$ with the internuclear
distances and assume it to be constant (it is simply the 6s to 6p dipole
transition). We will primarily use two approximations:

$\displaystyle A(E)$ $\displaystyle\propto$ $\displaystyle\int P_{\rm
g}(\bm{Q})\delta[E-(V_{\rm e}(\bm{Q})-V_{\rm g}(\bm{Q}))]{\rm d}\bm{Q}$ (4)
$\displaystyle A(E)$ $\displaystyle\propto$
$\displaystyle\sum_{i}P_{i}\int\left|\Psi_{i}(\bm{Q})\right|^{2}\delta[E-(V_{\rm
e}(\bm{Q})-E_{i})]{\rm d}\bm{Q}$ (5)

The first one (equation (4)) is called the Mulliken approximation. The second
one (equation (5)) is referred to as the Reflection approximation. $V_{\rm
e}(\bm{Q})$ is the excited electron (6p) potential energy curve, and $V_{\rm
g}(\bm{Q})$ is the ground state one.

$P_{\rm g}(\bm{Q})$ can be chosen by classical statistics $P_{\rm
g}(\bm{Q})\propto e^{-V_{\rm g}(\bm{Q})/k_{B}T}$, and the formulae are thus
called classical. However, it is clearly better to choose the quantum-
statistical mechanical probability distribution $P_{\rm
g}(\bm{Q})=\sum_{i}P_{i}\left|\Psi_{i}(\bm{Q})\right|^{2}$, in which case we
call them semi-classical.

Depending on the derivation chosen, the classical Mulliken approximation can
be derived from the Reflection one with an additional approximation, as in
Ref. lax1952franck . Conversely, the reflection approximation can be derived
from the Mulliken one, as in Ref. heller1978quantum . Hence, it is not clear
which formula is the most accurate. However, proper derivation (see Appendix
D) tends to favor the (semi-)classical Mulliken approximation. Indeed, the
Mulliken difference potential $E-(V_{\rm e}(\bm{Q})-V_{\rm g}(\bm{Q}))$ is
known to be a quite good approximation for a Franck-Condon factor
mulliken1971role ; tellinghuisen1987franck  as it favors transitions where the
kinetic energy term is identical between ground and excited state. This is
where the phase of the ground and excited wavefunctions matches best, thereby
favoring their overlap, as also highlighted by the one-dimensional WKB
approximation.

Thus, we start with the semi-classical Mulliken approximation.

We will use Monte Carlo integration to estimate the integrals numerically.
Thanks to the spin-orbit interaction and the interaction matrix we have, in
the Cs(6p1/2,6p3/2) manifold, six eigenvalues $V_{i,\rm exc}(\bm{Q})$ (the
degeneracy of which depends on the symmetry of the trapping sites). The
spectra are simply given by summing the histogram of all eigenvalue cases.

####  III.4.3 Phonons thermal effects

At temperature $T$, the average number of phonons with frequency $\omega$ is
given by the Bose-Einstein distribution
$\bar{n}_{\omega}=\frac{1}{e^{\hbar\omega/k_{B}T}-1}$. The probability of
exciting this phonon mode with $n$ quanta (of energy
$E_{n}=(n+1/2)\hbar\omega$) is given by the geometric distribution
$P_{n}=e^{-n\hbar\omega/k_{B}T}\left[1-e^{-\hbar\omega/k_{B}T}\right]$. In the
semi-classical approximation, we use the true (quantum) probability
distribution function of a harmonic phonon mode of normal mode coordinate
$Q_{\omega}$ and angular oscillation frequency $\omega$, as given by
barragan2018harmonic :
$P_{\omega}(T)=\frac{1}{\sqrt{2\pi\hbar[\bar{n}_{\omega}+1/2]/\omega}}e^{-\frac{\omega^{2}Q_{\omega}^{2}}{2\hbar\omega[\bar{n}_{\omega}+1/2]}}$.
Thus, if we perform the classical Mulliken calculation with a Boltzmann
distribution at a temperature $T^{\prime}=T^{\prime}(\omega)$ such that

$\frac{T^{\prime}(\omega)}{T}=\frac{\hbar\omega}{k_{B}T}[\bar{n}_{\omega}+1/2]=\frac{\hbar\omega}{2k_{B}T}\left[\tanh{\frac{\hbar\omega}{2k_{B}T}}\right]^{-1},$
(6)

we will get the correct quantum spectral band contour in the approximation of
harmonic motion under the $\frac{1}{2}\omega^{2}Q_{\omega}^{2}$ potential, as
if we had performed the full real quantum calculation at $T$. Thus, a
classical Mulliken simulation at $T^{\prime}$ is equivalent to a semi-
classical Mulliken simulation at $T$. This scaled temperature is a way to
include phenomenologically the position oscillation due to the zero point
energy in classical simulations bergsma1984electronic . The scaled temperature
was used, for instance, to study Na in Ar using Molecular Dynamic simulation
in Ref. ryan2010investigations , where $T=$12\text{\,}\mathrm{K}$$ and so
$T^{\prime}=$45\text{\,}\mathrm{K}$$ was used for the pure argon subsystem
that thermalizes the Na deposition on other Ar layers.

In argon, the Debye frequency is
$\omega_{D}=k_{B}\times$93.3\text{\,}\mathrm{K}$/\hbar$, leading to an
effective temperature of $T^{\prime}\approx$44\text{\,}\mathrm{K}$$ at low
temperatures. Because we have several frequencies (one per mode), we used the
scaled temperature $T^{\prime}(\omega)$ for each of them and chose the mode
coordinate $Q_{\omega}$ according to the distribution
$P_{\omega}(T^{\prime}(\omega))\propto
e^{-\frac{\omega^{2}Q_{\omega}^{2}}{2k_{B}T^{\prime}(\omega)}}$. The energy
levels are then calculated for the chosen nucleus position, as done previously
for the equilibrium positions. A histogram for the energy difference between
excited levels and the ground state one finally gives the absorption
coefficient at each energy according to Eq. (4). For this simulation, only Ar
atoms within a specific cut-off distance $R_{c}$ from the Cs atom are
considered. For all simulations presented here, convergence is obtained for
$R_{c}=\qty{3}{nm}$. This corresponds to approximately 3000 Ar atoms, the
positions of which have been optimized to relax the ground state energy.
However, already $R_{c}=\qty{1}{nm}$ leads to reasonable convergence.

To benchmark our method, we compared it with Na in Ar results from
ryan2010investigations  and indeed found very similar absorption spectra for
the $T_{d}$ site.

For our Cs-Ar system, calculations are done with a bare two-body Cs-Ar
potential. We found no effect on the shape when adding third order corrections
or when using other potentials such as from blank2012m+ .

Figure 8: Comparison between experimental absorption spectra and theoretical
semi-classical Mulliken formula (eq.(4)) at different temperatures
$T=\qtylist{6;8;10;12;14;20}{K}$ and also for $T=\qty{0}{K}$ for the theory.
The theoretical spectra are created from a kernel density estimation of the
Monte Carlo simulation results. The results for the Td and Oh site have been
shifted in energy by a global offset (that are
$500\text{\,}\mathrm{c}\mathrm{m}^{-1}$ and
$-380\text{\,}\mathrm{c}\mathrm{m}^{-1}$ for Oh and Td respectively) to enable
easier comparison with the experimental data. For better visual clarity,
vertical offsets have been applied to both the experimental and the
theoretical spectra.

The results for 10000 simulated samples of the semi-classical Mulliken
simulations are given in Fig. 8 for the Oh and Td cases. However, as discussed
previously, the exact position cannot be perfectly reproduced for all
presented simulations. Therefore, we shift the theoretical spectra for each Oh
and Td cases. The experimentally observed line broadening is well reproduced,
but the triplets are not visible because the splitting of the 6p3/2 doublet is
too low to be resolved. The line positions barely evolve with temperature and
the widths only increase by less than $\qty{20}{cm^{-1}}$ for the considered
temperature range. The reason that almost no temperature effect arises is that
we have a phonon angular frequencies spectral density very similar to the one
for a pure fcc matrix, so ranging from $\sim$0\text{\,}\mathrm{K}$$ to the
Debye frequency $\sim$94\text{\,}\mathrm{K}$$. Thus, whatever the experimental
temperature $T=$6\text{\,}\mathrm{K}14\text{\,}\mathrm{K}$$, the energy
spanned $V_{\rm g}(\bm{Q})$ by the ground state motion does not vary much with
the temperature because the scaled temperature $T^{\prime}(\omega)$ is always
$T^{\prime}(\omega)\approx\frac{\hbar\omega}{2k_{B}}$ and that way it is not
very sensitive to the actual experimental temperature.

It is interesting to note that a purely classical model using $T$, not
$T^{\prime}(\omega)$, as done for Na in Ar in Ref. boatz1994monte , is able to
resolve the peak degeneracy and to show temperature effects when $T$ varies.
This is an example of a faulty model that gives erroneous results but that
might, at first glance, look reasonable.

Another intriguing aspect concerns the fact that a simple Crystal Field model
used in Ref. Kanda1971  was able to reproduce the triplet structure quite
well. Thus, we now study this Crystal Field model based only on first-order
electron-lattice interaction in order to understand why a simplification of
the model might, strangely, give better results than the more complete one we
just used.

####  III.4.4 Crystal Field: First order electron-lattice interaction

The crystal field model is based on the natural expectation that, for better
physical insight, it is worth linking the motional modes to the symmetry group
of a given trapping site using the so-called crystal field theory
pavarini2012crystal .

The normal modes $\bm{Q}^{\rm n}$ coordinates have been chosen to obtain a
diagonal Hamiltonian with ground state potential energy given by $\langle
6s|\hat{V}|6s\rangle=V_{0}+\sum_{k}\frac{1}{2}\omega_{k}^{2}{Q^{\rm
n}_{k}}^{2}$. However, within the Cs(6p) manifold, these coordinates
$\bm{Q}^{\rm n}$ are not necessarily the most appropriate anymore. Indeed, if
we restrict ourselves to the first order series in nuclear coordinates, each
of the 9 coefficients $V_{m^{\prime},m}=\langle
6pm^{\prime}|\hat{V}|6pm\rangle$ of the interaction matrix (2), when summed
over all Ar positions, contains only linear combinations of the $\bm{Q}^{\rm
n}$. Thus, at most 9 different so-called interaction modes matter, each of
them being a linear combination of the normal modes. These interaction modes
are much more appropriate than the normal modes to give a physical insight
into the Cs(6p) interaction within the Ar matrix.

They can be calculated when the interactions, so the potential curves, are
known. But, in the crystal field theory, the $Q_{\Gamma\gamma}$ interaction
mode coordinates are simply predicted using the symmetries of the trapping
site. A given trapping site is an element of a symmetry group with irreducible
representations $\Gamma$ (and its row $\gamma$). From it, we find the
$Q_{\Gamma\gamma}$, which are the linear combination of the normal modes
invariant under the $\Gamma$ symmetries. The symmetries are also reflected on
the electronic wavefunction that should have an adapted basis set
$|\Gamma\gamma\rangle$. Finally, the electron-lattice interaction can be
written at first order (so-called crystal-field approximation) as $H_{\rm
CF}=V^{(0)}({\bm{r}})+\sum_{\Gamma\gamma}V_{\Gamma\gamma}^{(1)}({\bm{r}})Q_{\Gamma\gamma}.$

Focusing on the $6p$ manifold, the symmetries of the zero-order
$V^{(0)}({\bm{r}})$ potential in the cubic Oh case suggest $\Gamma=T_{\rm
1u}$. The adapted electronic basis $|\Gamma\gamma\rangle$ is the $6p$ real
(tesseral) spherical harmonics $|x\rangle,|y\rangle,|z\rangle$. Similarly,
$\Gamma=\rm T_{\rm 1}$ applies in the $\rm T_{d}$ case.

To evaluate the first-order term
$\sum_{\Gamma\gamma}V_{\Gamma\gamma}^{(1)}({\bm{r}})Q_{\Gamma\gamma}$ within
this basis, we use the Wigner-Eckart theorem tsukerblat2006group ;
dresselhaus2007group ; gtpack2 :

$\langle\Gamma_{1}\gamma_{1}|V_{\Gamma\gamma}^{(1)}|\Gamma_{2}\gamma_{2}\rangle=\langle\Gamma_{1}||\hat{V}_{\Gamma}^{(1)}||\Gamma_{2}\rangle\langle\Gamma_{1}\gamma_{1}|\Gamma_{2}\gamma_{2},\Gamma\gamma\rangle$
(7)

This theorem delivers the selection rules and predicts the non-zero
coefficients
$\langle\Gamma_{1}\gamma_{1}|V_{\Gamma\gamma}^{(1)}|\Gamma_{2}\gamma_{2}\rangle$
that only arise for some $Q_{\Gamma\gamma}$ van1939jahn ;
toyozawa1966dynamical ; cho1968optical . The calculations’ specifics are
presented in appendix A for multiple symmetry groups in the
$|x\rangle,|y\rangle,|z\rangle$ basis. To better align with prior
calculations, we also offer the results in the $|m\rangle=|-1,0,+1\rangle$
basis. Here, the Hamiltonian matrix $M_{\rm CF}$ for the Td case (refer to Eq.
(13)) is given by:

$\begin{pmatrix}V_{A_{1}}Q_{A_{1}}-\frac{V_{E}Q_{E,2}}{\sqrt{3}}&V_{T_{2}}\frac{iQ_{T_{2},1}+Q_{T_{2},2}}{\sqrt{2}}&-V_{E}Q_{E,1}-iV_{T_{2}}Q_{T_{2},3}\cr
V_{T_{2}}\frac{-iQ_{T_{2},1}+Q_{T_{2},2}}{\sqrt{2}}&V_{A_{1}}Q_{A_{1}}+2\frac{V_{E}Q_{E,2}}{\sqrt{3}}&V_{T_{2}}\frac{-iQ_{T_{2},1}-Q_{T_{2},2}}{\sqrt{2}}\cr-
V_{E}Q_{E,1}+iV_{T_{2}}Q_{T_{2},3}&V_{T_{2}}\frac{iQ_{T_{2},1}-Q_{T_{2},2}}{\sqrt{2}}&V_{A_{1}}Q_{A_{1}}-\frac{V_{E}Q_{E,2}}{\sqrt{3}}\end{pmatrix}$
(8)

For Oh, the results are analogous, with the symmetries’ names slightly
altered. For simplicity, we maintain these notations for both Td and Oh in the
following discussions. Comparing with the matrix (2) provides physical
insight.

The $A_{1}$ mode, denoted as $Q_{A_{1}}$, which varies as $X^{2}+Y^{2}+Z^{2}$,
is a symmetrical radial mode that preserves symmetry. The $E_{g}$ mode is a
tetragonal distortion, represented by $Q_{E,1}$ varying as $X^{2}-Y^{2}$ and
$Q_{E,2}$ as $2Z^{2}-X^{2}-Y^{2}$. It transforms Oh symmetry into D4h. The
$T_{2}$ mode bends the crystal to D3d or C3v, with $Q_{T,1}$, $Q_{T,2}$, and
$Q_{T,3}$ varying as $XZ$, $YZ$, and $XY$ respectively.

Note that even if the matrix elements are linear in the interaction mode
coordinates, the eigenvalues, which are the potential curves, are not linear.
Hence, to calculate the Cs(6s-6p) absorption spectra, a numerical integration
method such as a Monte Carlo simulation is necessary. In the Mulliken formula
(4), we have $V_{\rm g}(\bm{Q})=V_{\rm Cs(6s)-Ar}(\bm{Q})+V_{\rm Ar-
Ar}(\bm{Q})$ and $V_{\rm e}(\bm{Q})=V_{\rm Cs(6p)-Ar}(\bm{Q})+V_{\rm Ar-
Ar}(\bm{Q})$ leading to

$V_{\rm e}(\bm{Q})-V_{\rm g}(\bm{Q})=V_{\rm Cs(6p)-Ar}(\bm{Q})-V_{\rm
Cs(6s)-Ar}(\bm{Q}).$ (9)

In the linear crystal field approximation, $V_{\rm Cs(6s)-Ar}(\bm{Q})$ is null
(or a constant that we offset), hence Eq. (4) simplifies to

$A(E)\propto\int\sum_{i=1}^{6}\delta[E-X_{i}(\bm{Q}^{\Upgamma})]e^{-\sum_{\Gamma\gamma}\omega_{\Gamma}^{2}Q_{\Gamma\gamma}^{2}/2k_{B}T}{\rm
d}\bm{Q}^{\Upgamma},$ (10)

where $X_{i}(\bm{Q}^{\Upgamma})$ are the eigenvalues of $M_{\rm CF}$ plus the
spin-orbit matrix (14), $A^{\rm SO}{\bm{L}}\cdot{\bm{S}}$. This can be
estimated using a Monte Carlo method where each
$\frac{\omega_{\Gamma}Q_{\Gamma\gamma}}{\sqrt{k_{B}T}}$ follows a standard
unit normal Gaussian distribution cho1968optical .

Assuming in the Boltzmann distribution that the ground state energy can be
written as $V_{\rm
g}(\bm{Q})=\frac{1}{2}\sum_{\Gamma\gamma}\omega_{\Gamma}^{2}Q_{\Gamma\gamma}^{2}$,
which isn’t usually the case as $Q_{\Gamma\gamma}$ are not the normal mode
coordinates, we found that in our study of the $6s-6p$ transition, if
projected to the sole ”active” coordinates
$Q_{A_{1}},Q_{E,1},Q_{E,2},Q_{T,1},Q_{T,2},Q_{T,3}$, the ground state
Hamiltonian indeed takes the desired quadratic diagonal form.

In this case, as used in Ref. toyozawa1966dynamical ; cho1968optical  (see
discussion in Appendix C.3), the oscillation frequencies $\omega^{\Upgamma}$
are not needed for the Mulliken semi-classical case, because we can include
them in the definition of new scaled coordinates
$\tilde{Q}_{\Gamma\gamma}=\frac{\omega_{\Upgamma}}{\sqrt{2}}Q$ with new
interaction coefficients
$\tilde{V}_{\Gamma\gamma}=\frac{\sqrt{2}}{\omega_{\Upgamma}}V$. This leads to
$V_{\Gamma\gamma}Q_{\Gamma\gamma}=\tilde{V}_{\Gamma\gamma}\tilde{Q}_{\Gamma\gamma}$,
thus Eq. (10) becomes

$A(E)\propto\int\sum_{i=1}^{6}\delta[E-X_{i}({\tilde{\bm{Q}}}^{\Upgamma})]e^{-\sum_{\Gamma\gamma}{\tilde{Q}}_{\Gamma\gamma}^{2}/k_{B}T}{\rm
d}{\tilde{\bm{Q}}}^{\Upgamma}.$ (11)

Usually, crystal field theory is used as an effective theory; the parameters
$\tilde{V}_{\Gamma\gamma}$ are adjusted to fit the data according to the
trapping site symmetry. This fitting process, along with the 6s offset for the
line position and potential fine structure splitting, typically yields a
strong agreement between experiment and theory due to the large number of free
parameters.

Compared to the Monte Carlo or molecular dynamics methods, the crystal field
method is significantly faster. It only requires a few parameters (in our
case, $\tilde{V}_{A_{1}},\tilde{V}_{E},\tilde{V}_{T_{2}}$) and Gaussian random
distributions
($\tilde{Q}_{A_{1}},\tilde{Q}_{E,1},\tilde{Q}_{E,2},\tilde{Q}_{T,1},\tilde{Q}_{T,2},\tilde{Q}_{T,3}$
in our case). Furthermore, the interaction is given by the analytical matrix
(8). It is therefore a very attractive method. As previously mentioned in Ref.
Kanda1971 , the observed triplets, similar to our experimental data, were
attributed to cubic symmetry sites and were well-fitted using the matrix (8)
with parameters listed in table 2. The fit was performed at a temperature of
$T=\qty{1.8}{K}$ (classically, without any temperature scaling). Using these
crystal field parameters, the calculated line shapes $A(E)$ from Eq. (11) well
match the observed shape, and the triplet structure can be resolved. By
varying $\tilde{V}_{A_{1}},\tilde{V}_{E},\tilde{V}_{T_{2}}$ for the two sites,
an almost perfect match can be achieved.

Site | $O_{h}$ | $T_{d}$
---|---|---
modes | $A_{1g}$ | $E_{g}$ | $T_{2g}$ | $A_{1}$ | $E$ | $T_{2}$
$\hbar\omega_{\Gamma}/k_{B}$ | 60 | 75 | 65 | 64 | 68 | 72
$\frac{(V_{\Gamma}/\omega_{\Gamma})^{2}}{k_{B}}$ | 50 | 120 | 140 | 110 | 130 | 150
from Kanda1971  $\frac{(V_{\Gamma}/\omega_{\Gamma})^{2}}{k_{B}}$ | 13600 | 84700 | 12600 | 13600 | 84700 | 12600
Table 2: Crystal field parameters, given in Kelvin, for the two types of
sites. Oscillation frequencies are such that the ground state potential is
given by $\langle
6s|\hat{V}|6s\rangle=V_{0}+\sum_{\Gamma\gamma}\frac{1}{2}\omega_{\Gamma}^{2}Q_{\Gamma\gamma}^{2}$,
while the parameters $V_{A_{1}},V_{E},V_{T_{2}}$ serve as coupling parameters
for the excited state matrix (8) incorporating linear terms of the
$V_{\Gamma}Q_{\Gamma\gamma}$ type. For comparison, values from the
experimental spectrum of Ref. Kanda1971  are included.

However, the coefficients $V_{A_{1}},V_{E},V_{T_{2}}$ can also be computed by
comparing the matrix (8) to the gradient of the interaction matrix calculated
via the summation of Cs-Ar potentials (refer to Appendix C and formulas (19)
and (20) for details). The parameter $(V_{\Gamma}/\omega_{\Gamma})^{2}$
represents a Jahn-Teller energy shift for the excited potential, evident from
the formula
$V_{\Gamma}Q_{\Gamma\gamma}+\frac{1}{2}\omega_{\Gamma}^{2}Q_{\Gamma\gamma}^{2}=\frac{1}{2}\omega_{\Gamma}^{2}\left(Q_{\Gamma\gamma}+\frac{V_{\Gamma}}{\omega_{\Gamma}^{2}}\right)^{2}-\frac{1}{2}\frac{V_{\Gamma}^{2}}{\omega_{\Gamma}^{2}}$.
Results are presented in Table 2, where frequencies $\omega_{\Gamma}$ and
interaction coefficients $V_{\Gamma}$ are displayed in temperature units using
$\hbar\omega_{\Gamma}/k_{B}$ and $(V_{\Gamma}/\omega_{\Gamma})^{2}/k_{B}$.

It is noteworthy that the $V_{\Gamma}$ parameters derived from realistic
interaction potentials markedly differ from those fitted in Ref. Kanda1971 .
This shows the risk of using solely crystal field parameters for fitting,
which may not align with physical reality, leading to potential inaccuracies
and misleading predictions for subsequent studies lund1984magnetic . It is
therefore crucial to validate them with an independent study, such as magnetic
circular dichroism (MCD) lund1984magnetic , temperature dependence of
absorption lines, or theoretical predictions ammeter1973electronic ;
bersuker1975jahn ; forstmann1977analysis ; kolb1978analysis .

Despite these findings, it is important to emphasize that using our derived
$V_{A_{1}},V_{E},V_{T_{2}}$ coefficients, and limiting our absorption band
shape Mulliken semi-classical simulation to the first order in the
$\bm{Q}^{\Upgamma}$ coordinates for the excited state interaction matrix, we
achieve crystal field results that are nearly identical to those obtained from
the full Monte Carlo classical simulation (shown in Fig. 8). Here, excited
states are calculated from the full potentials. Indeed, for our temperature
range, the first-order matrix element values in the Monte Carlo simulation
deviate less than 10% from the actual ones.

In conclusion, both the full potential approach and the crystal field first-
order approximation for the excited state matrix coefficients generate almost
identical results. However, they fail to reproduce the observed triplet
structure and its temperature evolution. The problem might originate from the
inadequacy of the classical model for ground state motion at low temperatures
where the vibrational energies are quantized.

####  III.4.5 Reflection approximation and quantization of the energy

Classical simulations, despite temperature scaling via the appropriate spatial
Gaussian distribution $\Psi_{i}(\bm{Q})$, inherently fall short as they use a
continuous energy $P_{\rm g}(\bm{Q})$ in Eq. (4). In contrast, real systems
possess quantized ground state energy $E_{i}$. For instance, at low
temperatures with only a single zero phonon mode occupation, the ground state
exhibits a single energy value, whereas classical simulation permits all
energies from $0$ to this zero-point energy.

Such a discrepancy suggests that classical simulations can induce artificial
line broadening on the order of half the Debye oscillation frequency. In
argon, this equates to roughly $\sim$30\text{\,}\mathrm{c}\mathrm{m}^{-1}$$
per vibrational (phonon) mode. While this constitutes an extreme case, given
that the majority of lattice vibrations between ground and excited states are
similar.

This indicates that the semi-classical Mulliken approximation, often employed
to simulate atomic spectra in rare gas matrices may not always be perfectly
adequate o1971jahn ; boatz1994monte ; ryan2010investigations ;
jacquet2011spectroscopic ; tao2015heat ; kleshchina2019stable . Previous
studies have shown that this approximation cannot adequately account for
observed line shapes, especially at low temperatures, for the $A\rightarrow T$
spectral band shape ($6s\rightarrow 6p$ is $A_{1_{g}}\rightarrow T_{\rm 1u}$
for $O_{h}$ or $A_{1}\rightarrow T_{\rm 1}$ for $T_{d}$) o1971jahn ;
stockmann1984influence ; bersuker2013jahn .

In most cryogenic rare gas matrix experiments, a majority of phonon modes
exhibit energies significantly larger than the thermal energy, $\hbar\omega\gg
k_{B}T$. As a result, vibronic transitions often occur from zero-phonon
occupation, a phenomenon observed in heavy atoms like Eu, Sm, or Yb embedded
in Ar jakob1976zero ; jakob1977observation ; luchner1979phonon ;
lambo2012electronic ; lambo2021high . Such cases, confined to a single $A_{1}$
mode, lead to vibrational quantization that produces well-resolved Huang-Rhys-
Pekar type peaks in the excited states due to the displacement of the harmonic
potential in the ground state wagner1964structural ; keil1965shapes ;
shi2019huang ; zhang2019applications . These peaks can only be elucidated via
a full quantum treatment of kinetic energy, as opposed to a semi-classical
Franck-Condon-Mulliken approximation.

We are dealing with a particularly complex case: the so-called
$T\otimes(a+e+t_{2})$ coupling, where a 6p triplet $T$ couples with the
lattice through several modes plus spin-orbit interaction. This situation
bears a strong similarity to the extensively studied Jahn-Teller
$T\otimes(e+t_{2})$ coupling, which has been investigated using a variety of
approximation techniques in full quantum mechanical treatments STURGE1968a ;
bersuker1975jahn ; bates1978jahn ; bates1987analysis ; bersuker2006jahn ;
bersuker2013jahn . Even with some justifiable approximations for our case (see
Table 2), of a $p$ state equally coupled to $E$ and $T$ vibrations o1971jahn ,
the problem remains intricate. Nevertheless, full quantum treatments seem to
agree with classical approaches o1980calculation ; cho1968optical .

The optimal path forward to derive a definitive answer would be to perform a
full quantum calculation of nuclear motion in the calculation of probability
transfers. However, this is outside the scope of this article, as our aim is
to use semi-classical expressions to calculate electronic absorption spectra
without the need to compute molecular vibrational eigenstates.

In our scenario, most phonon energies $\hbar\omega$ and Jahn-Teller energy
shifts $\hbar(V_{\Gamma}/\omega_{\Gamma})^{2}$ exceed the kinetic energy
temperature $k_{B}T$. Therefore, the reflection approximation in Eq. (5)
should provide greater accuracy than the Franck-Condon Mulliken approximation,
as it treats the excited-state oscillator classically while the ground state
is quantized lax1952franck ; keil1965shapes .

Regrettably, utilizing the reflection approximation gives rise to numerous
complexities. This is particularly true when employing the harmonic
approximation, as shown by: $V_{\rm g}(\bm{Q})=V_{\rm
Cs(6s)-Ar}(\bm{Q})+V_{\rm Ar-
Ar}(\bm{Q})=\sum_{k}\frac{1}{2}\omega_{k}^{2}{Q^{\rm n}_{k}}^{2}$. From this,
we must calculate the following: $V_{\rm e}(\bm{Q})-E_{i}=V_{\rm
Cs(6p)-Ar}(\bm{Q})-V_{\rm
Cs(6s)-Ar}(\bm{Q})+\sum_{k}\frac{1}{2}\omega_{k}^{2}{Q^{\rm
n}_{k}}^{2}-E_{i}$. In the above equation, the total quantized ground state
energy $E_{i}=\sum_{k}\hbar\omega_{k}(n_{k}+1/2)$ and
$\sum_{k}\frac{1}{2}\omega_{k}^{2}{Q^{\rm n}_{k}}^{2}$ both increase linearly
with the number of modes. We can bypass this nonphysical divergence by
categorizing the modes as either ”active” or ”passive”, depending on whether
they impact $V_{\rm Cs(6p)-Ar}(\bm{Q})-V_{\rm Cs(6s)-Ar}(\bm{Q})$. This can be
done using symmetry-adapted coordinates, in which only the 6 active
coordinates significantly affect the Cs interaction with the Ar crystal.

Figure 9: Comparison of observed absorption with the simulated one using the
temperature scaled first order crystal fields model. The model incorporates
the modified reflection approximation (Eq. (12)), using parameters from Table
2. The theoretical spectra are generated from a density kernel estimation of
the Monte Carlo simulation results. Results for the Td and Oh sites are
shifted in energy by a global offset for easier comparison with the
experimental data (that are $1650\text{\,}\mathrm{c}\mathrm{m}^{-1}$ and
$85\text{\,}\mathrm{c}\mathrm{m}^{-1}$ respectively for Oh and Td). Vertical
offsets have been applied to both experimental and theoretical spectra to
enhance visual clarity.

Unfortunately, the ground state Hamiltonian no longer remains diagonal in the
interaction modes, making the energy ill-defined. One possible solution to
this problem is to use the fact that, in the low temperature regime, the soft
normal coordinates phonon modes with low frequency ($\omega\gtrsim
k_{B}T/\hbar$) are populated. This population results in a different spatial
distribution for these normal modes and, thus, for their linear combinations
present in the six interaction mode coordinates. For each temperature $T$, we
scale each normal mode frequency, transforming $\omega$ into
$\omega^{\prime}$, where $\omega^{\prime 2}/T=\omega^{2}/T^{\prime}(\omega)$,
to maintain the correct spatial distribution for the normal modes. The
resulting orthonormal transformation produces new diagonal
$\omega^{\prime}_{\Gamma\gamma}(T)$ oscillation frequencies for the six
interaction mode coordinates. At this temperature $T$, we simulate a spectrum
using the classical Boltzmann distribution, $\propto
e^{-\frac{1}{2}\omega^{\prime}_{\Gamma\gamma}(T)^{2}Q_{\Gamma\gamma}^{2}/k_{B}T}$,
for each $Q_{\Gamma\gamma}$ and a classical energy
$E_{i}=\sum_{\Gamma\gamma}\frac{1}{2}\omega^{\prime}_{\Gamma\gamma}(T)^{2}Q_{\Gamma\gamma}^{2}$.
The formula used, with $\bm{Q}$ representing the interaction mode coordinates,
is as follows:

$\displaystyle
A(E)\propto\sum_{i}P_{i}\int\left|\Psi_{i}(\bm{Q})\right|^{2}\delta[E-(V_{\rm
Cs(6p)-Ar}(\bm{Q})$
$\displaystyle-\sum_{k}\frac{1}{2}{\omega^{\prime}}_{k}^{2}{Q^{\rm
n}_{k}}^{2})]{\rm d}\bm{Q}$ (12)

The results are presented in Fig. 9. This modified reflection approximation
restores the triplet shapes, but the splittings do not match closely with the
observed ones at the lowest temperature of $T=\qty{6}{K}$.

Furthermore, the simulated spectra show some temperature-dependent evolution
that does not fully align with observed data. Experimentally, the red trapping
site (attributed to Td), and the blue trapping site (attributed to Oh), behave
differently with temperature. However, the simulation predicts similar
behaviors for Oh or Td sites, as the Cs$(6p)$ interaction with the Ar matrix
is governed by similar couplings (as evidenced by the similar values for the
coupling coefficients in Table 2). Theoretically we observe, for a Gaussian
fit, a red shift of all lines (compare to blue shift of some lines in the
experiment) of nearly $\qty{30}{cm^{-1}}$ for the temperature evolution up to
$\qty{20}{K}$. We also see a quite strong modification (almost 50%) of the
theoretical width, where the experimental evolution is smaller in relative
value. Also care must be taken with experimental temperature calibration, as
the temperature of a matrix sample may differ from a temperature sensor nearby
rose1986magnetic . There are, however, more critical discrepancies that need
to be addressed. The most notable is the mismatch in the splitting between the
peaks in the experimental and theoretical data.

It remains unclear whether utilizing more sophisticated approximations for the
absorption shape would help alleviate these discrepancies. As an example, the
improved crystal field absorption line shape formula proposed in Ref.
stockmann1984influence  could be considered. Here, S-P transitions are
calculated for the most general case, where the coupling of all three
vibration modes $A_{1},E_{g},T_{2}$, as well as spin-orbit interaction, are
taken into account. However, because the $\omega_{\Gamma}$ frequencies are
almost identical for the Td and Oh case, this formula may also not produce
different temperature evolutions for the two sites.

Further, as highlighted in Ref. hupper1998uniform , semiclassical formulas
often require correction when the energy involved is comparable to the
quantization energy spacing. However, the proposed corrected terms often have
a similar order of magnitude, making the convergence unclear.

Hence, it may be more promising to explore other effects to resolve this
discrepancy, such as non-Born-Oppenheimer effects like nonadiabatic coupling.

##  IV Conclusion

In conclusion, we have observed Cs absorption spectra in an Ar matrix
containing six major resonances. These were assigned to transitions from the
$6s_{1/2}$ Cs ground state to the $6p_{1/2}$ and $6p_{3/2}$ excited states,
with lifted $|m|$-degeneracy and two different trapping environments. Pair-
wise Ar-Ar and Cs-Ar potentials facilitated a stability study, indicating that
the two observed triplets are likely due to a 4 vacancy tetrahedral Td
trapping site and a 6 vacancy cubic Oh one in an fcc Ar matrix.

To qualitatively reproduce the line locations using the pair-wise approach, we
derived effective two-body potential curves. These included zero-point energy
and third-order terms. We found that zero-point energy has a minor effect, but
adding third-order many-body corrections significantly modifies the
potentials. Our derivation of the effective two-body potential from the long-
range part of the triple dipolar interaction, based on an exact derivation of
the long-range part, is not very accurate. More work is needed on effects such
as the Pauli repulsion, exchange interaction, and three-body dispersion
interactions for the rare gases - alkali atoms, especially in excited states
and for the short-range part anatole2010two .

Our study reasonably assigns the nearly unshifted (relative to the gas phase)
triplet to the 4 vacancy $T_{\mathrm{d}}$ site, and the blue-shifted one to
the 6 vacancy cubic $O_{\mathrm{h}}$ trapping site. Using the computationally
demanding core polarization pseudopotentials will be an interesting approach
to confirm this fact and to achieve better agreement between theory and
experiment for the line separations gross1998pseudopotential ;
jacquet2011spectroscopic .

Within our current approach, we demonstrate that a triplet shape can only be
reproduced by using quantized energy for the ground state vibrational motions.
Indeed, the classical simulation may produce artificial broadening of the
lines, on the order of half the Debye oscillation frequency. Therefore, in Ar
this would be on the order of $\sim 30\,cm^{-1}$. Hence, our semi-classical
Mulliken simulation fails to produce a triplet structure, instead producing a
broad doublet. A simulation using the reflection approximation, on the other
hand, only yields a triplet at $0\text{\,}\mathrm{K}$, but otherwise diverges
with increasing numbers of Ar atoms due to unphysical treatment of the
phononic modes.

Finally, a quantized simulation based on a modified reflection approximation
reproduces a triplet structure. We show that modifying the population of the
soft (low-frequency) phonon modes can also alter the absorption shapes. Using
only the matrix gradient (linear, first-order approximation) for the excited
state does not alter the quality of the results but leads to significantly
faster calculations. Moreover, at low temperatures when only ground state
vibrational (zero-phonon) levels are populated, the simplest crystal field
model suffices to reproduce the observed line broadening caused by the dynamic
Jahn-Teller thermal effect. This requires only three electron-lattice
parameters given in Table 2 (together with Eq. (10) and diagonalization of
spin-orbit and (8) matrices) cho1968optical . The reduced number of modes is
advantageous as it leads to faster calculations than the full Monte Carlo
method and captures the essential aspects of the electrostatic interaction,
providing a powerful starting point for further studies.

However, our work was unable to reproduce the triplet splitting and did not
provide accurate red and blue line separations. The reason for this
discrepancy is unclear, but could be due to our neglect of the dipole
transition dependence on the $Q$ coordinates, lack of coupling with other
excited states (like Cs$(5d)$), or an inappropriate inclusion of the third-
order non-additive terms in the potentials. The role of nonadiabatic couplings
should also be studied.

Following ryan2010investigations , modifying the potential curves, such as the
B$\Upsigma$ state (recently suspected to be inaccurate hewitt2023csar ), might
lead to better agreement between experiment and theory. However, our
preliminary attempts suggest that substantial modification of the position and
the slope of the curve in the region near the first nearest neighbors (see
Fig. 6) may be necessary. This should be further investigated, along with
experimental and theoretical studies of other possibilities such as trapping
in other crystal defects of lower symmetry (for instance, 5 vacancy of D3h
symmetry in the hcp phase davis2018investigation ; ozerov2021generic ), which
may explain the discrepancy between experiment and simulation.

Experimentally, many tools are available such as magnetic circular dichroism,
fluorescence, electron (paramagnetic) spin resonance, oscillator strength
determination, time-resolved emission spectra, nonradiative branching ratios,
optical pumping, or bleaching studies. Comparing with the rubidium case and
using other rare gases, or even mixed ones, would be interesting
corbin1994optical , to confirm that the matrix shift seen in the S-P
transitions is approximately linear with the polarizability of the matrix host
laursen1991multiplicity ; ryan2010investigations ; tao2015heat .

We hope that our work will be useful for further studies of similar systems
used for fundamental physics experiments in rare gas solid crystals
guarise2022particle ; lambo2022calculation ; koyanagi2022accurate .

##  V Acknowledgements

T.B. and S.L. contributed equally to this work.

This research was financed in whole or in part by Agence Nationale de la
Recherche (ANR) under the project ANR-21-CE30 -0028-01.

A CC-BY public copyright license has been applied by the authors to the
present document and will be applied to all subsequent versions up to the
Author Accepted Manuscript arising from this submission, in accordance with
the grant’s open access conditions.

We acknowledge W. Chin, B. Darquier, O. Dulieu, B. Gervais, H. Lignier, Ch.
Malbrunot, C. Smai, B. Viaris for fruitful discussions, and B. Vivan and L.
Marriaux for the design and mechanical realization.

##  Appendix A Perturbation matrices for different symmetries

A given trapping site is defined by the near neighbor environment of the
impurity atom and by its overall crystal orientation in the laboratory frame.
We detail here the crystal field method that we use to simulate the absorption
spectrum.

Using vector notation for the cartesian coordinate of each position of the Cs
and the $(N-n)$ Ar nuclei $\bm{R}=\\{\bm{R}_{1},\ldots,\bm{R}_{3(N-n+1)}\\}$)
and $\bm{r}$ as the Cs valence electron coordinates, the electron-lattice
interaction can be written as $V=V(\bm{r},\bm{R})$, where for simplicity, we
use, in this article, the same notation for operators and their
representations in coordinate basis.

###  A.1 Crystal field interpretation

A given single trapping site, at the equilibrium position $\bm{R}_{0}$,
possesses a certain symmetry group with irreducible representation $\Gamma$
(and its row $\gamma$). So it is more convenient to use coordinates,
$Q_{\Gamma\gamma}$, that transforms well under such a representation. They are
linear combinations of the $R_{i}-{R_{0}}_{i}$ ones, such that at a linear
(first order) approximation, in the nuclear coordinate $\bm{R}$ or
$\bm{Q}^{\Upgamma}$, the electron-lattice interaction can be written as
$\hat{V}=V^{(0)}({\bm{r}})+\sum_{\Gamma\gamma}V_{\Gamma\gamma}^{(1)}({\bm{r}})Q_{\Gamma\gamma}$.
We refer to this expression as the ”crystal-field approximation”.

At zero order, the electron lattice $V^{(0)}({\bm{r}})=V(\bm{r},\bm{R}_{0})$
modifies the Cs electron wavefunction throughout the $\langle
n_{1}l_{1}m_{1}|V^{(0)}|n_{2}l_{2}m_{2}\rangle$ couplings (we use here the
non-relativistic wavefunctions, without fine structure that depend on the
spherical harmonics $Y_{l}^{m}$). The trapping site equilibrium position
$\bm{R}_{0}$ determines the symmetry group for the nuclear environment that
produces the electron lattice interaction $V_{0}(\bm{r})$. For simplicity, we
keep the same notation, as for the nuclear coordinate representation $\Gamma$,
and name the eigenstates of the electronic Hamiltonian that are linear
combinations of the $|nlm\rangle$ states, $|\Gamma\gamma\rangle$. They are
given by the projection operators for the irreducible representation $\Gamma$.

Once the Cs electronical eigenstates $|\Gamma\gamma\rangle$, that are adapted
for the ground state equilibrium position, are found, we can use them to
calculate the next order of the electron lattice interaction
$\sum_{\Gamma\gamma}V_{\Gamma\gamma}^{(1)}({\bm{r}})Q_{\Gamma\gamma}$.

The key result tsukerblat2006group ; dresselhaus2007group  is that the
$\langle\Gamma_{1}\gamma_{1}|V_{\Gamma\gamma}^{(1)}|\Gamma_{2}\gamma_{2}\rangle$
coupling terms can be calculated using the Wigner-Eckart’s theorem
$\langle\Gamma_{1}\gamma_{1}|V_{\Gamma\gamma}^{(1)}|\Gamma_{2}\gamma_{2}\rangle=\langle\Gamma_{1}||\hat{V}_{\Gamma}^{(1)}||\Gamma_{2}\rangle\langle\Gamma_{1}\gamma_{1}|\Gamma_{2}\gamma_{2},\Gamma\gamma\rangle$.
Non zero Clebch-Gordan coefficient
$\langle\Gamma_{1}\gamma_{1}|\Gamma_{2}\gamma_{2},\Gamma\gamma\rangle$ occurs
only if the $\Gamma$ representation is part of the
$\Gamma_{1}\times\Gamma_{2}$ representation tsukerblat2006group ;
dresselhaus2007group . Since $V$ is real, we only need to consider the
symmetric part of the product: $[\Gamma_{1}\times\Gamma_{2}]$.

Thus, group theory automatically leads to the proper basis
$|\Gamma\gamma\rangle$ and provides directly the useful representations
$\Gamma$ with non-zero matrix elements
$V_{\Gamma}=\langle\Gamma_{1}||\hat{V}_{\Gamma}^{(1)}||\Gamma_{2}\rangle$ and
the Clebch-Gordan coupling coefficients.

We will now illustrate this crystal field procedure for several possible
trapping site symmetries that might occur for atoms in matrices
cotterill1966morse ; damask1968defects ; johnson1973empirical ;
lam1985multiple ; lam1986calculations ; sabochick1988atomistic ;
nemirovich2007vacancies ; wang2013defect ; zenia2016stability ;
lounis2016stability ; lara2019effects ; peng2020solute . For this we use the
free Mathematica group theory package GTPack gtpack1 ; gtpack2 .

###  A.2 Cubic: Oh or Td symmetry

In fcc crystals, the 1 vacancy substitution and the 6 vacancies site have both
Oh symmetry (same symmetry as the SF6 molecule) while the 4 vacancies site has
a Td symmetry (as CH4).

For such symmetries, the electron-lattice potential at equilibrium
$V^{(0)}({\bm{r}})$ does not lift the 6s nor the 6p degeneracy. This can be
understood in a simple manner by decomposing it in spherical harmonics:
$V^{(0)}({\bm{r}})=\sum_{l,m}r^{l}a_{l}^{m}Y_{l}^{m}(\theta,\varphi)$. We now
only have to take into account that this interaction potential has necessarily
the same Oh or Td symmetry as the nuclear position that creates it. Thus most
of the terms are zero, and the potential becomes
$V^{(0)}({\bm{r}})=a_{0}Y_{0}^{0}+r^{4}a_{4}(Y_{4}^{-4}+\sqrt{14/5}Y_{4}^{0}+Y_{4}^{4})+\cdots$
($a_{0}$ is not to be confused with the Bohr radius). These leading order
terms formed the so-called Devonshire’s potential devonshire1936rotation . The
triangle inequality $|l_{1}-l_{2}|\leq l\leq l_{1}+l_{2}$ for the Clebsch-
Gordan coefficients appearing when calculating $\langle
n_{1}l_{1}m_{1}|V^{(0)}=\sum_{l,m}r^{l}a_{l}^{m}Y_{l}^{m}|n_{2}l_{2}m_{2}\rangle$
indicates that the $l\geq 4$ terms cannot couple to $6s$ ($l=0$) nor to $6p$
($l=1$) states but would lift the degeneracy of the 5d ($l=2$) level for
instance.

For completeness, we indicate that the same calculation can be done elegantly
using directly the quantum states $|nlm\rangle$, without using spherical
harmonics, by directly writing the potential $\hat{V}(\bm{r},\bm{R}_{0})$
using angular momentum operators. For this, we simply have to replace, in the
Devonshire’s potential expression, $Y_{0}^{0}\rightarrow 1$, $Y_{4}^{\pm
4}\rightarrow L_{\pm}^{4}$ and $Y_{4}^{0}\rightarrow
35L_{z}^{4}-(30L(L+1)-25)L_{z}^{2}+2L^{2}(L+1)^{2}-6L(L+1)$ where $\bm{L}$ is
the orbital angular momentum operator misra2011spin ; ryabov2009operator .

It would be possible to use the $|6sm=0\rangle$ and $|6pm\rangle=|m\rangle$
basis but for cubic symmetry the so-called real (or tesseral) spherical
harmonics:
$|x\rangle=(|m=-1\rangle-|m=+1\rangle)/\sqrt{2},\;|y\rangle=i(|m=-1\rangle+|m=+1\rangle)/\sqrt{2},\rangle,\;|z\rangle=|m=0\rangle$
is better. Indeed, the projection operator technique for the angular momentum
basis $|6pm\rangle$ under the symmetry group Oh or Td and their
representations $\Gamma$ leads to the following: for the 6s electronic state
$\Gamma=A_{1_{g}}$ (or $\Gamma=A_{1}$ for $T_{d}$) and $|6sm=0\rangle$ forms a
good basis. For the 6p electronics state $\Gamma=T_{\rm 1u}$ (or
$\Gamma=T_{\rm 1}$ in the $\rm T_{d}$ case) and the basis
$|\Gamma\gamma\rangle$ is just given by $|x\rangle,\;|y\rangle,\;|z\rangle$.

Now that the zero-order symmetry $|\Gamma\gamma\rangle$ basis adapted to the
nuclear equilibrium configuration is established, we can use it to calculate
the first-order correction
$\sum_{\Gamma\gamma}\langle\Gamma_{1}\gamma_{1}|V_{\Gamma\gamma}^{(1)}({\bm{r}})|\Gamma_{2}\gamma_{2}\rangle
Q_{\Gamma\gamma}$. Only the non zero $\langle 6p=T_{\rm
1u}||V^{(1)}||6p=T_{\rm 1u}\rangle$ terms are interesting. Then because
$[T_{1u}\times T_{1u}]=A_{1g}\oplus E_{g}\oplus T_{2g}$ (or $[T_{1}\times
T_{1}]=A_{1}\oplus E\oplus T_{2}$ in the Td case ), the
$\langle\Gamma_{1}\gamma_{1}|V_{\Gamma\gamma}|\Gamma_{2}\gamma_{2}\rangle
Q_{\Gamma\gamma}$ calculation leads to the matrix (with the $u,g$ label
removed in the $T_{d}$ case), in the 6p manifold, for the total
$\langle\Gamma_{1}\gamma_{1}|\hat{V}=V^{(0)}+V^{(1)}|\Gamma_{2}\gamma_{2}\rangle
Q_{\Gamma\gamma}$ crystal field (CF) interaction:

$M_{\rm
CF}=\left(\frac{a_{0}}{\sqrt{4\pi}}+V_{A_{1g}}Q_{A_{1g}}\right)\bm{I}_{3}+\bordermatrix{&|x\rangle&|y\rangle&|z\rangle\cr&V_{E_{g}}\left(Q_{E_{g},1}-\frac{1}{\sqrt{3}}Q_{E_{g},2}\right)&V_{T_{2g}}Q_{T_{2g},3}&V_{T_{2g}}Q_{T_{2g},2}\cr&V_{T_{2g}}Q_{T_{2g},3}&-V_{E_{g}}\left(Q_{E_{g},1}+\frac{1}{\sqrt{3}}Q_{E_{g},2}\right)&V_{T_{2g}}Q_{T_{2g},1}\cr&V_{T_{2g}}Q_{T_{2g},2}&V_{T_{2g}}Q_{T_{2g},1}&V_{E_{g}}\frac{2}{\sqrt{3}}Q_{E_{g},2}\cr}$
(13)

Where we have used
$V_{A_{1}}=\langle\Gamma_{1}||\hat{V}_{A_{1}}^{(1)}||\Gamma_{2}\rangle/\sqrt{3}$,
$V_{E_{g}}=\langle\Gamma_{1}||\hat{V}_{E_{g}}^{(1)}||\Gamma_{2}\rangle/\sqrt{2}$,
$V_{T_{2g}}=\langle\Gamma_{1}||\hat{V}_{T_{2g}}^{(1)}||\Gamma_{2}\rangle/\sqrt{2}$
and performed an orthonormal transformation for the
$\\{Q_{E_{g},1},Q_{E_{g},2}\\}$ coordinates to obtain the same matrix as the
one given in Ref. cho1968optical .

This matrix highlights again the fact that, with no thermal effect, so with
$Q_{\Gamma\gamma}=0$, the 6p levels are not lifted by the crystal field, only
a ($a_{0}$) global shift of the free transition $6s\rightarrow 6p$ will occur.

In order to complete the calculation, we add the spin degree of freedom with
the spin $1/2$ basis sets $|\pm\rangle$.

The spin-orbit matrix is

$H_{\rm SO}=\frac{A^{\rm
SO}}{2}\bordermatrix{&|x+\rangle&|y+\rangle&|z+\rangle&|x-\rangle&|y-\rangle&|z-\rangle\cr&0&-i&0&0&0&1\cr&i&0&0&0&0&-i\cr&0&0&0&-1&i&0\cr&0&0&-1&0&i&0\cr&0&0&-i&-i&0&0\cr&1&i&0&0&0&0\cr}$
(14)

And the crystal field matrix takes a simple block form:

$H_{\rm CF}=\begin{pmatrix}M_{\rm CF}&0\\\ 0&M_{\rm CF}\end{pmatrix}$

###  A.3 Hexagonal: D3h, C3v

####  A.3.1 Hexagonal: D3h

The D3h symmetry (as C2H6 Eclipsed Ethane), can arise from 5 vacancy in hcp
phase, for instance. Its basis vectors are
$|E^{\prime},1\rangle,\;|E^{\prime},2\rangle$, and
$|A^{\prime\prime}_{2}\rangle=|z\rangle$. The crystal field is, at the lowest
order, $V^{(0)}({\bm{r}})=a_{0}^{0}Y_{0}^{0}+r^{2}a_{2}Y_{2}^{0}+\cdots$ (with
$\hat{Y}_{2}^{0}=3L_{z}^{2}-L(L+1)$).

And with $[E^{\prime}\times E^{\prime}]=A^{\prime}_{1}\oplus E^{\prime}$ we
have

$M_{\rm
CF}=\left(\frac{a_{0}}{\sqrt{4\pi}}\right)\bm{I}_{3}+\bordermatrix{&|E^{\prime},1\rangle&|E^{\prime},2\rangle&|z\rangle\cr&-\frac{\langle
r^{2}\rangle
a_{2}}{\sqrt{20\pi}}+V_{E^{\prime}}Q_{E^{\prime},1}&V_{A^{\prime}_{1}}Q_{A^{\prime}_{1}}&0\cr&V_{A^{\prime}_{1}}Q_{A^{\prime}_{1}}&-\frac{\langle
r^{2}\rangle
a_{2}}{\sqrt{20\pi}}+V_{E^{\prime}}Q_{E^{\prime},2}&0\cr&0&0&\frac{2\langle
r^{2}\rangle a_{2}}{\sqrt{20\pi}}\cr}$

where $\langle r^{2}\rangle=\langle 6p|r^{2}|6p\rangle$ is the square of the
size of a Cs atom in its 6p level.

Here the $6p$ ($l=1$) degeneracy is partially lifted into two levels by the
crystal field and only the spin-orbit will create the triplet structure that
is then broadened and modified by the $A^{\prime}_{1}$ and $E^{\prime}$ terms.
This symmetry would likely produce a triplet structure.

####  A.3.2 Hexagonal: C3v

The C3v symmetry (as Ammonia NH3 molecule), can arise in 4 vacancy in an hcp
phase. It is very similar to the D3h symmetry case and can be obtained from it
through the substitutions $E^{\prime}\rightarrow E$,
$A^{\prime}_{1}\rightarrow A_{1}$, and $A^{\prime\prime}_{2}\rightarrow
A_{1}$. Therefore, the basis becomes
$|E,1\rangle,|E,2\rangle,|A_{1}\rangle=|z\rangle$ and the crystal field that
at lowest order is given by:
$V^{(0)}({\bm{r}})=a_{0}^{0}Y_{0}^{0}+r^{2}a_{2}Y_{2}^{0}+\cdots$. With
$[E\times E]=A_{1}\oplus E$ we have

$M_{\rm
CF}=\left(\frac{a_{0}}{\sqrt{4\pi}}\right)\bm{I}_{3}+\bordermatrix{&|E,1\rangle&|E,2\rangle&|z\rangle\cr&-\frac{\langle
r^{2}\rangle
a_{2}}{\sqrt{20\pi}}+V_{E}Q_{E,1}&V_{A_{1}}Q_{A_{1}}&0\cr&V_{A_{1}}Q_{A_{1}}&-\frac{\langle
r^{2}\rangle a_{2}}{\sqrt{20\pi}}+V_{E}Q_{E,2}&0\cr&0&0&\frac{2\langle
r^{2}\rangle a_{2}}{\sqrt{20\pi}}\cr}$

This is identical to the previous case, so this symmetry would likely produce
a triplet structure.

###  A.4 Tetragonal: C4v

The C4v symmetry ( as SF5Cl or IF5), can arise in the 10 vacancy case. Its
basis is $|E,1\rangle=|x\rangle$, $|E,2\rangle=|y\rangle$, and
$|A_{1}\rangle=|z\rangle$. The crystal field is at the lowest order:
$V^{(0)}({\bm{r}})=a_{0}^{0}Y_{0}^{0}+r^{2}a_{2}Y_{2}^{0}+\cdots$. And with
$[E\times E]=A_{1}\oplus B_{2}\oplus B_{1}$ we get

$M_{\rm
CF}=\left(\frac{a_{0}}{\sqrt{4\pi}}+V_{A_{1g}}Q_{A_{1g}}\right)\bm{I}_{3}+\bordermatrix{&|x\rangle&|y\rangle&|z\rangle\cr&-\frac{\langle
r^{2}\rangle
a_{2}}{\sqrt{20\pi}}-V_{B_{2}}Q_{B_{2}}&V_{B_{1}}Q_{B_{1}}&0\cr&V_{B_{1}}Q_{B_{1}}&-\frac{\langle
r^{2}\rangle a_{2}}{\sqrt{20\pi}}+V_{B_{2}}Q_{B_{2}}&0\cr&0&0&\frac{2\langle
r^{2}\rangle a_{2}}{\sqrt{20\pi}}\cr}$

This is very similar to the previous cases, and this symmetry would likely
produce a triplet structure.

###  A.5 Lower symmetry: C2v

The C2v symmetry (as H2O water molecule), can arise in the 8 vacancy case and
has the basis
$|B_{2}\rangle=|x\rangle,|B_{1}\rangle=|y\rangle,|A_{1}\rangle=|z\rangle$. The
crystal field is at the lowest order:
$V^{(0)}({\bm{r}})=a_{0}^{0}Y_{0}^{0}+r^{2}[a_{2}^{0}Y_{2}^{0}+a_{2}^{2}(Y_{2}^{-2}+Y_{2}^{2})]+\cdots$.
(with $\hat{Y}_{2}^{\pm 2}=L_{\pm}^{2}$) and its matrix takes a diagonal form:

$M_{\rm
CF}=\left(\frac{a_{0}}{\sqrt{4\pi}}+V_{A_{1g}}Q_{A_{1g}}\right)\bm{I}_{3}+\frac{\langle
r^{2}\rangle}{\sqrt{20\pi}}\bordermatrix{&|x\rangle&|y\rangle&|z\rangle\cr&-a_{2}^{0}+\sqrt{6}a_{2}^{2}&0&0\cr&0&-a_{2}^{0}-\sqrt{6}a_{2}^{2}&0\cr&0&0&2a_{2}^{0}\cr}$

Here the $6p$ ($l=1$) degeneracy is fully lifted in the three levels by the
crystal field. The spin-orbit will also shift the levels and the triplet
structure will be broadened in a symmetric way by the $A_{1g}$ terms.

##  Appendix B Potential curves calculation

We discuss here the Cs-Ar potential for Cs in its 6p excited level, the origin
of the energy is taken at the 6p level in gas phase.

###  B.1 Hund’s case (a) and (c) curves

Without spin-orbit interaction the Cs$(6p)$-Ar interaction is diagonal in the
basis $|\Uppi,\pm\rangle=|L=1,M=\pm 1\rangle$,
$|\Upsigma\rangle=|L=1,M=0\rangle$ with $V_{\Uppi}(R),V_{\Upsigma}(R)$ Hund’s
case (a) potential curves, where $R$ is the Cs-Ar separation.

When adding the spin-orbit $A^{\rm SO}(R)\bm{L}.\bm{S}$, the Hamiltonian
becomes

$\bordermatrix{&|\Uppi,+\rangle|m_{s}=-1/2\rangle&|\Upsigma\rangle|m_{s}=1/2\rangle&|\Uppi,+\rangle|m_{s}=1/2\rangle\cr&V_{\Uppi}-\frac{A^{\rm
SO}}{2}&\frac{A^{\rm SO}}{\sqrt{2}}&0\cr&\frac{A^{\rm
SO}}{\sqrt{2}}&V_{\Upsigma}&0\cr&0&0&V_{\Uppi}+\frac{A^{\rm SO}}{2}\cr}$

and the matrix will be identical for the basis with opposite signs for all
values of the projection of angular momenta:
$|\Uppi,-\rangle|m_{s}=1/2\rangle,|\Upsigma\rangle|m_{s}=-1/2\rangle,|\Uppi,-\rangle|m_{s}=-1/2\rangle$.
The spin-orbit parameter $A^{\rm SO}(R)$ is not in general constant and
deviates from the asymptotic atomic value as shown in Fig. 3
lewis2017theoretical . Taking the eigenvalues leads to the potential curves
with spin-orbit
$V_{\Uppi_{\frac{1}{2}}},V_{\Uppi_{\frac{3}{2}}},V_{\Upsigma_{\frac{1}{2}}}$
with the correspondence:

$\displaystyle V_{\Upsigma}$ $\displaystyle=$
$\displaystyle\frac{1}{3}(2V_{\Uppi_{\frac{1}{2}}}-V_{\Uppi_{\frac{3}{2}}}+2V_{\Upsigma_{\frac{1}{2}}}+\sqrt{V_{\Uppi_{\frac{1}{2}}}^{2}+2V_{\Uppi_{\frac{1}{2}}}V_{\Uppi_{\frac{3}{2}}}-2V_{\Uppi_{\frac{3}{2}}}^{2}-4V_{\Uppi_{\frac{1}{2}}}V_{\Upsigma_{\frac{1}{2}}}+2V_{\Uppi_{\frac{3}{2}}}V_{\Upsigma_{\frac{1}{2}}}+V_{\Upsigma_{\frac{1}{2}}}^{2}})$
(15) $\displaystyle V_{\Uppi}$ $\displaystyle=$
$\displaystyle\frac{1}{6}(V_{\Uppi_{\frac{1}{2}}}+4V_{\Uppi_{\frac{3}{2}}}+V_{\Upsigma_{\frac{1}{2}}}-\sqrt{V_{\Uppi_{\frac{1}{2}}}^{2}+2V_{\Uppi_{\frac{1}{2}}}V_{\Uppi_{\frac{3}{2}}}-2V_{\Uppi_{\frac{3}{2}}}^{2}-4V_{\Uppi_{\frac{1}{2}}}V_{\Upsigma_{\frac{1}{2}}}+2V_{\Uppi_{\frac{3}{2}}}V_{\Upsigma_{\frac{1}{2}}}+V_{\Upsigma_{\frac{1}{2}}}^{2}})$
$\displaystyle A^{\rm SO}$ $\displaystyle=$
$\displaystyle\frac{1}{3}(-V_{\Uppi_{\frac{1}{2}}}+2V_{\Uppi_{\frac{3}{2}}}-V_{\Upsigma_{\frac{1}{2}}}+\sqrt{V_{\Uppi_{\frac{1}{2}}}^{2}+2V_{\Uppi_{\frac{1}{2}}}V_{\Uppi_{\frac{3}{2}}}-2V_{\Uppi_{\frac{3}{2}}}^{2}-4V_{\Uppi_{\frac{1}{2}}}V_{\Upsigma_{\frac{1}{2}}}+2V_{\Uppi_{\frac{3}{2}}}V_{\Upsigma_{\frac{1}{2}}}+V_{\Upsigma_{\frac{1}{2}}}^{2}})$

From this correspondence we have extracted the $V_{\Upsigma}(R)$ and
$V_{\Uppi}(R)$ (and $A^{\rm SO}(R)$) potential curves from the
$V_{\Uppi_{\frac{1}{2}}},V_{\Uppi_{\frac{3}{2}}},V_{\Upsigma_{\frac{1}{2}}}$
ones. The calculation is done using the most recent $V_{\rm
X\,^{2}\Sigma_{1/2}^{+}}$, $V_{\rm A\,\Uppi_{1/2}}$, $V_{\rm A\,\Uppi_{3/2}}$,
$V_{\rm B\,^{2}\Sigma_{1/2}^{+}}$ potentials (in Hund’s case (a) notation, but
calculated with fine structure included) with spin-orbit interactions from
kobayashi2016ab  and with a cubic interpolation between points and above
$2\text{\,}\mathrm{n}\mathrm{m}$ the long-range part is taken from blank2012m+
.

###  B.2 Cs-Ar sum

Using the primitive cell and lattice vectors, we sum the interaction $\langle
L^{\prime}M^{\prime}|\hat{V}_{Cs,Ar}(\bm{R}_{Cs,Ar})|LM\rangle$ between a Cs
atom and all Ar ones. Where $|LM\rangle=|LM\rangle_{z}$ is quantized along the
fixed $z$ axis.

We then use the fixed right-handed frame, right-hand screw counterclockwise
rule active interpretation (extrinsic) rotation z y z convention ${\cal
R}(\alpha,\beta,\gamma)=e^{-i\alpha L_{z}}e^{-i\beta L_{y}}e^{-i\gamma L_{z}}$
varshalovich1988quantum ; man2016wigner  such that ${\cal R}$ rotates the
$|LM\rangle_{z}$ state to the $|LM\rangle_{{\cal R}(z)}$ where the new axis
$z^{\prime}={\cal R}(z)$ is along $\bm{R}_{Cs,Ar}$. Thus $\bm{L}_{{\cal
R}(z)}={\cal R}\bm{L}_{z}{\cal R}^{\dagger}$. The calculation is done using
active rotation of the state vector and Wigner D-matrix with the convention:
${}_{z}\langle LM^{\prime}|{\cal
R}(\alpha,\beta,\gamma)|LM\rangle_{z}=D^{L}_{M^{\prime}M}(\alpha,\beta,\gamma)=e^{-i\alpha
m^{\prime}}{}_{z}\langle L^{\prime}M^{\prime}|e^{-i\beta
L_{y}}|LM\rangle_{z}e^{-i\gamma m}={\rm
WignerD}[\\{L,M^{\prime},M\\},-\alpha,-\beta,-\gamma]$ (the last notation
being adapted for the Mathematica software). Using the spherical coordinate,
polar angle $\theta$ and azimuthal angle $\varphi$, for the vector
$\bm{R}_{Cs,Ar}$ (from Cs to Ar) we have
$\alpha=\varphi,\beta=\theta,\gamma=0$.

Finally, we find the following matrix for the $\langle
L^{\prime}M^{\prime}|\hat{V}_{Cs,Ar}(\bm{R}_{Cs,Ar})|LM\rangle$ interaction
boatz1994monte ; kenney2005theory :

$\bordermatrix{&|\Uppi,-\rangle&|\Upsigma\rangle&|\Uppi,+\rangle\cr&\frac{3V_{\Uppi}+V_{\Upsigma}+(V_{\Uppi}-V_{\Upsigma})\cos(2\theta)}{4}&\frac{e^{i\varphi}(V_{\Upsigma}-V_{\Uppi})\cos\theta\sin\theta}{\sqrt{2}}&\frac{e^{2i\varphi}(V_{\Uppi}-V_{\Upsigma})\sin^{2}\theta}{2}\cr&\frac{e^{-i\varphi}(V_{\Upsigma}-V_{\Uppi})\cos\theta\sin\theta}{\sqrt{2}}&V_{\Upsigma}\cos^{2}\theta+V_{\Uppi}\sin^{2}\theta&\frac{e^{i\varphi}(V_{\Uppi}-V_{\Upsigma})\cos\theta\sin\theta}{\sqrt{2}}\cr&\frac{e^{-2i\varphi}(V_{\Uppi}-V_{\Upsigma})\sin^{2}\theta}{2}&\frac{e^{-i\varphi}(V_{\Uppi}-V_{\Upsigma})\cos\theta\sin\theta}{\sqrt{2}}&\frac{3V_{\Uppi}+V_{\Upsigma}+(V_{\Uppi}-V_{\Upsigma})\cos(2\theta)}{4}\cr}$
(16)

for a vector $\bm{R}$ between the Cs atoms and the Ar atoms with $X,Y,Z$
cartesian coordinate where $\cos(\theta)=Z/R,e^{\pm
i\varphi}\sin(\theta)=(X\pm iY)/R$. Its cartesian form is given by Eq. 2.

We then sum over all Ar atoms, using $V_{\Upsigma}$ and $V_{\Uppi}$ potentials
without spin-orbit interactions in Eq. 2, and after we add a constant spin-
orbit interaction that is

$H_{\rm SO}=A^{\rm
SO}\bordermatrix{&|\Uppi_{-}\downarrow\rangle&|\Upsigma\downarrow\rangle&|\Uppi_{+}\downarrow\rangle&|\Uppi_{-}\uparrow\rangle&|\Upsigma\uparrow\rangle&|\Uppi_{+}\uparrow\rangle\cr&\frac{1}{2}&0&0&0&0&0\cr&0&0&0&\frac{1}{\sqrt{2}}&0&0\cr&0&0&-\frac{1}{2}&0&\frac{1}{\sqrt{2}}&0\cr&0&\frac{1}{\sqrt{2}}&0&-\frac{1}{2}&0&0\cr&0&0&\frac{1}{\sqrt{2}}&0&0&0\cr&0&0&0&0&0&\frac{1}{2}\cr}$
(17)

##  Appendix C Vibrational modes

In order to study the vibronic interactions hergenhahn2004vibrational  and the
Jahn-Teller effect STURGE1968a ; bersuker2006jahn ; bersuker2013jahn  we use
several different mode coordinates. We will detail here all coordinates we
used:

  * •

cartesian $\bm{R}=\\{R_{j}\\}$.

  * •

mass-weighted $\bm{Q}^{\rm m}=\\{Q^{\rm
m}_{j}\\}={\bm{m}}^{1/2}(\bm{R}-\bm{R}_{0})$

  * •

normal modes $\bm{Q}^{\rm n}=\\{Q^{\rm n}_{j}\\}={\bm{O}_{\rm m}^{\rm
n}}^{\dagger}\bm{Q}^{\rm m}$ with oscillations $\bm{\omega}_{\rm n}$ simply
noted $\bm{\omega}$.

  * •

interaction modes $\bm{Q}^{\rm int}=\\{Q^{\rm
int}_{j}\\}=\bm{V}_{r}^{\dagger}{\bm{Q}^{\rm n}}$

  * •

crystal field symmetry adapted interaction modes
$\bm{Q}^{\Upgamma}=\\{Q_{\Upgamma\gamma}\\}={\bm{O}^{\Upgamma}_{\rm
n}}^{\dagger}{\bm{Q}^{\rm n}}$ with oscillations $\bm{\omega}^{\Upgamma}$.

  * •

frequency scaled coordinates ${\tilde{\bm{Q}}}^{\rm
n}=\frac{1}{\sqrt{2}}\bm{\omega}\bm{Q}^{\rm n}$ or
${\tilde{\bm{Q}}}^{\Upgamma}=\frac{1}{\sqrt{2}}\bm{\omega}^{\Upgamma}\bm{Q}^{\rm\Upgamma}$
with unit oscillations frequencies.

and their corresponding momenta conjugate respectively
$\hat{\bm{P}},{\hat{\bm{P}}}_{Q^{\rm m}},{\hat{\bm{P}}}_{Q^{\rm
n}},{\hat{\bm{P}}}_{Q^{\rm
int}},{\hat{\bm{P}}}_{Q^{\Upgamma}},{\hat{\bm{P}}}_{\tilde{Q^{\rm n}}}$,
${\hat{\bm{P}}}_{\tilde{Q^{\Upgamma}}}$ .

###  C.1 Quantization in the normal modes

We consider the matrix elements $V(\bm{R})=\langle
n^{\prime}l^{\prime}m^{\prime}|\hat{V}|nlm\rangle$ for the Cs(6s) and 6p
levels.

Neglecting the anharmonic terms, the electron-lattice crystal field
$V=V_{0}+{\bm{\nabla}}_{\\!\\!\bm{R}}V^{\dagger}(\bm{R}-\bm{R}_{0})+\dfrac{1}{2}(\bm{R}-\bm{R}_{0})^{\dagger}{\bm{\nabla}}_{\\!\\!\bm{R}}^{2}V(\bm{R}-\bm{R}_{0})$
is numerically evaluated by finite difference:
${\bm{\nabla}}_{\\!\\!\bm{R}}V_{k}=\frac{\partial V}{\partial
R_{k}}(\bm{R}^{0})\approx\frac{V(\bm{R}^{0}+\delta R_{k})-V(\bm{R}^{0}-\delta
R_{k})}{2\delta R_{k}}$ and
${\bm{\nabla}}_{\\!\\!\bm{R}}^{2}V_{kl}=\frac{\partial^{2}V}{\partial
R_{k}\partial R_{l}}(\bm{R}^{0})\approx\frac{\frac{\partial V}{\partial
R_{l}}(\bm{R}^{0}+\delta R_{k})-\frac{\partial V}{\partial
R_{l}}(\bm{R}^{0}-\delta R_{k})}{2\delta R_{k}}$ Where we have used standard
row vector
$\bm{R}^{\dagger}=(R_{1},R_{2},\cdots),{\bm{\nabla}}_{\\!\\!\bm{R}}V^{\dagger}=(\frac{\partial
V}{\partial R_{1}},\frac{\partial V}{\partial R_{2}},\cdots)$ († being the
Hermitian adjoint, that is the conjugate transpose) and matrix
${\bm{\nabla}}_{\\!\\!\bm{R}}^{2}V$ notations.

In order to treat all kinetic energy terms $\frac{{\hat{P}}_{k}^{2}}{2m_{k}}$,
(${\hat{P}}_{k}=-i\hbar\frac{\partial}{\partial R_{k}}$ being the quantized
momentum conjugate coordinate) in a similar manner, we use the $3(N-n+1)$
mass-weighted Cartesian coordinates $Q^{\rm
m}_{j}=(R_{j}-{R_{0}}_{j})\sqrt{m_{j}}$ where $m_{j}$ is the mass of the atom
on which the $j^{\rm th}$ coordinate resides. We define the diagonal mass
matrix $\bm{m}$ with $m_{jj}=m_{j}$ and ${\bm{Q}^{\rm m}}^{\dagger}=(Q^{\rm
m}_{1},Q^{\rm m}_{2},\cdots),{{\bm{\nabla}}_{\\!\\!{\bm{Q}}^{\rm
m}}}^{\dagger}=(\frac{\partial V}{\partial Q^{\rm m}_{1}},\frac{\partial
V}{\partial Q^{\rm m}_{2}},\cdots)$.

So the nuclear Hamiltonian in the Born-Oppenheimer approximation
$\hat{H}=\dfrac{1}{2}{\hat{\bm{P}}}^{\dagger}{\bm{m}}^{-1}{\hat{\bm{P}}}+V_{0}+{\bm{\nabla}}_{\\!\\!\bm{R}}V^{\dagger}(\bm{R}-\bm{R}_{0})+\dfrac{1}{2}(\bm{R}-\bm{R}_{0})^{\dagger}{\bm{\nabla}}_{\\!\\!\bm{R}}^{2}V(\bm{R}-\bm{R}_{0})$
becomes simply $\hat{H}=\dfrac{1}{2}{\hat{\bm{P}}_{\rm Q^{\rm
m}}}^{\dagger}{\hat{\bm{P}}_{\rm Q^{\rm
m}}}+V_{0}+{\bm{\nabla}}_{\\!\\!\bm{Q}^{\rm m}}V^{\dagger}\bm{Q}^{\rm
m}+\dfrac{1}{2}{\bm{Q}^{\rm m}}^{\dagger}{\bm{\nabla}}_{\\!\\!\bm{Q}^{\rm
m}}^{2}V\bm{Q}^{\rm m}$, where we have used ${\hat{\bm{P}}_{\rm Q^{\rm
m}}}={\bm{m}}^{-1/2}\hat{\bm{P}}$ vector notation for the momentum conjugate
with the $\bm{Q}^{\rm m}$ coordinates, that verifies the canonical commutation
relation $[\bm{Q}^{\rm m},{\hat{\bm{P}}_{\rm Q^{\rm
m}}}]=[\bm{R},\hat{\bm{P}}]=\bm{R}{\hat{\bm{P}}}^{\dagger}-{\hat{\bm{P}}}{\bm{R}}^{\dagger}=i\hbar\bm{1}$.

The eigendecomposition of the Hessian matrix ${\bm{\nabla}}_{\\!\\!\bm{Q}^{\rm
m}}^{2}V=\bm{O}_{\rm m}^{\rm n}{\bm{\omega}}^{2}{\bm{O}_{\rm m}^{\rm
n}}^{\dagger}$ creates the the normal vibrational modes $\bm{Q}^{\rm
n}={\bm{O}_{\rm m}^{\rm n}}^{\dagger}\bm{Q}^{\rm m}$ with ${\bm{O}_{\rm
m}^{\rm n}}$ the change-of-basis orthonormal matrix and ${\bm{\omega}}$ a
diagonal matrix with the k-th element being the (mass scaled) oscillation
angular frequency $\omega_{k}$. The Hamiltonian is thus written as

$\displaystyle\hat{H}$ $\displaystyle=$
$\displaystyle\dfrac{1}{2}{\hat{\bm{P}}_{\rm Q^{\rm
n}}}^{\dagger}{\hat{\bm{P}}_{\rm Q^{\rm
n}}}+V_{0}+{\bm{\nabla}}_{\\!\\!\bm{Q}^{\rm n}}V^{\dagger}\bm{Q}^{\rm
n}+\dfrac{1}{2}{\bm{Q}^{\rm n}}^{\dagger}\bm{\omega}^{2}\bm{Q}^{\rm n}$
$\displaystyle=$
$\displaystyle\sum_{k}-\hbar^{2}\frac{\partial^{2}}{\partial^{2}Q^{\rm
n}_{k}}+V_{0}+({\bm{\nabla}}_{\\!\\!\bm{Q}^{\rm n}}V^{*})_{k}Q^{\rm
n}_{k}+\frac{1}{2}\omega_{k}^{2}{Q^{\rm n}_{k}}^{2}$

with ${\bm{\nabla}}_{\\!\\!\bm{Q}^{\rm n}}V={\bm{O}_{\rm m}^{\rm
n}}^{\dagger}{\bm{\nabla}}_{\\!\\!\bm{Q}^{\rm m}}V$ and ${\hat{\bm{P}_{\rm
Q^{\rm n}}}}={\bm{O}_{\rm m}^{\rm n}}^{\dagger}{\hat{\bm{P}}_{\rm Q^{\rm m}}}$
is the conjugate momentum with $[\bm{Q}^{\rm n},{\hat{\bm{P}}_{\rm Q^{\rm
n}}}]=i\hbar\bm{1}$.

###  C.2 First order: crystal field coefficient

We will now use the decomposition (C.1) for $V(\bm{R})=\langle
6pm^{\prime}|\hat{V}|6pm\rangle$.

####  C.2.1 Crystal field coefficients

In order to determine the (real) crystal fields coefficients $V_{\Gamma}$,
that are the
$({\bm{\nabla}}_{\\!\\!\bm{Q}^{\Upgamma}}V^{\dagger})_{m^{\prime}m}=\langle
m^{\prime}|\frac{\partial V}{\partial Q_{\Gamma\gamma}}|m\rangle$ coefficients
present in the 3$\times$3 crystal field matrix $M_{\rm CF}$, we restrict
ourselves to only the first order series in nuclear coordinates. Thus
$V=\langle 6pm^{\prime}|\hat{V}|6pm\rangle$ contains only
$(V_{0})_{m^{\prime}m}+({\bm{\nabla}}_{\\!\\!\bm{Q}^{\rm
n}}V^{\dagger})_{m^{\prime}m}\bm{Q}^{\rm n}$ type of terms. We will put the
$({\bm{\nabla}}_{\\!\\!\bm{Q}^{\rm n}}V^{\dagger})_{m^{\prime}m}$ values in a
3$\times$3=9 rows matrix, with each row linked to a given $m^{\prime},m$.

The first order potential can be calculated using the normal mode basis
$\bm{Q}^{\rm n}$ or expressed in the expected $M_{\rm CF}$ crystal field ones.
Equalizing the two expressions leads to the equation
${\bm{\nabla}}_{\\!\\!\bm{Q}^{\Upgamma}}V^{\dagger}\bm{Q}^{\Upgamma}={\bm{\nabla}}_{\\!\\!\bm{Q}^{\rm
n}}V^{\dagger}\bm{Q}^{\rm n}$.

Without knowing the orthonormal transformation,
$\bm{Q}^{\Upgamma}={\bm{O}^{\Upgamma}_{\rm n}}^{\dagger}{\bm{Q}^{\rm n}}$ that
links the crystal field symmetry adapted interaction mode coordinates
$Q_{\Gamma\gamma}$ to the normal mode coordinates $Q^{\rm n}_{k}$, the use of
the fact that the transformation is orthonormal leads to a 9$\times$9 matrix
equation

${\bm{\nabla}}_{\\!\\!\bm{Q}^{\Upgamma}}V^{\dagger}{\bm{\nabla}}_{\\!\\!\bm{Q}^{\Upgamma}}V={\bm{\nabla}}_{\\!\\!\bm{Q}^{\rm
n}}V^{\dagger}{\bm{\nabla}}_{\\!\\!\bm{Q}^{\rm n}}V$ (19)

that, when solved, gives the crystal field coefficients
${\bm{\nabla}}_{\\!\\!\bm{Q}_{\Gamma}}V$.

####  C.2.2 Interaction mode coordinates

In order to perform Monte Carlo simulations, it might be important to know how
the crystal field symmetry adapted interaction mode coordinates
$\bm{Q}^{\Upgamma}$ can be calculated.

They can be found using the singular value decomposition of the 9 rows matrix
${\bm{\nabla}}_{\\!\\!\bm{Q}^{\rm n}}V^{\dagger}=\bm{U}\bm{W}\bm{V}^{\dagger}$
of rank $r$ (the matrix $\bm{V}$ not being confused with the potential $V$).
Equation
${\bm{\nabla}}_{\\!\\!\bm{Q}^{\Upgamma}}V^{\dagger}\bm{Q}^{\Upgamma}={\bm{\nabla}}_{\\!\\!\bm{Q}^{\rm
n}}V^{\dagger}\bm{Q}^{\rm n}$ leads to
${\bm{\nabla}}_{\\!\\!\bm{Q}_{\Gamma}}V^{\dagger}\bm{Q}^{\Upgamma}=(\bm{U}\bm{W})(\bm{V}^{\dagger}{\bm{Q}^{\rm
n}})$. We cannot simply equalize
${\bm{\nabla}}_{\\!\\!\bm{Q}_{\Gamma}}V^{\dagger}$ with $\bm{U}\bm{W}$ and
$\bm{Q}^{\Upgamma}$ with $\bm{V}^{\dagger}{\bm{Q}^{\rm n}}$ because
$\bm{Q}^{\rm int}=(\bm{V}^{\dagger}{\bm{Q}^{\rm n}})$ is only one choice for
the interaction mode coordinates, that is not necessarily the same basis (for
instance not in the same order, same sign, …) as the
$\bm{Q}^{\Upgamma}=\\{Q_{\Gamma\gamma}\\}$ used to determine from symmetry
consideration the crystal field matrix $M_{\rm CF}$. However, using the only
useful first $r$ rows of the $r$ blocks matrix restriction (noted with index
r) of
$(\bm{W}_{r})^{-1}\bm{U}_{r}^{\dagger}\left({\bm{\nabla}}_{\\!\\!\bm{Q}^{\Upgamma}}V^{\dagger}\cdot\bm{Q}^{\Upgamma}\right)=\bm{Q}^{\rm
int}=(\bm{V}_{r}^{\dagger}{\bm{Q}^{\rm n}})$ leads to the transformation
$\bm{Q}^{\Upgamma}={\bm{O}^{\Upgamma}_{\rm n}}^{\dagger}{\bm{Q}^{\rm n}}$,
where, for simplicity, we have kept the same notation $\bm{Q}^{\Upgamma}$ or
$\bm{Q}^{\rm int}$ for the first $r$ interaction modes coordinates than for
the $3N$ ones.

The important fact is that the diagonal matrix $\bm{W}$ has only $r$ non-zero
singular values on the diagonal (so $\bm{W}_{r}$ is square diagonal of size
$r$), and thus we will often be interested only in the $r$ relevant ($r$
first) interaction modes coordinates.

With these new coordinates it turns out that the $r\times r$ matrix
$(\bm{\omega}_{\Upgamma})^{2}={\bm{O}^{\Upgamma}_{\rm
n}}^{\dagger}\bm{\omega}^{2}{\bm{O}^{\Upgamma}_{\rm n}}$ matrix is diagonal.
We thus have

$\langle
6s|\hat{V}|6s\rangle=V_{0}+\dfrac{1}{2}{\bm{Q}^{\Upgamma}}^{\dagger}(\bm{\omega}^{\Upgamma})^{2}\bm{Q}^{\Upgamma}=V_{0}+\sum_{\Gamma\gamma=1}^{r}\frac{1}{2}\omega_{\Gamma}^{2}Q_{\Gamma\gamma}^{2}$
(20)

Once again, we see the advantage of using the interaction mode coordinates to
reduce drastically the number of modes to be calculated. The Monte Carlo
simulation becomes very simple because we have to calculate only few ($r\leq
9$) interaction mode coordinates, such as
$\bm{Q}^{\Upgamma}=\\{Q_{A_{1}},Q_{E,1},Q_{E,2},Q_{T,1},Q_{T,2},Q_{T,3}\\}$ to
be compared to $\sim 3N$ modes in standard molecular dynamic simulation
(almost 10000 in our case for the $\sim 3000$ movable atoms).

###  C.3 Frequency scaled coordinates

Following toyozawa1966dynamical ; cho1968optical , it might be interesting to
use the oscillation frequency scaled coordinates
${\tilde{\bm{Q}}}^{n}=\frac{1}{\sqrt{2}}\bm{\omega}\bm{Q}^{\rm n}$ and
${\tilde{\bm{Q}}}^{\Upgamma}=\frac{1}{\sqrt{2}}\bm{\omega}^{\Upgamma}\bm{Q}^{\rm\Upgamma}$
such that Eq. (20) become

$\displaystyle\langle 6s|\hat{V}|6s\rangle$ $\displaystyle=$ $\displaystyle
V_{0}+{{\tilde{\bm{Q}}}^{\rm n}}{}^{\dagger}{\tilde{\bm{Q}}}^{\rm
n}=V_{0}+\sum_{k}{{\tilde{Q}}_{k}^{\rm n}}{}^{2}$ $\displaystyle=$
$\displaystyle
V_{0}+\sum_{\Gamma\gamma=1}^{r}{{\tilde{Q}}_{\Gamma\gamma}}^{2}$

This form is interesting, because in a pure classical molecular dynamical
simulation, the Boltzmann statistics indicates that ${\tilde{\bm{Q}}}^{\rm n}$
follows a multivariate normal distribution $\propto
e^{-\frac{1}{2}{\tilde{\bm{Q}}}^{\rm
n}{}^{\dagger}{\tilde{\bm{\Sigma}}}^{-1}{\tilde{\bm{Q}}}^{\rm n}}$ with a
variance ${\tilde{\bm{\Sigma}}}^{\rm n}=\frac{1}{2}k_{B}T\bm{I}$. That way the
interaction mode coordinate vector
${\tilde{\bm{Q}}}^{\Upgamma}={\bm{O}^{\Upgamma}_{\rm
n}}^{\dagger}{\tilde{\bm{Q}}}^{n}$ follows a multivariate normal distribution
with a variance that is also diagonal:
${\tilde{\bm{\Sigma}}}^{\Upgamma}={\bm{O}^{\Upgamma}_{\rm
n}}^{\dagger}{\tilde{\bm{\Sigma}}}^{\rm n}{\bm{O}^{\Upgamma}_{\rm
n}}=\frac{1}{2}k_{B}T\bm{I}$. Thus, in a standard classical Monte Carlo
simulation at temperature $T$, each
$\frac{{{\tilde{Q}}_{\Gamma\gamma}}}{\sqrt{k_{B}T/2}}$ coordinate has the same
distribution: a standard unit normal (Gaussian).

It is even possible to keep the normalized coordinates in a classical
simulation that includes the temperature scales $T^{\prime}$ by having
${\tilde{\bm{Q}}}^{\rm n}$ following a multivariate normal distribution with a
new diagonal variance ${\tilde{\bm{\Sigma}}}{}^{\rm n}{}^{\prime}$ with
diagonal element given by $\frac{1}{2}k_{B}T^{\prime}(\omega_{k})$ according
to Eq. (6), namely $\frac{{\tilde{Q}}^{\rm
n}_{k}}{\sqrt{k_{B}T^{\prime}(\omega_{k})/2}}$. In this case,
${\tilde{\bm{Q}}}^{\Upgamma}$ has simply to be chosen with a Gaussian variance
${\tilde{\bm{\Sigma}}}^{\Upgamma}{}^{\prime}={\bm{O}^{\Upgamma}_{\rm
n}}^{\dagger}{\tilde{\bm{\Sigma}}}^{\rm n}{}^{\prime}{\bm{O}^{\Upgamma}_{\rm
n}}$. Another strategy, that leads to the same final result, is to keep the
unit normal distribution for $\frac{\tilde{Q}_{k}^{\rm n}}{\sqrt{k_{B}T/2}}$
but change each $\omega_{k}$ into $\omega_{k}^{\prime}$ such that
$\omega_{k}^{\prime 2}/T=\omega_{k}^{2}/T^{\prime}(\omega_{k})$.

So, as used in Ref. toyozawa1966dynamical ; cho1968optical , such frequency
scaled coordinates are useful in classical simulations because the
oscillations frequencies $\omega^{\Upgamma}$ (or the scaled ones
$\omega^{\prime\Upgamma}$) of the modes are not needed, if we include them in
the definition of new crystal field parameters
$\tilde{V}_{\Upgamma}=\sqrt{2}V_{\Upgamma}/\omega_{\Gamma}$ . Then all terms
in the interaction matrix can be written using
$V_{\Upgamma}Q_{\Upgamma\gamma}=\tilde{V}_{\Upgamma}\tilde{Q}_{\Gamma\gamma}$
with now Eq. (10) becoming Eq. (11) that is

$A(E)\propto\int\sum_{i=1}^{6}\delta[E-X_{i}({\tilde{\bm{Q}}}^{\Upgamma})]e^{-\sum_{\Gamma\gamma}{\tilde{Q}}_{\Gamma\gamma}^{2}/k_{B}T}{\rm
d}{\tilde{\bm{Q}}}^{\Upgamma}$ (21)

with the $\frac{{{\tilde{Q}}_{\Gamma\gamma}}}{\sqrt{k_{B}T/2}}$ coordinates
having a standard unit normal (Gaussian) distribution.

To illustrate this, we mention that the values given in Ref. Kanda1971  for
the crystal field in Table 2 are

$A=\frac{\sqrt{k_{B}T/2}}{A^{SO}}\tilde{V}_{A_{1}}=\frac{\sqrt{k_{B}T/2}}{A^{SO}}\frac{\sqrt{2}}{\omega_{A_{1}}}V_{A_{1}}=0.294$,

$B=\frac{\sqrt{k_{B}T/3}}{A^{SO}}\tilde{V}_{E}=\frac{\sqrt{k_{B}T/3}}{A^{SO}}\frac{\sqrt{2}}{\omega_{E}}V_{E}=0.6$,

$C=\frac{\sqrt{k_{B}T/4}}{A^{SO}}\tilde{V}_{T_{2}}=\frac{\sqrt{k_{B}T/4}}{A^{SO}}\frac{\sqrt{2}}{\omega_{T_{2}}}V_{T_{2}}=0.2$.

Indeed the linear part of the potential is

$V_{\Gamma\gamma}Q_{\Gamma\gamma}=A^{SO}\left(\frac{\sqrt{k_{B}T/x_{\Gamma}}}{A^{SO}}\frac{\sqrt{2}}{\omega_{\Gamma}}V_{\Gamma_{\gamma}}\right)\left(\frac{1}{\sqrt{k_{B}T/x_{\Gamma}}}\frac{\omega_{\Gamma}}{\sqrt{2}}Q_{\Gamma_{\gamma}}\right)$

($x_{\Gamma}$ is 2 for the case of $\Gamma=A_{1}$, 3 for $E_{2}$ and 4 for
$T_{2}$cho1968optical ).

In the above notation, the ground state potential reduces to
$\sum_{\Gamma\gamma}\frac{1}{2}\omega_{\Gamma}^{2}Q_{\Gamma\gamma}^{2}=\sum_{\Gamma\gamma}{{\tilde{Q}}_{\Gamma\gamma}}^{2}$
.

For the quantized simulation, this frequency scaled coordinates procedure does
not work anymore because the quantized energy levels are intrinsically linked
to the frequencies, that thus cannot be simply scaled. We could think of using
the frequency scaled coordinates with unity frequencies such that all modes
look similar and the number of different modes would appear to matter less.
However this approach is futile because the kinetic energy of the lattice
(from Eq. (C.1)) would not be diagonal in the momenta conjugate to the scaled
interaction mode anymore, as already noticed in Ref. toyozawa1966dynamical .
This is not important in the classical Franck-Condon Mulliken approximation
where the kinetic energy plays no role because there, we only need to find the
positions that are determined by the Boltzmann distribution involving only the
potential energy. It, however, becomes crucial for quantized energies, such as
in the reflection approximation.

###  C.4 Coordinate distribution

Let us mention that at $T=\qty{0}{K}$ each normal mode
$|\Psi_{i}(\bm{Q}_{k})|^{2}$ distribution is gaussian (ground state of the
harmonic oscillator) and so the multivariate normal distribution $\propto
e^{-\frac{1}{2}{{\bm{Q}}}^{\rm
n}{}^{\dagger}{{\bm{\Sigma}}}^{-1}{{\bm{Q}}}^{\rm n}}$ of the $\bm{Q}^{\rm n}$
normal modes is transformed in a multivariate normal distribution with a
variance ${{\bm{\Sigma}}}^{\Upgamma}={\bm{O}^{\Upgamma}_{\rm
n}}^{\dagger}{{\bm{\Sigma}}}^{\rm n}{\bm{O}^{\Upgamma}_{\rm n}}$ for the
interaction mode. It turns out that the marginal distribution for the 6 mode
coordinates is indeed a variance that is diagonal. This is a consequence of
the fact that the projection on the 6 mode coordinates on the ground state
Hessian matrix is diagonal with the frequencies $\omega^{\Upgamma}$ given in
Table 2. Together with the fact that the first order approximation is quite
accurate for the excited state interaction this leads to the fact that the
$T=\qty{0}{K}$ spectra can be calculated in a pure crystal field model with
only 6 coordinates simply using frequencies $\omega^{\Upgamma}$ to determine
the gaussian $Q^{\Upgamma}$ coordinates using temperature scaling or the
ground state wavefunction of the
$\frac{1}{2}\omega_{\Gamma\gamma}^{2}Q_{\Gamma\gamma}^{2}$ harmonic
oscillator. The result for $T=\qty{0}{K}$ gives the same spectra as the one in
Fig. 9. Indeed at $T=\qty{0}{K}$ each $E_{i}$ is constant and is the zero
point energy of the $i$-th modes. By offsetting to this value, we have at
$T=\qty{0}{K}$, $E_{i}=0$ and so the reflection approximation is
$A(E)\propto\sum_{i=1}^{6}P_{i}\int\left|\Psi_{i}(\bm{Q})\right|^{2}\delta[E-(V_{\rm
Cs(6p)-Ar}(\bm{Q})-E_{i})]{\rm d}\bm{Q}$ which can be compared to the formula
(12).

##  Appendix D Semi-classical approximations

To study laser excitation of the $(6s)$ cesium atom toward the $6p$ manifold
we start (thanks to the Beer-Lambert-Bouguer’s law) by using the fact that the
spectral density optical absorption coefficient $A(E)$ for a photon of energy
$E$, is given by the sum over all initial vibronic levels $|i\rangle$,
populated with probability $P_{i}$ with wavefunction $\Psi_{i}(\bm{Q})$ of
energy $E_{i}$, towards all possible final ones $\Psi_{f}(\bm{Q})$ of energy
$E_{f}$:

$A(E)\propto\sum_{if}P_{i}\delta[E-(E_{f}-E_{i})]\left|\int\Psi_{f}^{\dagger}(\bm{Q})d_{if}(\bm{Q})\Psi_{i}(\bm{Q}){\rm
d}\bm{Q}\right|^{2}$ (22)

Thus $A(E)dE$ is the absorption coefficient for a photon in the energy band
$E,E+dE$.

###  D.1 Semi-classical Franck-Condon approximation

In our case, the variation of the dipole transition strength from X to
$\Upsigma$ or $\Uppi$ states is on the order of 10% blank2012m+ . So we can
reasonably consider variation of the dipole $d_{if}(\bm{Q})$ with the
internuclear distances $\bm{Q}$ as negligible (this is the ”Condon”
approximation lax1952franck ). If needed, a better approximation would be the
Herzberg-Teller centroid one with $d_{if}$ taken at the point $\bar{\bm{Q}}$,
which cancels the first order evolution in $\bm{Q}$. In the following, for
simplicity we will assume $d_{if}(\bm{Q})\approx d_{if}(\bar{\bm{Q}})$ to be
constant (it is simply the 6s to 6p dipole transition and can thus be left out
of the equation). Furthermore, we can also justify an isotropic nature of the
dipole in our case because of the poly-crystalline structure of our sample
leads to random orientations of the crystal axes.

In our case, after laser excitation of the Cs(6s) atom from its equilibrium
position, the Cs(6p) atom is far from being in its equilibrium position and
the excited electron potentials are mostly linear and not quadratic. Thus the
motional excited states will be states with large vibrational quantum numbers
or even quasi continuum one, which, according to the Bohr correspondence
principle, can be treated as approximately classical. Thus, following
lax1952franck  we can assume transitions near the classical turning points and
for the sum over $f$ we can replace $E_{f}$ in the delta-function by a mean
value independent of $f$ that is the electron excited (6p) potential energy
curve $V_{\rm e}(\bm{Q})$ so:

$A(E)\propto\sum_{i}P_{i}\int\delta[E-(V_{\rm
e}(\bm{Q})-E_{i})]\left|\Psi_{i}(\bm{Q})\right|^{2}{\rm d}\bm{Q}$ (23)

This expression is also called the reflection approximation heller1978quantum
.

If many initial states (vibrational for instance) $i$ are involved, it is
convenient to use a mean-value approximation using for the ground state
energies, the potential energy $V_{\rm g}(\bm{Q})$ leading to the so-called
semiclassical Franck-Condon formula

$A(E)\propto\int\delta[E-(V_{\rm e}(\bm{Q})-V_{\rm g}(\bm{Q}))]P_{\rm
g}(\bm{Q}){\rm d}\bm{Q}$ (24)

with $P_{\rm g}(\bm{Q})=\sum_{i}P_{i}\left|\Psi_{i}(\bm{Q})\right|^{2}$ is the
quantum-statistical mechanical probability distribution. In the case where
$P_{\rm g}(\bm{Q})$ is given by classical statistics, the expression becomes
the standard classical Franck-Condon formula lax1952franck .

###  D.2 Semiclassical transition in energy or in phase space

It is useful to present the spectral density transition probability
$P(E_{0})$, for a transition energy $E_{0}$ given by Eq. (22), in a state
picture

$P(E_{0})=\sum_{f,i}P_{i}\left|\langle\Psi_{f}|\hat{d}|\Psi_{i}\rangle\right|^{2}\delta(E_{0}-(E_{f}-E_{i}))$
(25)

because $2\pi P(E_{0})/\hbar$ can now be interpreted as a rate given by
Fermi’s Golden rule siebrand1967radiationless ; freed1975nuclear ;
stockmann1984influence .

This can be written simply as (in the Condon or centroid approximation where
we forget about $d^{2}$ from now on), for a pure state ($P_{i}=1$ to simplify
the notation) as

$\displaystyle P(E_{0})$ $\displaystyle=$
$\displaystyle\sum_{f,i}\langle\psi_{f}|\psi_{i}\rangle\langle\psi_{i}|\psi_{f}\rangle\langle\psi_{f}|\delta(E_{0}-(\hat{H}_{\rm
f}-E_{i}))|\psi_{f}\rangle$ $\displaystyle=$
$\displaystyle\sum_{f,i}\langle\psi_{f}|\psi_{i}\rangle\langle\psi_{i}|\sum_{f^{\prime}}|\psi_{f^{\prime}}\rangle\langle\psi_{f^{\prime}}|\delta(E_{0}-(\hat{H}_{\rm
f}-E_{i}))|\psi_{f}\rangle$ $\displaystyle=$
$\displaystyle\sum_{f}\langle\psi_{f}|\Big{(}\sum_{i}|\psi_{i}\rangle\langle\psi_{i}|\delta(E_{0}-(\hat{H}_{\rm
f}-E_{i}))\Big{)}|\psi_{f}\rangle$ $\displaystyle=$ $\displaystyle
Tr\Big{[}\sum_{i}|\psi_{i}\rangle\langle\psi_{i}|\delta(E_{0}-(\hat{H}_{\rm
f}-E_{i}))\Big{]}$ $\displaystyle=$ $\displaystyle Tr[\hat{\rho}_{\rm
i}\delta(E_{0}-(\hat{H}_{\rm f}-E_{i}))]$ $\displaystyle=$
$\displaystyle\langle\delta(E_{0}-(\hat{H}_{\rm f}-E_{i}))\rangle$
$\displaystyle=$ $\displaystyle\langle\delta(E_{0}-(\hat{V}_{\rm
f}-\hat{V}_{\rm i})-(\hat{H}_{\rm i}-E_{i}))\rangle$

where $\hat{H}_{\rm i}=\sum_{k}{\hat{P}}_{k}^{2}/(2m_{k})+V_{\rm
i}(\hat{\bm{Q}})$ and $\hat{H}_{\rm
f}=\sum_{k}{\hat{P}}_{k}^{2}/(2m_{k})+V_{\rm f}(\hat{\bm{Q}})$ are the
Hamilton operators respectively for the initial and final electronic state and
$\hat{\rho}_{\rm i}=\sum_{i}|\Psi_{i}\rangle\langle\Psi_{i}|$ is the density
matrix operators for the initial state. The formula can also be written, in
the case of a pure state as $P(E_{0})=\frac{1}{2\pi\hbar}\int
dte^{i(E_{0}-\hat{E}_{\rm i})t/\hbar}\langle\Psi_{i}|\Psi_{i}(t)\rangle$where
$|\Psi_{i}(t)\rangle=e^{-i\hat{H}_{\rm f}t/\hbar}|\Psi_{i}\rangle$. This time
(and Fourier transform) picture opens another way to treat the transition
using the (wave packet) time evolution of the system to find $|\psi(t)\rangle$
heller1976wigner ; mukamel1982semiclassical ; shao1998quantum ;
hupper1998uniform ; ansari2022instanton .

One of the best approximations for the formula is clearly to use $\hat{H}_{\rm
i}-E_{i}$ as a perturbation, especially when projected on $|\psi_{i}\rangle$
because $\langle\psi_{i}|\hat{H}_{\rm i}-E_{i}|\psi_{i}\rangle=0$. Therefore
the simplest natural approximation is the zero-order one for a $d$ dimensional
space:

$\displaystyle P(E_{0})$ $\displaystyle=$
$\displaystyle\langle\delta(E_{0}-(\hat{V}_{\rm f}-\hat{V}_{\rm i}))\rangle$
$\displaystyle=$ $\displaystyle(2\pi\hbar)^{-d}\sum_{i}\int
d\bm{Q}|\psi_{i}(\bm{Q})|^{2}\delta((E-{V_{\rm f}}(\bm{Q})+{V_{\rm
i}}(\bm{Q})))$

This is exactly the same as the formula (24). But, note that this derivation
differs considerably from the previous traditional derivation despite the
final results being the same. For instance, it gives a very clear physical
understanding that the reflection approximation of Eq. (5) is reproduced under
the extra approximation of $\delta(E-({H_{\rm f}}_{\rm cl}-{H_{\rm i}}_{\rm
cl})$ replaced by $\delta(E-(V_{\rm f}(\bm{Q})-E_{i}))$.

A very similar derivation has been done by moving to the Wigner representation
of quantum mechanics 1932Wigner ; heller1976wigner ; tatarskiui1983wigner ;
hillery1984distribution ; case2008wigner ; ferry2018wigner  where the exact
formula for $P(E_{0})$ becomes, in $d$ dimension, simply
$P(E_{0})=(2\pi\hbar)^{-d}\int d\bm{P}d\bm{Q}{\rho_{\rm
i}}_{W}\delta(E_{0}-(\hat{H}_{\rm f}-E_{i}))_{W}$ where W designing the Wigner
transform. The leading order in $\hbar$ is simply the semi-classical limit
with $\hat{H}_{\rm f,i}$ replaced by their classical counterpart $H_{\rm
f,i}=\sum_{k}{\hat{P}}_{k}/(2m_{k})+V_{\rm f,i}(\bm{Q})$ and ${\rho_{\rm
i}}_{W}$ replaced by the classical phase space density $\rho_{\rm cl}\propto
e^{-H_{\rm cl}/k_{B}T}$. It is obviously slightly better to keep the true
Wigner function ${\rho_{\rm i}}_{W}$ (also because it is known analytically
for the harmonic potential case) to have segev2003fermi :
$P(E)=(2\pi\hbar)^{-d}\int d\bm{P}d\bm{Q}{\rho_{\rm i}}_{W}\delta(E-({H_{\rm
f}}_{\rm cl}-{H_{\rm i}}_{\rm cl}))$. This is again the same formula (24) if
using the fact that integration of the Wigner function over $\bm{P}$ gives
exactly the wavefunction probability distribution $|\psi_{i}(\bm{Q})|^{2}$
tatarskiui1983wigner . This derivation has the advantage of justifying the use
of the true quantum (or Wigner) phase-space initial distribution and not the
classical one. Thus the probability $P_{i}$ has to be chosen accordingly to
the actual quantum distribution.

Another advantage of this semi-classical phase space picture is that it allows
for systematic series expansion, typically in power of $\hbar$ order
corrections 1932Wigner ; heller1976wigner ; heller1978quantum ;
shao1998quantum ; japha2002semiclassical ; segev2003fermi ;
sergeev2003semiclassical ; pollak2023thermal . Finally, the method can be
generalized for non-radiative transitions, such as surface hoping non-
adiabatic effects or electron transfer with instanton theory or other non
Born–Oppenheimer effects kallush2002surface ; segev2003fermi ;
segev2003dominant ; heller2020semiclassical ; agostini2022chemistry .

##  Appendix E Third order many-body term

We will detail some calculations to determine the effective two-body potential
from the sum of the third-order terms. For this, we will first derive the
general expression for the third-order terms and calculate it for the long-
range dipolar case and in a simple two-level approximation. This will allow us
to create a mean field of long-range third-order terms that we compare to the
long-range part of the two-body case. We finally extend this comparison for
all internuclear distances to produce an effective two-body potential that
includes the third-order effect.

###  E.1 Perturbation theory up to third order