# Optical phonon contribution to the thermal conductivity of a quantum
paraelectric

Pankaj Bhalla  Beijing Computational Science Research Center, Beijing, 100193,
China  Nabyendu Das  Department of Physics, The LNM-Institute of Information
Technology, Jaipur 302031, India

###### Abstract

Motivated by recent experimental findings, we study the contribution of a
quantum critical optical phonon branch to the thermal conductivity of a
paraelectric system. We consider the proximity of the optical phonon branch to
transverse acoustic phonon branch and calculate its contribution to the
thermal conductivity within the Kubo formalism. We find a low temperature
power law dependence of the thermal conductivity as $T^{\alpha}$, with
$1<\alpha<2$, (lower than $T^{3}$ behavior) due to optical phonons near the
quantum critical point. This result is in accord with the experimental
findings and indicates the importance of quantum fluctuations in the thermal
conduction in these materials.

##  I Introduction

In the last few decades, the study of quantum paraelectrics such as SrTiO3 and
KTaO3 has become a topic of considerable interest [1, 2, 3, 4, 5, 6]. In
general, interest in the quantum fluctuation induced novel states is on rise
as it can potentially lead to technological advancements. The main feature of
these materials is a suppressed long range ferroelectric order down to the
zero temperature due to quantum fluctuations. It is found that the dielectric
constant saturates at a very high value on lowering the temperature. Thus
these materials in their pristine state may be considered to be on the verge
of a phase transition near $T=0$. In principle, the ferroelectric transition
can be induced by some non thermal parameters such as pressure, chemical
composition, etc. at a zero temperature. A continuous zero temperature
transition without any discontinuity in the order-parameter variations are
termed as quantum critical point. The “zero temperature transition” or a
“quantum critical point” are not just theoretical concepts with academic
interests only. Once the non-thermal parameters in a system are tuned to value
that lead to a zero temperature continuous transition, they determine the
finite temperature properties of a system in terms of quantum critical
fluctuations of the associated degrees of freedom. In such cases, the finite
temperature properties shows some universal behaviour which are independent of
many details of the microscopic interations [7]. A set of dielectric materials
such as SrTiO3, KTaO3 associated with such phenomena as a result of quantum
critical polarization fluctuations are coined as quantum paraelectrics. In
these materials (suppressed) ferroelectric transition is associated with a
nearly soft zone center optical mode and in the low temperature regime, many
of the thermodynamic and transport properties can be described by quantum
fluctuations of the associated optic branch and its interactions with other
degrees of freedom. Such novel dielectric behavior was first reported by
Müller et al. [8] in strontium titanate compounds. In recent times, pressure
induced quantum phase transitions and associated quantum critical behavior
have been experimentally confirmed [2]. A lot of theoretical and experimental
activities in this direction to describe the dielectric behavior of these
systems have been successfully taken place in the recent times [3, 9, 10, 11,
12, 13, 14, 15, 16, 17, 18].
To explore consequences of quantum fluctuations in these materials, the study
of heat transport properties is also equally important like their dielectric
behavior. This may shed further light on the behavior of quantum paraelectrics
at or near the quantum critical point. In an insulator heat conduction is
mainly driven by the flow of phonons. Optical phonons have higher energies
$\sim$ Debye temperature and they often lead to novel thermal conductivity
behavior at the said temperature regime [19, 20]. However, at low temperature
thermal conductivity is dominated by the flow of the acoustic phonons. In case
of a quantum paraelectric, zone center transverse optic mode is nearly soft
and the phonons in the respective transverse optical (TO) branch have lowered
energy compared to the usual insulators. This is to mention that LO phonons in
this materials still have higher energy owing to the long range dipolar
interaction [21].This scenario is observed in experiments [22, 23] and
schematically shown in Fig. 1. Thus it is plausible that optical phonons in a
quantum paraelectric shape the low temperature behavior in certain systems,
such as KTaO3.

Figure 1: A schematic description of proximity of the transverse acoustic (TA)
and transverse optical (TO) phonon modes associated with (incipient)
ferroelectric transition in a typical quantum paraelectric. Here
$\Delta=\Delta_{0}+aT^{2}$ is the temperature dependent gap between the TA and
TO phonons, where $T$ is the temperature, $a$ is material dependent parameter,
and $\Delta_{0}$ is the phonon energy gap at $T=0$.

Being high electronic energy gap semiconductors, electronic exciations are
gapped out at a low temperature and transport properties in these materials
are governed by the interactions between lowered energy optical phonons with
the acoustic phonons [24, 25, 26] and other degrees of freedom like
impurities, structural domains[16], etc.
The basic understanding of the phonon thermal conductivity is as follows. In a
classical kinetic theory estimate, the thermal conductivity $\kappa\sim
v_{\text{ph}}c_{V}\lambda_{\text{ph}}$, where $v_{\text{ph}}$, $c_{V}$, and
$\lambda_{\text{ph}}$ are the phonon velocity, phonon specific heat and the
phonon mean free path respectively [27, 28, 29]. Thermal conductivity being a
transport coefficient, depends on the relaxation of the phonons due to its
interactions with other degrees of freedom and or anharmonicities. The same is
encoded in the $\lambda_{\text{ph}}=v_{\text{ph}}\tau_{ph}$, $\tau_{ph}$ is
the relaxation time. At low temperature and for acoustic phonon dominated
scenario, $c_{V}\sim T^{3}$ and $\lambda_{\text{ph}}$ is determined by the
scattering with the impurities and is temperature independent. Thus
$\kappa\sim T^{3}$. Experimental findings on this materials also support such
$T^{3}$ behavior of the specific heat [30]. Thus any departure from the said
$T^{3}$ behavior of the thermal conductivity, especially in a regime where
some other interacting degrees of freedom become active, indicates the
presence of some new scattering mechanism which results in frequancy and
temperature dependence of the mean free path or the scattering rate.
A recent experiment shows a few such novel behavior of the thermal
conductivity in quantum paraelectrics. For SrTiO3 there is a stronger than
$T^{3}$ increase of the thermal conductivity at low temperature. Poiseuille
flow of the phonons has been argued as a possible explanation of this
behavior. On the other hand it has been reported that in case of KTaO3, there
is a weaker than $T^{3}$ behavior for the thermal conductivity [24]. The later
behavior is typically attributed to the mutual interactions between more than
one low lying modes [31]. These experimental results suggest that quantum
paraelectrics are wonderful playgrounds to explore them. We put forward such
an explanation in terms of the thermal conduction of the optical phonons and
its interaction with the acoustic phonons which are the relevant degrees of
freedom in this system and in the temperature regime of our interest.
Our main result can be summarized as follows. Within a minimal model with
coupled TA-TO phonons, it is shown that near the quantum critical point, the
temperature variation of the thermal conductivity due the nearly soft optical
phonon follows a power law $T^{\alpha}$, $2>\alpha>1$ and its numerical value
can increase an order of magnitude. The same is shown in Fig. 2.

Figure 2: Power law behavior of the optical phonon contribution to the thermal
conductivity near the quantum critical point as obtained in the present
formalism.

This paper is organized as follows. First, we provide the mathematical
description of the anharmonic effects in Sec. II by introducing the model
Hamiltonian terms of the order of cubic and quartic in phonon displacement
functions. Then, we derive the expressions of the decay constants and thermal
conductivity in Sec. III. In Sec. IV, we present the numerical results. In
Sec. V, we discuss our results and conclude.

##  II Phonon decay due to anharmonicity

In a perfect crystal, harmonic phonons do not decay due to their infinitely
long lived nature and lead to infinite thermal conductivity [32]. We need to
introduce symmetry permitted anharmonic interactions to study the finite
lifetime of the phonons and hence the associated transport porperties [33, 34]
In the present work, we also focus on the study of anharmonic effects on the
thermal conductivity of quantum paraelectrics.
For a system in which atoms are displaced from their equilibrium positions,
the total Hamiltonian can be written as

$\displaystyle H$ $\displaystyle=$ $\displaystyle
H_{\text{har}}+H_{\text{anhar}},$ (1)

where $H_{\text{har}}$ corresponds to a free phonon part of the Hamiltonian
and $H_{\text{anhar}}$ represents the corrections to the Harmonic
approximation. The Harmonic part in terms of the phonon coordinates can be
expressed like

$\displaystyle H_{\text{har}}$ $\displaystyle=$
$\displaystyle\sum_{\alpha\mathbf{q}}\frac{1}{2}\omega_{\mathbf{\alpha
q}}A_{\alpha\mathbf{q}}A_{\alpha\mathbf{-q}}.$ (2)

Here
$A_{\alpha\mathbf{q}}=a^{\dagger}_{\alpha\mathbf{q}}+a_{\alpha\mathbf{q}}$ is
the sum of the creation $(a^{\dagger}_{\alpha\mathbf{q}})$ and the
annihilation $(a_{\alpha\mathbf{q}})$ phonon operators,
$\omega_{\alpha\mathbf{q}}$ is the phonon dispersion for the $\alpha$-th
branch. The latter quantity for small wavelengths can have form

$\omega_{\alpha q}=\sqrt{{\Delta_{\alpha 0}}^{2}+c_{\alpha}^{2}q^{2}},$ (3)

having $\Delta_{\alpha 0}$ the phonon energy gap, $c_{\alpha}$ a sound
velocity and $q$ a phonon wave vector. The dispersion correspond to the phonon
branch can be associated with the gap. For the case $\Delta_{0}=0$, the
dispersion corresponds to acoustic phonon and for $\Delta_{0}\neq 0$ it
represents optical phonon. Since we are considering a minimal model with one
TA and one TO branch, we drop the branch index $\alpha$ for further
discussions and $\Delta_{0}$ represents the energy gap for the zone center TO
phonon at $T=0$. Further, the ferroelectric instability in materials having
anharmonic effects is determined with $q=0$ optical phonon or
$\omega_{q}=\Delta_{0}$. The ferroelectric quantum critical point, it is
determined by $\Delta_{0}=0$ at $T=0$ for the optic branch. At this quantum
critical point, the phonon energy gap for the optical branch shows quadratic
temperature dependent phonon dispersion at finite temperature such as
$\Delta=\Delta_{0}+aT^{2}$ having $a$ as a material dependent parameter [2, 1]
as depicted in Fig. 1.

Figure 3: a): Lowest order three phonon process due to the cubic anharmonic
term where the quantum fluctuations of the relevant optic branch contribute
the most and b) Second order four phonon process due to the quartic anharmonic
term. Here the dotted and the dashed curves represent two different kind of
phonons.

In the anharmonic part $H_{\text{anhar}}$ of the total Hamiltonian, we
consider two lowest order terms of the following forms.

$\displaystyle H_{3}$ $\displaystyle=$
$\displaystyle\sum_{\alpha\beta\gamma\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}V_{\alpha\beta\gamma}(\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3})A_{\alpha,\mathbf{k}_{1}}A_{\beta,\mathbf{k}_{2}}A_{\beta,\mathbf{k}_{3}}$
$\displaystyle H_{4}$ $\displaystyle=$
$\displaystyle\sum_{\alpha\beta\gamma\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}}V_{\alpha\beta\gamma\delta}(\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3},\mathbf{k}_{4})A_{\alpha,\mathbf{k}_{1}}A_{\beta,\mathbf{k}_{2}}A_{\beta,\mathbf{k}_{3}}A_{\delta,\mathbf{k}_{4}}$

Here $V$ in both terms represents the anharmonic coefficients and the relative
weights of these coefficients depend on the symmetry of system. The
Hamiltonian $H_{3}$ corresponds to the three phonon process and $H_{4}$ to the
four phonon process. The lowest order processes form both types of
anharmonicities which involve momentum transfer between two kind of phonons
are shown in Fig. 3. Here the first order effect from the fourth order
anharmonicity does not associated with any net momentum transfer. It
contributes to the quadratic temperature dependence to the phonon energy gap
as shown in Fig. 4.

Figure 4: Mass enhancement due to the lowest order quartic anharmonic effect
which leads to Eq. 3 with renormalized gap $\Delta_{0}\sim T^{2}.$

###  II.1 Three Phonon process

For the three phonon process, decay rate or the inverse lifetime of a phonon
at the phonon frequency due to the cubic anharmonicity and at a temperature
$T$ is given as:

$\Gamma^{3}(\omega,T)=18\pi\sum_{\mathbf{k}_{1}\mathbf{k}_{2}\lambda_{1}\lambda_{2}}|V_{\alpha\beta}(-\mathbf{k}\lambda;\mathbf{k}_{1}\lambda_{1};\mathbf{k}_{2}\lambda_{2})|^{2}\bigg{\\{}(n_{\omega_{1}}-n_{-\omega_{2}})\delta_{12}^{s}+(n_{\omega_{1}}-n_{\omega_{2}})\delta_{12}^{d}\bigg{\\}}.$
(5)

Here $\omega_{i}\equiv\left\\{\omega(\mathbf{k}_{i}\lambda_{i})\right\\}$,
$n_{\omega_{i}}=(e^{\beta\omega_{i}}-1)^{-1}$ is the Bose-Einstein
distribution function having $\beta$ a inverse of the temperature, and the
dirac delta functions denote
$\delta_{12}^{s}=\delta_{\omega-\omega_{1}-\omega_{2}}-\delta_{\omega+\omega_{1}+\omega_{2}}$
and
$\delta_{12}^{d}=\delta_{\omega+\omega_{1}-\omega_{2}}-\delta_{\omega-\omega_{1}+\omega_{2}}$.
Considering the anharmonic coefficient
$|V(-\mathbf{k}\lambda;\mathbf{k}_{1}\lambda_{1};\mathbf{k}_{2}\lambda_{2})|^{2}=\frac{V_{0}^{2}}{\omega(-\mathbf{k}\lambda)\omega(\mathbf{k}_{1}\lambda_{1})\omega(\mathbf{k}_{2}\lambda_{2})}$[35]
and then converting the summations into integrals, we have

$\Gamma^{3}(\omega,T)=18\pi\frac{V_{0}^{2}}{\omega}\int
d\omega_{1}\rho_{1}^{\omega_{1}}\bigg{\\{}\rho_{2,\omega}^{\omega_{1}}(n_{\omega_{1}}-n_{-\omega+\omega_{1}})+\rho_{2,\omega}^{-\omega_{1}}(n_{\omega_{1}}-n_{\omega+\omega_{1}})-\text{terms
with}(\omega\rightarrow-\omega)\bigg{\\}},$ (6)

where
$\rho_{i,\omega}^{\pm\omega_{1}}=\rho_{i}(\omega\mp\omega_{1})/(\omega\mp\omega_{1})$
having $\rho_{i}(\omega)$ the density of states. For a certain phonon branch,
it takes a form $\rho_{i}({\omega})\sim$ $\omega\sqrt{\omega^{2}-\Delta^{2}}$
for $\omega>\Delta$ and zero in other cases. Here we consider
$\Delta=\Delta_{0}+aT^{2}$ with $a$ eV${}^{-1})$ as a constant for optic
phonon which is not explicitly calculated here and taken from other references
[1, 2] and shown in Fig. 1.
To simplify the expression of decay rate, we have used $\rho_{1}(\omega)$ for
optic phonon and $\rho_{2}(\omega)$ for acoustic phonon. This yields
$\Gamma^{3}(\omega,T)$

$\Gamma^{3}(\omega,T)=\frac{18\pi}{4\pi^{4}}\frac{V_{0}^{2}}{\omega}\int
d\omega_{1}\sqrt{\omega_{1}^{2}-\Delta^{2}}\bigg{\\{}(\omega-\omega_{1})\big{[}n_{\omega_{1}}-n_{-\omega+\omega_{1}}\big{]}+(\omega+\omega_{1})\big{[}n_{\omega_{1}}-n_{\omega+\omega_{1}}\big{]}\bigg{\\}}.$
(7)

This is an expression for the decay rate of the optical phonon due to three
phonon process.

###  II.2 Four Phonon process

Due to the fourth order anharmonicity, the phonon decay rate can be expressed
as

$\displaystyle\Gamma^{4}(\omega,T)=$ $\displaystyle
96\pi\sum_{\mathbf{k}_{1}\lambda_{1}}\sum_{\mathbf{k}_{2}\lambda_{2}}\sum_{\mathbf{k}_{3}\lambda_{3}}|V(-\mathbf{k}\lambda;\mathbf{k}_{1}\lambda_{1};\mathbf{k}_{2}\lambda_{2};\mathbf{k}_{3}\lambda_{3})|^{2}$
$\displaystyle\quad\bigg{\\{}(-n_{-\omega_{1}}n_{-\omega_{2}}n_{-\omega_{3}}-n_{\omega_{1}}n_{\omega_{2}}n_{\omega_{3}})\delta_{123}^{s}+3(n_{\omega_{1}}n_{-\omega_{2}}n_{-\omega_{3}}+n_{-\omega_{1}}n_{\omega_{2}}n_{\omega_{3}})\delta_{123}^{d}\bigg{\\}}$
, (8)

where $\delta_{123}$ represents the dirac delta function due to the
combination of different phonon frequencies. The corresponding factors with
superscript $s$ and $d$ are
$\delta_{\omega-\omega_{1}-\omega_{2}-\omega_{3}}-\delta_{\omega+\omega_{1}+\omega_{2}+\omega_{3}}$,
and
$\delta_{\omega+\omega_{1}-\omega_{2}-\omega_{3}}-\delta_{\omega-\omega_{1}+\omega_{2}+\omega_{3}}$
respectively.
Assuming the anharmonic coefficient for this process in a similar fashion to
the three phonon process
$|V(-\mathbf{k}\lambda;\mathbf{k}_{1}\lambda_{1};\mathbf{k}_{2}\lambda_{2};\mathbf{k}_{3}\lambda_{3})|^{2}=\frac{V_{1}^{2}}{\omega(-\mathbf{k}\lambda)\omega(\mathbf{k}_{1}\lambda_{1})\omega(\mathbf{k}_{2}\lambda_{2})\omega(\mathbf{k}_{3}\lambda_{3})}$
and simplifying the above equation, the the decay rate for four phonon process
can been written as

$\displaystyle\Gamma^{4}(\omega,T)=$ $\displaystyle
96\pi\frac{V_{1}^{2}}{\omega}\int d\omega_{1}\rho_{1}^{\omega_{1}}\int
d\omega_{2}\rho_{2}^{\omega_{2}}\bigg{\\{}\rho_{2}^{-\omega_{1}-\omega_{2}}(-n_{-\omega_{1}}n_{-\omega_{2}}n_{-\omega+\omega_{1}+\omega_{2}}-n_{\omega_{1}}n_{\omega_{2}}n_{\omega-\omega_{1}-\omega_{2}})$
$\displaystyle+3\rho_{2}^{\omega_{1}-\omega_{2}}(n_{\omega_{1}}n_{-\omega_{2}}n_{-\omega-\omega_{1}+\omega_{2}}+n_{-\omega_{1}}n_{\omega_{2}}n_{\omega+\omega_{1}-\omega_{2}})-\text{terms
with}(\omega\rightarrow-\omega)\bigg{\\}}$

Here we consider restrict us to the simplest model which captures the
essential aspects of the quantum fluctuations associated with the nearly soft
optic mode and consider $\rho_{1}(\omega)$ for optic phonon and other density
of states for the case of acoustic phonon branches.

##  III Thermal conductivity

In order to compute the thermal conductivity, we start with its expression
based on the Kubo approach [36] and is given by

$\displaystyle\kappa(\omega,T)$ $\displaystyle=$
$\displaystyle\frac{\kappa_{0}}{T}\sum_{\mathbf{k}}\omega^{2}(\mathbf{k})v^{2}_{\mathbf{k}}\int_{-\infty}^{\infty}d\omega^{\prime}f_{T}(\omega,\omega^{\prime})\mathcal{S}(\omega,\omega^{\prime}).$

Here
$f_{T}(\omega,\omega^{\prime})=\frac{n_{\omega^{\prime}}-n_{\omega+\omega^{\prime}}}{\omega}$
is the thermal weight and the spectral weight function
$\mathcal{S}(\omega,\omega^{\prime})$ is defined as
$\mathcal{S}(\omega,\omega^{\prime})=\mathcal{A}(\mathbf{k},\omega^{\prime})\mathcal{A}(\mathbf{k},\omega+\omega^{\prime})$.
Further in zero frequency limit, the above equation can be written as

$\kappa(\omega=0,T)=\frac{\kappa_{0}}{T}\sum_{\mathbf{k}}\omega^{2}(\mathbf{k})v^{2}_{\mathbf{k}}\int_{-\infty}^{\infty}d\omega^{\prime}\frac{\partial
n(\omega^{\prime})}{\partial\omega^{\prime}}\mathcal{A}^{2}(\mathbf{k},\omega^{\prime}).$
(11)

Here the spectral function $\mathcal{A}(\mathbf{k},\omega)$ is defined as the
imaginary part of the phonon propagator and multiplied by a factor $1/\pi$
i.e.,

$\displaystyle\mathcal{A}(\mathbf{k},\omega)$ $\displaystyle=$
$\displaystyle-\frac{1}{\pi}Im\mathcal{D}(\omega,\mathbf{k})$ (12)
$\displaystyle=$
$\displaystyle-\frac{1}{\pi}\frac{2\omega(\mathbf{k})\Gamma(\mathbf{k},\omega)}{\left(\omega^{2}-\omega^{2}(\mathbf{k})\right)^{2}+4\omega(\mathbf{k})^{2}\Gamma(\mathbf{k},\omega)^{2}}.$

Using the expressions for decay rate Eqs. (7) and (LABEL:eqn:DR4) and the
phonon dispersion in the above equation, the thermal conductivity Eq. (11)can
be calculated. This calculation is within linear response theory by beyond
relaxation time approximation. We calculate these expressions numerically and
discuss the behavior of relevant quantities in the next section.

(a) noonleline

(b) noonleline

Figure 5: (a). Normalized temperature dependent thermal decay rate as a
function of $T/\Theta_{D}$ having $\Theta_{D}$ a Debye’s temperature
corresponds to three (T) and four phonon processes. Here the solid blue curve
refers to the three phonon contribution and dotted black curve to the three
and four phonon (T+F) combination.(b). Decay rate at a fixed temperature as a
function of frequency.

##  IV Results

In Fig. 5(a), we plot the normalized thermal decay rate $\Gamma(T)/\Gamma_{0}$
using Eqs. 7 and LABEL:eqn:DR4 as a function of $T/\Theta_{D}$ having
$\Theta_{D}$ a Debye’s temperature at a gap $\Delta_{0}=0$. Here, we observe
that the scattering rate increases with temperature. Both three and phonon
processes contribute to rise in $\Gamma(T)/\Gamma_{0}$. However, the
contribution from four phonon process is quite small as shown by black dotted
curve in Fig. 5(b). Further it shows nonlinear behavior at temperature greater
than $0.1$ Debye temperature. In frequency dependent case, it has been
observed that the decay rate varies nonlinearly with the freqeuncy, shown in
Fig. 5(b).

(a) noonleline

(b) noonleline

Figure 6: Plot for the normalized thermal decay rate of the optical phonon as
a function of (a): scaled temperature $T/\Theta_{D}$, (b): frequency at
different gap values $\Delta_{0}$ such as $0.0$, $0.1$, and $1.0$ meV. Figure
7: Thermal conductivity of a paraelectric due to three (T) phonon process and
the combination of three and four (T+F) processes. The curve has been plotted
at the quantum critical point $\Delta_{0}=0$.

Refering Fig. 6(a) for the temperature dependent normalized scattering rate at
different gap values, we have following findings. Firstly, the increase in the
gap $\Delta_{0}$ suppresses $\Gamma(T)/\Gamma_{0}$. Physically, this arises
due to the less number of excitations of phonon from both acoustic and optical
branches. Secondly, on moving away from quantum critical point $\Delta_{0}=0$,
the contribution made by acoustic phonons decreases which results the less
decay rate. In the frequency dependent case Fig. 6(b), the variation of
scattering rate at different gap is not clearly visible due to the less
magnitude of gap as compared to frequency scale.
The thermal conductivity of a paraelectric material is presented in Fig. 7 for
different phonon processes. At low temperatures, four phonon process does not
show appreciable contribution to the thermal conductivity using fitting
parameters $a=0.01$ eV-1 and $V_{1}$ and $V_{0}=1000$eV-1. Our choices of
parameter is such that the maximum value of the optical phonon decay rates
which $\sim\Gamma_{0}\sim$ a few meV. This choice is in accord with
experimental findings on the spectral width of the optical phonons due to TA-
TO interactions [37, 22] and is well justified for the temperature regime of
our interest. Even with same or slightly higher values of $V_{1}$, the fourth
order anharmonicity gives a little contribution which indicates that third
order ahnarmonicity effects play the most significant role in the thermal
conduction as expected.

This is to mention that usually four phonon process is considered to be
important near the Debye temperature [38]. However we had an anticipation that
owing to the presence of a nearly soft optic phonon which has significant
amount of excitation even at low temperature, it might turnout to be
important. Our calculation rules out that possibility.

Here the plot of $\kappa(T)/\kappa_{0}$ vs $T/\Theta_{D}$ shows a peak around
a temperature

Figure 8: Plot for normalized thermal conductivity with $T/\Theta_{D}$ at
different values of gap. The blue curve corresponds to the quantum critical
scenario.

$T\approx 0.075\Theta_{D}$ and then decays towards the zero as the temperature
is lowered. Further, the broadening of the curve is decided by the decay rate
$\Gamma(T)$ whose magnitude depends on the strength of the anharmonic
coefficients. To look at its behavior away from a quantum critical point, we
plot thermal conductivity at different gap values in Fig. 8 as a function of
$T/\Theta_{D}$. We see that at points away from the quantum critical point,
the optical phonon does not conduct, resulting vanishing contribution to the
thermal conductivity at low temperatures. Also the peak shifts towards the
high temperatures as shown in Fig. 8 corresponding to $\Delta_{0}=0.1$, $1$
meV. These results are qualitatively in agreement with experimental results
[24] as discussed.

##  V Discussion and Conclusion

Our results can be understood as follows. Like acoustic phonons, near quantum
critical point optical phonon density is also strongly temperature dependent.
Thus the decay rate of the optical phonon is strongly temperature dependent as
plotted in Fig. 6(a). We may describe the temperature dependence as $\sim
T^{1+x}$ where $0<x<1$. A finite temperature scaling description should give
$c_{v}\sim T^{3}$ for a quantum critical optical phonon which is same as the
acoustic phonons. Considering the above and following the classical
expression, at low temperature, $\kappa\sim$ $c_{V}/\Gamma\sim$ $T^{3-1-x}$
$\sim T^{\alpha}$ $\,\,2>\alpha=(2-x)>1$. The same has been obtained in our
study and shown in the log plot for thermal conductivity in Fig. 2. Though we
do not have a closed analytic expression for the thermal conductivity, the
power law obtained here is universal in a sense that is determined by only
available energy scale $k_{B}T$, spatial dimension, symmetry of the order
parameter and the dynamical critical exponent $z=1$ and do not depend on the
specific values of the interaction parameters. Similar behavior has been
obtained in various model calculations on different quantum critical systems:
e.g. the loop contribution in 2D superconductor-diffusive metal with z=2,
$\kappa\sim\ln T$ [39].
Other important finding form the numerical evaluation of the expression for
thermal conductivity is near the critical point, thermal conductivity due to
optic phonons enhanced by an order of magnitude than when it is away from a
quantum critical point. Thus we propose that TA-Quantum critical TO scattering
introduces an additional channel for thermal conduction which may dominate
over other contributions in a certain temperature regime. As the related
experimental results suggests, the results seem to be relevant for KTaO3 [24],
where any structural transition is absent. Our analysis paves a way for
addressing quantum critical aspects of thermal conductivity and thermal hall
effects in quantum paraelectrics [26, 40] and may be extended to bicritical
scenarios in a certain multiferroic sytems [4, 41, 10].

## Acknowledgments

PB acknowledges the National Key Research and Development Program of China
(grant No. 2017YFA0303400), China postdoctoral science foundation (grant no.
2019M650461) and NSFC grant no. U1930402 for financial support.

## References

  * Das and Mishra [2009] N. Das and S. G. Mishra, Fluctuations and criticality in quantum paraelectrics, Journal of Physics: Condensed Matter 21, 095901 (2009).
  * Rowley  _et al._ [2014] S. E. Rowley, L. J. Spalek, R. P. Smith, M. P. M. Dean, M. Itoh, J. F. Scott, G. G. Lonzarich, and S. S. Saxena, Ferroelectric quantum criticality, Nature Physics 10, 367 (2014).
  * Pálová  _et al._ [2009] L. Pálová, P. Chandra, and P. Coleman, Quantum critical paraelectrics and the casimir effect in time, Phys. Rev. B 79, 075101 (2009).
  * Das [2012] N. Das, Quantum critical behavior of a magnetic quantum paraelectric, Physics Letters A 376, 2683 (2012).
  * Das [2013] N. Das, Effects of strain coupling and marginal dimensionality in the nature of phase transition in quantum paraelectrics, International Journal of Modern Physics B 27, 1350028 (2013).
  * Das [2014] N. Das, On the possibility of mixed phases in disordered quantum paraelectrics, Modern Physics Letters B 28, 1450167 (2014).
  * Sachdev [1999] S. Sachdev,  _Quantum Phase Transitions_ (Cambridge University Press, UK, 1999).
  * Müller and Burkard [1979] K. A. Müller and H. Burkard, SrTiO3: An intrinsic quantum paraelectric below 4 k, Phys. Rev. B 19, 3593 (1979).
  * Conduit and Simons [2010] G. J. Conduit and B. D. Simons, Theory of quantum paraelectrics and the metaelectric transition, Phys. Rev. B 81, 024102 (2010).
  * Narayan  _et al._ [2019] A. Narayan, A. Cano, A. V. Balatsky, and N. A. Spaldin, Multiferroic quantum criticality, Nature Materials 18, 223 (2019).
  * Li  _et al._ [2019] X. Li, T. Qiu, J. Zhang, E. Baldini, J. Lu, A. M. Rappe, and K. A. Nelson, Terahertz field–induced ferroelectricity in quantum paraelectric SrTiO3, Science 364, 1079 (2019).
  * Ang  _et al._ [2000] C. Ang, Z. Yu, and Z. Jing, Impurity-induced ferroelectric relaxor behavior in quantum paraelectric SrTiO3 and ferroelectric BaTiO3, Phys. Rev. B 61, 957 (2000).
  * Khaetskii  _et al._ [2020] A. Khaetskii, V. Juričič, and A. V. Balatsky, Thermal magnetic fluctuations of a ferroelectric quantum critical point, Journal of Physics: Condensed Matter 33, 04LT01 (2020).
  * Geirhos  _et al._ [2020] K. Geirhos, P. Lunkenheimer, M. Blankenhorn, R. Claus, Y. Matsumoto, K. Kitagawa, T. Takayama, H. Takagi, I. Kézsmárki, and A. Loidl, Quantum paraelectricity in the kitaev quantum spin liquid candidates H3LiIr2O6 and D3LiIr2O6, Phys. Rev. B 101, 184410 (2020).
  * He  _et al._ [2020] X. He, D. Bansal, B. Winn, S. Chi, L. Boatner, and O. Delaire, Anharmonic eigenvectors and acoustic phonon disappearance in quantum paraelectric SrTiO3, Phys. Rev. Lett. 124, 145901 (2020).
  * Kustov  _et al._ [2020] S. Kustov, I. Liubimova, and E. K. H. Salje, Domain dynamics in quantum-paraelectric SrTiO3, Phys. Rev. Lett. 124, 016801 (2020).
  * Nataf  _et al._ [2020] G. F. Nataf, M. Guennou, J. M. Gregg, D. Meier, J. Hlinka, E. K. H. Salje, and J. Kreisel, Domain-wall engineering and topological defects in ferroelectric and ferroelastic materials, Nature Reviews Physics 2, 634 (2020).
  * Sim  _et al._ [2021] S. Sim, H. Yang, H.-L. Kim, M. J. Coak, M. Itoh, Y. Noda, and J.-G. Park, Sizable suppression of thermal hall effect upon isotopic substitution in srtio3, Phys. Rev. Lett. 126, 015901 (2021).
  * Ye  _et al._ [2015] Z.-Q. Ye, B.-Y. Cao, W.-J. Yao, T. Feng, and X. Ruan, Spectral phonon thermal properties in graphene nanoribbons, Carbon 93, 915 (2015).
  * Zou  _et al._ [2019] J.-H. Zou, X.-T. Xu, and B.-Y. Cao, Size-dependent mode contributions to the thermal transport of suspended and supported graphene, Applied Physics Letters 115, 123105 (2019).
  * Coak  _et al._ [2019] M. J. Coak, C. R. S. Haines, C. Liu, G. G. Guzmán-Verri, and S. S. Saxena, Pressure dependence of ferroelectric quantum critical fluctuations, Phys. Rev. B 100, 214111 (2019).
  * Shirane  _et al._ [1967] G. Shirane, R. Nathans, and V. J. Minkiewicz, Temperature dependence of the soft ferroelectric mode in kta${\mathrm{o}}_{3}$, Phys. Rev. 157, 396 (1967).
  * Farhi  _et al._ [2000] E. Farhi, A. Tagantsev, R. Currat, B. Hehlen, E. Courtens, and L. Boatner, Low energy phonon spectrum and its parameterization in pure ktao3 below 80 k, EPJB 15, 615 (2000).
  * Martelli  _et al._ [2018] V. Martelli, J. L. Jiménez, M. Continentino, E. Baggio-Saitovitch, and K. Behnia, Thermal transport and phonon hydrodynamics in strontium titanate, Phys. Rev. Lett. 120, 125901 (2018).
  * Yang  _et al._ [2020] Y.-f. Yang, G.-M. Zhang, and F.-C. Zhang, Universal behavior of the thermal hall conductivity, Phys. Rev. Lett. 124, 186602 (2020).
  * Li  _et al._ [2020] X. Li, B. Fauqué, Z. Zhu, and K. Behnia, Phonon thermal hall effect in strontium titanate, Phys. Rev. Lett. 124, 105901 (2020).
  * Ziman [1960] J. M. Ziman,  _Electrons and Phonons_ (Clarendon Oxford, 1960).
  * Tritt [2005] T. M. Tritt,  _Conductivity Theory, Properties, and Applications_ (Springer, Berlin, 2005).
  * Berman [1976] R. Berman,  _Thermal Conduction in Solids_ (Clarendon Press, Oxford, 1976).
  * Bourgeal  _et al._ [1988] S. Bourgeal, R. Villar, S. Vieira, and V. A. Trepakov, Low temperature specific heat of ktao3 from 0.3k, Ferroelectrics 79, 237 (1988).
  * Gurzhi [1968] R. N. Gurzhi, Hydrodynamic effects in solids at low temperature, Phys. Usp. 11, 255 (1968).
  * Ashcroft and Mermin [1976] N. W. Ashcroft and N. D. Mermin,  _Solid state physics_ (Science: Physics (Saunders College), 1976).
  * Klemens [1966] P. G. Klemens, Anharmonic decay of optical phonons, Phys. Rev. 148, 845 (1966).
  * Chang and Jones [1962] G. K. Chang and R. E. Jones, Low-temperature thermal conductivity of amorphous solids, Phys. Rev. 126, 2055 (1962).
  * Peierls [1974] R. Peierls,  _Quantum Theory of Solids_ (Oxford University Press, 1974).
  * Mahan [1990] G. D. Mahan,  _Many-Particle Physics_ (Plenum, New York and London, 2nd. Ed., 1990).
  * Axe  _et al._ [1970] J. D. Axe, J. Harada, and G. Shirane, Anomalous acoustic dispersion in centrosymmetric crystals with soft optic phonons, Phys. Rev. B 1, 1227 (1970).
  * Glassbrenner and Slack [1964] C. J. Glassbrenner and G. A. Slack, Thermal conductivity of silicon and germanium from 3k to the melting point, Phys. Rev. 134, A1058 (1964).
  * Podolsky  _et al._ [2007] D. Podolsky, A. Vishwanath, J. Moore, and S. Sachdev, Thermoelectric transport near pair breaking quantum phase transition out of $d$-wave superconductivity, Phys. Rev. B 75, 014520 (2007).
  * Chen  _et al._ [2020] J.-Y. Chen, S. A. Kivelson, and X.-Q. Sun, Enhanced thermal hall effect in nearly ferroelectric insulators, Phys. Rev. Lett. 124, 167601 (2020).
  * Morice  _et al._ [2017] C. Morice, P. Chandra, S. E. Rowley, G. Lonzarich, and S. S. Saxena, Hidden fluctuations close to a quantum bicritical point, Phys. Rev. B 96, 245104 (2017).


# Multidimensional King-plot analysis for accurate extraction of Cd nuclear
charge radii: a challenge for nuclear structure theory

J. Z. Han  State Key Laboratory of Precision Measurement Technology and
Instruments, Key Laboratory of Photon Measurement and Control Technology of
Ministry of Education, Department of Precision Instrument, Tsinghua
University, Beijing 100084, China  C. Pan  State Key Laboratory of Nuclear
Physics and Technology, School of Physics, Peking University, Beijing 100871,
China  K. Y. Zhang  State Key Laboratory of Nuclear Physics and Technology,
School of Physics, Peking University, Beijing 100871, China  X. F. Yang  State
Key Laboratory of Nuclear Physics and Technology, School of Physics, Peking
University, Beijing 100871, China  S. Q. Zhang  State Key Laboratory of
Nuclear Physics and Technology, School of Physics, Peking University, Beijing
100871, China  J. C. Berengut  School of Physics, University of New South
Wales, Sydney 2052, Australia  S. Goriely  Institut d′Astronomie et
d′Astrophysique, CP-226, Universite Libré de Bruxelles, 1050 Brussels, Belgium
H. Wang  School of Economics and Management, Tsinghua University, Beijing
100084, China  Y. M. Yu<EMAIL_ADDRESS>Beijing National Laboratory for
Condensed Matter Physics, Institute of Physics, Chinese Academy of Sciences,
Beijing 100190, China  J. Meng<EMAIL_ADDRESS>State Key Laboratory of
Nuclear Physics and Technology, School of Physics, Peking University, Beijing
100871, China  J. W. Zhang<EMAIL_ADDRESS>State Key Laboratory of
Precision Measurement Technology and Instruments, Key Laboratory of Photon
Measurement and Control Technology of Ministry of Education, Department of
Precision Instrument, Tsinghua University, Beijing 100084, China  L. J. Wang
<EMAIL_ADDRESS>State Key Laboratory of Precision Measurement
Technology and Instruments, Key Laboratory of Photon Measurement and Control
Technology of Ministry of Education, Department of Precision Instrument,
Tsinghua University, Beijing 100084, China

###### Abstract

High-accuracy determination of isotope shift factors, which plays a crucial
role in accurate extraction of nuclear charge radius, is well-known to be
challenging experimentally and theoretically. Nonetheless, based on an
accurate measurement of the isotope shifts for the Cd+
$5s~{}^{2}S_{1/2}-5p~{}^{2}P_{3/2}$ and the Cd
$5s^{2}~{}^{1}S_{0}-5s5p~{}^{1}P_{1}$ transition, a multidimensional King-plot
analysis is performed to extract the atomic field shift and mass shift factors
accurately. The results are further cross-checked against results from atomic
structure calculations using a high-accuracy configuration interaction and
many-body perturbation theory. Combined with previous isotope shifts, nuclear
charge radii for 100-130Cd of the highest accuracy are reported. For the
neutron-rich region, accuracies for the charge radii are improved by nearly
one order of magnitude. This work provides a coherent and systematic
extraction of Cd nuclear charge radii from isotope shifts. The results offer
stringent testing and possible challenges to the latest advances in nuclear
theory and excellent benchmarking to the atomic structure calculations.

Introduction.— The nuclear charge radius is a fundamental property of an
atomic nucleus. It serves as an important probe in various nuclear phenomena,
such as the occurrence of magic numbers Koszorús  _et al._ (2021); Gorges  _et
al._ (2019); Goodacre  _et al._ (2021), shape staggering Marsh  _et al._
(2018) and evolution De Witte  _et al._ (2007); Yordanov  _et al._ (2012),
shape coexistence Yang  _et al._ (2016), proton Geithner  _et al._ (2008) and
neutron Sanchez  _et al._ (2006); Noertershaeuser  _et al._ (2009) halos, and
neutron skins Hagen  _et al._ (2016). Moreover, the subtle evolution of charge
radii along isotopic chains provides a stringent test and challenge for
nuclear models. Therefore, for both experimental and theoretical nuclear
physics, a highly accurate extraction of nuclear charge radii is paramount.
With the development of laser spectroscopy technology, the relative changes in
the nuclear charge radii may be precisely determined from isotope shifts (ISs)
Gebert  _et al._ (2015); Hammen  _et al._ (2018); Koszorús  _et al._ (2021).
Most radio-nuclide measurements currently rely on laser spectroscopy. However,
the extraction of the nuclear charge radius relies on atomic field shift
factors $F$ and mass shift factors $K$. At present, the accuracy of atomic
factors has become the most critical limiting factor in extracting the nuclear
charge radius and highlights the importance of determining them with high
accuracy and reliability.

Cd nuclei have attracted wide attention because of their special attributes:
i) Cd nuclei contain 48 protons, two proton holes below the $Z=50$ proton
magic number—nuclei near the shell having richer properties; ii) a radioactive
ion beam facility can provide Cd isotopes with mass numbers ranging from 100
to 130 Hammen  _et al._ (2018)— many isotopes of an atom provide considerable
sample information regarding the changes in nuclear properties with neutron
number; iii) in atomic physics, Cd/Cd+ systems have many spectral transitions
suitable for optical measurements taken using a Cd+ microwave clock Zhang  _et
al._ (2012); Wang  _et al._ (2013); Miao  _et al._ (2015) and sympathetic
cooling technology Zuo  _et al._ (2019); Han  _et al._ (2021)—here, the
measurement accuracy of the Cd+ ISs may be further improved. In the past, many
groups have measured the optical ISs of Cd isotopes and extracted nuclear
charge radii Gillespie  _et al._ (1975); Wenz  _et al._ (1981); Buchinger  _et
al._ (1987); Angeli and Marinova (2013); Fricke  _et al._ (2004); Yordanov
_et al._ (2013); Hammen  _et al._ (2018). However, the results extracted from
independent groups are not consistent (e.g. Refs. Fricke  _et al._ (2004);
Hammen  _et al._ (2018); Angeli and Marinova (2013)). For Cd nuclei, with
nucleon numbers near the magic numbers, $50$ and $82$, the accuracy of the
nuclear charge radii is relatively low Hammen  _et al._ (2018).

Here, for 100-130Cd nuclear charge radii, we present our latest results which
currently have the highest accuracy. The atomic $F$ and $K$ factors used to
extract the nuclear radii are obtained experimentally from a multidimensional
King-plot analysis (MKA) and were cross-checked using theoretical atomic
structure calculations. In experiments, laser-induced fluorescence (LIF) from
sympathetically-cooled large Cd+ ion crystal and saturated absorption
spectroscopy are used to measure the ISs in the Cd+
$5s~{}^{2}S_{1/2}\rightarrow 5p~{}^{2}P_{3/2}$ transition line (214.5 nm) and
the Cd $5s^{2}~{}^{1}S_{0}\rightarrow 5s5p~{}^{1}P_{1}$ transition line (228.8
nm). Combined with the IS results of the Cd
$5s5p~{}^{3}P_{2}-5s6s~{}^{3}S_{1}$ transition line (508.6 nm) and the nuclear
size parameter $\lambda$, MKA plots were used to extract the atomic factors.
From theory, high-accuracy atomic structure calculations were performed using
the configuration interaction and many-body perturbation theory (CI+MBPT)
approach to obtain the atomic factors for the typical Cd/Cd+ transitions.
These atomic factor results ensure the accuracy and reliability of our latest
calculations of the 100-130Cd nuclear charge radii, especially for those
nuclei with neutron numbers approaching magic numbers, $N=50$ and $82$. This
Letter reports on a more consistent way to determine the atomic factors and
extract the nuclear charge radii through IS calculations. High-accuracy data
of nuclear charge radii will prompt further testing and provide guidance of
nuclear theory.

Experiment.— The MKA starts from high-accuracy IS measurements of the
transitions for two Cd systems specifically, the Cd+ 214.5 nm and Cd 228.8 nm
transition lines. The ISs of the Cd+ 214.5 nm line were measured using LIF
spectroscopy. In a linear Paul trap, we prepared two-species crystals
consisting of $2\times 10^{5}$ Ca+ and $4\times 10^{5}$ Cd+ ions (see Ref. Han
_et al._ (2021) for details). The coolant ions, Ca+ ions, were Doppler-cooled
by lasers, and the Cd+ ions were sympathetically cooled to less than 0.5 K
through the Coulomb interaction between ions. The ISs of the transition line
for cooled Cd+ were measured by scanning the frequencies in the weak 214.5-nm
probe laser beam. The frequencies of the probe laser beam spans 5 GHz thus
permitting the measurement of all frequency resonances of the Cd+ isotopes of
even-numbered neutrons. The frequencies of the laser beam were measured using
a high-precision wavemeter, which was calibrated against an optical frequency
comb referenced to a hydrogen maser. The measured LIF spectrum is shown in
Fig. 1(a). The ISs of the Cd 228.8 nm line were measured through SAS using a
hollow cathode lamp. The 228.8-nm laser beam was split into two beams, a pump
and a probe, which were subsequently configured to overlap at the center of
the lamp. From the measured SAS data taken with a sample with a natural
abundance of Cd [Fig. 1(b)], the source of uncertainty of the IS results is
mainly statistical uncertainty. With many of the systematic effects being
common to all the neutron-even isotopes and the potential systematic shifts
obtained from differences of line frequencies, many potential systematic
shifts are common and cancel to a high degree. The measurement results of the
ISs of Cd+ 214.5 nm line and Cd 228.8 nm line are listed in Table 1. For the
former, our results are consistent with those of Hammen and colleagues Hammen
_et al._ (2018), although the accuracy is improved by 2–6 times. For the Cd
228.8 nm line, our results are to our best knowledge the first high-precision
IS measurements, being 30 times more accurate than earlier experimental
results Kelly and Tomchuk (1961).

Figure 1: Measurement results for (a) the Cd+ $5s~{}^{2}S_{1/2}\rightarrow 5p~{}^{2}P_{3/2}$ transition line (214.5 nm) fitted using Voigt profiles, and (b) the Cd $5s^{2}~{}^{1}S_{0}\rightarrow 5s5p~{}^{1}P_{1}$ transition line (228.8 nm) fitted using the Lorentz profile (signal) and the Doppler profile (background, not shown). For both plots, each measurement point is the average of three measurements, and each colored line represents the fitted curve of an isotope, and the red solid line represents the total fitted curve. Table 1: Measured values of the isotope shifts of the Cd+ $5s~{}^{2}S_{1/2}-5p~{}^{2}P_{3/2}$ (214.5 nm) and Cd $5s^{2}~{}^{1}S_{0}-5s5p~{}^{1}P_{1}$ (228.8 nm) transition lines compared with Ref. Hammen  _et al._ (2018) and Kelly and Tomchuk (1961).  $A^{\prime}$ |  $\delta\nu^{114,A^{\prime}}_{214.5}$ (MHz) |  $\delta\nu^{114,A^{\prime}}_{228.8}$ (MHz)
---|---|---
| This work | Hammen  _et al._ (2018) | This work | Kelly and Tomchuk (1961)
116 | -521.2(1.2) | -526.5(2.7) | -310.1(1.2) |
112 | 675.4(1.2) | 674.6(2.7) | 411.2(1.2) | 396(30)
110 | 1429.8(1.2) | 1432.2(3.8) | 874.2(1.2) | 906(35)
108 | 2187.8(1.2) | 2194.0(5.0) | 1339.0(2.1) |
106 | 2984.9(1.2) | 2991.1(6.6) | 1826.6(1.2) |

Multidimensional King-plot Analysis.— According to IS theory BreÍt (1958);
King (2013), the relationship between isotope shift
$\delta\nu^{A,A^{\prime}}=\nu^{A}-\nu^{A^{\prime}}$ and nuclear size $\lambda$
is

$\mu\delta\nu^{A,A^{\prime}}=F\cdot\mu\lambda^{A,A^{\prime}}+K,$ (1)

where $\mu=m_{A}\cdot m_{A^{\prime}}/(m_{A}-m_{A^{\prime}})$. The nuclear size
data are obtained from muonic atom spectroscopy and electron scattering
experiments. The KA procedure introduces the $F$ and $K$ factors determined
from the slope and intercept of the fitted curve of the
$\mu\delta\nu^{A,A^{\prime}}$ and $\mu\lambda^{A,A^{\prime}}$ data. This
procedure mitigates inaccuracies of the atomic model calculation. However, it
strongly relies on the precision of $\lambda^{A,A^{\prime}}$ obtained through
muonic atom spectroscopy and electron scattering experiments. Previously,
based on the ISs of the Cd+ 214.5 nm line and the $\lambda$ input, Hammen and
colleagues performed a two-dimensional King-plot analysis (2D KA) and obtained
values of $K$ and $F$ for the Cd+ 214.5-nm transition line, but with large
uncertainties Hammen  _et al._ (2018). These uncertainties have become the
main source of uncertainty in Cd nuclear charge radii. In the MKA method, we
used the ISs of the Cd+ 214.5-nm line, of the Cd 228.8-nm line, of Cd 508.6-nm
line Hammen  _et al._ (2018), and values of the nuclear size $\lambda$ Fricke
_et al._ (2004) as input to reduce the uncertainties of the extracted $K$ and
$F$ factors. In the MKA method, the $K$ and $F$ factors are constrained by the
assumption underlying Eq. (1). There are four experimental data points
corresponding to isotope pairs
$(A,A^{\prime})=(110,108),(112,110),(114,112),(116,114)$, which should lie on
a straight line in the four-dimensional isotope space formed from
$\mu\delta\nu^{A,A^{\prime}}_{214.5}$, $\mu\delta\nu^{A,A^{\prime}}_{228.8}$,
$\mu\delta\nu^{A,A^{\prime}}_{508.6}$, and $\mu\lambda^{A,A^{\prime}}$. The
objective of the MKA method is to ensure all correlations and uncertainties
are made correctly thereby enabling the $K$ and $F$ factors to be extracted
more accurately and reliably.

The procedures of the MKA method are as follows: i) The uncertainties of the
experimental data points $D_{1},\cdots,D_{k}$ in the multidimensional space
are assumed to be independent of each other and comply with the multivariate
normal distribution. The covariance matrix of the multivariate normal
distribution are the corresponding measurement uncertainties Gebert  _et al._
(2015); ii) Two arbitrary Monte Carlo points are generated within the
1$\sigma$ range of two selected experimental points (e.g. $D_{1}$ and
$D_{k}$), and a multidimensional straight line $l$ is generated using these
two Monte Carlo points; iii) Projection points (labeled $P_{1},\cdots,P_{k}$)
for each data point ($D_{1},\cdots,D_{k}$) on the straight line $l$ are
generated. The probability density of each data point at the corresponding
projection points $f_{D_{1}}(P_{1}),\cdots,f_{D_{k}}(P_{k})$ is obtained; iv)
The product of the probability densities of each projection point,
$f=f_{D_{1}}(P_{1})\cdot f_{D_{2}}(P_{2})\cdot\cdots\cdot f_{D_{k}}(P_{k})$,
is used as the selection criterion. We generate $N$ sets of Monte Carlo
samples (in this work, we chose $N=3\times 10^{7}$) and selected $M$ sets of
samples that meet the condition $f/f_{max}>1\sigma$, where $f_{max}$ is the
maximum probability density of all $N$ sets of Monte Carlo samples; v) The
selected samples generate $M$ sets of multidimensional straight lines $l$ that
were then used to extract the $K$ and $F$ factors. These factors are thus
obtained from the probability density of the data points and have
uncertainties with 1$\sigma$ range.

The MKA provides model-independent atomic factors for the Cd+ 214.5-nm, Cd
228.8-nm and Cd 508.6-nm transition lines (see Table 2). Comparing the results
of the MKA and 2D KA, and with the addition of two other dimensions, the
accuracy of the $K$ and $F$ factors was improved by about 3.5 times, thereby
confirming that the MKA method functions as claimed. The MKA method is
expandible to higher dimensions and may be applied to elements with at least
three measured nuclear size parameters $\lambda$. More high-accuracy IS
measurements of other transition lines in the Cd ions may further improve the
accuracy of both $K$ and $F$ in the MKA. Comparing the 2D KA results with
those in Ref. Hammen  _et al._ (2018), we find that using the more accurate IS
results of the Cd+ 214.5-nm transition line, the extracted atomic factors
become twice as accurate, highlighting the importance of further improving the
accuracy of the IS measurements for stable nuclei. Comparing the MKA results
with those of Hammen and colleagues, the accuracy of the atomic factors is
improved six-fold, ensuring that we can extract Cd nuclear charge radii of
greater accuracy. The MKA results also provide excellent reference data for
benchmarking and cross-checking atomic structure calculations.

Table 2: Field shift ($F$) and mass shift ($K$) factors for Cd+ and Cd obtained by the MKA method, 2D KA method, and the CI+MBPT calculations. The atomic factors for several typical transition lines of the Cd/Cd+ systems obtained from CI+MBPT calculations are also listed. $F$ (MHz/fm2) |  $K$ (GHz u) | Source
---|---|---
Cd+ $5s~{}^{2}S_{1/2}-5p~{}^{2}P_{3/2}$ (214.5 nm)
-6192(320) | 1851(315) | MKA (This work)
-6176(1107) | 1831(1093) | 2D KA (This work)
-6144(300) | 1667(300) | CI+MBPT (This work)
-6260(1860) | 1860(1920) | 2D KA Hammen  _et al._ (2018)
Cd $5s^{2}~{}^{1}S_{0}-5s5p~{}^{1}P_{1}$ (228.8 nm)
-4140(214) | 1491(211) | MKA (This work)
-4129(741) | 1477(731) | 2D KA (This work)
-4024(200) | 1428(220) | CI+MBPT (This work)
Cd $5s5p~{}^{3}P_{2}-5s6s~{}^{3}S_{1}$ (508.6 nm)
1194(64) | -20(63) | MKA (This work)
1216(226) | -38(222) | 2D KA (This work)
1228(60) | -63(400) | CI+MBPT(This work)
Cd+ $5s~{}^{2}S_{1/2}-5p~{}^{2}P_{1/2}$ (226.5 nm)
-6067(300) | 1770(300) | CI+MBPT (This work)
-6077(734) | 1468(734) | semi-empirical Bishop and King (1971)
-6174(778) | 558(385) | semi-empirical Bauche  _et al._ (1985)
Cd $5s^{2}~{}^{1}S_{0}-5s5p~{}^{3}P_{1}$ (326.1 nm)
-4559(230) | 1865(400) | CI+MBPT (This work)
-4420 | 1717(330) | 2D KA Fricke  _et al._ (2004)
-3910(460) | 876(230) | semi-empirical Angeli and Marinova (2013)
-3900(460) | 809(407) | semi-empirical Buchinger  _et al._ (1987)

Atomic structure calculation.— The calculations of the $F$ and $K$ factors set
strong challenges for many-body atomic theory. The mass shift, in particular,
has proven very difficult to calculate even for single-valence-electron Cd+
isotopes. We used the CI+MBPT method, implemented in the AMBiT software Kahl
and Berengut (2019), to calculate the atomic factors. The CI+MBPT calculations
start from a self-consistent Dirac–Fock–Breit (DFB) calculation employing the
potential of a closed-shell Pd-like core. Next, a large valence basis of one-
particle orbitals is generated by diagonalizing a set of B-splines over the
one-electron DFB operator. Then, a set of many-electron configurations for the
CI expansion are constructed with these orbitals. For neutral Cd, the CI
configurations include all valence excitations up to level $20spdf$. For Cd+,
we also need to access the hole states, so we employ a particle–hole CI+MBPT
Berengut (2016), for which the core $4d$ shell is unfrozen and included in the
CI configurations. The configuration set is generated from all single and
double excitations up to $10spdf$ from leading configurations $5s$, $5p$,
$5d$, and $4d^{-1}\,5s^{2}$. Correlations with core shells that are not
included in the CI expansion are accounted for using the MBPT for both Cd and
Cd+. This includes excitations to virtual orbitals up to $30spdfg$.

The finite-field approach is implemented in the $F$ and $K$ factors
calculation of the AMBiT software Berengut  _et al._ (2003). Therein, the
operators for $F$ and $K$ are added to the Coulomb potential from the very
beginning of the DFB calculation. An IS scaling factor $\Gamma$ ahead of the
$F$ and $K$ operators is varied, and the level energies $E$ are calculated as
a function of $\Gamma$. The gradient of an $E$ versus $\Gamma$ graph then
gives us the field and mass shift factors. Values of the $F$ and $K$ factors
for Cd+ and Cd obtained by the CI+MBPT calculations are listed in Table 2. For
the Cd+ 214.5-nm, Cd 228.8-nm, and Cd 508.6-nm transition lines, the CI+MBPT
results show excellent consistency with the MKA results, thus highlighting the
high accuracy and reliability of this method. The CI+MBPT method has also
checked for accuracy and reliability in calculations of the atomic factors of
other neutral atoms and monovalent ions systems (see for example Ca+ Gebert
_et al._ (2015); Solaro  _et al._ (2020); Müller  _et al._ (2020) and Yb+
Counts  _et al._ (2020)). Cross-checking of the atomic factors obtained using
the CI+MBPT and MKA further confirms the reliability of our results.

Table 3: Extracted changes in the mean-square charge radii $\delta\langle r^{2}\rangle$ (fm2) relative to 114Cd and the nuclear charge radius $R_{ch}$ (fm) of 100-130Cd. We used the Cd+ 214.5-nm IS data of this work for those isotopes tagged by a superscripted $*$. For those isotopes tagged by a $\dagger$ and those untagged isotopes, we used respectively the Cd 508.6-nm IS data and the Cd+ 214.5-nm IS data of Hammen and colleagues Hammen  _et al._ (2018). The various uncertainties associated with the $\delta\langle r^{2}\rangle$ and $R_{ch}$ values are given in parentheses: the first relate to uncertainties arising from the IS measurements, the second relate to systematic errors stemming from the calculation of the atomic $F$ and $M$ factors, and the third for the $R_{ch}$ values are the uncertainties originating from the 114Cd measurements. The isotope-dependent linear transformation parameters $\alpha$ (fm2 u) and the corresponding isotope-dependent mass shift factors $K_{\alpha}$ (GHz u) are also listed.  N | $\alpha$ | $K_{\alpha}$ | $\delta\langle r^{2}\rangle$ | $R_{ch}$
---|---|---|---|---
52 | 1061 | -4719(34) | -1.4350(20)(85) | 4.45372(22)(96)(100)
53 | 1061 | -4719(34) | -1.3196(18)(78) | 4.46665(20)(88)(100)
54 | 1036 | -4564(29) | -1.1533(16)(57) | 4.48522(18)(63)(100)
55 | 1040 | -4589(30) | -1.0551(15)(54) | 4.49616(17)(60)(100)
56 | 1017 | -4447(26) | -0.9105(14)(39) | 4.51221(15)(43)(100)
57 | 1015 | -4434(26) | -0.8302(12)(40) | 4.52110(13)(44)(100)
58∗ | 1006 | -4378(25) | -0.6991(2)(28) | 4.53558(3)(31)(100)
59 | 1026 | -4502(27) | -0.6298(10)(29) | 4.54321(11)(32)(100)
60∗ | 1004 | -4366(24) | -0.5131(2)(20) | 4.55604(3)(22)(100)
61 | 1034 | -4552(28) | -0.4488(7)(22) | 4.56308(8)(24)(100)
62∗ | 1003 | -4360(24) | -0.3355(2)(13) | 4.57548(3)(14)(100)
63 | 1089 | -4892(41) | -0.2911(5)(21) | 4.58033(6)(23)(100)
64∗ | 990 | -4279(24) | -0.1603(2)(6) | 4.59459(3)(7)(100)
65 | 1219 | -5697(79) | -0.1160(4)(14) | 4.59941(4)(15)(100)
66 | Reference | 0(0)(0) | 4.61200(0)(0)(100)
67 | 758 | -2843(76) | 0.0418(5)(13) | 4.61653(5)(14)(100)
68∗ | 917 | -3827(32) | 0.1331(2)(9) | 4.62640(3)(10)(100)
69† | 855 | 1000(9) | 0.1700(20)(23) | 4.63040(22)(25)(100)
70† | 884 | 1035(8) | 0.2424(21)(25) | 4.63821(23)(27)(100)
71† | 858 | 1004(9) | 0.2811(29)(38) | 4.64238(31)(41)(100)
72† | 865 | 1012(9) | 0.3405(20)(43) | 4.64876(22)(46)(100)
73 | 847 | -3394(50) | 0.3737(11)(55) | 4.65233(11)(59)(100)
74† | 849 | 993(10) | 0.4276(25)(63) | 4.65813(27)(67)(100)
75 | 834 | -3313(53) | 0.4553(12)(76) | 4.66110(13)(81)(100)
76 | 838 | -3338(52) | 0.5081(12)(81) | 4.66676(13)(87)(100)
77 | 822 | -3239(57) | 0.5295(14)(97) | 4.66906(15)(104)(100)
78 | 828 | -3276(55) | 0.5823(15)(102) | 4.67471(16)(109)(100)
79 | 811 | -3171(60) | 0.5941(16)(121) | 4.67597(18)(129)(100)
80 | 821 | -3233(57) | 0.6561(17)(122) | 4.68259(18)(130)(100)
81 | 790 | -3041(66) | 0.6316(20)(153) | 4.67997(21)(163)(100)
82 | 804 | -3128(62) | 0.6992(20)(151) | 4.68719(21)(161)(100)

Nuclear Charge Radii of 110-130Cd.— Using the MKA method and our measured IS
of the Cd+ 214.5-nm and Cd 228.8-nm transition lines, the 100-130Cd nuclear
charge radii can be extracted with greater accuracy. An isotope-dependent
linear transformation
$\mu\lambda_{\rm{Cd}}\rightarrow(\mu\lambda_{\rm{Cd}}-\alpha)$ is used to
reduce the extracted uncertainties Hammen  _et al._ (2018). For each isotope,
the nuclear size factor $\mu\lambda_{\rm{Cd}}$ is translated to the position
at which the slope and intercept of the fitted curve have the least
correlation. Then, the MKA method is used to extract the mass shift factor
$K_{\alpha}$ from the set of data points
($\mu\delta\nu_{214.5},~{}\mu\delta\nu_{228.8},~{}\mu\delta\nu_{508.6},~{}\mu\lambda_{\rm{Cd}}-\alpha$).
The changes in the nuclear charge radii $\delta\langle r^{2}\rangle$ are
subsequently extracted using

$\displaystyle\lambda^{A,A^{\prime}}$ $\displaystyle=$
$\displaystyle\frac{\delta\nu_{IS}^{A,A^{{}^{\prime}}}-K_{\alpha}/\mu^{A,A^{{}^{\prime}}}}{F}+\frac{\alpha}{\mu^{A,A^{{}^{\prime}}}}$
(2) $\displaystyle=$ $\displaystyle\delta\langle
r^{2}\rangle^{A,A^{{}^{\prime}}}+\frac{C_{2}}{C_{1}}\delta\langle
r^{4}\rangle^{A,A^{{}^{\prime}}}+\frac{C_{3}}{C_{1}}\delta\langle
r^{6}\rangle^{A,A^{{}^{\prime}}},$

where $C_{1,2,3}$ denote the coefficients associated with the contributions of
the IS from higher radial moments Seltzer (1969). For the radioactive
isotopes, a contribution of $-2.7\%$ is assumed Fricke  _et al._ (2004).

For the naturally stable isotopes with $A=$110, 112, 114, 116, and 118, the
new measurements of the IS for the Cd+ 214.5-nm transition line have improved
accuracies; for the other isotopes of Cd, the measurements by Hammen and
colleagues give IS values of sufficient precision Hammen  _et al._ (2018).
These results lead to a new determination of the 100-130Cd nuclear charge
radii based on our $F$ and $K$ values obtained by the MKA method. The
recommended values of $\delta\langle r^{2}\rangle$ and nuclear charge radii
$R_{ch}$ are listed in Table 3. The isotope-dependent linear transformation
used in extracting the nuclear charge radius from the IS data eliminates the
quadratic trend in the uncertainty, and the more accurate values of $F$ and
$K$ yield more precise determination of $R_{ch}$, especially for isotopes from
the neutron-deficient ($N<60$) and neutron-rich ($N>70$) regions (Fig. 2). The
results reported in this Letter agree with those of previous work lying well
within the errors of the latter and with much improved accuracy. For the
neutron-rich region, the charge-radii accuracy has been improved by nearly one
order of magnitude. The improved accuracy highlights the odd-even staggering
(OES) of the charge radii in Cd isotopes for $N>70$, the phase of which was
uncertain in other work.

Figure 2: Comparison of the accuracy for the 100-130Cd nuclear charge radii of
this work and the results in Ref. Hammen  _et al._ (2018). Theoretical
predictions obtained using Fayans functionals Fy($\Delta r$) and the
Skyrme–Hartree–Fock–Bogoliubov model HFB-24 are shown for comparison.

The present high-precision results of the Cd charge radii set a formidable
benchmark for the latest advances in nuclear theory. In Fig. 2, we include two
modern theoretical predictions for comparison, the Fayans functional
Fy($\Delta r$) Reinhard and Nazarewicz (2017); Hammen  _et al._ (2018) and one
of best fits of the Skyrme–Hartree–Fock–Bogoliubov mass formulas, HFB-24
Goriely  _et al._ (2013). Note that other predictions from different density
functional calculations have been compared in Ref. Hammen  _et al._ (2018),
but most fail to reproduce the isotopic trend as a whole and the OES in
detail. The Fayans functional, which includes gradient terms in the surface
term and the pairing functional Fayans  _et al._ (2000); Reinhard and
Nazarewicz (2017) reproduce satisfyingly the experimental data. For this
parametrization, which is tuned with the Ca chain, the functional performs
well for the Cd chain, capturing not only the isotopic trend but also the
correct phase of the OES Hammen  _et al._ (2018).

Nevertheless, the Fy($\Delta r$) results for the neutron-rich Cd nuclei near
the magic number $N=82$, which were within the errors of previous experimental
data, are now beyond the high-accuracy results obtained in this work. As for
the K isotopes, the amplitudes of the OES in general are also overestimated by
Fy($\Delta r$), Koszorús  _et al._ (2021). The HFB-24, in which the strength
in contact pairing for neutrons (protons) depends on the neutron and proton
densities, provides a remarkable quantitative description of the experimental
data. Nevertheless, some apparent discrepancies exist in the OES, e.g., the
opposite phase near $N=80$. Other best Skyrme HFB-n mass formulas, such as
HFB-25, HFB-29, and HFB-31 Goriely  _et al._ (2013); Goriely (2015); Goriely
_et al._ (2016), show similar performances as HFB-24 and are thus not shown
explicitly in Fig. 2. Therefore, the present Cd charge radii with its improved
accuracy, particularly for the neutron-rich region, are set to provide
stringent tests and challenges for advanced nuclear models.

Conclusion.— High accuracy measurements of the IS of the Cd+
$5s~{}^{2}S_{1/2}-5p~{}^{2}P_{3/2}$ (214.5 nm) and the Cd
$5s^{2}~{}^{1}S_{0}-5s5p~{}^{1}P_{1}$ (228.8 nm) transition lines were
reported. Combined with the nuclear size parameter $\lambda$ and the ISs of
the Cd $5s5p~{}^{3}P_{2}-5s6s~{}^{3}S_{1}$ (508.6 nm) transition line, a MKA
was performed to extract the $F$ and $K$ factors of these transition lines.
The MKA method placed strong constraints on linear fits, and the uncertainties
for the $F$ and $K$ factors are well-determined by the multivariate normal
distribution. For the Cd+ $5s~{}^{2}S_{1/2}-5p~{}^{2}P_{3/2}$ transition line,
the accuracy of the extracted $F$ and $K$ factors from the MKA were improved
six-fold compared with previous results and furthermore were cross-checked in
high-accuracy CI+MBPT calculations. From the improved values for $F$ and $K$,
more accurate values of the nuclear charge radii of 100-130Cd were extracted.
The results are consistent with previous work, being well within previous
errors but with a greater accuracy, especially for isotopes of the neutron-
rich and neutron-deficient regions. The high-precision nuclear charge radii of
100-130Cd reported in this Letter set stringent constraints and pose
challenges to current nuclear theory, in which state-of-the-art density
functional calculations, HFB-24 and Fy($\Delta r$), fail to reproduce
perfectly currently-accepted results. In addition, accurate nuclear charge
radii may help in future investigations of Cd shell closure when extending the
neutron number up to magic number $N=82$ if new neutron-rich Cd isotopes are
fabricated.

Acknowledgments.— We thank V. A. Dzuba, W. Nörtershäuser, D. T. Yordanov, P.
O. Schmidt and T. Li for helpful discussions. This work is supported by the
National Key R&D Program of China (Nos. 2016YFA0302100, 2017YFE0116700,
2018YFA0404400), National Natural Science Foundation of China (Nos. 12073015,
11935003, 11875075) and Beijing Natural Science Foundation (1202011).

## References

  * Koszorús  _et al._ (2021) Á. Koszorús, X. Yang, W. Jiang, S. Novario, S. Bai, J. Billowes, C. Binnersley, M. Bissell, T. Cocolios, B. Cooper, _et al._ , Nature Physics 17, 439 (2021).
  * Gorges  _et al._ (2019) C. Gorges, L. Rodríguez, D. Balabanski, M. Bissell, K. Blaum, B. Cheal, R. G. Ruiz, G. Georgiev, W. Gins, H. Heylen, _et al._ , Physical review letters 122, 192502 (2019).
  * Goodacre  _et al._ (2021) T. D. Goodacre, A. Afanasjev, A. Barzakh, B. Marsh, S. Sels, P. Ring, H. Nakada, A. Andreyev, P. Van Duppen, N. Althubiti, _et al._ , Physical Review Letters 126, 032502 (2021).
  * Marsh  _et al._ (2018) B. Marsh, T. D. Goodacre, S. Sels, Y. Tsunoda, B. Andel, A. Andreyev, N. Althubiti, D. Atanasov, A. Barzakh, J. Billowes, _et al._ , Nature Physics 14, 1163 (2018).
  * De Witte  _et al._ (2007) H. De Witte, A. Andreyev, N. Barré, M. Bender, T. E. Cocolios, S. Dean, D. Fedorov, V. Fedoseyev, L. M. Fraile, S. Franchoo, _et al._ , Physical review letters 98, 112502 (2007).
  * Yordanov  _et al._ (2012) D. Yordanov, M. Bissell, K. Blaum, M. De Rydt, C. Geppert, M. Kowalska, J. Krämer, K. Kreim, A. Krieger, P. Lievens, _et al._ , Physical review letters 108, 042504 (2012).
  * Yang  _et al._ (2016) X. Yang, C. Wraith, L. Xie, C. Babcock, J. Billowes, M. Bissell, K. Blaum, B. Cheal, K. Flanagan, R. G. Ruiz, _et al._ , Physical review letters 116, 182502 (2016).
  * Geithner  _et al._ (2008) W. Geithner, T. Neff, G. Audi, K. Blaum, P. Delahaye, H. Feldmeier, S. George, C. Guénaut, F. Herfurth, A. Herlert, _et al._ , Physical review letters 101, 252502 (2008).
  * Sanchez  _et al._ (2006) R. Sanchez, W. Nörtershäuser, G. Ewald, D. Albers, J. Behr, P. Bricault, B. A. Bushaw, A. Dax, J. Dilling, M. Dombsky, _et al._ , Physical review letters 96, 033002 (2006).
  * Noertershaeuser  _et al._ (2009) W. Noertershaeuser, D. Tiedemann, M. Zakova, Z. Andjelkovic, R. Cazan, J. Kraemer, A. Krieger, R. Neugart, K. Blaum, D. Yordanov, _et al._ , Physical Review Letters 102 (2009).
  * Hagen  _et al._ (2016) G. Hagen, A. Ekström, C. Forssén, G. Jansen, W. Nazarewicz, T. Papenbrock, K. Wendt, S. Bacca, N. Barnea, B. Carlsson, _et al._ , Nature Physics 12, 186 (2016).
  * Gebert  _et al._ (2015) F. Gebert, Y. Wan, F. Wolf, C. N. Angstmann, J. C. Berengut, and P. O. Schmidt, Physical review letters 115, 053003 (2015).
  * Hammen  _et al._ (2018) M. Hammen, W. Nörtershäuser, D. Balabanski, M. Bissell, K. Blaum, I. Budinčević, B. Cheal, K. Flanagan, N. Frömmgen, G. Georgiev, _et al._ , Physical review letters 121, 102501 (2018).
  * Zhang  _et al._ (2012) J. Zhang, Z. Wang, S. Wang, K. Miao, B. Wang, and L. Wang, Physical Review A 86, 022523 (2012).
  * Wang  _et al._ (2013) S. Wang, J. Zhang, K. Miao, Z. Wang, and L. Wang, Optics express 21, 12434 (2013).
  * Miao  _et al._ (2015) K. Miao, J. Zhang, X. Sun, S. Wang, A. Zhang, K. Liang, and L. Wang, Optics letters 40, 4249 (2015).
  * Zuo  _et al._ (2019) Y. Zuo, J. Han, J. Zhang, and L. Wang, Applied Physics Letters 115, 061103 (2019).
  * Han  _et al._ (2021) J. Han, H. Qin, N. Xin, Y. Yu, V. Dzuba, J. Zhang, and L. Wang, Applied Physics Letters 118, 101103 (2021).
  * Gillespie  _et al._ (1975) W. Gillespie, S. Brain, A. Johnston, E. Lees, R. Singhal, A. Slight, M. Brimicombe, D. Stacey, V. Stacey, and H. Huhnermann, Journal of Physics G: Nuclear Physics 1, L6 (1975).
  * Wenz  _et al._ (1981) R. Wenz, A. Timmermann, and E. Matthias, Zeitschrift für Physik A Atoms and Nuclei 303, 87 (1981).
  * Buchinger  _et al._ (1987) F. Buchinger, P. Dabkiewicz, H.-J. Kluge, A. Mueller, and E.-W. Otten, Nuclear Physics A 462, 305 (1987).
  * Angeli and Marinova (2013) I. Angeli and K. P. Marinova, Atomic Data and Nuclear Data Tables 99, 69 (2013).
  * Fricke  _et al._ (2004) G. Fricke, K. Heilig, and H. F. Schopper, (2004).
  * Yordanov  _et al._ (2013) D. T. Yordanov, D. Balabanski, J. Bieroń, M. Bissell, K. Blaum, I. Budinčević, S. Fritzsche, N. Frömmgen, G. Georgiev, C. Geppert, _et al._ , Physical review letters 110, 192501 (2013).
  * Kelly and Tomchuk (1961) F. Kelly and E. Tomchuk, Proceedings of the Physical Society (1958-1967) 78, 1304 (1961).
  * BreÍt (1958) G. BreÍt, Reviews of modern physics 30, 507 (1958).
  * King (2013) W. H. King,  _Isotope shifts in atomic spectra_ (Springer Science & Business Media, 2013).
  * Bishop and King (1971) D. Bishop and W. King, Journal of Physics B: Atomic and Molecular Physics 4, 1798 (1971).
  * Bauche  _et al._ (1985) J. Bauche, H. Hühnermann, D. Stacey, V. Stacey, and M. Wilson, Zeitschrift für Physik A Atoms and Nuclei 320, 157 (1985).
  * Kahl and Berengut (2019) E. V. Kahl and J. C. Berengut, Computer Physics Communications 238, 232 (2019).
  * Berengut (2016) J. C. Berengut, Phys. Rev. A 94, 012502 (2016).
  * Berengut  _et al._ (2003) J. Berengut, V. Dzuba, and V. Flambaum, Physical Review A 68, 022502 (2003).
  * Solaro  _et al._ (2020) C. Solaro, S. Meyer, K. Fisher, J. C. Berengut, E. Fuchs, and M. Drewsen, Physical Review Letters 125, 123003 (2020).
  * Müller  _et al._ (2020) P. Müller, K. König, P. Imgram, J. Krämer, and W. Nörtershäuser, Physical Review Research 2, 043351 (2020).
  * Counts  _et al._ (2020) I. Counts, J. Hur, D. P. A. Craik, H. Jeon, C. Leung, J. C. Berengut, A. Geddes, A. Kawasaki, W. Jhe, and V. Vuletić, Physical Review Letters 125, 123002 (2020).
  * Seltzer (1969) E. Seltzer, Physical Review 188, 1916 (1969).
  * Reinhard and Nazarewicz (2017) P.-G. Reinhard and W. Nazarewicz, Physical Review C 95, 064328 (2017).
  * Goriely  _et al._ (2013) S. Goriely, N. Chamel, and J. Pearson, Physical Review C 88, 024308 (2013).
  * Fayans  _et al._ (2000) S. Fayans, S. Tolokonnikov, E. Trykov, and D. Zawischa, Nuclear Physics A 676, 49 (2000).
  * Goriely (2015) S. Goriely, Nuclear Physics A 933, 68 (2015).
  * Goriely  _et al._ (2016) S. Goriely, N. Chamel, and J. Pearson, Physical Review C 93, 034337 (2016).