# Shift of nuclear clock transition frequency in 229Th ions due to hyperfine
interaction.

V. A. Dzuba and V. V. Flambaum  School of Physics, University of New South
Wales, Sydney 2052, Australia

###### Abstract

We calculate hyperfine structure of 229Th and its ions (Th IV, Th III, Th II,
Th I) to reveal the dependence of the nuclear clock frequency on the hyperfine
interaction (hfi). We calculate first and second-order hfi shifts and
demonstrate that due to the differences in the hyperfine structure for
different ions and for the ground and isomeric nuclear states, the nuclear
frequencies in Th IV, Th III, Th II and Th I are also different. The first-
order shift of frequency is large for a particular hyperfine component, but it
vanishes after averaging over all hyperfine states. The second-order shift is
small but it does not vanish after averaging. It is three to six orders of
magnitude smaller than the shift of the nuclear frequencies due to the Coulomb
electron-nucleus interaction considered in our previous work (V. A. Dzuba and
V. V. Flambaum, arXiv:2309.11176 (2023)). However, it is six to eight orders
of magnitude larger than the projected accuracy of the nuclear clock
($10^{-19}$).

##  I Introduction

The 229Th isotope is considered for building nuclear clocks of exceptional
accuracy (see, e.g. PeikTamm ; Thm ). This is because the 229Th nucleus has a
unique feature of having very low-energy excitation connected to the nuclear
ground state by the magnetic dipole (M1) transition (see, e.g. reviews Rev ;
Rev1  and references therein). The relative uncertainty of the proposed clock
is expected to reach $10^{-19}$ Th3+ . There are strong arguments that this
nuclear clock would be very sensitive to physics beyond standard model
including space-time variation of the fundamental constants, violation of the
Lorentz invariance and Einstein equivalence principle, and search for scalar
and axion dark matter fields vara ; Lorentz ; varq ; vara1 ; vara2 ; vara3 ;
vara4 ; Arvanitaki ; Stadnik . There is good progress in recent years of
measurements of the frequency of this transition wn1 ; wn2 ; wn3 ; wn4 ; wn5 .
The latest, most precise measurements, give the value of 8.338(24) eV wn .
There are plans to use Th ions of different ionisation degree Th3+ ; Th+ ; dr
and solid-state Th nuclear clock Hudson ; ThSS ; ThSS1 . In our previous paper
previous  we demonstrated that in all these systems the frequency of the
nuclear clock will be different. This is due to the Coulomb interaction of
atomic electrons with the nucleus, leading to the significant electronic shift
of the nuclear transition frequency. The effect is similar to the field
isotopic shifts and isomeric shifts of atomic transition frequencies but here
electrons affect the nuclear transition. The shift is caused by the difference
in electronic structure between Th ions and the difference in the nuclear
radii for the ground and isomeric nuclear states.

In the present paper we consider another effect which leads to the difference
in nuclear frequency for different Th systems. It comes from the hyperfine
interaction (hfi) of atomic electrons with the nucleus. The difference in
frequencies is caused by the difference in electronic structure of Th ions and
the difference in nuclear parameters (nuclear spin $I$, nuclear magnetic
moment $\mu$ and nuclear electric quadrupole moment $Q$) for the ground and
isomeric nuclear states.

The hyperfine structure (hfs) of Th and its ions has beed studied for many
years both experimentally and theoretically Th+old ; ThQ ; Th-mu-Q ; Beloy ;
Th+hfs ; ThIIIhfs . However, the data is still incomplete. There are no
experimental data for the hfs of the ground states of Th I and Th III. To the
best of our knowledge, there are no calculations for Th II and Th I hfs. We
address these shortcoming by performing the hfs calculations for Th IV, Th
III, Th II and Th I. For the purpose of the present work we need only hfs of
the ground states of all these systems. However, we calculate hfs constants
for some excited states too for the convenience of comparing with other data
either experimental or theoretical. This gives us some understanding of the
accuracy of our calculations.

First-order hfi leads only to splitting of atomic levels. The total shift,
averaged over hfs components, is zero. In contrast, the second-order hfi does
lead to the $F$-dependent shifts of the levels which do not vanish after
averaging. Here $F$ is the value of the total atomic angular momentum, which
is the sum of electronic angular momentum $J$ and nuclear spin $I$
($\mathbf{F}=\mathbf{J}+\mathbf{I}$). In the present paper we calculate both,
the first-order hfs and the second-order hfi shift.

##  II First-order hfs shift

First-order hfs leads to splitting of the atomic levels. Total shift averaged
over hfs components with statistical weight $(2F+1)$ is zero. However, it was
suggested in Ref. Th3+  to use the stretched hyperfine states ($F=I+J$,
$F_{z}=F$ and $F_{z}=-F$) because in such states atomic wave function is a
simple product of electronic and nuclear wave functions. This absence of
entanglement between electron and nuclear variables leads to strong
suppression of systematic effects in the nuclear transition since external
fields mainly interact with electrons. For a stretch state the energy shift
caused by the magnetic dipole hyperfine interaction is

$\Delta E_{A}=AIJ,$ (1)

where $A$ is the magnetic dipole hfs constant, $I$ is the nuclear spin, $J$ is
the total electron angular momentum. The difference in shifts for the ground
and isomeric nuclear states comes from different values of nuclear magnetic
moments $\mu$ of the states. It is convenient to separate electronic and
nuclear variables by writing $A=A_{0}\mu/I$. Then the difference is given by

$\delta(\Delta E_{A})=A_{0}J\Delta\mu.$ (2)

For the electric quadrupole hyperfine shift the general formula reads

$\displaystyle\Delta
E_{B}=\frac{B}{2}\left[C(C+1)-\frac{4}{3}J(J+1)I(I+1)\right],$ $\displaystyle
C=F(F+1)-J(J+1)-I(I+1).$

To get the result for the stretched state we should substitute $F=I+J$. This
gives

$\Delta E_{B}=\frac{B}{3}IJ(2I-1)(2J-1).$ (4)

Separating electronic and nuclear variables by introducing $B=B_{0}Q$, where
$Q$ is nuclear electric quadrupole moment, we see that the difference between
ground and isomeric nuclear states comes from the different values of $Q$ and
$I$. To calculate frequency differences caused by hfi we need to know
electronic structure factors $A_{0}$, $B_{0}$ for the ground states of Th IV,
Th III, Th II, and Th I, as well as nuclear parameters $\mu$, $I$ and $Q$ for
the ground and isomeric nuclear states of 229Th. Nuclear parameters are known
Thm ; Th-mu-Q ; ThQold , while experimental data on the hfs of the ground
states are available for Th IV ThQ ; Th-mu-Q  and Th II Th+old ; Th+hfs .

Table 1: Nuclear parameters of the ground and isomeric states of 229Th taken from Refs. Thm ; Th-mu-Q ; ThQold . Nucleus | $I$ |  $\mu$ [$\mu_{N}$] |  $Q$ [barn]
---|---|---|---
229Th | 5/2 | 0.360(7)111Ref. Th-mu-Q . | 3.11(6)111Ref. Th-mu-Q .; 3.15(3)222Ref. ThQold .
229mTh | 3/2 | -0.37(6)333Ref. Thm . | 1.73(10)333Ref. Thm .

###  II.1 Hyperfine structure of Th IV

Table 2: Experimental and theoretical values of the magnetic dipole and electric quadrupole hfs constants $A$ and $B$ (in MHz) for four low states of 229Th IV. It is assumes in the calculations that $\mu=0.360\mu_{N}$, $I=5/2$ Thm , and $Q=3.11b$ Th-mu-Q . State | $A_{\rm expt}$ | $A_{\rm theor}$ | $B_{\rm expt}$ | $B_{\rm theor}$
---|---|---|---|---
| ThQ  | Th-mu-Q  | Beloy  | This | ThQ  | Th-mu-Q  | Beloy  | This
|  |  |  | work |  |  |  | work
$5f_{5/2}$ | 82.2(6) | 82.5 | 91.9 | 81.0 | 2269(6) | 2254 | 2258 | 2264
$5f_{7/2}$ | 31.4(7) | 31.0 | 36.6 | 31.0 | 2550(12) | 2515 | 2630 | 2569
$6d_{3/2}$ | 155.3(12) | 155.3 |  | 160 | 2265(9) | 2295 |  | 2276
$6d_{5/2}$ | -12.6(7) | -13.2 |  | -27 | 2694(7) | 2716 |  | 2724

Hyperfine structure of the ground state of the Th IV is known experimentally
ThQ . However, it is instructive to perform the calculations to check the
accuracy of our approach.

We start the calculations from the relativistic Hartree-Fock (RHF) procedure
for the closed-shell Rn-like core. The single-electron basis states for
valence space are calculated in the field of frozen core using the B-spline
technique B-spline . These basis states are used to calculate the single-
electron correlation operator $\hat{\Sigma}^{(2)}_{1}$ (correlation potential
CPM ). The same basis is used for the configuration interaction (CI)
calculations in next subsection. We use the second-order many-body
perturbation theory to calculate $\hat{\Sigma}^{(2)}_{1}$.

The single-electron valence states of Th IV are calculated using this
correlation potential $\hat{\Sigma}^{(2)}_{1}$,

$(\hat{H}^{\rm RHF}+\lambda\hat{\Sigma}^{(2)}_{1}-\epsilon_{v})\psi^{\rm
BO}_{v}=0.$ (5)

Here $\hat{H}^{\rm RHF}$ is the RHF Hamiltonian, index $v$ numerates valence
states, $\lambda$ is rescaling coefficient, its value is chosen to fit
experimental energies. We use $\lambda=0.75$ for $d$-states and $\lambda=0.8$
for $f$-states. The fitting with $\lambda$ imitates the effect of higher-order
correlations. The solutions of (5) are usually called Brueckner orbitals (BO).

To include the interaction of atomic electrons with the magnetic dipole and
electric quadrupole fields of the nucleus (i.e., to calculate hfs) we use the
version of the random phase approximation (RPA) developed in Ref. CPM . The
RPA equations have the form

$(\hat{H}^{\rm RHF}-\epsilon_{c})\delta\psi_{c}=-(\hat{F}+\delta V^{F}_{\rm
core})\psi_{c}.$ (6)

Here $F$ is the operator of external field (magnetic dipole or electric
quadrupole field of the nucleus), index $c$ numerates core states,
$\delta\psi_{c}$ is the correction to core orbital $\psi_{c}$ caused by
external field, $\delta V^{F}_{\rm core}$ is the correction to the self-
consistent RHF potential of the core caused by the change of all core
orbitals. After the RPA equations are solved, the energy shift for valence
state $v$ is given by

$\delta\epsilon_{v}=\langle\psi^{\rm BO}_{v}|\hat{F}+\delta V^{F}_{\rm
core}+\delta\Sigma^{(2)}_{1}|\psi^{\rm BO}_{v}\rangle.$ (7)

Here $\delta\Sigma^{(2)}_{1}$ is the change of the correlation potential due
to the change of all basis states in the external field. Corresponding term is
called structure radiation. It was not included in our previous calculations
for the electric quadrupole hfs constant $B$ vara2 . Inclusion of the
structure radiation leads to significant improvement of the results for both
hfs constants $A$ and $B$.

Energy shifts (7) are used to calculate hfs constants $A$ and $B$. The results
are presented in Table 2. Comparison with experimental data from Ref. ThQ  and
earlier calculations of Ref. Th-mu-Q ; Beloy  shows good agreement in all
cases except $A$ constant for the $6d_{5/2}$ state. Here the value of $A$ is
small and dominated by the many-body corrections which change the sign of $A$
.

For the purpose of the present work we need only hfs constants $A$ and $B$ for
the ground state. We can use experimental values for Th IV. However, the
calculations are useful since they serve as a starting point for more
complicated calculations for Th III, Th II and Th I (see next subsection)
where experimental data are incomplete.

###  II.2 Hyperfine structure of Th I, Th II and Th III.

Table 3: Experimental and theoretical values of the magnetic dipole and electric quadrupole hfs constants $A$ and $B$ (in MHz) for some states of 229Th III, 229Th II, and 229Th I. It is assumes in calculations that $\mu=0.360\mu_{N}$, $I=5/2$ Thm , and $Q=3.11b$ Th-mu-Q . Ion/ | State | Energy |  $A$ [MHz] |  $B$ [MHz]
---|---|---|---|---
Atom | Conf. | $J^{p}$ | Ref. NIST  | CI+MBPT | MCDF | This | Expt. | CI+MBPT | MCDF | This | Expt.
|  |  |  | Ref. ThIIIhfs  | Ref. ThIIIhfs  | work | Ref. Thm  | Ref. ThIIIhfs  | Ref. ThIIIhfs  | work | Ref. Thm
Th III | $5f6d$ | $4^{-}$ | 0 | 64(17) | 81(4) | 76.3 |  | 3287(630) | 3008(260) | 3828 |
| $6d^{2}$ | $2^{+}$ | 63 | 143(47) | 162(8) | 146 | 151(8) | 68(23) | 71(7) | -35 | 73(27)
| $5f^{2}$ | $4^{+}$ | 15148 | 38(3) | 72(3) | 60.2 |  | 1221(390) | 1910(200) | 2469 |
| $5f6d$ | $1^{-}$ | 20711 | 109(36) | 90(4) | 111 | 88(5) | 839(220) | 689(110) | 672 | 901(18)
| $5f^{2}$ | $4^{+}$ | 21784 | 8(36) | 26(2) | 4.9 |  | 65(21) | 39(45) | -600 |
Th II | $6d^{7}s$ | $3/2^{+}$ | 0 |  |  | -472 | -444.2(3.4)111Ref. Th+hfs . |  |  | 507 | 308(13)111Ref. Th+hfs .
|  |  |  |  |  |  | -444.2(1.9)222Ref. Th+old . |  |  |  | 303(6)222Ref. Th+old .
Th I | $6d^{2}7s^{2}$ | $2^{+}$ | 0 |  |  | 41 |  |  |  | 320 |

To calculate hfs of atoms with two, three and four valence electrons we use
the all-orders SD+CI SD+CI  ( single-double coupled cluster method combined
with the configuration interaction technique). The SD method produces the all-
orders single electron correlation operator $\Sigma^{\infty}_{1}$, similar to
the second-order operator $\Sigma^{(2)}_{1}$ used in the previous section. It
also produces the all-orders two-electron correlation operator
$\Sigma^{\infty}_{2}$. The use of the all-orders correlation operators usually
leads to more accurate results. The effective CI Hamiltonian has the form

$\hat{H}^{\rm CI}=\sum_{i=1}^{N_{v}}\left(\hat{H}^{\rm
RHF}+\hat{\Sigma}^{\infty}_{1}\right)_{i}+\sum_{i<j}\left(\frac{e^{2}}{r_{ij}}+\hat{\Sigma}^{\infty}_{2ij}\right).$
(8)

Here summation goes over valence electrons, $N_{v}$ is the number of valence
electrons (2,3 or 4). The same RPA technique as in previous section is used to
include external field. The energy shift for valence state $v$ is given by

$\delta E_{v}=\langle v|\sum_{i}^{N_{v}}(\hat{F}+\delta V^{F}_{\rm
core})_{i}|v\rangle.$ (9)

The results of calculations are presented in Table 3 and compared with
available experimental data and previous calculations. Comparison shows that
the uncertainty of the present calculations for the magnetic hfs constant $A$
is probably within 10%. The uncertainty for $B$ is larger, it is about 30% for
large $B$ constants and even higher for small constants. There are no
experimental data for the hfs of the ground states of Th III and Th I. There
are also no other calculations for Th I.

###  II.3 Frequency shifts

Table 4: Frequency shifts in the stretched states caused by the first order hfi. The numbers in square brackets stand for powers of ten. Ion/ | $J_{0}$ | $\delta(\Delta E_{A})$ | $\delta(\Delta E_{B})$ | $\Delta\omega_{N}/\omega_{N}$
---|---|---|---|---
Atom |  | [MHz] | [MHz] |
Th IV | 5/2 | -1.049[+3] | -4.667[+1] | -5.436[-7]
Th III | 4 | -1.558[+3] | -1.307[+2] | -8.377[-7]
Th II | 3/2 | 3.399[+3] | -1.400[+1] | 1.679[-6]
Th I | 2 | -4.145[+2] | -2.800[+1] | -2.195[-7]

To calculate the first-order hfs shift between ground and excited nuclear
states in the stretched atomic state we use formulae (2) and (4), nuclear
parameters from Table 1, and the hyperfine structure constants $A$ and $B$ for
the ground states of Th IV, Th III, Th II, and Th I from Tables 2 and 3. For
Th IV and Th I ions we use experimental values. For Th III ion we use weighted
theoretical values from Table 3. We assume theoretical uncertainty of 10% for
the hfs constant $A$ and 30% for $B$ in our results to use in this weighting
procedure. The results for the first order energy shifts in the stretched
states are presented in Table 4.

##  III Second-oder hfs shift.

Table 5: Second-order hfs shift (MHz). See Eq. (11) for definitions of $\Delta E_{a}$, $\Delta E_{ab}$, and $\Delta E_{b}$. $\langle\Delta E\rangle$ is the total shift averaged over all hfs components. The numbers in square brackets stand for powers of ten. | Ground nuclear state, $I=5/2$. | Isomeric nuclear state, $I=3/2$.
---|---|---
$F$ | $\Delta E_{a}$ | $\Delta E_{ab}$ | $\Delta E_{b}$ | Sum | $\Delta E_{a}$ | $\Delta E_{ab}$ | $\Delta E_{b}$ | Sum
Th IV, $5f\ J_{0}=5/2$
1.0 | 2.87[-5] | -2.84[-4] | 7.01[-4] | 4.46[-4] |  |  |  |
2.0 | 7.76[-5] | -5.96[-4] | 1.15[-3] | 6.27[-4] | 1.04[-4] | 5.49[-4] | 7.24[-4] | 1.38[-3]
3.0 | 1.29[-4] | -5.68[-4] | 6.24[-4] | 1.85[-4] | 1.95[-4] | 3.43[-4] | 1.51[-4] | 6.88[-4]
4.0 | 1.58[-4] |  |  | 1.58[-4] | 1.95[-4] | -5.71[-4] | 4.19[-4] | 4.23[-5]
5.0 | 1.29[-4] | 7.10[-4] | 9.74[-4] | 1.81[-3] |  |  |  |
| $\langle\Delta E\rangle$ | 7.54[-4] | $\langle\Delta E\rangle$ | 5.03[-4]
Th III, $5f6d\ J_{0}=4$
1.5 | 1.24[-3] | 1.03[-3] | 9.77[-4] | 3.25[-3] |  |  |  |
2.5 | 1.76[-3] | 9.25[-4] | 2.90[-3] | 5.58[-3] | 1.68[-3] | -6.00[-4] | 2.15[-3] | 3.23[-3]
3.5 | 2.33[-3] | 7.67[-4] | 5.55[-3] | 8.64[-3] | 2.96[-3] | -6.74[-4] | 3.90[-3] | 6.18[-3]
4.5 | 2.76[-3] | 4.88[-4] | 5.52[-3] | 8.77[-3] | 3.71[-3] | -3.66[-4] | 1.21[-4] | 3.47[-3]
5.5 | 2.84[-3] | -1.04[-4] | 1.42[-4] | 2.88[-3] | 3.30[-3] | 1.05[-3] | 1.08[-4] | 4.46[-3]
6.5 | 2.31[-3] | -1.39[-3] | 2.70[-4] | 1.19[-3] |  |  |  |
| $\langle\Delta E\rangle$ | 4.72[-3] | $\langle\Delta E\rangle$ | 4.36[-3]
Th II, $6d7s^{2}\ J_{0}=3/2$
0.0 |  |  |  |  | 5.26[-2] | -1.11[-2] | 8.12[-4] | 4.23[-2]
1.0 | 8.41[-2] | 1.03[-2] | 9.51[-4] | 9.53[-2] | 1.25[-1] | -3.31[-3] | 4.29[-4] | 1.22[-1]
2.0 | 1.41[-1] | 1.56[-3] | 6.40[-4] | 1.44[-1] | 2.28[-1] | -1.32[-4] | 4.30[-4] | 2.28[-1]
3.0 | 1.93[-1] | 3.74[-4] | 7.84[-4] | 1.94[-1] | 2.76[-1] | 3.10[-3] | 1.58[-4] | 2.80[-1]
4.0 | 1.95[-1] | -4.59[-3] | 3.33[-4] | 1.91[-1] |  |  |  |
| $\langle\Delta E\rangle$ | 1.70[-1] | $\langle\Delta E\rangle$ | 2.19[-1]
Th I, $6d^{2}7s^{2}\ J_{0}=2$
0.5 | 4.42[-4] | 3.27[-4] | 1.52[-4] | 9.21[-4] | 4.91[-4] | -2.42[-4] | 2.11[-4] | 4.59[-4]
1.5 | 6.25[-4] | 4.23[-4] | 3.87[-4] | 1.44[-3] | 8.61[-4] | -3.42[-4] | 1.24[-3] | 1.75[-3]
2.5 | 8.50[-4] | 3.40[-4] | 2.11[-3] | 3.30[-3] | 1.24[-3] | -9.87[-5] | 1.11[-4] | 1.25[-3]
3.5 | 9.97[-4] | -4.45[-5] | 2.22[-4] | 1.17[-3] | 1.26[-3] | 3.06[-4] | 2.42[-5] | 1.59[-3]
4.5 | 8.97[-4] | -4.03[-4] | 6.03[-5] | 5.54[-4] |  |  |  |
| $\langle\Delta E\rangle$ | 1.41[-3] | $\langle\Delta E\rangle$ | 1.41[-3]
Table 6: The difference in nuclear frequencies between Th ions due to the second-order hfs shifts. The average shifts $\langle\Delta E\rangle$ from Table 5 are used to obtain the numbers. The numbers in square brackets stand for powers of ten. Ions | $\Delta\omega_{N}$ | $\Delta\omega_{N}/\omega_{N}$
---|---|---
|  |  | (MHz) | (eV) |
Th I | $-$ | Th II | -4.916[-2] | -2.033[-10] | -2.438[-11]
Th I | $-$ | Th III | 3.495[-4] | 1.446[-12] | 1.734[-13]
Th I | $-$ | Th IV | 2.459[-4] | 1.017[-12] | 1.220[-13]
Th II | $-$ | Th III | 4.951[-2] | 2.048[-10] | 2.456[-11]
Th II | $-$ | Th IV | 4.940[-2] | 2.043[-10] | 2.450[-11]
Th III | $-$ | Th IV | -1.037[-4] | -4.287[-13] | -5.142[-14]

The second-order hfs shift is given by

$\displaystyle\Delta E_{v}=\sum_{n}\left[\frac{\langle
v|\hat{H}_{A}|n\rangle^{2}}{E_{v}-E_{n}}+\right.$ (10)
$\displaystyle\left.2\frac{\langle v|\hat{H}_{A}|n\rangle\langle
n|\hat{H}_{B}|v\rangle}{E_{v}-E_{n}}+\frac{\langle
v|\hat{H}_{B}|n\rangle^{2}}{E_{v}-E_{n}}\right]$ $\displaystyle\equiv\Delta
E_{a}+\Delta E_{ab}+\Delta E_{b}.$ (11)

Here $\hat{H}_{A}$ and $\hat{H}_{B}$ are the magnetic dipole and electric
quadrupole hfi operators. Corresponding reduced matrix elements are given by
(see, e.g. Ref. hfs2 )

$\displaystyle\langle
J_{v},F|\hat{H}_{A}|J_{n},F\rangle=(-1)^{I+F+J_{n}}\mu\times$ (12)
$\displaystyle\sqrt{\frac{(2I+1)(I+1)}{I}}\left\\{\begin{array}[]{ccc}I&J_{n}&F\\\
J_{v}&I&1\end{array}\right\\}\langle J_{v}||\hat{H}^{e}_{A}||J_{n}\rangle,$
(15)

and

$\displaystyle\langle
J_{v},F|\hat{H}_{B}|J_{v},F\rangle=(-1)^{I+F+J_{n}}\frac{Q}{2}\times$ (16)
$\displaystyle\sqrt{\frac{(2I+1)(I+1)(2I+3)}{I(2I-1)}}\left\\{\begin{array}[]{ccc}I&J_{n}&F\\\
J_{v}&I&2\end{array}\right\\}\langle J_{v}||\hat{H}^{e}_{B}||J_{n}\rangle.$
(19)

We include only limited number of low-energy electron states in the summation
over $n$ in (10). They dominate over the rest of the sum due to small energy
denominators. Saturation of the summation has been observed on the level of
few per cent.

Nuclear excited states $|n\rangle$ also give a contribution to the sum in the
expression for $\Delta E_{v}$ in Eq. (10). However, a minimal energy
denominator $E_{v}-E_{n}$ in this case is the interval between the ground and
isomeric nuclear states equal to 8.3 eV. This interval is 1-2 orders of
magnitude bigger than the interval between electron states (fine structure
intervals in this case). As a result, the nuclear excitations contribution may
be neglected in comparison with the electron excitations contribution.

As in case of the first-order hfs shifts, the second-order shift (10) is
different for the ground and isomeric nuclear states due to different nuclear
parameters $\mu$, $I$ and $Q$. The results of calculations are presented in
Tables 5 and 6. The second order hfi shift is much smaller than the first
order shift. However, the first-order shift vanishes after averaging over hfs
components while the second-order shift does not. Note that the second-order
hfi shift is three to six orders of magnitude smaller than the shift due to
the electrostatic interaction between atomic electrons and the nucleus
considered in our previous work previous . However, it is six to eight orders
of magnitude larger than the projected relative uncertainty of the nuclear
clock ($10^{-19}$).

###### Acknowledgements.

This work was supported by the Australian Research Council Grants No.
DP230101058 and DP200100150.

## References

  * (1) E. Peik and Chr. Tamm, Nuclear laser spectroscopy of the 3.5 eV transition in Th-229, Europhys. Lett. 61, 181 (2003).
  * (2) J. Thielking, M. V. Okhapkin, P. Glowacki, D. M. Meier, L. von der Wense, B. Seiferle, C. E. Düllmann, P. G. Thirolf, and E. Peik, Laser spectroscopic characterization of the nuclear-clock isomer 229mTh, Nature 556, 321 (2018). https://doi.org/10.1038/s41586-018-0011-8
  * (3) E. Peik, T. Schumm, M. S. Safronova, A. Pálffy, J. Weitenberg, and P. G. Thirolf, Nuclear clocks for testing fundamental physics, Quantum Sci. Technol. 6, 034002 (2021).
  * (4) K. Beeks, T. Sikorsky, T. Schumm, J. Thielking, M. V. Okhapkin and E. Peik, The thorium-229 low-energy isomer and the nuclear clock, Nat. Rev. Phys. 3, 238 (2021).
  * (5) C. J. Campbell, A. G. Radnaev, A. Kuzmich, V. A. Dzuba, V. V. Flambaum, and A. Derevianko, Single-Ion Nuclear Clock for Metrology at the 19th Decimal Place, Phys. Rev. Lett. 108, 120802 (2012).
  * (6) V. V. Flambaum, Enhanced Effect of Temporal Variation of the Fine Structure Constant and the Strong Interaction in 229Th, Phys. Rev. Lett. 97, 092502 (2006).
  * (7) V. V. Flambaum, Enhancing the effect of Lorentz invariance and Einstein’s equivalence principle violation in nuclei and atoms, Phys. Rev. Lett. 117, 072501, (2016).
  * (8) V.V. Flambaum, R.B. Wiringa, Enhanced effect of quark mass variation in 229Th and limits from Oklo data, Phys. Rev. C 79, 034302 (2009).
  * (9) Pavel Fadeev, Julian C. Berengut, and Victor V. Flambaum, Sensitivity of 229Th nuclear clock transition to variation of the fine-structure constant, Phys. Rev. A 102, 052833 (2020).
  * (10) J.C. Berengut, V.A. Dzuba, V.V. Flambaum, S.G. Porsev, Proposed experimental method to determine the alpha sensitivity of splitting between ground and 7.6 eV isomeric states in 229Th, Phys. Rev. Lett. 102, 210801 (2009).
  * (11) E. Litvinova, H. Feldmeier, J. Dobaczewski, V.V. Flambaum, Nuclear structure of lowest 229Th states and time-dependent fundamental constants, Phys. Rev. C 79, 064303 (2009).
  * (12) P. Fadeev, J. C. Berengut, V. V. Flambaum, Effects of variation of the fine structure constant $\alpha$ and quark mass $m_{q}$ in Mössbauer nuclear transitions, Phys. Rev. C 105, L051303 (2022).
  * (13) A. Arvanitaki, J. Huang, and K. Van Tilburg, Searching for dilaton dark matter with atomic clock, Phys. Rev. D 91, 015015 (2015).
  * (14) Y. V. Stadnik, V. V. Flambaum, Can dark matter induce cosmological evolution of the fundamental constants of Nature? Phys. Rev. Lett. 115, 201301, (2015).
  * (15) L. von der Wense, B. Seiferle, M. Laatiaoui, J. B. Neumayr, H.-J. Maier, H.-F. Wirth, C. Mokry, J. Runke, K. Eberhardt, C. E. Düllmann et al, Direct detection of the 229Th nuclear clock transition, Nature 533, 47 (2016).
  * (16) B. Seiferle, L. von der Wense, and P. G. Thirolf, Lifetime measurement of the 229mTh nuclear isomer, Phys. Rev. Lett. 118, 042501 (2017).
  * (17) B. Seiferle, L. von der Wense, P. V. Bilous, I. Amersdorffer, C. Lemell, F. Libisch, S. Stellmer, T. Schumm, C. E. Düllmann, A. Pálffy, and P. G. Thirolf, Energy of the 229Th nuclear clock transition, Nature 573, 243 (2019).
  * (18) A. Yamaguchi, H. Muramatsu, T. Hayashi, N. Yuasa, K. Nakamura, M. Takimoto, H. Haba, K. Konashi, M. Watanabe, H. Kikunaga, K. Maehata, N.Y. Yamasaki,and K. Mitsuda, Energy of the 229Th nuclear clock isomer determined by absolute $\gamma$-ray energy difference, Phys. Rev. Lett. 123, 222501 (2019).
  * (19) T. Sikorsky, J. Geist, D. Hengstler, S. Kempf, L. Gastaldo, C. Enss, C. Mokry, J. Runke, C. E. Düllmann, P. Wobrauschek, et al, Measurement of the 229Th isomer energy with a magnetic microcalorimeter, Phys. Rev. Lett. 125, 142503 (2020).
  * (20) S. Kraemer, J. Moens, M. Athanasakis-Kaklamanakis, S. Bara, K. Beeks, P. Chhetri, K. Chrysalidis, A. Claessens, T. E. Cocolios, J. G. M. Correia, et al, Observation of the radiative decay of the 229Th nuclear clock isomer, Nature 617, 706 (2023).
  * (21) A. Claessens, F. Ivandikov, S. Bara, P. Chhetri, A. Dragoun, Ch.E. Düllmann, Y. Elskens, R. Ferrer, S. Kraemer, Yu. Kudryavtsev, et al, Laser ionization scheme development for in-gas-jet spectroscopy studies of Th+, Nuclear Inst. and Methods in Physics Research B 540, 224 (2023).
  * (22) M. S. Safronova, S. G. Porsev, M. G. Kozlov, J. Thielking, M. V. Okhapkin,P. Glowacki, D. M. Meier, and E. Peik, Nuclear Charge Radii of 229Th from Isotope and Isomer Shifts, Phys. Rev. Lett. 121, 213001 (2018).
  * (23) Wade G. Rellergert, D. DeMille, R. R. Greco, M. P. Hehlen, J. R. Torgerson, and Eric R. Hudson, Constraining the evolution of the fundamental constants with a solid-state optical frequency reference based on the 229Th nucleus, Phys. Rev. Lett. 104, 200802 (2010).
  * (24) G. A. Kazakov, A. N. Litvinov, V. I. Romanenko, L. P. Yatsenko, A. V. Romanenko, M. Schreitl, G. Winkler, and T. Schumm, Performance of a 229Thorium solid-state nuclear clock, New J. Phys. 14, 083019 (2012).
  * (25) K. Beeks, T. Sikorsky, V. Rosecker, M. Pressler, F. Schaden, D. Werban, N. Hosseini, L. Rudischer, F. Schneider, P. Berwian, et al, Growth and characterization of thorium-doped calcium fluoride single crystals, Sci Rep 13, 3897 (2023). https://doi.org/10.1038/s41598-023-31045-5
  * (26) V. A. Dzuba and V. V. Flambaum, Effects of electrons on nuclear clock transition frequency in 229Th ions, arXiv:2309.11176 (2023).
  * (27) W. Kälber, J. Rink, K. Bekk, W. Faubel, S. Göring, G. Meisel, H. Rebel, and R. C. Thompson, Nuclear Radii of Stored Ions of Thorium Isotopes from Laser Spectroscopy, Z. Phys. A - Atomic Nuclei 334, 103 (1989).
  * (28) C. J. Campbell, A. G. Radnaev, and A. Kuzmich, Wigner Crystals of 229Th for Optical Excitation of the Nuclear Isomer, Phys. Rev. Lett. 106, 223001 (2011).
  * (29) M. S. Safronova, U. I. Safronova, A. G. Radnaev, C. J. Campbell, and A. Kuzmich, Magnetic dipole and electric quadrupole moments of the 229Th nucleus, Phys. Rev. A 88, 060501(R) (2013).
  * (30) C. E. Bemis, F. K. McGowan, J. L. C. Ford Jr, W. T. Milner, R. L. Robinson, P. H. Stelson, G. A. Leander and C. W. Reich, Coulomb excitation of states in 229Th, Phys. Scr. 38, 657 (1988).
  * (31) K. Beloy, Hyperfine Structure in 229gTh3+ as a Probe of the 229gTh $\rightarrow$ 229mTh Nuclear Excitation Energy, Phys. Rev. Lett. 112, 062503 (2014).
  * (32) M. V. Okhapkin, D. M. Meier, E. Peik, M. S. Safronova, M. G. Kozlov, and S. G. Porsev, Observation of an unexpected negative isotope shift in 229Th+ and its theoretical explanation, Phys. Rev. A 92, 020503(R) (2015).
  * (33) R. A. Müller, A. V. Maiorova, S. Fritzsche, A. V. Volotka, R. Beerwerth, P. Glowacki, J. Thielking, D.-M. Meier, M. Okhapkin, E. Peik, and A. Surzhykov, Hyperfine interaction with the 229Th nucleus and its low-lying isomeric state, Phys. Rev. A 98, 020503(R) (2018).
  * (34) W. R. Johnson, and J. Sapirstein, Computation of Second-Order Many-Body Corrections in Relativistic Atomic Systems, Phys. Rev. Lett. 57, 1126 (1986).
  * (35) V. A. Dzuba, V. V. Flambaum, P. G. Silvestrov, O. P. Sushkov, Correlation potential method for the calculation of energy levels, hyperfine structure and E1 transition amplitudes in atoms with one unpaired electron, J. Phys. B: At. Mol. Phys., 20, 1399-1412 (1987).
  * (36) V. A. Dzuba, Combination of the single-double coupled cluster and the configuration interaction methods: application to barium, lutetium and their ions, Phys. Rev. A 90, 012517 (2014).
  * (37) A. Kramida, Yu. Ralchenko, J. Reader, and NIST ASD Team (2022). NIST Atomic Spectra Database (ver. 5.10), [Online]. Available: https://physics.nist.gov/asd [2023, October 8]. National Institute of Standards and Technology, Gaithersburg, MD. DOI: https://doi.org/10.18434/T4W30F
  * (38) V. A. Dzuba and V. V. Flambaum, Hyperfine-induced electric dipole contributions to the electric octupole and magnetic quadrupole atomic clock transitions, Phys. Rev. A 93, 052517 (2016).


# Supplementary Material for “Irreversible Monte Carlo algorithms for hard
disk glasses: From event-chain to collective swaps”

Federico Ghimenti  Laboratoire Matière et Systèmes Complexes (MSC), Université
Paris Cité & CNRS (UMR 7057), 75013 Paris, France  Ludovic Berthier
Laboratoire Charles Coulomb (L2C), Université de Montpellier & CNRS (UMR
5221), 34095 Montpellier, France  Gulliver, UMR CNRS 7083, ESPCI Paris, PSL
Research University, 75005 Paris, France  Frédéric van Wijland  Laboratoire
Matière et Systèmes Complexes (MSC), Université Paris Cité & CNRS (UMR 7057),
75013 Paris, France

##  I A model with a different polydispersity

To show that the efficiency of the cSwap algorithm does not depend on the
specific polydispersity of the model studied in the main text, we study the
relaxation dynamics of a second model. It consists of $N=1024$ hard disks with
a diameter distribution following a power law,
$\pi(\sigma)\propto\sigma^{-3}$, with
$\sigma_{\text{min}}\leq\sigma\leq\sigma_{\text{max}}$, as studied earlier
[1]. We recall that the units of length are chosen so that
$\overline{\sigma}=1$. We have chosen the support of the diameter distribution
such that the polydispersity is $\Delta=23\%$. We ran $NVT$ simulations using
the six different algorithms studied in the main text at a large packing
fraction, $\phi=0.86$, where we measured the time decay of $C_{\psi}$. The
results are shown in Fig. 1. The hierarchy of speedups obtained by the
different algorithm, and their relative values, is comparable to the results
reported in the main text for a similar value of $\phi$.

Figure 1: Time dependence of the hexatic correlation function for a system of
$N=1024$ particles with a power law distribution of the diameters. Time is
measured in units of $Nt_{\text{move}}$, and the efficiency of the algorithms
that involve swap moves is comparable to the results shown in Fig. 2 of the
main text near $\phi\approx 0.85$.

##  II Diameter distribution of the hard disk system

Figure 2 shows the probability distribution function of the diameters $\sigma$
for the polydisperse mixture of hard disks studied in the main text. The
polydipersity of the system, defined as
$\Delta\equiv\frac{\sqrt{\overline{\sigma^{2}}-\overline{\sigma}^{2}}}{\overline{\sigma}}\approx
25\%$. This diameter distribution results from a gradient descent protocol in
an extended space composed of particles positions and radii, which will be
detailed in [2], which represents an improvement over earlier versions [3].

Figure 2: Particle size distribution $\pi(\sigma)$ for the system investigated
in this work.

##  III Event-Chain Monte Carlo

We recall the implementation of the Event-chain Monte Carlo (ECMC) dynamics.
We resort to the so-called ‘straight’ version, as described in [4, 5].

The state of the system is described by a configuration
$\mathcal{C}=\\{\text{\bf r}^{N},\sigma^{N},i,\text{\bf v}\\}$ where we store
the position of the $N$ particles $\text{\bf r}^{N}$, their diameters
$\sigma^{N}$, the label of the active particle $i$ and the direction of self
propulsion v, a two-dimensional vector of unit norm. In the mathematical
literature on Markov chains, $i$ and v are called lifting degrees of freedom,
and they govern the nonequilibrium dynamics of the system. An ECMC move
consists in moving the active particle $i$ along the direction v, until an
event, i.e. a collision with another particle, whose label is denoted as $j$,
occurs. The distance $\delta\ell_{ij}$ traveled by $i$ during a move is thus
determined by the equation

$\delta\ell_{ij}=\text{\bf r}_{ji}\cdot\text{\bf
v}-\sqrt{\sigma_{ij}^{2}-\left(\text{\bf r}_{ji}\cdot\text{\bf
v}_{\perp}\right)^{2}},$ (1)

where $\text{\bf r}_{ji}=\text{\bf r}_{j}-\text{\bf r}_{i}$ is the vector
joining particle $i$ with particle $j$,
$\sigma_{ij}\equiv\frac{\sigma_{i}+\sigma_{j}}{2}$ is the effective diameter,
and $\text{\bf v}_{\perp}=(v_{y},-v_{x})$ is the direction orthogonal to v. In
practice, in the simulations the inter-particle distance $\text{\bf r}_{ij}$
is computed according to the minimum image convention [6] to keep into account
periodic boundary conditions, and the particle $j$ is identified through an
event-driven scheme, by minimizing $\delta\ell_{ik}$ among all possible target
particles $k$, i.e. $j=\arg\min_{k}{\delta\ell_{ik}}$. After an event occurs,
a lifting move is performed: the activity label changes from $i$ to $j$. In
the next ECMC move particle $j$ will perform directed motion along the
direction v. The composition of several ECMC moves builds up a chain of
particles performing directed motion. When the displacement performed by the
active particles add up to a fixed parameter $\ell$, which fixes the length of
the chain, the activity label and the self propulsion direction are resampled
uniformly in their domain of definition, which are respectively the set
$\\{1,2,\ldots,N\\}$ for the activity label and the set $\\{\text{\bf
e}_{x},\text{\bf e}_{y}\\}$ for the self-propulsion, for the case of straight
ECMC. As the density of the system increases, more and more particles
participate to the chain. This is demonstrated in Figure 3, which shows
$n_{\ell}$, the average number of particles in a chain as a function of the
packing fraction $\phi$ for a fixed value of $\ell$.

Figure 3: Average number of particles $n_{\ell}$ participating in a chain of
fixed length $\ell$ as a function of the packing fraction $\phi$.

##  IV Hexatic correlation function

We recall the definition of the time correlations for the hexatic order
parameter, which have been extensively used in the past to monitor the
relaxation of two-dimensional glassy systems [7]. The hexatic correlation
function $C_{\psi}(t)$ is defined as

$C_{\psi}\equiv\frac{\left\langle\psi^{*}(t)\psi(0)\right\rangle}{\left\langle\lvert\psi(0)\rvert^{2}\right\rangle},$
(2)

where $\psi(t)$ is the global hexatic order parameter. It is defined as a sum
of local terms $\psi(t)\equiv\frac{1}{N}\sum_{i=1}^{N}\phi_{6,i}(t)$ where the
local hexatic order parameter for particle $i$, $\phi_{6,i}(t)$ is defined as

$\phi_{6,i}(t)\equiv\frac{1}{n_{i}}\sum_{j=1}^{n_{i}}e^{6i\theta_{ij}},$ (3)

with the sum running over the $n_{i}$ neighbors of particle $i$, defined
through a Voronoi tessellation, constructed using the freud libary [8], and
$\theta_{ij}$ is the angle between the vector $\text{\bf
r}_{ij}\equiv\text{\bf r}_{i}-\text{\bf r}_{j}$ and the vector $\text{\bf
e}_{x}$. From the time decay of $C_{\psi}$ we define the relaxation time
$\tau_{\alpha}$, such that $C_{\psi}(\tau_{\alpha})=e^{-1}$.

##  V Ergodicity of cSwap in an $N=4$ system

Here we prove the ergodicity of the cSwap dynamics in the case of a small
polydisperse system at low densities. The discrete-time dynamics we consider
is made up only by cSwap moves, the positions of the disks $\text{\bf r}^{N}$
being fixed at all times. The state of the system
$\mathcal{C}=\\{\sigma^{N},i\\}$ is specified by assigning a diameter to each
of the $N$ particles –the resulting permutation of diameters is denoted by
$\sigma^{N}$– and by the lifting degree of freedom $i$, i.e. the label of the
active particle. The system can access a subset of $N!\times N$
configurations. The configurations in the subset satisfy the non-overlapping
hard disks condition. The probability $\pi_{t}(\mathcal{C})$ for the system to
be in state $\mathcal{C}$ at time $t$ obeys the discrete time Markov dynamics

$\pi_{t+1}(\mathcal{C})=\sum_{\mathcal{C}^{\prime}}P(\mathcal{C}^{\prime}\to\mathcal{C})\pi_{t}(\mathcal{C}^{\prime}),$
(4)

with $P$ the transition matrix encoding the probability to jump from one
configuration to another during a discrete time-step. The convergence
properties of the dynamics are encoded in the spectrum of $P$. Proving that
the cSwap dynamics is ergodic amounts to showing that the spectrum of $P$ has
a unique, nondegenerate, eigenvalue $\lambda=1$ lying on the unit circle, and
that all the other eigenvalues have a norm strictly smaller than $1$ [9].

We consider the case where $N=4$ and the particles are far away from each
other, so that no overlap between them can be generated upon permutation of
the radii. To write down the transition matrix describing the cSwap dynamics,
we decompose the configuration space into a tensor product $S_{4}\otimes i$,
where $S_{4}$ is an element of the permutation group of 4 elements, and $i$ is
the lifting degree of freedom labeling the active, expanding particle. The
permutations are labeled by natural numbers in the following way:

$\begin{array}[]{cccc}1\equiv\\{1,2,3,4\\}&2\equiv\\{2,3,4,1\\}&3\equiv\\{3,4,1,2\\}&4\equiv\\{4,1,2,3\\}\\\
5\equiv\\{1,3,4,2\\}&6\equiv\\{3,4,2,1\\}&7\equiv\\{4,2,1,3\\}&8\equiv\\{2,1,3,4\\}\\\
9\equiv\\{1,2,4,3\\}&10\equiv\\{2,4,3,1\\}&11\equiv\\{4,3,1,2\\}&12\equiv\\{3,1,2,4\\}\\\
13\equiv\\{2,4,1,3\\}&14\equiv\\{4,1,3,2\\}&15\equiv\\{1,3,2,4\\}&16\equiv\\{3,2,4,1\\}\\\
17\equiv\\{2,3,1,4\\}&18\equiv\\{3,1,4,2\\}&19\equiv\\{1,4,2,3\\}&20\equiv\\{4,2,3,1\\}\\\
21\equiv\\{2,1,4,3\\}&22\equiv\\{1,4,3,2\\}&23\equiv\\{4,3,2,1\\}&24\equiv\\{3,2,1,4\\}\\\
\end{array}$ (5)

Here the numbers inside the brackets denote the index of the particles,
ordered according to their diameters, from the smallest to the largest. For
example, $\\{2,3,4,1\\}$ is the configurations where particle with label $1$
has the largest diameter, particle with label $4$ has the second largest
diameter, and so on. If the system is in a configuration where particle $3$ is
the active one, we would have $\mathcal{C}=\\{2,3,4,1\\}\otimes 3$.

With this notation, the transition matrix can be written as a $24\times 24$
block matrix, each block being made by a $4\times 4$ sub-block. Transitions
between blocks represent changes in the diameter assignments for the different
particles, while transitions within a block represent change in the active
degree of freedom. The transition matrix describing the cSwap dynamics is

$P=\left[\begin{array}[]{*{24}c}D_{4}&B_{1}&\cdot&\cdot&B_{2}&\cdot&\cdot&\cdot&B_{3}&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot\\\
\cdot&D_{4}&B_{1}&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&B_{2}&\cdot&\cdot&\cdot&B_{3}&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot\\\
\cdot&\cdot&D_{4}&B_{1}&\cdot&B_{3}&\cdot&\cdot&\cdot&\cdot&\cdot&B_{2}&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot\\\
B_{1}&\cdot&\cdot&D_{4}&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&B_{3}&\cdot&\cdot&\cdot&\cdot&\cdot&B_{2}&\cdot&\cdot&\cdot&\cdot\\\
\cdot&\cdot&\cdot&\cdot&D_{4}&B_{1}&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&B_{3}&\cdot&\cdot&\cdot&B_{2}&\cdot&\cdot&\cdot&\cdot&\cdot\\\
\cdot&\cdot&B_{3}&\cdot&\cdot&D_{4}&B_{1}&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&B_{2}\\\
\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&D_{4}&B_{1}&\cdot&\cdot&\cdot&\cdot&\cdot&B_{2}&\cdot&\cdot&\cdot&\cdot&\cdot&B_{3}&\cdot&\cdot&\cdot&\cdot\\\
\cdot&B_{2}&\cdot&\cdot&B_{1}&\cdot&\cdot&D_{4}&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&B_{3}&\cdot&\cdot&\cdot\\\
B_{3}&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&D_{4}&B_{1}&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&B_{2}&\cdot&\cdot\\\
\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&D_{4}&B_{1}&\cdot&B_{3}&\cdot&\cdot&\cdot&B_{2}&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot\\\
\cdot&\cdot&\cdot&B_{2}&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&D_{4}&B_{1}&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&B_{3}&\cdot\\\
\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&B_{1}&\cdot&\cdot&D_{4}&\cdot&\cdot&\cdot&B_{2}&\cdot&B_{3}&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot\\\
\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&B_{2}&\cdot&B_{3}&\cdot&\cdot&D_{4}&B_{1}&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot\\\
\cdot&\cdot&\cdot&B_{3}&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&D_{4}&B_{1}&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&B_{2}&\cdot\\\
\cdot&\cdot&\cdot&\cdot&B_{3}&\cdot&\cdot&\cdot&B_{2}&\cdot&\cdot&\cdot&\cdot&\cdot&D_{4}&B_{1}&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot\\\
\cdot&\cdot&B_{2}&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&B_{1}&\cdot&\cdot&D_{4}&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&B_{3}\\\
\cdot&B_{3}&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&D_{4}&B_{1}&\cdot&\cdot&B_{2}&\cdot&\cdot&\cdot\\\
\cdot&\cdot&\cdot&\cdot&\cdot&B_{2}&\cdot&\cdot&\cdot&\cdot&\cdot&B_{3}&\cdot&\cdot&\cdot&\cdot&\cdot&D_{4}&B_{1}&\cdot&\cdot&\cdot&\cdot&\cdot\\\
B_{2}&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&D_{4}&B_{1}&\cdot&B_{3}&\cdot&\cdot\\\
\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&B_{3}&\cdot&\cdot&\cdot&B_{2}&\cdot&\cdot&\cdot&\cdot&\cdot&B_{1}&\cdot&\cdot&D_{4}&\cdot&\cdot&\cdot&\cdot\\\
\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&B_{3}&\cdot&B_{2}&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&D_{4}&B_{1}&\cdot&\cdot\\\
\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&B_{2}&\cdot&\cdot&\cdot&B_{3}&\cdot&\cdot&D_{4}&B_{1}&\cdot\\\
\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&B_{2}&\cdot&\cdot&\cdot&B_{3}&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&D_{4}&B_{1}\\\
\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&\cdot&B_{3}&\cdot&B_{2}&\cdot&\cdot&B_{1}&\cdot&\cdot&D_{4}\\\
\end{array}\right].$ (6)

The dot symbol $\cdot$ denote $4\times 4$ matrices whose entries are all
zeros. The diagonal block is a sum of two terms, $D_{4}\equiv A_{4}+B_{4}$.
The block of type $A_{4}$ encodes transitions involving only the active degree
of freedom,

$A_{4}\equiv\begin{bmatrix}\alpha&\alpha&\alpha&\alpha\\\
\alpha&\alpha&\alpha&\alpha\\\ \alpha&\alpha&\alpha&\alpha\\\
\alpha&\alpha&\alpha&\alpha\\\ \end{bmatrix},$ (7)

with $\alpha=1/4^{2}=1/16$. They are associated to random resampling of the
activity label, and they are used here to ensure that $P$ is ergodic. The
blocks $B_{i}$ encode transitions generated by the cSwap moves, involving an
inflation of the active particle and the cascade of swaps. They are given by:

$\begin{array}[]{cccc}B_{1}\equiv\begin{bmatrix}\circ&\circ&\circ&1-\frac{1}{4}\\\
\circ&\circ&\circ&\circ\\\ \circ&\circ&\circ&\circ\\\
\circ&\circ&\circ&\circ\end{bmatrix}&B_{2}\equiv\begin{bmatrix}\circ&\circ&\circ&\circ\\\
1-\frac{1}{4}&\circ&\circ&\circ\\\ \circ&\circ&\circ&\circ\\\
\circ&\circ&\circ&\circ\end{bmatrix}&B_{3}\equiv\begin{bmatrix}\circ&\circ&\circ&\circ\\\
\circ&\circ&\circ&\circ\\\ \circ&1-\frac{1}{4}&\circ&\circ\\\
\circ&\circ&\circ&\circ\end{bmatrix}&B_{4}\equiv\begin{bmatrix}\circ&\circ&\circ&\circ\\\
\circ&\circ&\circ&\circ\\\ \circ&\circ&\circ&\circ\\\
\circ&\circ&1-\frac{1}{4}&\circ\end{bmatrix},\end{array}$ (8)

with the symbol ‘$\circ$’ denoting the entries with value $0$.

The matrix $P$ is doubly stochastic, i.e. the sum of all the elements
belonging to a fixed row, or fixed column, is one. A stationary solution of
the stochastic process associated with $P$ is given by

$\pi_{ss}\equiv\frac{1}{4\times
4!}\bigoplus_{i=1}^{24}\left[1,1,1,1\right]^{T},$ (9)

which is the tensor product of the Boltzmann distribution for hard disks times
a uniform distribution of the active degree of freedom among the $4$
particles. To show ergodicity, we inspect the spectrum of $P$. Its eigenvalues
$\lambda$ are shown in the complex plane in Fig. 4. The only eigenvalue lying
on the unit circle is $\lambda=1$, thus proving the ergodicity of the cSwap
Markov chain. We observe that the other eigenvalues tend to accumulate at the
vertices of an octagon inside the unit circle. This is a consequence of the
introduction of a refreshment probability for the label of the active
particle. Without such a refreshment dynamics, the Markov chain would be
periodic, and there would be $2N$ eigenvalues lying on the unit circle. The
introduction of a refreshment rate pushes the eigenvalues inside the unit
circle, making the Markov Chain aperiodic.

Figure 4: Representation in the complex plane of the eigenvalues $\lambda$ of
the transition matrix $P$ given by Eq. (6) for the cSwap dynamics in the case
of a system of $N=4$ particles. The dotted line is the unit circle.

##  VI Additional numerical evidence of ergodicity for cSwap dynamics

We present numerical tests supporting the ergodicity of the cSwap and ECMC
algorithms. The first test is the numerical calculation of the equation of
state $Z(\phi)$ of the polydisperse system considered in the main text. The
equation of state relates the reduced pressure $Z=\frac{\beta P}{\rho}$ to the
packing fraction $\phi$. Here $P$ is the pressure applied to the system,
$\rho$ is the number density and $\beta^{-1}=k_{\text{\tiny{B}}}T$. It can be
obtained from simulations in the $NPT$ ensemble, where one has access to the
running averages of $Z$ and $\phi$ at fixed applied pressure $P$. The results
are shown in Fig. 5(a),

Figure 5: Numerical test of ergodicity for the cSwap and ECMC algorithms. (a)
Equation of state $Z(\phi)$ for the polydisperse hard disks system. (b)
Rescaled radial distribution function $g(r/\sigma-1)$ as a function of the
distance from its first peak $r/\sigma-1$. (c) Probability distribution
function of the packing fractions $\pi(\phi)$ explored by the system during
$NPT$ simulations at $\beta P=16.3\overline{\sigma^{2}}$.

where they are compared with an extension to polydisperse systems of the
empirical Henderson formula [10, 11] which reads

$\begin{split}Z_{\rm
eH}&=\frac{1-\left(1-\frac{\overline{\sigma}^{2}}{\overline{\sigma^{2}}}\right)\phi+(b-3)\frac{\overline{\sigma}^{2}}{\overline{\sigma^{2}}}\phi^{2}}{\left(1-\phi\right)^{2}},\\\
b&\equiv\frac{16}{3}-\frac{4\sqrt{3}}{\pi}.\end{split}$ (10)

The agreement with the empirical formula and with the conventional Metropolis
and Swap algorithms is excellent, and this serves as a guide to detect
deviations from one algorithm to another. We find that all algorithms agree
with each other.

We next compare the rescaled radial distribution function $g(r/\sigma)$ in the
NPT ensemble in Fig. 5(b). Its expression is given by

$\begin{split}g(x)&\equiv\sum_{\begin{subarray}{c}i,j\\\
i<j\end{subarray}}C_{ij}\int_{x_{b}}^{x_{b}+\Delta
x}\delta\left(x^{\prime}-\frac{r_{ij}}{\sigma_{ij}}\right)dx^{\prime}\\\
C_{ij}&\equiv\frac{1}{2\pi(x_{b}+\Delta x/2)\Delta
x\sigma_{ij}^{2}\rho(N-1)}\end{split}$ (11)

where $\sigma_{ij}\equiv\frac{1}{2}\left(\sigma_{i}+\sigma_{j}\right)$,
$\rho=N/L^{2}$ is the number density of the system, and we collect the
rescaled interparticle distances in a histogram with bin width $\Delta x$.
$x_{b}=b\Delta_{x}$ is the coordinate of the $b$-th bin, with $b$ chosen so
that $x_{b}\leq x<x_{b}+\Delta x$. Again, the curves obtained with the
different algorithms superimpose on each other.

Finally, we report in Fig. 5(c) the histogram of the packing fractions
explored during an NPT simulation, $\pi(\phi)$, for a fixed pressure $P$.
Again all algorithms explore the same fluctuations, including in the tails of
the distribution, showing that the same Boltzmann distribution is indeed
properly sampled in all our algorithms.

##  VII Choice of time units for the different Monte Carlo algorithms

In the main text a unit of time $t_{\text{move}}$ was chosen based upon the
elementary transitions involved in each algorithm. Here, we clarify its
relation with the CPU time required by the algorithms studied in the main
text. We first observe that in all the algorithms investigated, during a unit
of $t_{\text{move}}$ only one evaluation of the overlap of a given disk with
its neighbors is performed. This is usually the most computationally demanding
task in hard disks simulations. We thus expect the trend of the relaxation
time in units of CPU time to resemble the one obtained using units of
$t_{\text{move}}$. This is confirmed in Fig. 6, where we show the relaxation
time $\tau_{\alpha}$, now measured in units of CPU time, as a function of
$\phi$. The trend of the curves is virtually identical to the one shown in
Fig. 1 of the main text.

Figure 6: Equilibrium relaxation times of the six algorithms investigated in
the main text, in units of CPU time. Here $t_{\text{CPU}}=1$ second.

When using the ECMC algorithm, it is also possible to measure times using the
number of directed chains of particles that have moved [4]. When comparing
ECMC with other types of algorithms, however, this choice is a poor indicator
of its efficiency, as it hides the number of particles–and hence of event
determinations–that are involved in each chain. This is demonstrated in Fig.
7(a), where the time relaxation of $C_{\psi}$ is shown as a function the
number of chain displaced during an ECMC in the $NVT$ ensemble, for different
values of the chain length $\ell$. When times are measured according to the
number of chains displaced, longer chains have a stronger impact on the system
relaxation. However, if one measures times in units of $Nt_{\text{move}}$,
thus counting individual particle displacements, we see that all the
relaxation curves now collapse. This collapse implies that the efficiency of
ECMC (in CPU time) is nearly independent of $\ell$.

Figure 7: Equilibrium correlation function for the ECMC dynamics, using
different values of the chain length $\ell$. Time is measured in units of (a)
the number of chains displaced, and (b) $Nt_{\text{move}}$.

## References

  * Berthier  _et al._ [2019] L. Berthier, P. Charbonneau, A. Ninarello, M. Ozawa, and S. Yaida, Zero-temperature glass transition in two dimensions, Nature communications 10, 1508 (2019).
  * [2] V. M. Bolton-Lum, R. C. Dennis, P. K. Morse and E. I. Corwin, in preparation.
  * Hagh  _et al._ [2022] V. F. Hagh, S. R. Nagel, A. J. Liu, M. L. Manning, and E. I. Corwin, Transient learning degrees of freedom for introducing function in materials, Proceedings of the National Academy of Sciences 119, e2117622119 (2022).
  * Bernard  _et al._ [2009] E. P. Bernard, W. Krauth, and D. B. Wilson, Event-chain monte carlo algorithms for hard-sphere systems, Phys. Rev. E 80, 056704 (2009).
  * Bernard [2011] E. Bernard,  _Algorithms and applications of the Monte Carlo method: Two-dimensional melting and perfect sampling_ , Ph.D. thesis, Université Pierre et Marie Curie-Paris VI (2011).
  * Frenkel and Smit [2001] D. Frenkel and B. Smit,  _Understanding molecular simulation: from algorithms to applications_ , Vol. 1 (Elsevier, 2001).
  * Flenner and Szamel [2015] E. Flenner and G. Szamel, Fundamental differences between glassy dynamics in two and three dimensions, Nature communications 6, 7392 (2015).
  * Ramasubramani  _et al._ [2020] V. Ramasubramani, B. D. Dice, E. S. Harper, M. P. Spellings, J. A. Anderson, and S. C. Glotzer, freud: A software suite for high throughput analysis of particle simulation data, Computer Physics Communications 254, 107275 (2020).
  * Newman and Barkema [1999] M. E. Newman and G. T. Barkema,  _Monte Carlo methods in statistical physics_ (Clarendon Press, 1999).
  * Henderson [1975] D. Henderson, A simple equation of state for hard discs, Molecular Physics 30, 971 (1975).
  * Santos  _et al._ [1999] A. Santos, S. B. Yuste, and M. L. De Haro, Equation of state of a multicomponent d-dimensional hard-sphere fluid, Molecular Physics 96, 1 (1999).