The energy shift $\Delta E$ of the energy $E$ of the full system of a given
atom A (Cs in our case) in a given state $|lm\rangle$ within the crystal
formed by many other atoms (Ar in our case) can be estimated by the
perturbation theory from the full hamiltonian energy
$H=\sum_{i}H_{0}(i)+\frac{1}{2}\sum_{i\neq j}V_{ij}$ where $H_{0}(i)$ is the
single $i^{\rm th}$ atom hamiltonian and
$V_{ij}(\bm{r}_{i},\bm{r}_{j},\bm{R}_{i},\bm{R}_{j})$ contains all
electrostatic interactions between the $i^{\rm th}$ atom and the $j^{\rm th}$
with electron coordinates $\bm{r}_{i},\bm{r}_{j}$ and nuclear ones
$\bm{R}_{i},\bm{R}_{j}$.

We note the eigenstates of the single atom Hamiltonian $H_{0}(i)$ as
$|0_{i}\rangle$ for the initial (not necessarily the ground state in the Cs
case) state of atom $i$ of energy $E_{0_{i}}$ and $|m_{i}\rangle$ the other
states of energy $E_{m_{i}}$.

We calculate the energy shift $\Delta E(\bm{R})$ that only depends on the
nuclear coordinates $\bm{R}$ by perturbation theory. Up to the third-order,
the full crystal energy is
$\frac{1}{2}\sum_{AB}E_{AB}^{(2)}+\frac{1}{6}\sum_{ABC}E_{ABC}^{(3)}=\sum_{A<B}E_{AB}^{(2)}+\sum_{A<B<C}E_{ABC}^{(3)}$.
The first order is zero and the second order leads to a shift depending on
atom A given by $\sum_{B}E_{AB}^{(2)}$ where
$E_{12}^{(2)}=\sum_{m_{1},m_{2}}-\frac{|\langle
0_{1}0_{2}|V_{12}|m_{1}m_{2}\rangle|^{2}}{\Delta_{m_{1}}(1)+\Delta_{m_{2}}(2)}$
with $\Delta_{m_{i}}(i)=E_{m_{i}}-E_{0_{i}}$.

The third-order perturbation theory terms depend on three atoms A, B, and C.
Using the symmetry group $S_{3}$ of permutations of the 3 atoms A=1, B=2, C=3,
and using $V_{ij}=V_{ji}$, $E_{ABC}^{(3)}=E_{123}^{(3)}$ is the sum of six
terms (see bell1970multipolar  (7b)):

$E_{123}^{(3)}=\sum_{\sigma\in S_{3},m_{1},m_{2},m_{3}}\frac{\langle
0_{\sigma(1)}0_{\sigma(2)}|V_{\sigma(1)\sigma(2)}|m_{\sigma(1)}m_{\sigma(2)}\rangle\langle
m_{\sigma(2)}0_{\sigma(3)}|V_{\sigma(2)\sigma(3)}|0_{\sigma(2)}m_{\sigma(3)}\rangle\langle
m_{\sigma(3)}m_{\sigma(1)}|V_{\sigma(3)\sigma(1)}|0_{\sigma(3)}0_{\sigma(1)}\rangle}{(\Delta_{m_{\sigma(1)}}(\sigma(1))+\Delta_{m_{\sigma(2)}}(\sigma(2)))(\Delta_{m_{\sigma(1)}}(\sigma(1))+\Delta_{m_{\sigma(3)}}(\sigma(3)))}$
(27)

###  E.2 dipole dipole dipole long-range interaction

We will calculate these terms in the dipolar long-range interaction case.

We thus assume a dipolar type of interaction (with a fixed quantization axis
named z). For two atoms labeled 1 and 2 with internuclear separation
$\bm{R}_{12}=R_{12}\bm{e}_{12}$ where $\bm{e}_{12}$ is a unit vector, the
dipole-dipole interaction is
$V_{12}=\frac{e^{2}}{4\pi\epsilon_{0}R_{12}^{3}}[\bm{r}(1).\bm{r}(2)-3(\bm{r}(1).\bm{e}_{12})(\bm{r}(2).\bm{e}_{12})]=\frac{e^{2}}{4\pi\epsilon_{0}R_{12}^{3}}\bm{r}(1).(1-3\bm{e}_{12}\bm{e}_{12}).\bm{r}(2)$.

That can be written using the irreducible tensors notation,
$r_{q}=r\sqrt{\frac{4\pi}{3}}Y_{1q}(\bm{r})$, as

$V_{12}=-\frac{e^{2}}{4\pi\epsilon_{0}R^{3}}\sqrt{6}C(\bm{e}_{12}).\\{\bm{r}(1)\otimes\bm{r}(2)\\}_{2}$
where $C_{M}(\bm{e}_{12})=\sqrt{\frac{4\pi}{5}}Y_{2M}(\bm{e}_{12})$ and
$\\{\bm{r}(1)\otimes\bm{r}(2)\\}_{2M}=\sum_{q,q^{\prime}}C_{1q1q^{\prime}}^{2M}r_{q}(1)r_{q^{\prime}}(2)$.
So

$V_{12}=-\frac{\sqrt{6}e^{2}}{4\pi\epsilon_{0}R_{12}^{3}}\sum_{M,q_{1},q_{2}}(-1)^{M}C_{-M}(\bm{e}_{12})C_{1q_{1}1q_{2}}^{2M}r_{q_{1}}(1)r_{q_{2}}(2)$

We will often use atomic units where $\frac{e^{2}}{4\pi\epsilon_{0}}=1$. This
can be simply realized by changing $r(i)$ by
$r(i)\sqrt{\frac{e^{2}}{4\pi\epsilon_{0}}}$.

####  E.2.1 Two level approximation

One of the simplest approximations is to consider a single effective atomic
energy transition, in the pure case, a $S\leftrightarrow P$ transition for
each atom (energy difference $E_{SP}(Cs)$ for Cs and $E_{SP}(Ar)$ for Ar). In
our case, all Ar atoms are in the ground state $S$, but the Cs can be either
in ground $S$ or in excited state $P$. Therefore, $S$ or $P$ degeneracy means
that, in Eq. (27), all energy terms are independent of $m_{i}$. We write
$\Delta_{m_{i}}(\sigma(i))$ as $\Delta(\sigma(i))$ that can be factorized out
of the sum in Eq. (27).

Using Cs as atom 1 and Ar for 2 and 3, we have for $i=2,3$
$|0_{i}\rangle=|l=0m=0\rangle$ and $|m_{i}\rangle=|l=1m=m_{i}\rangle$ states.

We then follow the elegant derivation done in ref. midzuno1956non . We first
look on terms depending on $m_{2}$ so for $i$ label for which $\sigma(i)=2$,
in Eq. (27) we have terms $\langle 0_{2}|V(j2)|m_{2}\rangle$ and $\langle
m_{2}|V(2k)|0_{2}\rangle$. So using the Wigner-Eckart theorem $\langle
l^{\prime}m^{\prime}|r_{2}^{(q)}|lm\rangle=\frac{C_{lm1q}^{l^{\prime}m^{\prime}}}{\sqrt{2l^{\prime}+1}}\langle
l^{\prime}||r_{2}||l\rangle$ the $\sum_{m_{2}}\langle
0_{2}|V(j2)|m_{2}\rangle....\langle m_{2}|V(2k)|0_{2}\rangle$ sum (the $....$
are simply here to indicate that the terms are not necessarily neighbors in
Eq. (27)) contains $\sum_{m_{2}}\langle
0_{2}|r_{q_{2}}(2)|m_{2}\rangle....\langle
m_{2}|r_{q_{2}^{\prime}}(2)|0_{2}\rangle=-\sum_{m_{2}}C_{1m_{2}1q_{2}}^{00}...\frac{C_{001q_{2}^{\prime}}^{1m_{2}}}{\sqrt{3}}\langle
1||r(2)||0\rangle^{2}$, where we have used $\langle 0||r(2)||1\rangle=-\langle
1||r(2)||0\rangle$. Then, $\sum_{m_{2}=-1}^{m_{2}=1}$ leads to $\frac{\langle
1||r(2)||0\rangle^{2}}{3}\delta_{q_{2},-q_{2}^{\prime}}(-1)^{q_{2}}$. Similar
results arise for atoms 3, for $j$ label for which $\sigma(j)=3$.

####  E.2.2 Interaction between ground state atoms

If the Cs atom is in its ground state, the levels are
$|0_{1}\rangle=|l=0\,m=0\rangle$ and $|m_{1}\rangle=|l=1\,m_{1}\rangle$.

The second order leads to

$\displaystyle E_{12}^{(2)}$ $\displaystyle=-\frac{C_{6}}{R_{12}^{6}}$ (28)
$\displaystyle C_{6}$ $\displaystyle=$
$\displaystyle\frac{2}{3}\frac{r_{SP}^{2}(Cs)r_{SP}^{2}(Ar)}{E_{SP}(Cs)+E_{SP}(Ar)}$

with $r_{SP}(Cs)=\langle 1||r(1)||0\rangle$ and $r_{SP}(Ar)=\langle
1||r(2)||0\rangle$

For the third order calculation, the results we just derived, in section
E.2.1, indicate that, in the sum of Eq. (27), one numerator becomes

$\displaystyle E$ $\displaystyle=$
$\displaystyle-\frac{\sqrt{6}}{R_{12}^{3}}\frac{\sqrt{6}}{R_{13}^{3}}\frac{\sqrt{6}}{R_{23}^{3}}\frac{\langle
1||r(1)||0\rangle^{2}}{3}\frac{\langle 1||r(2)||0\rangle^{2}}{3}\frac{\langle
1||r(3)||0\rangle^{2}}{3}\times$
$\displaystyle\sum_{q_{1},q_{2},q_{3},q^{\prime}_{1},q^{\prime}_{2},q^{\prime}_{3}}(-1)^{q_{1}+q_{2}+q_{3}}\delta_{q_{1},-q_{1}^{\prime}}\delta_{q_{2},-q_{2}^{\prime}}\delta_{q_{3},-q_{3}^{\prime}}\times$
$\displaystyle\sum_{M,M^{\prime},M^{\prime\prime}}(-1)^{M+M^{\prime}+M^{\prime\prime}}C_{-M}(\bm{e}_{12})C_{-M^{\prime}}(\bm{e}_{13})C_{-M^{\prime\prime}}(\bm{e}_{23})\times$
$\displaystyle
C_{1q_{1}1q_{2}}^{2M}C_{1q^{\prime}_{1}1q_{3}}^{2M^{\prime}}C_{1q^{\prime}_{2}1q^{\prime}_{3}}^{2M^{\prime\prime}}$

but using (8.4(10) and 8.7(15) from varshalovich1988quantum ) we have
$\sum_{q_{1},q_{2},q_{3}}(-1)^{q_{1}+q_{2}+q_{3}}C_{1q_{1}1q_{2}}^{2M}C_{1-q_{1}1q_{3}}^{2M^{\prime}}C_{1-q_{2}1-q_{3}}^{2M^{\prime\prime}}=-(-1)^{-M^{\prime\prime}}5C_{2M2M^{\prime}}^{2-M^{\prime\prime}}\begin{Bmatrix}1&1&2\\\
2&2&1\end{Bmatrix}=-(-1)^{-M^{\prime\prime}}C_{2M2M^{\prime}}^{2-M^{\prime\prime}}\sqrt{7/12}$.
So

$\displaystyle E$ $\displaystyle=$
$\displaystyle\frac{\sqrt{14}}{9}\frac{\langle
1||r(1)||0\rangle^{2}}{R_{12}^{3}}\frac{\langle
1||r(2)||0\rangle^{2}}{R_{13}^{3}}\frac{\langle
1||r(3)||0\rangle^{2}}{R_{23}^{3}}\times$
$\displaystyle\sum_{M,M^{\prime},M^{\prime\prime}}(-1)^{M+M^{\prime}}C_{-M}(\bm{e}_{12})C_{-M^{\prime}}(\bm{e}_{13})C_{-M^{\prime\prime}}(\bm{e}_{23})C_{2M2M^{\prime}}^{2-M^{\prime\prime}}$

and finally using $z$ axis along $\bm{e}_{12}$ (so $M=0$ in the previous sum)
and $\theta_{i}$ being the inner angles of the 123 triangle with
$\cos\theta_{1}=\bm{e}_{12}.\bm{e}_{13}$ and
$\cos(\pi-\theta_{2})=\bm{e}_{12}.\bm{e}_{23}$. We get

$\displaystyle E$ $\displaystyle=$ $\displaystyle\frac{\langle
1||r(1)||0\rangle^{2}\langle 1||r(2)||0\rangle^{2}\langle
1||r(3)||0\rangle^{2}}{36R_{12}^{3}R_{13}^{3}R_{23}^{3}}\times$
$\displaystyle(1-3(\cos(2\theta_{1})+\cos(2\theta_{2})+\cos(2\theta_{3})))$

So when summing the 6 terms, this finally leads to

$\displaystyle E_{123}^{(3)}$ $\displaystyle=$ $\displaystyle\frac{4\langle
1||r(1)||0\rangle^{2}\langle 1||r(2)||0\rangle^{2}\langle
1||r(3)||0\rangle^{2}}{36R_{12}^{3}R_{13}^{3}R_{23}^{3}}\times$
$\displaystyle(1-3(\cos(2\theta_{1})+\cos(2\theta_{2})+\cos(2\theta_{3})))\times$
$\displaystyle\frac{\Delta(1)+\Delta(2)+\Delta(3)}{(\Delta(1)+\Delta(2))(\Delta(2)+\Delta(3))(\Delta(1)+\Delta(3))}$

So in summary

$\displaystyle E_{123}^{(3)}$ $\displaystyle=$
$\displaystyle\frac{C_{9}}{R_{12}^{3}R_{23}^{3}R_{31}^{3}}\frac{1-3(\cos(2\theta_{1})+\cos(2\theta_{2})+\cos(2\theta_{3}))}{4}$
(29) $\displaystyle C_{9}$ $\displaystyle=$
$\displaystyle\frac{4}{9}r_{SP}^{2}(Cs)r_{SP}^{4}(Ar)\frac{E_{SP}(Cs)+2E_{SP}(Ar)}{2(E_{SP}(Cs)+E_{SP}(Ar))^{2}E_{SP}(Ar)}$

We wrote the $C_{9}$ coefficient such that we can restore the historical
(Axilrod-Tenner-Mutô axilrod1943interaction ; muto1943force  form using
$1-3(\cos(2\theta_{1})+\cos(2\theta_{2})+\cos(2\theta_{3}))=4(1+3\cos(\theta_{1})\cos(\theta_{2})\cos(\theta_{3}))$.

These calculations are technical but similar in a way to the case of two atoms
under the effect of a static external electric field lepers2011long , the
physical picture of the dipole of one atom (1 for instance) interacting with
the vector electric fields produced by the other dipoles (2) and (3) helps to
understand that the final result depends on the angles and distance between
atoms.

####  E.2.3 Interaction with Cs in an excited state $|l=1\,m\rangle$

To our knowledge, no general simple formula has been derived for the dipole-
dipole-dipole interaction with one atom in the excited state (see however Ref.
yan2021long  and reference therein). We derive one formula here. For this, we
start by the fact that in this case, $|0_{1}\rangle=|l=1\,m\rangle$ and
$|m_{1}\rangle=|l=0\,m=0\rangle$. So

$\displaystyle\sum_{m_{1}}\langle 0_{1}|r_{q_{1}}(1)|m_{1}\rangle....\langle
m_{1}|r_{q_{1}^{\prime}}(1)|0_{1}\rangle$ $\displaystyle=\langle
1m|r_{q_{1}}(1)|00\rangle....\langle 00|r_{q_{1}^{\prime}}(1)|1m\rangle$
$\displaystyle=-C_{1m1q^{\prime}_{1}}^{00}\frac{C_{001q_{1}}^{1m}}{\sqrt{3}}\langle
1||r(1)||0\rangle^{2}$
$\displaystyle=\delta_{mq_{1}}\delta_{q_{1},-q_{1}^{\prime}}(-1)^{q_{1}}\frac{\langle
1||r(1)||0\rangle^{2}}{3}$

The second order leads to $C_{6}^{*}=\frac{2}{3}\frac{\langle
1||r(1)||0\rangle^{2}\langle 1||r(2)||0\rangle^{2}}{\Delta(1)+\Delta(2)}$ with

$\displaystyle E_{12}^{(2)}$ $\displaystyle=$
$\displaystyle-\frac{2C_{6}^{*}}{3R_{12}^{6}}=-\frac{C_{6}^{*}(\Uppi)}{R_{12}^{6}}\mbox{
for }m=0$ (30) $\displaystyle E_{12}^{(2)}$ $\displaystyle=$
$\displaystyle-\frac{C_{6}^{*}}{6R_{12}^{6}}=-\frac{C_{6}^{*}(\Sigma)}{R_{12}^{6}}\mbox{
for }m=\pm 1$ $\displaystyle C_{6}^{*}$ $\displaystyle=$
$\displaystyle\frac{2}{3}\frac{r_{SP}^{2}(Cs)r_{SP}^{2}(Ar)}{E_{SP}(Ar)-E_{SP}(Cs)}$

so as expected $C_{6}^{*}(\Upsigma)+2C_{6}^{*}(\Pi)=C_{6}^{*}$.

For the third order, we have for a given numerator term:

$\displaystyle E$ $\displaystyle=$
$\displaystyle-\frac{\sqrt{6}}{R_{12}^{3}}\frac{\sqrt{6}}{R_{13}^{3}}\frac{\sqrt{6}}{R_{23}^{3}}\frac{\langle
1||r(1)||0\rangle^{2}}{3}\frac{\langle 1||r(2)||0\rangle^{2}}{3}\frac{\langle
1||r(3)||0\rangle^{2}}{3}\times$
$\displaystyle\sum_{q_{1},q_{2},q^{\prime}_{1},q^{\prime}_{2},q_{3},q^{\prime}_{3}}\delta_{mq_{1}}\delta_{q_{1},-q_{1}^{\prime}}(-1)^{q_{1}+q_{2}+q_{3}}\delta_{q_{2},-q_{2}^{\prime}}\delta_{q_{3},-q_{3}^{\prime}}\times$
$\displaystyle\sum_{M,M^{\prime},M^{\prime\prime}}(-1)^{M+M^{\prime}+M^{\prime}}C_{-M}(\bm{e}_{12})C_{-M^{\prime}}(\bm{e}_{13})C_{-M^{\prime\prime}}(\bm{e}_{23})\times$
$\displaystyle
C_{1q_{1}1q_{2}}^{2M}C_{1q^{\prime}_{1}1q_{3}}^{2M^{\prime}}C_{1q^{\prime}_{2}1q^{\prime}_{3}}^{2M^{\prime\prime}}$

but using $z$ axis along $\bm{e}_{12}$ (so $M=0$) we have

$\displaystyle E$ $\displaystyle=$
$\displaystyle-\sqrt{\frac{8}{27}}\frac{\langle
1||r(1)||0\rangle^{2}}{R_{12}^{3}}\frac{\langle
1||r(2)||0\rangle^{2}}{R_{13}^{3}}\frac{\langle
1||r(3)||0\rangle^{2}}{R_{23}^{3}}\times$
$\displaystyle\sum_{M^{\prime},M^{\prime\prime}}(-1)^{M^{\prime}+M^{\prime\prime}}C_{0}(\bm{e}_{12})C_{-M^{\prime}}(\bm{e}_{13})C_{-M^{\prime\prime}}(\bm{e}_{23})\times$
$\displaystyle
C_{1m1-m}^{20}\sum_{q_{3}}(-1)^{q_{3}}C_{1-m1q_{3}}^{2M^{\prime}}C_{1m1-q_{3}}^{2M^{\prime\prime}}$

So, $E=\frac{\langle 1||r(1)||0\rangle^{2}\langle 1||r(2)||0\rangle^{2}\langle
1||r(3)||0\rangle^{2}}{54R_{12}^{3}R_{13}^{3}R_{23}^{3}}f(\theta_{1},\theta_{2},\theta_{3},m)$
with

$\displaystyle f(\theta_{1},\theta_{2},\theta_{3},m)$ $\displaystyle=$
$\displaystyle-(1+3\cos(2\theta_{1})+3\cos(2\theta_{2})+9\cos(2\theta_{3}))$
(31) $\displaystyle\mbox{for}\ m=0$ $\displaystyle
f(\theta_{1},\theta_{2},\theta_{3},m)$ $\displaystyle=$
$\displaystyle(5-3\cos(2\theta_{1})-3\cos(2\theta_{2})+9\cos(2\theta_{3}))/4$
(32) $\displaystyle\mbox{for}\ m=\pm 1$

Then, summing the 6 terms (see midzuno1956non ) finally leads to

$\displaystyle E_{123}^{(3)}$ $\displaystyle=$ $\displaystyle\frac{2\langle
1||r(1)||0\rangle^{2}\langle 1||r(2)||0\rangle^{2}\langle
1||r(3)||0\rangle^{2}}{27R_{12}^{3}R_{13}^{3}R_{23}^{3}}\times$ $\displaystyle
f(\theta_{1},\theta_{2},\theta_{3},m)\times$
$\displaystyle\frac{\Delta(1)+\Delta(2)+\Delta(3)}{(\Delta(1)+\Delta(2))(\Delta(2)+\Delta(3))(\Delta(1)+\Delta(3))}$

so,

$\displaystyle E_{123}^{(3)}$ $\displaystyle=$
$\displaystyle\frac{C_{9}^{*}}{R_{12}^{3}R_{13}^{3}R_{23}^{3}}\frac{f(\theta_{1},\theta_{2},\theta_{3},m)}{6}$
(33) $\displaystyle C_{9}^{*}$ $\displaystyle=$
$\displaystyle\frac{4}{9}r_{SP}^{2}(Cs)r_{SP}^{4}(Ar)\frac{2E_{SP}(Ar)-E_{SP}(Cs)}{2(E_{SP}(Ar)-E_{SP}(Cs))^{2}E_{SP}(Ar)}$

###  E.3 Mean field effect and effective two-body potentials

The most accurate way to include this third-order correction for the energy of
atoms A is to sum over all B,C pairs of atoms. But due to the lack of
information on the short-range part and due to the crude (two-level)
estimation made up to now, we can simplify the problem further and use the
mean-field approach done in Ref. stenschke1994effective  (see also
dridi2022development ; muser2022interatomic ) by creating a mean-field
potential $\bar{E}_{AB}^{(3)}=\sum_{C}E_{ABC}^{(3)}$ such that the full
crystal energy
$\frac{1}{2}\left(\sum_{AB}E_{AB}^{(2)}+\frac{1}{3}\sum_{ABC}E_{ABC}^{(3)}\right)$
can be written as sum of two body terms $E_{AB}^{\rm
eff}=(E_{AB}^{(2)}+\frac{1}{3}\bar{E}_{AB}^{(3)}$.

We can thus stay at a two-body level simply by modifying the two-body
interaction between atoms by adding $\frac{1}{3}\bar{E}_{AB}^{(3)}$ to the
$E_{AB}^{(2)}$ potential interaction curves we already have. This mean-field
approach has been proven to be quite accurate for pure crystals, as shown by
the so-called Marcelli-Wang-Sadus potential grimme2016dispersion ; xu2020many
; muser2022interatomic ; deiters2021interatomic ; stroker2022thermodynamic ;
dridi2022development  and as we demonstrate in Fig. 3 in the Ar-Ar case.

####  E.3.1 Mean field homogeneous assumption

Strangely enough, almost no mathematical derivation of the mean field approach
exists though it has been performed as back as in the original paper by Mutô
muto1943force  (whose journal reference is almost always wrong) and by
Stenschkestenschke1994effective . Following these works, we replace the sum
$\bar{E}_{AB}^{(3)}=\sum_{C}E_{ABC}^{(3)}$ by an integral assuming Ar atoms
uniformly distributed with density $\rho=N/V$ (for solid argon the density
$\rho=$280\text{\,}\mathrm{a}_{\mathrm{0}}^{-3}$$ where $a_{0}$ is the Bohr
radius). We thus sum over the C atoms (the number 3 above); that is, we
average the potential $V=\sum_{j}E_{12j}^{(3)}$ as $\bar{V}=\rho\int
Vr^{2}\sin\theta{\rm d}r{\rm d}\theta{\rm d}\varphi=2\pi\rho\int
Vr^{2}\sin\theta{\rm d}r{\rm d}\theta$ in spherical coordinate with
$\theta=\theta_{1}$ and $r=R=R_{12}$. We found simpler to use
$r_{1}=R_{13},r_{2}=R_{23}$ coordinates, as in Ref. stenschke1994effective ,
such that $\bar{V}=2\pi\rho\int V\frac{r_{1}r_{2}}{R}{\rm d}r_{1}{\rm
d}r_{2}$, with an exclusion sphere of radius $\sigma$ around the atoms A=1 and
B=2. Equations (29) and (33) lead to

$\displaystyle\bar{E}_{12}^{(3)}$ $\displaystyle=$ $\displaystyle
2\pi\rho\frac{C_{9}}{R^{6}}\frac{4}{3}$ (34)
$\displaystyle\bar{E^{*}_{12}}^{(3)}$ $\displaystyle=$ $\displaystyle
2\pi\rho\frac{C_{9}^{*}}{R^{6}}\frac{4}{9}(4-3\log[(R^{2}-\sigma^{2})/\sigma^{2}])\mbox{
for }m=0$ $\displaystyle\bar{E^{*}_{12}}^{(3)}$ $\displaystyle=$
$\displaystyle
2\pi\rho\frac{C_{9}^{*}}{R^{6}}\frac{2}{9}(-1+3\log[(R^{2}-\sigma^{2})/\sigma^{2}])\mbox{
for }m=\pm 1$

Therefore, the triple dipole interaction does not depend on $\sigma$ only when
the state is isotropic, such as the S ground state, or when summing over
$m=-1,0,1$ for the P state.

####  E.3.2 Effective potentials

All results are summarized by Eqs (28-34) from where we can extract the
effective potentials for the internuclear distance $R$:

$\displaystyle E_{AB}^{\rm eff}(R)$ $\displaystyle=$ $\displaystyle
E_{AB}^{(2)}(R)\left(1-\frac{8\pi\rho}{9}\frac{C_{9}}{C_{6}}\right)$ (35)
$\displaystyle E_{AB}^{*,\rm eff}(\Upsigma)$ $\displaystyle=$ $\displaystyle
E_{AB}^{(2)}(\Upsigma)\left(1-\frac{4\pi\rho}{9}\frac{C_{9}^{*}}{C_{6}^{*}}\left(4-3\log\left[\frac{R^{2}-\sigma^{2}}{\sigma^{2}}\right]\right)\right)$
$\displaystyle E_{AB}^{*,\rm eff}(\Uppi)$ $\displaystyle=$ $\displaystyle
E_{AB}^{(2)}(\Uppi)\left(1-\frac{8\pi\rho}{9}\frac{C_{9}^{*}}{C_{6}^{*}}\left(-1+3\log\left[\frac{R^{2}-\sigma^{2}}{\sigma^{2}}\right]\right)\right)$

This form, derived from the long-range triple dipole interaction, can however
naturally be extended to the full range of the potential by simply keeping the
same formula for all $R$. The only requirement is the choice of the cut-off
because $R>2\sigma$. Because the short-range part of the potential is clearly
not a dipole-dipole one, we smoothly reduce the $C_{9}^{*}$ part for
$R<2\sigma$ by multiplying it in Eq. (35), for all $R$, by a sharp cut-off
function $1-e^{-(R/2\sigma)^{10}}$. A reasonable value for $\sigma$ might be
the nearest neighbor distance in the fcc crystal, the lattice size or the
first nearest neighbor Cs-Ar atoms in a given trapping site (see Fig. 6). This
may depend on the chosen trapping site to be studied. Another cut-off, in
principle independent, should be of the order of the LeRoy radius that
corresponds to the change in the multipole expansion at small $R$. In that
case the theoretical curves are too different from the long-range part, which
we used to derive the formula le1974long ; kielkopf1974semiempirical . From
all this, we found that a reasonable value can be
$\sigma\sim$8\text{\,}\mathrm{a}_{\mathrm{0}}$$. However, all these parameters
(cut-offs and the power 10 in the sharp function) are arbitrary choices that
can be optimized. We have modified them a bit (typically by a factor of 2) to
produce the uncertainty in the line position shown in Fig. 5.

####  E.3.3 Approximate calculation of $C_{6}$ and $C_{9}$ coefficients

Our previous formulae for second Eq. (28) and (30), and third-order
perturbations Eq. (29) and (33), contain the van der Waals interaction
coefficients $C_{6},C_{6}^{*},C_{9},C_{9}^{*}$, which we now need to evaluate.
For this, we will use an effective two-level approximation, known to be 25%
accurate for ground state interactions between rare-gas and alkali atoms
tang1969dynamic ; diaz1980combination , and that we are going to use for the
excited states as well. We will also compare the value found by fitting them
to the long-range part of our potential curves to existing results to assess
the method accuracy.

From now on, we will use atomic units for this calculation.

The atomic unit value for the $C_{6}$ coefficient for Ar-Ar, found by fitting
the theoretical Ar-Ar potential curve by $-C_{6}/R^{6}$, is 67; and 380 for
Ar-Cs. Using Eq. (28) with the well-known experimental values,
$E_{SP}(Cs)\approx 0.0524$ and $E_{SP}(Ar)\approx 0.43$ for the first excited
$P$ levels, leads to $r_{SP}(Cs)=5.4$ and $r_{SP}(Ar)=3.0$.

These values are acceptable because they are reasonable: the $C_{6}$ Cs-Cs
coefficient is 5500 (compare to the real value of nearly 6700 gould2016c ) ;
$C_{9}(Ar-Ar-Ar)\approx 650$ (compared to $525$ in Refs. bell1966van ;
tang1969dynamic ); and $C_{9}(Cs-Cs-Cs)\approx 1500000$ (compared to
$C_{9}(Cs-Cs-Cs)\approx 2200000$ in Ref. diaz1980combination ). However, these
values for $r_{SP}(Cs),r_{SP}(Ar),E_{SP}(Cs)$ and $E_{SP}(Ar)$ lead to quite
wrong values for the excited Cs∗-Ar coefficient of
$C_{6}^{*}(\Uppi)=\frac{2}{3}C_{6}^{*}=310$ compared to 600 in our potential
curves and $C_{6}^{*}(\Upsigma)=\frac{1}{6}C_{6}^{*}=80$ compared to 270.
Unfortunately, it is not possible to solve this discrepancy for all excited
states at once simply because Eq. (30) gives a ratio 4 between
$C_{6}^{*}(\Pi)$ and $C_{6}^{*}(\Sigma)$ while in our case, the theoretical
potential curves give only a factor $\sim 2$. This factor 2 is indeed found in
many other cases bussery1986calculated ; zhang2007long .

This discrepancy probably comes from the fact that we have neglected other
excited states (like Cs$(7s)$) that are quite near the Cs excited (6p) level.
But, because our intention is to simulate the Cs spectrum in argon, we will
focus more on the Cs-Ar and Cs∗-Ar interaction curves than on the Ar-Ar one in
order to try to reduce the discrepancy. For this, we choose an intermediate
strategy and allow modification of the dipole transitions $r_{SP}(Cs)$,
$r_{SP}(Ar)$ and the energy level $E_{SP}(Ar)$, that become effective values.
This allows us to better reproduce the long-range part of the theoretical
potential curve for Cs-Ar and Cs∗-Ar. For this, we start with fitting the
theoretical potential curve that leads to $\frac{C_{6}^{*}}{C_{6}}\approx 3$.
Then Eq. (28) and (30) give
$\frac{C_{6}^{*}}{C_{6}}=\frac{E_{SP}(Ar)+E_{SP}(Cs)}{E_{SP}(Ar)-E_{SP}(Cs)}$.
This leads to $E_{SP}(Ar)=2E_{SP}(Cs)=0.10$. Then, using expression Eq. (28)
for the $C_{6}$ coefficient for Ar-Ar (67) and Ar-Cs (380) leads to (in atomic
units):

$\displaystyle E_{SP}(Cs)$ $\displaystyle=$ $\displaystyle 0.052$ (36)
$\displaystyle E_{SP}(Ar)$ $\displaystyle\approx$ $\displaystyle 0.10$
$\displaystyle r_{SP}(Cs)$ $\displaystyle\approx$ $\displaystyle 4.4$
$\displaystyle r_{SP}(Ar)$ $\displaystyle\approx$ $\displaystyle 2.1$
$\displaystyle C_{6}$ $\displaystyle\approx$ $\displaystyle 370$ (37)
$\displaystyle C_{6}^{*}$ $\displaystyle\approx$ $\displaystyle 1200$
$\displaystyle C_{9}$ $\displaystyle\approx$ $\displaystyle 9100$
$\displaystyle C_{9}^{*}$ $\displaystyle\approx$ $\displaystyle 54000$

Figure 10: Cs(6p)-Ar potential curves compared to the one with the effective
third-order effect included.

The long-range parts of some potential curves are shown in Fig. 10. These
parameter values, and the modification by a factor $\sim 2$ for the cut-off
and the power of the sharp function, have been used in order to produce the
uncertainty in the line positions shown in Fig. 5.

## References

  * (1) A. Barnes, W. Orville-Thomas, A. Müller, and R. Gaufrès,  Matrix isolation spectroscopy (Springer Dordrecht, 1981).
  * (2) M. Almond and A. Downs,  Spectroscopy of Matrix Isolated Species (Wiley, Chichester, 1989).
  * (3) V. E. Bondybey, A. M. Smith, and J. Agreiter,  Chemical reviews 96, 2113 (1996).
  * (4) C. Crepin-Gilbert and A. Tramer,  International Reviews in Physical Chemistry 18, 485 (1999).
  * (5) B. M. Davis, B. Gervais, and J. G. McCaffrey,  The Journal of Chemical Physics 148, 124308 (2018).
  * (6) G. K. Ozerov, D. S. Bezrukov, and A. A. Buchachenko,  Physical Review B 103, 184110 (2021).
  * (7) C. Pryor and F. Wilczek,  Physics Letters B 194, 137 (1987).
  * (8) M. Arndt, S. Kanorsky, A. Weis, and T. Hänsch,  Physics Letters A 174, 298 (1993).
  * (9) M. Kozlov and A. Derevianko,  Physical review letters 97, 063001 (2006).
  * (10) A. C. Vutha, M. Horbatsch, and E. A. Hessels,  Atoms 6 (2018).
  * (11) S. Upadhyay et al.,  Physical Review A 100, 063419 (2019).
  * (12) C. Braggio et al.,  Applied Sciences 12, 6492 (2022).
  * (13) D. Budker et al.,  arXiv preprint arXiv:2203.09488 (2022).
  * (14) M. Safronova et al.,  Reviews of Modern Physics 90, 025008 (2018).
  * (15) C. Ien, V. Bowers, E. Cochran, and S. Foner,  Physical Review 126, 1749 (1962).
  * (16) J. Goldsborough and T. Koehler,  Physical Review 133, A135 (1964).
  * (17) W. Weyhmann and F. Pipkin,  Physical Review 137, A490 (1965).
  * (18) T. Kanda and T. Ebisu,  Journal of the Physical Society of Japan 31, 957 (1971).
  * (19) L. Balling and J. Wright,  The Journal of Chemical Physics 78, 592 (1983).
  * (20) J. L. Hollenberg and D. A. Dows,  The Journal of chemical physics 37, 1300 (1962).
  * (21) P. Groner, I. Stolkin, and H. H. Gunthard,  Journal of Physics E: Scientific Instruments 6, 122 (1973).
  * (22) G. J. Jiang, W. B. Person, and K. G. Brown,  The Journal of Chemical Physics 62, 1201 (1975).
  * (23) A. Sinnock and B. Smith,  Physical Review 181, 1297 (1969).
  * (24) R. Hahn et al.,  Review of Scientific Instruments 93, 043302 (2022).
  * (25) D. A. Steck,  Cesium D line data (2019).
  * (26) S. A. Nepijko, I. Rabin, and W. Schulze,  ChemPhysChem 6, 235 (2005).
  * (27) J. D. Patterson and B. C. Bailey,  Solid-state physics: introduction to the theory (third edition) (Springer Science & Business Media, 2018).
  * (28) B. Smith,  Contemporary Physics 11, 125 (1970).
  * (29) J. Venables and B. Smith,  (1977).
  * (30) S. Kovalenko and N. Bagrov,  Soviet Physics Solid State, USSR 11, 2207 (1970).
  * (31) J. ShakhsEmampour, R. Pyzalski, M. Vala, and J.-C. Rivoal,  Journal de Physique 45, 953 (1984).
  * (32) R. Rudman,  Acta Crystallographica Section A: Crystal Physics, Diffraction, Theoretical and General Crystallography 34, 639 (1978).
  * (33) K. S. Song and R. T. Williams,  Self-trapped excitons (Springer Science & Business Media, 2013).
  * (34) W. Schulze and D. Kolb,  Journal of the Chemistry Society, Faraday Transactions 2 70, 1098 (1974).
  * (35) G. L. Pollack,  Reviews of Modern Physics 36, 748 (1964).
  * (36) Y. K. Tovbin, S. Titov, and V. Komarov,  Physics of the Solid State 57, 360 (2015).
  * (37) Y. K. Tovbin, V. Komarov, and E. Gvozdeva,  Protection of Metals and Physical Chemistry of Surfaces 53, 591 (2017).
  * (38) N. Tishchenko,  Physica Status Solidi. A, Applied Research 73, 279 (1982).
  * (39) N. Tishchenko,  physica status solidi (a) 83, 513 (1984).
  * (40) I. V. Leibin, I. S. Kalinina, D. S. Bezrukov, and A. A. Buchachenko,  The Journal of Chemical Physics 154, 044305 (2021).
  * (41) L. Blank, D. E. Weeks, and G. S. Kedziora,  The Journal of Chemical Physics 136, 124315 (2012).
  * (42) T. Kobayashi, K. Yuki, and L. Matsuoka,  Chemistry Letters 45, 1400 (2016).
  * (43) A. A. Medvedev, V. V. Meshkov, A. V. Stolyarov, and M. C. Heaven,  Physical Chemistry Chemical Physics 20, 25974 (2018).
  * (44) G. K. Ozerov, D. S. Bezrukov, and A. A. Buchachenko,  Low Temperature Physics 45, 301 (2019).
  * (45) N. N. Kleshchina, I. S. Kalinina, I. V. Leibin, D. S. Bezrukov, and A. A. Buchachenko,  The Journal of Chemical Physics 151, 121104 (2019).
  * (46) U. K. Deiters and R. J. Sadus,  The Journal of Chemical Physics 150, 134504 (2019).
  * (47) V. F. Lotrich and K. Szalewicz,  Physical review letters 79, 1301 (1997).
  * (48) P. Schwerdtfeger, R. Tonner, G. E. Moyano, and E. Pahl,  Angewandte Chemie International Edition 55, 12200 (2016).
  * (49) D. Bezrukov, N. Kleshchina, I. Kalinina, and A. Buchachenko,  Russian Journal of Physical Chemistry A 93, 1505 (2019).
  * (50) M. Lara-Moreno, J. Alvarez-Hernández, H. Negrín-Yuvero, J. G. McCaffrey, and G. Rojas-Lorenzo,  Low Temperature Physics 45, 697 (2019).
  * (51) N. N. Kleshchina, K. A. Korchagina, D. S. Bezrukov, and A. A. Buchachenko,  The Journal of Physical Chemistry A 121, 2429 (2017).
  * (52) L.-G. Tao et al.,  The Journal of Chemical Physics 143, 174306 (2015).
  * (53) N. Lam, N. Doan, and L. Dagens,  Journal of Physics F: Metal Physics 15, 799 (1985).
  * (54) N. Lam and L. Dagens,  Journal of Physics F: Metal Physics 16, 1373 (1986).
  * (55) M. Sabochick, S. Yip, and N. Lam,  Journal of Physics F: Metal Physics 18, 349 (1988).
  * (56) H. Zenia, K. Lounis, E. Megchiche, and C. Mijoule,  Computational Materials Science 124, 428 (2016).
  * (57) K. Lounis, H. Zenia, E. Megchiche, and C. Mijoule,  Computational Materials Science 118, 279 (2016).
  * (58) R. Cotterill and M. Doyama,  Calculation of the Properties of Vacancies and Interstitials: Proceedings 13, 47 (1966).
  * (59) A. Damask,  Science 162, 448 (1968).
  * (60) R. Johnson,  Journal of Physics F: Metal Physics 3, 295 (1973).
  * (61) L. Y. Nemirovich-Danchenko, A. Lipnitskiui, and S. Kul’kova,  Physics of the Solid State 49, 1079 (2007).
  * (62) H. Wang, D. Rodney, D. Xu, R. Yang, and P. Veyssière,  Philosophical Magazine 93, 186 (2013).
  * (63) J. Peng, S. Bahl, A. Shyam, J. A. Haynes, and D. Shin,  Acta Materialia 196, 747 (2020).
  * (64) We use the following articles cotterill1966morse ; damask1968defects ; johnson1973empirical ; lam1985multiple ; lam1986calculations ; sabochick1988atomistic ; nemirovich2007vacancies ; wang2013defect ; zenia2016stability ; lounis2016stability ; lara2019effects ; peng2020solute  and using their notations, this include for instance for 0 vacancy octahedral, tetrahedral, crowdion and $<$100$>$, $<$110$>$ and $<$111$>$ interstitial sites; for 1 vacancy cornered or centered substitutional sites; for 2 vacancies corner to: center (2v-a) to corner (2v-b) or to tetrahedral (2v-c) sites; for 3 vacancies corner to center-center (3v-a) to corner-center (3v-b) to center-corner (3v-c: linear) sites; for 4 vacancies tetrahedral, rhombic or squared; the 5,7 and 8 vacancies are typically taken from the 6 vacancy plus (or minus) centered, tetraedral or corner sites.
  * (65) D. P. Kingma and J. Ba,  arXiv preprint arXiv:1412.6980, Published as a conference paper at the 3rd International Conference for Learning Representations, San Diego, 2015 (2014).
  * (66) T. Dozat,  (2016).
  * (67) S. Ruder,  arXiv preprint arXiv:1609.04747 (2016).
  * (68) S. H. Haji and A. M. Abdulazeez,  PalArch’s Journal of Archaeology of Egypt/Egyptology 18, 2715 (2021).
  * (69) J. Simons,  An introduction to theoretical chemistry (Cambridge University Press, 2003).
  * (70) M. Heaven and A. Stolyarov,  Potential energy curves for alkali metal-rare gas exciplex lasers,  in 41st Plasmadynamics and Lasers Conference, p. 4877, 2010.
  * (71) M. Ryan, M. Collier, P. d. Pujo, C. Crépin, and J. G. McCaffrey,  The Journal of Physical Chemistry A 114, 3011 (2010).
  * (72) E. Jacquet, D. Zanuttini, J. Douady, E. Giglio, and B. Gervais,  The Journal of chemical physics 135, 174503 (2011).
  * (73) M. Diaz Peña, C. Pando, and J. Renuncio,  The Journal of Chemical Physics 73, 1750 (1980).
  * (74) S. Dridi, M. B. Amar, M. Abderraba, and J.-P. Passarello,  Fluid Phase Equilibria 562, 113563 (2022).
  * (75) H. Stenschke,  The Journal of chemical physics 100, 4704 (1994).
  * (76) B. Akhouri and J. Solana,  Physica A: Statistical Mechanics and its Applications , 128280 (2022).
  * (77) H. Berendsen,  Simulating the physical world (Cambridge University Press Cambridge, 2007).
  * (78) P. Ströker, R. Hellmann, and K. Meier,  Physical Review E 105, 064129 (2022).
  * (79) M. Lax,  The Journal of chemical physics 20, 1752 (1952).
  * (80) E. J. Heller,  The Journal of Chemical Physics 68, 2066 (1978).
  * (81) R. S. Mulliken,  The Journal of Chemical Physics 55, 309 (1971).
  * (82) J. Tellinghuisen,  The franck-condon principle,  in Photons and Continuum States of Atoms and Molecules: Proceedings of a Workshop Cortona, Italy, June 16–20, 1986, pp. 149–156, Springer, 1987.
  * (83) L. Barragán-Gil and R. Walser,  American Journal of Physics 86, 22 (2018).
  * (84) J. P. Bergsma, P. H. Berens, K. R. Wilson, D. R. Fredkin, and E. J. Heller,  The Journal of Physical Chemistry 88, 612 (1984).
  * (85) J. A. Boatz and M. E. Fajardo,  The Journal of chemical physics 101, 3472 (1994).
  * (86) E. Pavarini,  Correlated electrons: from models to materials 2, 6 (2012).
  * (87) B. S. Tsukerblat,  Group theory in chemistry and spectroscopy: a simple guide to advanced usage (Courier Corporation, 2006).
  * (88) M. S. Dresselhaus, G. Dresselhaus, and A. Jorio,  Group theory: application to the physics of condensed matter (Springer Science & Business Media, 2007).
  * (89) W. Hergert and R. M. Geilhufe,  Group Theory in Solid State Physics and Photonics: Problem Solving with Mathematica (Wiley-VCH, 2018),  isbn: 978-3-527-41133-7.
  * (90) J. Van Vleck,  The Journal of Chemical Physics 7, 72 (1939).
  * (91) Y. Toyozawa and M. Inoue,  Journal of the Physical Society of Japan 21, 1663 (1966).
  * (92) K. Cho,  Journal of the Physical Society of Japan 25, 1372 (1968).
  * (93) P. Lund, D. Smith, S. Jacobs, and P. Schatz,  The Journal of Physical Chemistry 88, 31 (1984).
  * (94) J. Ammeter and D. Schlosnagle,  The Journal of Chemical Physics 59, 4784 (1973).
  * (95) I. Bersuker,  Coordination Chemistry Reviews 14, 357 (1975).
  * (96) F. Forstmann, D. Kolb, D. Leutloff, and W. Schulze,  The Journal of Chemical Physics 66, 2806 (1977).
  * (97) D. Kolb, D. Leutloff, W. Schulze, and F. Forstmann,  Berichte der Bunsengesellschaft für physikalische Chemie 82, 33 (1978).
  * (98) M. O’brien,  Journal of Physics C: Solid State Physics 4, 2524 (1971).
  * (99) H.-J. Stöckmann,  Zeitschrift für Physik B Condensed Matter 54, 229 (1984).
  * (100) I. Bersuker,  The Jahn-Teller effect and vibronic interactions in modern chemistry (Springer Science & Business Media, 2013).
  * (101) M. Jakob, H. Micklitz, and K. Luchner,  Physics Letters A 57, 67 (1976).
  * (102) M. Jakob, H. Micklitz, and K. Luchner,  Physics Letters A 61, 265 (1977).
  * (103) K. Luchner and H. Micklitz,  Journal of Luminescence 18, 882 (1979).
  * (104) R. Lambo et al.,  The Journal of chemical physics 137, 204315 (2012).
  * (105) R. Lambo et al.,  Physical Review A 104, 062809 (2021).
  * (106) M. Wagner,  The Journal of Chemical Physics 41, 3939 (1964).
  * (107) T. H. Keil,  Physical Review 140, A601 (1965).
  * (108) X. Shi-Jie,  Acta Phys. Sin 68 (2019).
  * (109) Y. Zhang,  Journal of Semiconductors 40, 091102 (2019).
  * (110) M. Sturge,  The jahn-teller effect in solids,  in The Jahn-Teller Effect in Solids, edited by F. Seitz, D. Turnbull, and H. Ehrenreich, , Solid State Physics Vol. 20, pp. 91–211, Academic Press, 1968.
  * (111) C. A. Bates,  Physics reports 35, 187 (1978).
  * (112) C. Bates, J. Dunn, and E. Sigmund,  Journal of Physics C: Solid State Physics 20, 1965 (1987).
  * (113) I. Bersuker,  The Jahn-Teller (Cambridge University Press, New York, 2006).
  * (114) M. O’Brien and S. Evangelou,  Journal of Physics C: Solid State Physics 13, 611 (1980).
  * (115) J. Rose, D. Smith, B. Williamson, P. Schatz, and M. O’Brien,  The Journal of Physical Chemistry 90, 2608 (1986).
  * (116) B. Hüpper and B. Eckhardt,  Physical Review A 57, 1536 (1998).
  * (117) O. Anatole von Lilienfeld and A. Tkatchenko,  The Journal of chemical physics 132, 234109 (2010).
  * (118) M. Gross and F. Spiegelmann,  The European Physical Journal D-Atomic, Molecular, Optical and Plasma Physics 4, 219 (1998).
  * (119) J. D. Hewitt et al.,  The Journal of Physical Chemistry A (2023).
  * (120) R. A. Corbin and M. E. Fajardo,  The Journal of chemical physics 101, 2678 (1994).
  * (121) S. L. Laursen and H. E. Cartland,  The Journal of chemical physics 95, 4751 (1991).
  * (122) M. Guarise,  The European Physical Journal Plus 137, 1 (2022).
  * (123) R. Lambo et al.,  arXiv preprint arXiv:2212.09232 (2022).
  * (124) G. Koyanagi et al.,  arXiv preprint arXiv:2211.14804 (2022).
  * (125) R. M. Geilhufe and W. Hergert,  Frontiers in Physics 6, 86 (2018).
  * (126) A. Devonshire,  Proceedings of the Royal Society of London. Series A-Mathematical and Physical Sciences 153, 601 (1936).
  * (127) S. K. Misra,  Multifrequency Electron Paramagnetic Resonance: Theory and Applications , 327 (2011).
  * (128) I. D. Ryabov,  Applied Magnetic Resonance 35, 481 (2009).
  * (129) C. D. Lewis and D. E. Weeks,  The Journal of Physical Chemistry A 121, 3340 (2017).
  * (130) D. Varshalovich, A. Moskalev, and V. Khersonskii,  Quantum theory of angular momentum (World Scientific, 1988).
  * (131) P. P. Man,  Concepts in Magnetic Resonance Part A 45, e21385 (2016).
  * (132) J. W. Kenney III, J. A. Boatz, and H. A. Terrill Vosbein,  International journal of quantum chemistry 103, 854 (2005).
  * (133) U. Hergenhahn,  Journal of Physics B: Atomic, Molecular and Optical Physics 37, R89 (2004).
  * (134) W. Siebrand,  The Journal of Chemical Physics 46, 440 (1967).
  * (135) K. F. Freed and S. Lin,  Chemical Physics 11, 409 (1975).
  * (136) E. J. Heller,  The Journal of Chemical Physics 65, 1289 (1976).
  * (137) S. Mukamel,  The Journal of Chemical Physics 77, 173 (1982).
  * (138) J. Shao, J.-L. Liao, and E. Pollak,  The Journal of chemical physics 108, 9711 (1998).
  * (139) I. M. Ansari, E. R. Heller, G. Trenins, and J. O. Richardson,  Philosophical Transactions of the Royal Society A 380, 20200378 (2022).
  * (140) E. Wigner,  Physical Review 40, 749 (1932).
  * (141) V. Tatarskiĭ,  Soviet Physics Uspekhi 26, 311 (1983).
  * (142) M. Hillery, R. F. O’Connell, M. O. Scully, and E. P. Wigner,  Physics reports 106, 121 (1984).
  * (143) W. B. Case,  American Journal of Physics 76, 937 (2008).
  * (144) D. K. Ferry and M. Nedjalkov,  The Wigner function in science and technology (IoP Publishing, 2018).
  * (145) B. Segev,  Journal of Optics B: Quantum and Semiclassical Optics 5, S381 (2003).
  * (146) Y. Japha and B. Segev,  Physical Review A 65, 063411 (2002).
  * (147) A. Sergeev and B. Segev,  The Journal of chemical physics 118, 5852 (2003).
  * (148) E. Pollak and J. Cao,  Physical Review A 107, 022203 (2023).
  * (149) S. Kallush, B. Segev, A. Sergeev, and E. Heller,  The Journal of Physical Chemistry A 106, 6006 (2002).
  * (150) B. Segev and A. Sergeev,  Chemical physics letters 367, 382 (2003).
  * (151) E. R. Heller and J. O. Richardson,  The Journal of Chemical Physics 152, 244117 (2020).
  * (152) F. Agostini and B. F. Curchod,  Chemistry without the born–oppenheimer approximation, 2022.
  * (153) R. Bell,  Journal of Physics B: Atomic and Molecular Physics (1968-1987) 3, 751 (1970).
  * (154) Y. Midzuno and T. Kihara,  Journal of the Physical Society of Japan 11, 1045 (1956).
  * (155) B. Axilrod and E. Teller,  The Journal of Chemical Physics 11, 299 (1943).
  * (156) Y. Muto,  Nippon Sugaku-Buturigakkwaishi 17, 629 (1943).
  * (157) M. Lepers and O. Dulieu,  Physical Chemistry Chemical Physics 13, 19106 (2011).
  * (158) P.-G. Yan, L.-Y. Tang, Z.-C. Yan, and J. F. Babb,  Physical Review A 104, 022807 (2021).
  * (159) M. H. Muser, S. V. Sukhomlinov, and L. Pastewka,  arXiv preprint arXiv:2204.09563 (2022).
  * (160) S. Grimme, A. Hansen, J. G. Brandenburg, and C. Bannwarth,  Chemical reviews 116, 5105 (2016).
  * (161) P. Xu, M. Alkan, and M. S. Gordon,  Chemical Reviews 120, 12343 (2020).
  * (162) U. K. Deiters and R. J. Sadus,  The Journal of Physical Chemistry B 125, 8522 (2021).
  * (163) R. J. Le Roy,  Canadian Journal of Physics 52, 246 (1974).
  * (164) J. Kielkopf,  The Journal of Chemical Physics 61, 4733 (1974).
  * (165) K. Tang,  Physical Review 177, 108 (1969).
  * (166) T. Gould and T. Bucko,  Journal of chemical theory and computation 12, 3603 (2016).
  * (167) R. Bell and A. Kingston,  Proceedings of the Physical Society (1958-1967) 88, 901 (1966).
  * (168) B. Bussery, M. Aubert-Frécon, and M. Saute,  Chemical physics 109, 39 (1986).
  * (169) J.-Y. Zhang and J. Mitroy,  Physical Review A 76, 022705 (2007).


# Fidelity and variability in the interlayer electronic structure of the
kagome superconductor CsV3Sb5

Aurland K. Watkins<EMAIL_ADDRESS>Materials Department, University of
California, Santa Barbara, California 93106  Dirk Johrendt  Department Chemie,
Ludwig-Maximilians-Universität München, Munich, Germany  Vojtech Vlcek
Department of Chemistry and Biochemistry, University of California, Santa
Barbara, California 93106  Materials Department, University of California,
Santa Barbara, California 93106  Stephen D. Wilson  Materials Department,
University of California, Santa Barbara, California 93106  Ram Seshadri
Materials Department, University of California, Santa Barbara, California
93106  Department of Chemistry and Biochemistry, University of California,
Santa Barbara, California 93106

###### Abstract

The _A_ V3Sb5 (_A_ = K, Rb, Cs) kagome materials host an interplay of emergent
phenomena including superconductivity, charge density wave states, and non-
trivial electronic structure topology. The band structures of these materials
exhibit a rich variety of features like Dirac crossings, saddle points
associated with van Hove singularities, and flat bands prompting significant
investigations into the in-plane electronic behavior. However, recent findings
including the charge density wave ordering and effects due to pressure or
chemical doping point to the importance of understanding interactions between
kagome layers. Probing this _c_ -axis electronic structure via experimental
methods remains challenging due to limitations of the crystals and, therefore,
rigorous computational approaches are necessary to study the interlayer
interactions. Here we use first-principles approaches to study the electronic
structure of CsV3Sb5 with emphasis on the _k z_ dispersion. We find that the
inclusion of nonlocal and dynamical many-body correlation has a substantial
impact on the interlayer band structure. We present new band behavior that
additionally supports the integration of symmetry in accurately plotting
electronic structures and influences further analysis like the calculation of
topological invariants.

##  I Introduction

First-principles electronic structure calculations are critical for the
identification and characterization of topological materials given the limited
and complex experimental tools designed to confirm non-trivial electronic
structure topology.[1, 2] While bulk band structure features like band
crossings or flat bands can point to non-trivial topology, confirmation
requires additional techniques. There are three common computational or
theoretical methods relying on the bulk band structure to classify topological
materials: 1) adiabatically evolving the Hamiltonian to match the band
structure of a material with known topology, 2) computing the surface
electronic structure to directly detect topological surface states, and 3)
calculating a topological invariant. This final method employs the underlying
symmetries and features of the bulk band structure to produce an index
describing the topology.[3]

A $\mathbb{Z}_{2}$ invariant ($\mathbb{\nu}$) is used to distinguish between
trivial ($\mathbb{\nu}$ = 0) and non-trivial ($\mathbb{\nu}$ = 1) time-
reversal invariant systems with a bulk electronic gap by specifying whether
topological surface states are expected to reside in this gap.[4, 5, 6] For
systems that additionally possess inversion symmetry, Fu and Kane developed a
protocol for calculating $\mathbb{Z}_{2}$ invariants that relies on the parity
(inversion) of the Bloch wave functions defined at specific _k_ -points that
are invariant under time-reversal symmetry. In a 3D Brillouin zone, there are
eight of these _k_ -points known as time-reversal invariant momenta (TRIM)
points. In the Fu-Kane protocol, the product of band parities across all bands
below the gap is calculated at individual TRIM points. Subsequently, a product
across all the TRIM point products is taken and the $\mathbb{Z}_{2}$ invariant
is calculated.[7]

While this method was originally developed for insulators or materials with
band gaps, this type of $\mathbb{Z}_{2}$ analysis does not depend on band
filling and, therefore, can be applied to metallic systems with partial band
occupancy.[8, 9, 10] However, to extend this calculation to metals, the band
structure must have a continuous gap throughout the Brillouin zone (BZ)
located near the Fermi level. This requirement ensures that the parity
products at each of the TRIM points consistently include the same bands. This
method was previously applied to metals featuring kagome nets to support their
non-trivial $\mathbb{Z}_{2}$ classification.[11, 12]

Figure 1: Crystal structure depictions of CsV3Sb5. (a) Quasi-2-dimensional
structure highlighting separation between kagome sheets. (b) and (c)
Coordination environment within and adjacent to kagome layers.

The kagome net of corner-sharing triangles (Figure 1) has drawn significant
interest as a platform for a variety of instabilities and emergent phenomena
like spin liquid states, superconductivity, charge density waves, and
topological surface states.[13, 14, 15] Recently, materials within the _A_
V3Sb5 (_A_ = K, Rb, Cs) kagome family were categorized as $\mathbb{Z}_{2}$
topological metals following the Fu-Kane method.[11, 12] The presence of two
continuous gaps throughout the BZ near _E F_ allowed for the calculation of
multiple $\mathbb{Z}_{2}$ invariants, both of which indicated non-trivial
topology suggesting that CsV3Sb5 exhibits topological surface states analogous
to a topological insulator. This original $\mathbb{Z}_{2}$ analysis relied on
a band structure that was obtained via a structural relaxation with the PBE
functional [16] and D3 correction with the addition of spin-orbit coupling
(SOC) in the self-consistent calculation. The D3 correction is an empirical
van der Waals total energy correction that was identified as a necessary
parameter to more closely match experimental lattice parameters (specifically
along the _c_ -axis) during relaxation.[17, 18]

Here, we present electronic structure analysis of CsV3Sb5, with special
consideration of anisotropic nonlocal correlation along the _k z_ or inter-
kagome-plane direction. Within the framework of density functional theory
(DFT), we find that including certain computational parameters during a
structural relaxation gives rise to a symmetry-allowed band crossing
previously not identified between $\Gamma$ and A as highlighted in Figure 2.
Yet, ultimately, this study suggests that it is critical to include nonlocal
and dynamical many-body correlation through the application of the _GW_
approximation to obtain accurate interlayer band behavior. The resulting
electronic structure shows significant shifts of the bands and the absence of
the symmetry-allowed crossing seen in standard DFT band structures. This
seemingly minor detail in the electronic structure impacts the calculation of
$\mathbb{Z}_{2}$ invariants in this material and points to a larger issue of
correct band identification or tracking when plotting band structures.

##  II Methods

First-principles calculations were performed using the Vienna Ab initio
Simulation Package (VASP) version 5.4.4. Pseudopotentials following the
projector-augmented wave method were selected with the following valence
configurations: Cs (5 _s_ 25 _p_ 66 _s_ 1), V (3 _s_ 23 _p_ 64 _s_ 23 _d_ 3),
Sb (5 _s_ 25 _p_ 3). An 11 $\times$ 11 $\times$ 5 $\Gamma$-centered _k_ -point
mesh was automatically generated by VASP and the plane-wave energy cutoff was
set to 500 eV. Geometric optimization of the experimental structure (ICSD
31841) was performed with all parameters in the description (eg. SOC, D3, TS,
etc.) unless otherwise stated and all degrees of freedom (atomic positions and
cell shape/volume) were allowed to relax until forces were converged within
10$-$7 eV/Å. All calculations had an energy convergence better than 10$-$8 eV.
A _k_ -point path for the band structure was generated using the AFLOW online
tool and the density of the path was set 50 _k_ -points per path segment (for
example between $\Gamma$ and A) to ensure high resolution of band features in
the final plot.[19] Irreducible representations were obtained using the Irvsp
program that interfaces with VASP outputs.[20]

_GW_ calculations implemented in VASP were performed on the experimental
(unrelaxed) structure and employed the PBE electronic structure as a starting
point. Both single-shot _G_ 0 _W_ 0 and partially self-consistent eigenvalue
_GW_ 0 (ev _GW_ 0) calculations were performed. _GW_ potentials were selected
based on the VASP recommendations with the following valence configurations:
Cs (5 _s_ 25 _p_ 66 _s_ 1), V (3 _s_ 23 _p_ 64 _s_ 23 _d_ 3), Sb (4 _d_ 105
_s_ 25 _p_ 3). The total number of bands was set to 4000 and the number of
frequency points was set to 100. Given the high computational cost of _GW_
calculations, a smaller 4 $\times$ 4 $\times$ 2 _k_ -point grid was used. This
setup yields quasiparticle energies converged to within 50 meV. For the _GW_ 0
calculation, seven steps were necessary for convergence, however, all results
detailed in this study were taken after ten steps. Since the orbitals are not
updated in the ev _GW_ 0 self-consistency scheme, the interpolation of bands
between $\Gamma$ and A in the _GW_ output should match that seen in standard
DFT even as the eigenvalues are updated. To visualize the resulting band
behavior, eigenvalues were matched to the ordering of states in DFT outputs.
The interpolation of bands between $\Gamma$ and A was preserved from the
unrelaxed PBE band structure, yet it was scaled such that the endpoint
energies (at $\Gamma$ and A) correspond to the updated eigenvalues from the
_GW_ calculations. Further details are provided in the Appendix.

##  III Results

Figure 2: PBE + SOC + D3 electronic structure of CsV3Sb5 along high-symmetry
_k_ -point path. A previously unidentified crossing along $\Gamma$-A
representing the _k z_ band dispersion is confirmed by irrep analysis.

###  III.1 General features of the band structure

The electronic structure of CsV3Sb5 hosts a variety of features associated
with unique physical phenomena. Many computational efforts have focused on the
Dirac crossings [21] and van Hove singularities [22] characteristic of the
kagome electronic structure. Specifically, there is significant interest in
the near-Fermi-level saddle points in the band structure that correspond to
van Hove singularities or divergences in the density of states (DOS). Due to
the high DOS, interactions can become pronounced when a saddle point becomes
populated and can nest across the Fermi surface. In CsV3Sb5, multiple saddle
points residing at the M-point give rise to competing instabilities depending
on band filling.[23, 22]

Yet, while the intra-layer interactions draw focus, the interlayer electronic
behavior is understudied.[24] Among kagome materials, the 135 family is highly
two-dimensional as seen in the stacking of the kagome layers (Figure 1). While
these materials are not considered traditional 2D materials given the presence
of _A_ -site ions between the layers, previous electronic structure work
concluded that van der Waals interactions were necessary to reproduce the
experimental _c_ -axis lattice parameter during structural relaxation. These
interactions are most relevant along the interlayer or _k z_ direction
represented along the $\Gamma$-A and K-H paths within the Brillouin zone. The
near-Fermi level bands along these paths are dominated by Sb pz-orbitals
unlike the V character bands found elsewhere in the BZ.[25]

###  III.2 Band crossing along $\Gamma$-A

When comparing published band structures for this material, there is a large
degree of variability with respect to the band behavior between $\Gamma$ and
A, specifically around 0.75 eV below the Fermi energy. These band structures
can be categorized into three main groups depending on the interaction of the
two bands in this energy range: 1) fully gapped bands [23, 26, 27], 2)
seemingly touching bands with a minor gap [12, 28], and 3) fully crossed bands
with crossings between $\Gamma$ and A and between A and L [29, 30]. We
similarly obtain these variations of electronic structures depending on the
parameters included in the calculation and the methods employed for plotting.
To understand the nature of and potentially validate these crossings, the
symmetry or the irreducible representations (irreps) of the bands at specific
_k_ -points needs to be analyzed. When bands possess the same symmetry or
irrep, the associated electronic states mix or hybridize and form a gap. Only
when bands have different irreps can a symmetry-allowed crossing occur.[31]

Irreps additionally change with the inclusion of SOC since a spinor
representation (as opposed to a vector representation) is required to capture
the spin symmetry. Within the Brillouin zone of CsV3Sb5, most of the paths
connecting _k_ -points have _C_ 2v symmetry which has only one spinor irrep.
Since bands cannot have differing irreps along these paths, there will not be
any symmetry-allowed crossings in these areas. However, the paths along
$\Gamma$-A and K-H have _C_ 6v symmetry with multiple spinor irreps allowing
for crossings. It is important to note that “symmetry-allowed crossing” does
not refer to symmetry-_enforced_ crossings arising from nonsymmorphic
symmetries. While the crossings of interest are dictated or “enforced” by the
differing symmetries of the bands, we avoid this terminology since it already
has a specific definition within topological band theory.[32, 33]

When band symmetry is included in the previously optimized calculations for
CsV3Sb5 (relaxed with PBE + D3 followed by PBE + D3 + SOC ground state
calculation in [12]), multiple band crossings are observed along $\Gamma$-A
and K-H. Since the band behavior between $\Gamma$ and A conflict with previous
reports, those crossings are the focus of this study. Based on symmetry
considerations, the two bands around 0.75 eV below the Fermi energy cross
between $\Gamma$ and A but are gapped beyond A since they possess the same
spinor irrep beyond that _k_ -point. Although these bands are significantly
below the Fermi level, they contribute to the $\mathbb{Z}_{2}$
characterization of this material and they provide insight into the interlayer
electronic behavior.

To identify the origin of this crossing, a series of calculations were
performed to track the evolution of the band structure and the relaxed
physical structure as summarized in Figure 3. From here on, when referring to
electronic structure calculations unless indicated otherwise, this includes a
structural relaxation, a self-consistent calculation, followed by band
structure and density of states (DOS) calculations in which the Fermi level
for the band structure is taken from the DOS calculation given the more
accurate sampling of _k_ -points (for more detailed discussion see methods
section). Each step in this procedure includes all listed parameters.

When the electronic structure is calculated without the D3 corrections, no
band crossing between $\Gamma$ and A is observed in the energy range of
interest. Specifically, the PBE band structure shows clearly gapped bands and
the PBE + SOC band structure shows a small gap consistent with irrep analysis.
These results indicate that some contribution from the D3 van der Waals
correction is responsible for the crossing. To deconvolute whether this
contribution is structural i.e. from the relaxation or electronic, the same
calculations were performed using the experimental structure without any
relaxation. The structurally unrelaxed calculations show the band crossing in
the absence of D3, suggesting that the contribution from this correction is
structural. In other words, the experimental geometry exhibits an electronic
structure with the crossing along $\Gamma$-A and the D3 correction maintains
this crossing since the relaxed structure most closely matches experimental
lattice parameters.

###  III.3 Comparison of band disperson corrections

It is well documented that DFT methods over-delocalize electrons leading to
overly disperse bands.[34, 35] These errors are potentially most pronounced
for band dispersions along the out-of-plane direction within quasi-two-
dimensional materials like the _A_ V3Sb5 kagomes compounds. Within the
Brillouin zone (Figure 2), this direction is represented by the $\Gamma$-A,
K-H, and M-L paths meaning that these paths are more likely affected by over-
dispersion. Since the crossing of interest is along $\Gamma$-A, any correction
of the band dispersion (i.e. decreasing the bandwidth) could open a gap
between the two bands.

Figure 3: Comparison of computationally-relaxed to experimental lattice
parameters. All calculations include SOC. Low deviation from experimental
values coincides with the presence of a symmetry-allowed band crossing along
$\Gamma$-A.

To probe the conditions under which a gap might open between the bands and,
thereby assess the robustness of this crossing, we compare a variety of
commonly applied DFT approaches that address anisotropic nonlocal exchange and
correlation at various levels. Again, all parameters or techniques were
benchmarked against experimental or unrelaxed results to parse the electronic
and structural contributions. Within electronic structure theory, there are a
variety of methods to deal with localized electrons including: 1) adding a
Hubbard _U_ correction to specific orbitals, 2) adding an energetic correction
term, 3) switching to a different exchange-correlation functional, and 4)
turning to more advanced methods that more accurately treat band dispersions.
While the following results do not constitute an exhaustive selection of all
methods impacting band dispersion, a representative set of calculations within
all of these categories were performed to provide generalized results as
summarized in Figure 3. Methods were selected based on their computational
accessibility and the specificity of the correction to the interlayer band
dispersion.

Hubbard _U_ corrections are commonly used because they provide a targeted
method for dealing with electron over-delocalization usually in d- and
f-orbitals in a tunable and computationally inexpensive way.[36, 37, 38, 39]
In the context of CsV3Sb5, primarily Sb p-orbitals contribute to the near _E
F_ bands at the Brillouin zone center. While p-orbitals typically do not
suffer from the same over-delocalization issues in DFT, adding a variable _U_
to these orbitals is a valuable exercise for tracking the behavior of the
crossing as a function of the bandwidths of the involved bands. As the value
of the _U_ is increased, the deviations along the _a_ \- and _c_ -axes become
progressively more positive signifying an expansion of the unit cell during
relaxation. In the presence of the D3 van der Waals correction, calculations
with a Hubbard _U_ from 1 eV to 6 eV exhibit a crossing. Without the D3
correction, these calculations show fully gapped bands with high lattice
parameter deviations along the _c_ -axis.

Since CsV3Sb5 is considered a quasi-two-dimensional material given the
separation of the kagome layers, van der Waals corrections beyond the D3 term
could provide more accurate band dispersions corresponding to localized
electrons along the _k z_ direction. These van der Waals terms are total
energy corrections that are added to the calculated Kohn-Sham DFT (KS-DFT)
energy and, therefore, do not impact the fundamental electronic structure. The
computational efficiency of adding these corrections and the relative success
in matching experimental lattice parameters make these terms a popular method
for addressing van der Waals systems.[40, 41] However, moving to less
empirical van der Waals corrections like the Tkatchenko-Sheffler (TS)
correction [42] does not yield an improvement in the relaxed lattice
parameters. In fact, the _c_ -axis lattice parameter is significantly shorter
compared to the experimental value.

Changing the exchange-correlation functional can be a more rigorous method for
incorporating dispersion changes within the electronic structure calculation.
For example, nonlocal van der Waals density functionals can approximate
dispersion interactions by including a nonlocal correlation functional.[43,
44, 45, 46] However, for CsV3Sb5 the selected van der Waals functionals (vdW-
DF, optB86b-vdW, SCAN + rVV10) do not perform better than PBE with the D3
total energy correction in reproducing the experimental lattice parameters.
Yet, using PBEsol [47], a modification of the PBE functional optimized for
extended solids, yielded a relaxed structure with lattice parameters within 2%
deviation from the experimental values without the addition of any explicit
van der Waals correction.

All of these corrections and methods were additionally compared to unrelaxed
calculations (not shown in Figure 3) to deconvolute the structural and
electronic contributions. Under all tested computational parameters, a
crossing is observed between the Sb bands along $\Gamma$-A, pointing to the
importance of the structural relaxation in controlling the presence of this
crossing.

Figure 4: Kohn-Sham (KS) DFT and _GW_ interlayer electronic structure
behavior. (a) Comparison of KS, single-shot _G_ 0 _W_ 0, and partially self-
consistent ev _G_ _W_ 0 band structures along $\Gamma$-A shows consistent band
dispersions yet significant energetic shifts with the opening of the crossing
between Sb bands in _GW_. (b) Band-decomposed charge densities corresponding
to bands in (a) ordered by degree of interlayer character. i) and ii)
represent out-of-kagome-plane and in-kagome-plane Sb pz charge density,
respectively, while iii) and iv) represent V d charge density with varying
interlayer intensity.

Moving outside of the framework of DFT or coupling DFT results with additional
methods can also provide a more accurate picture of the electronic structure.
For example, using the _GW_ approximation for the electronic self-energy can
appropriately narrow bands by more accurately capturing electron exchange and
correlation, however, at higher computational cost. The _GW_ approximation
relies on a quasiparticle (QP) picture of the system in which electrons are
“dressed” with positively-charged polarization clouds due to the Coulomb
repulsion of other electrons. These quasiparticles, which can be described
using single-particle Green’s functions (_G_), interact weakly via a
dynamically screened Coulomb potential (_W_). Typically, _GW_ calculations
provide good estimates of the band energies and dispersions when comparing
electronic structures from theory and experiment. Since _GW_ calculations
capture correlations stemming from dynamical and non-local density
fluctuations, these calculations are well-suited for systems in which van der
Waals forces play a role. For CsV3Sb5, these types of interactions likely
dominate the interlayer electronic behavior, meaning that the _GW_
approximation provides inherent advantages over static mean field methods like
KS-DFT.

_GW_ calculations were performed on the experimental structure without SOC and
compared to a range of DFT results using the LDA, PBE, SCAN, and nonlocal van
der Waals density functionals (vdW-DF, optB86b, SCAN + rVV10).[43, 44, 45, 46,
48] Single-shot _G_ 0 _W_ 0 and partially self-consistent eigenvalue _GW_ 0
yield similar band behavior as shown in Figure 4, indicating that single shot
calculations will give a relatively consistent picture of the band structure.
Based on the partially self-consistent eigenvalue _GW_ 0 results, the band
dispersions are fairly consistent with a maximum of 12% deviation in the
dispersions of the Sb bands. A similar investigation of V bands around 1.7 eV
below the Fermi level yield a significantly higher deviation of up to 48%.

While the dispersions of the revelant bands remain relatively consistent when
moving to _GW_ calculations, all bands experience significant energy shifts
and re-ordering especially below the Fermi level as indicated in Figure 4. In
particular, the two Sb bands shift downward by more than 1 eV, combining these
bands with lower-lying less-dispersive V bands that experience less than 1 eV
shifts downward. This shift induces a reordering of bands as some V bands now
lie between the two Sb bands. Quantifying the quasiparticle energy shifts at
$\Gamma$ and A shows larger shifts of around 1.6 eV and 1.3 eV for pz bands
associated with Sb outside the kagome layer and within the kagome layer,
respectively and smaller shifts of 0.80 eV and 0.65 eV for V dz2 and V
dxy/x2$-$y2 bands, respectively. Orbital contributions are attributed based on
band-decomposed charge densities depicted in Figure 4b and published
electronic structures with orbital projections.[25, 24] The relative energy
shifts of bands with varying degrees of interlayer character ultimately
results in bands that are fully gapped between $\Gamma$ and A.

##  IV Discussion

Within standard DFT, various band dispersion correction terms and methods lead
to a similar conclusion: when the structure matches or is comparable to the
experimental structure, the associated band structure will feature a symmetry-
allowed crossing along $\Gamma$-A. This conclusion is evident in calculations
where structural relaxations are performed; when there is a small deviation
between the relaxed and experimental lattice parameters, a crossing is
observed. Unrelaxed calculations reinforce this result by showing a crossing
that is robust to a variety of electronic corrections and terms. Additionally,
the dispersions of the Sb bands that introduce the crossing are maintained
across a variety of functionals and in _GW_ calculations. Yet, in moving to
quasiparticle energies, significant energetic shifts are observed. These
shifts are more pronounced for bands with higher degrees of interlayer orbital
character. For example, the purely intra-layer V dxy/x2$-$y2 bands experience
the smallest shift in quasiparticle versus Kohn-Sham energies, whereas the
out-of-plane Sb pz band experiences the largest shift. Since the in-plane Sb
band does not shift as much as the out-of-plane Sb, a gap is expected to open
between these bands (note the separation between bands i and ii in Figure 4).

Based on DFT calculations, the presence of this symmetry-allowed crossing and
its robustness with a variety of band dispersion corrections suggests a
reevaluation of the electronic structure of CsV3Sb5 specifically along the _k
z_ direction. _GW_ calculations provide further insights into the interlayer
electronic structure, revealing the absence of this crossing due to
significant energetic shifts of bands below the Fermi level. Additionally,
these results show interlayer interactions mediated through the V and out-of-
plane Sb atoms indicated by the almost degenerate i and iii bands shown in
Figure 4. The hybridization of these states suggests more interlayer phenomena
at play than indicated in previous electronic structures.

Given the two-dimensionality of this kagome system, it is particularly
important to study the interlayer electronic behavior as has been demonstrated
computationally and experimentally with effects due to pressure and chemical
doping and charge-density wave ordering.[26, 29, 30, 49] All of these
phenomena show dramatic changes along the _c_ \- or _k z_-axes pointing to the
critical role of understanding interlayer interactions. Rigorous electronic
structure calculations capturing the _k z_ dispersion are further necessary
since experimental verification of theses bands via ARPES remains challenging.
The highly two-dimensional single crystals make it difficult to obtain a clean
surface required for the technique along the relevant axis.

Additionally, the $\Gamma$-A crossing (or lack thereof) has greater
implications in the topological characterization of this system via
$\mathbb{Z}_{2}$ invariant calculations. A crossing implies the absence of a
continuous gap necessary for Fu-Kane $\mathbb{Z}_{2}$ calculations across a
band manifold. In the case of CsV3Sb5, there is an additional continuous gap
throughout the Brillouin zone that interacts with the Fermi level that yields
a non-trivial $\mathbb{Z}_{2}$ invariant. Since the Fu-Kane method involves
taking the product of parity eigenvalues, the ordering of bands below the
continuous gap does not affect the overall result meaning that the
$\mathbb{Z}_{2}$ calculation for this particular gap is unchanged from
previous reports even with the band re-ordering along $\Gamma$-A seen in the
_GW_ results.[12] Ultimately, the presence or absence of this symmetry-allowed
crossing does not change the non-trivial $\mathbb{Z}_{2}$ classification even
though it potentially reduces the number of continuous gaps. However, in other
systems that do not have additional continuous gaps, subtle crossings that are
not identified until irreps are considered could have more severe consequences
on topological invariant calculations.

This study motivates revisiting the protocol for defining $\mathbb{Z}_{2}$
invariants in metallic systems. Since the Fu-Kane method was originally
intended for insulators, the presence of a gap and, subsequently, a separable
band manifold is guaranteed by the system. [7] However, when moving to metals,
this condition of separable bands becomes a prerequisite step before
continuing with the calculation. From this study of CsV3Sb5, we find that
visual identification of a near-Fermi level gap that extends throughout the
entire Brillouin zone is not always sufficient fulfillment of this
prerequisite. Our results highlight the importance of confirming a gap via
irreps prior to performing $\mathbb{Z}_{2}$ analysis in metals.

Finally, broader lessons regarding the interpretation of band structures can
be learned from this study of CsV3Sb5. For example, recognizing that most band
structure outputs do not include band symmetries and can, therefore,
misidentify “continuous” bands advocates the integration of irrep analysis
when depicting and understanding band structures. This lack of proper band
tracking is reinforced if other computational parameters like _k_ -point
density along the Brillouin zone path are not high enough to resolve small
features. These issues are crucial when attempting to calculate topological
invariants based on bulk band structures but are also widely applicable to any
system in which an accurate picture of the electronic structure is necessary.

##  V Conclusion

Using first-principles calculations, we have identified new band behavior
along the interlayer direction in CsV3Sb5. Given the two-dimensional nature of
this system, capturing the electronic structure along the _c_ \- or _k z_-
direction is experimentally and computationally difficult. In this study, we
survey a variety of computational methods to address nonlocal exchange and
correlation at different scales. Density functional theory calculations show a
symmetry-allowed crossing along the _k z_-dispersion that arises due to the
structural contribution of dispersion corrections and methods employed during
relaxation to more closely match experimental lattice parameters. This
crossing is revealed upon irrep analysis of the bands. Yet, in turning towards
dynamical nonlocal correlation via _GW_ calculations, this crossing is absent
due to relative energy shifts of bands with varying interlayer character. The
resulting _k z_ electronic structure points to increased interactions
occurring out of the kagome plane and serves as a baseline band structure for
further studies of inter-kagome-layer phenomena.

##  VI Acknowledgments

We thank J. W. Harter, L. Balents, M. G. Vergniory, and I. Robredo for helpful
discussions and feedback. This work was supported by the National Science
Foundation (NSF) through Enabling Quantum Leap: Convergent Accelerated
Discovery Foundries for Quantum Materials Science, Engineering and Information
(Q-AMASE-i): Quantum Foundry at UC Santa Barbara (DMR 1906325). Additionally,
we thank NSF CAREER award DMR 1945098. We also acknowledge the use of
computational facilities administered by the Center for Scientific Computing
at UC Santa Barbara supported by NSF CNS 1725797 and NSF DMR 2308708.

*

##  Appendix A Plotting electronic structures

Since typical band structure data include the _k_ -points and energies of each
band, portions of this data need to be rearranged to reflect the true band
identity confirmed by irreps. For example, if two bands cross between _k_
-point 1 and _k_ -point 2, the output may incorrectly identify all lower
energies as a continuous band and similarly all higher energies as another
continuous band. In other words, the lower energy band before the crossing
will be spuriously connected to the lower energy band after the crossing.
However, looking at the irreps at _k_ -point 1 and _k_ -point 2 will indicate
there is a crossing because the irreps switch between the two _k_ -points.
This weakness in the output can be corrected by identifying the point at which
the crossing occurs and reconnecting the bands to match the irrep results.
This practice can be seen in Figure 2 where the bands between $\Gamma$ and A
have been reconnected and colored to match the irreps.

The _GW_ band structures in Figure 4 were constructed based on the unrelaxed
PBE band structure. Once the bands were reconnected based on irreps, each band
was rescaled and repositioned according to the renormalization computed in the
_G_ 0 _W_ 0 and _GW_ 0 calculations at $\Gamma$ and at A. As such, bands from
the _GW_ output may suffer from the same misidentification across crossings,
however since only the energies and not the orbitals are updated in the
employed partially self-consistent scheme, the _GW_ energies can be properly
matched based on the ordering of bands from the original PBE band structure.

## References

  * Xiao and Yan [2021] J. Xiao and B. Yan, First-principles calculations for topological quantum materials, Nat. Rev. Phys. 3, 283 (2021).
  * Bradlyn  _et al._ [2017] B. Bradlyn, L. Elcoro, J. Cano, M. G. Vergniory, Z. Wang, C. Felser, M. I. Aroyo, and B. A. Bernevig, Topological quantum chemistry, Nature 547, 298 (2017).
  * Bansil  _et al._ [2016] A. Bansil, H. Lin, and T. Das, Colloquium: Topological band theory, Rev. Mod. Phys. 88 (2016).
  * Kane and Mele [2005] C. L. Kane and E. J. Mele, Z2 Topological Order and the Quantum Spin Hall Effect, Phys. Rev. Lett. 95 (2005).
  * Fu  _et al._ [2007] L. Fu, C. L. Kane, and E. J. Mele, Topological Insulators in Three Dimensions, Phys. Rev. Lett. 98 (2007).
  * Moore and Balents [2007] J. E. Moore and L. Balents, Topological invariants of time-reversal-invariant band structures, Phys. Rev. B 75 (2007).
  * Fu and Kane [2007] L. Fu and C. L. Kane, Topological insulators with inversion symmetry, Phys. Rev. B 76, 045302 (2007).
  * Schoop  _et al._ [2015] L. M. Schoop, L. S. Xie, R. Chen, Q. D. Gibson, S. H. Lapidus, I. Kimchi, M. Hirschberger, N. Haldolaarachchige, M. N. Ali, C. A. Belvin, T. Liang, J. B. Neaton, N. P. Ong, A. Vishwanath, and R. J. Cava, Dirac metal to topological metal transition at a structural phase change in Au2Pb and prediction of Z2 topology for the superconductor, Phys. Rev. B 91, 214517 (2015).
  * Qian  _et al._ [2022] J. Qian, Z. Shen, X. Wei, and W. Li, Z2 nontrivial topology of rare-earth binary oxide superconductor LaO, Phys. Rev. B 105 (2022).
  * Nayak  _et al._ [2017] J. Nayak, S.-C. Wu, N. Kumar, C. Shekhar, S. Singh, J. Fink, E. E. D. Rienks, G. H. Fecher, S. S. P. Parkin, B. Yan, and C. Felser, Multiple dirac cones at the surface of the topological metal LaBi, Nat. Commun. 8 (2017).
  * Ortiz  _et al._ [2021] B. R. Ortiz, P. M. Sarte, E. M. Kenney, M. J. Graf, S. M. L. Teicher, R. Seshadri, and S. D. Wilson, Superconductivity in the Z2 kagome metal KV3Sb5, Phys. Rev. Mater. 5, 034801 (2021).
  * Ortiz  _et al._ [2020] B. R. Ortiz, S. M. Teicher, Y. Hu, J. L. Zuo, P. M. Sarte, E. C. Schueller, A. M. Abeykoon, M. J. Krogstad, S. Rosenkranz, R. Osborn, R. Seshadri, L. Balents, J. He, and S. D. Wilson, CsV3Sb5 : A Z2 Topological Kagome Metal with a Superconducting Ground State, Phys. Rev. Lett. 125, 247002 (2020).
  * Yan  _et al._ [2011] S. Yan, D. A. Huse, and S. R. White, Spin-Liquid Ground State of the S = 1/2 Kagome Heisenberg Antiferromagnet, Science 332, 1173 (2011).
  * Ko  _et al._ [2009] W.-H. Ko, P. A. Lee, and X.-G. Wen, Doped kagome system as exotic superconductor, Phys. Rev. B 79 (2009).
  * Guo and Franz [2009] H.-M. Guo and M. Franz, Topological insulator on the kagome lattice, Phys. Rev. B 80 (2009).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized Gradient Approximation Made Simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Grimme  _et al._ [2010] S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu, J. Chem. Phys. 132, 154104 (2010).
  * Grimme  _et al._ [2011] S. Grimme, S. Ehrlich, and L. Goerigk, Effect of the damping function in dispersion corrected density functional theory, J. Comput. Chem. 32, 1456 (2011).
  * Curtarolo  _et al._ [2012] S. Curtarolo, W. Setyawan, G. L. Hart, M. Jahnatek, R. V. Chepulskii, R. H. Taylor, S. Wang, J. Xue, K. Yang, O. Levy, M. J. Mehl, H. T. Stokes, D. O. Demchenko, and D. Morgan, AFLOW: An automatic framework for high-throughput materials discovery, Comput. Mater. Sci. 58, 218 (2012).
  * Gao  _et al._ [2021] J. Gao, Q. Wu, C. Persson, and Z. Wang, Irvsp: To obtain irreducible representations of electronic states in the VASP, Comput. Phys. Commun. 261, 107760 (2021).
  * Hao  _et al._ [2022] Z. Hao, Y. Cai, Y. Liu, Y. Wang, X. Sui, X.-M. Ma, Z. Shen, Z. Jiang, Y. Yang, W. Liu, Q. Jiang, Z. Liu, M. Ye, D. Shen, Y. Liu, S. Cui, J. Chen, L. Wang, C. Liu, J. Lin, J. Wang, B. Huang, J.-W. Mei, and C. Chen, Dirac nodal lines and nodal loops in the topological kagome superconductor CsV3Sb5, Phys. Rev. B 106 (2022).
  * Hu  _et al._ [2022] Y. Hu, S. M. Teicher, B. R. Ortiz, Y. Luo, S. Peng, L. Huai, J. Ma, N. C. Plumb, S. D. Wilson, J. He, and M. Shi, Topological surface states and flat bands in the kagome superconductor CsV3Sb5, Sci. Bull. 67, 495 (2022).
  * Kang  _et al._ [2022] M. Kang, S. Fang, J.-K. Kim, B. R. Ortiz, S. H. Ryu, J. Kim, J. Yoo, G. Sangiovanni, D. D. Sante, B.-G. Park, C. Jozwiak, A. Bostwick, E. Rotenberg, E. Kaxiras, S. D. Wilson, J.-H. Park, and R. Comin, Twofold van Hove singularity and origin of charge order in topological kagome superconductor CsV3Sb5, Nat. Phys. 18, 301 (2022).
  * Tsirlin  _et al._ [2022] A. Tsirlin, P. Fertey, B. R. Ortiz, B. Klis, V. Merkl, M. Dressel, S. Wilson, and E. Uykur, Role of Sb in the superconducting kagome metal CsV3Sb5 revealed by its anisotropic compression, SciPost Phys. 12 (2022).
  * Li  _et al._ [2022] C. Li, X. Wu, H. Liu, C. Polley, Q. Guo, Y. Wang, X. Han, M. Dendzik, M. H. Berntsen, B. Thiagarajan, Y. Shi, A. P. Schnyder, and O. Tjernberg, Coexistence of two intertwined charge density waves in a kagome system, Phys. Rev. Res. 4 (2022).
  * Oey  _et al._ [2022] Y. M. Oey, B. R. Ortiz, F. Kaboudvand, J. Frassineti, E. Garcia, R. Cong, S. Sanna, V. F. Mitrović, R. Seshadri, and S. D. Wilson, Fermi level tuning and double-dome superconductivity in the kagome metal CsV3Sb5-xSnx, Phys. Rev. Mater. 6 (2022).
  * Huang  _et al._ [2022] X. Huang, C. Guo, C. Putzke, M. Gutierrez-Amigo, Y. Sun, M. G. Vergniory, I. Errea, D. Chen, C. Felser, and P. J. W. Moll, Three-dimensional Fermi surfaces from charge order in layered CsV3Sb5, Phys. Rev. B 106 (2022).
  * Uykur  _et al._ [2021] E. Uykur, B. R. Ortiz, O. Iakutkina, M. Wenzel, S. D. Wilson, M. Dressel, and A. A. Tsirlin, Low-energy optical properties of the nonmagnetic kagome metal CsV3Sb5, Phys. Rev. B 104 (2021).
  * Chen  _et al._ [2021] K. Chen, N. Wang, Q. Yin, Y. Gu, K. Jiang, Z. Tu, C. Gong, Y. Uwatoko, J. Sun, H. Lei, J. Hu, and J.-G. Cheng, Double Superconducting Dome and Triple Enhancement of Tc in the Kagome Superconductor CsV3Sb5 under High Pressure, Phys. Rev. Lett. 126 (2021).
  * Yu  _et al._ [2022] F. Yu, X. Zhu, X. Wen, Z. Gui, Z. Li, Y. Han, T. Wu, Z. Wang, Z. Xiang, Z. Qiao, J. Ying, and X. Chen, Pressure-Induced Dimensional Crossover in a Kagome Superconductor, Phys. Rev. Lett. 128 (2022).
  * Dresselhaus  _et al._ [2008] M. S. Dresselhaus, G. Dresselhaus, and A. Jorio,  _Group Theory: Application to the Physics of Condensed Matter_ (Springer-Verlag Berlin Heidelberg, 2008) pp. 305–336.
  * Zhao and Schnyder [2016] Y. X. Zhao and A. P. Schnyder, Nonsymmorphic symmetry-required band crossings in topological semimetals, Phys. Rev. B 94 (2016).
  * Hirschmann  _et al._ [2021] M. M. Hirschmann, A. Leonhardt, B. Kilic, D. H. Fabini, and A. P. Schnyder, Symmetry-enforced band crossings in tetragonal materials: Dirac and Weyl degeneracies on points, lines, and planes, Phys. Rev. Mater. 5 (2021).
  * Sholl and Steckel [2009] D. S. Sholl and J. A. Steckel,  _Density Functional Theory: A Practical Introduction_ (Wiley, 2009) pp. 209–233.
  * Martin [2004] R. R. Martin,  _Electronic Structure: Best Theory and Practical Methods_ (Cambridge University Press, 2004) pp. 152–171.
  * Anisimov  _et al._ [1991] V. I. Anisimov, J. Zaanen, and O. K. Andersen, Band theory and Mott insulators: Hubbard U instead of Stoner I, Phys. Rev. B 44, 943 (1991).
  * Cococcioni and de Gironcoli [2005] M. Cococcioni and S. de Gironcoli, Linear response approach to the calculation of the effective interaction parameters in the LDA + U method, Phys. Rev. B 71 (2005).
  * Kulik [2015] H. J. Kulik, Perspective: Treating electron over-delocalization with the DFT+U method, J. Chem. Phys. 142, 240901 (2015).
  * Dudarev  _et al._ [1998] S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton, Electron-energy-loss spectra and the structural stability of nickel oxide: An LSDA+U study, Phys. Rev. B 57, 1505 (1998).
  * Reckien  _et al._ [2012] W. Reckien, F. Janetzko, M. F. Peintinger, and T. Bredow, Implementation of empirical dispersion corrections to density functional theory for periodic systems, J. Comput. Chem. 33, 2023 (2012).
  * Grimme  _et al._ [2016] S. Grimme, A. Hansen, J. G. Brandenburg, and C. Bannwarth, Dispersion-Corrected Mean-Field Electronic Structure Methods, Chem. Rev. 116, 5105 (2016).
  * Tkatchenko and Scheffler [2009] A. Tkatchenko and M. Scheffler, Accurate Molecular Van Der Waals Interactions from Ground-State Electron Density and Free-Atom Reference Data, Phys. Rev. Lett. 102 (2009).
  * Román-Pérez and Soler [2009] G. Román-Pérez and J. M. Soler, Efficient Implementation of a van der Waals Density Functional: Application to Double-Wall Carbon Nanotubes, Phys. Rev. Lett. 103 (2009).
  * Dion  _et al._ [2004] M. Dion, H. Rydberg, E. Schröder, D. C. Langreth, and B. I. Lundqvist, Van der Waals Density Functional for General Geometries, Phys. Rev. Lett. 92 (2004).
  * Peng  _et al._ [2016] H. Peng, Z.-H. Yang, J. P. Perdew, and J. Sun, Versatile van der Waals Density Functional Based on a Meta-Generalized Gradient Approximation, Phys. Rev. X 6 (2016).
  * Klimeš  _et al._ [2011] J. Klimeš, D. R. Bowler, and A. Michaelides, Van der Waals density functionals applied to solids, Phys. Rev. B 83 (2011).
  * Perdew  _et al._ [2008] J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke, Restoring the Density-Gradient Expansion for Exchange in Solids and Surfaces, Phys. Rev. Lett. 100 (2008).
  * Sun  _et al._ [2015] J. Sun, A. Ruzsinszky, and J. Perdew, Strongly Constrained and Appropriately Normed Semilocal Density Functional, Phys. Rev. Lett. 115 (2015).
  * Kautzsch  _et al._ [2023] L. Kautzsch, B. R. Ortiz, K. Mallayya, J. Plumb, G. Pokharel, J. P. C. Ruff, Z. Islam, E.-A. Kim, R. Seshadri, and S. D. Wilson, Structural evolution of the kagome superconductors AV3Sb5 (A = K, Rb, and Cs) through charge density wave order, Phys. Rev. Mater. 7 (2023).