# Effective light-induced Hamiltonian for atoms with large nuclear spin

D. Burba  Institute of Theoretical Physics and Astronomy, Vilnius University,
Saulėtekio 3, LT-10257, Vilnius, Lithuania  H. Dunikowski  Institute of
Physics PAS, Aleja Lotników 32/46, 02-668 Warszawa, Poland  M. Robert-de-
Saint-Vincent  Laboratoire de Physique des Lasers, CNRS UMR 7538, Université
Sorbonne Paris Nord, F-93430, Villetaneuse, France  E. Witkowska  Institute of
Physics PAS, Aleja Lotników 32/46, 02-668 Warszawa, Poland  G. Juzeliūnas
Institute of Theoretical Physics and Astronomy, Vilnius University, Saulėtekio
3, LT-10257, Vilnius, Lithuania

###### Abstract

Ultra-cold fermionic atoms, having two valence electrons, exhibit a
distinctive internal state structure, wherein the nuclear spin becomes
decoupled from the electronic degrees of freedom in the ground electronic
state. Consequently, the nuclear spin states are well isolated from the
environment, rendering these atomic systems an opportune platform for quantum
computation and quantum simulations. Coupling with off-resonance light is an
essential tool to selectively and coherently manipulate the nuclear spin
states. In this paper, we present a systematic derivation of the effective
Hamiltonian for the nuclear spin states of ultra-cold fermionic atoms due to
such an off-resonance light. We obtain compact expressions for the scalar,
vector and tensor light shifts taking into account both linear and quadratic
contributions to the hyperfine splitting. The analysis has been carried out
using the Green operator approach and solving the corresponding Dyson
equation. Finally, we analyze different scenarios of light configurations
which lead to the vector- and tensor-light shifts, as well as the pure spin-
orbit coupling for the nuclear spin.

##  I Introduction

Alkaline earth atoms stand as a prime platform for the realms of quantum
simulation and computation, as well as for high-performance optical clocks
Daley  _et al._ (2008); González-Cuadra  _et al._ (2023); Li  _et al._ (2023);
Gorshkov  _et al._ (2009); Lis  _et al._ (2023); Omanakuttan  _et al._ (2021);
Perlin  _et al._ (2022). These atoms have two valence electrons forming a
closed shell $ns^{2}$ with zero spin and angular momentum quantum numbers
($S=L=J=0$) in the electronic ground state ${}^{1}S_{0}$. They are also
characterized by forbidden optical transitions to electronic triplet states
${}^{3}P_{j}$, with $j=0,1,2$. Fermionic isotopes possess a nuclear spin that
is highly isolated from the environment rendering them a convenient platform
for the applications mentioned above Daley  _et al._ (2008); González-Cuadra
_et al._ (2023); Li  _et al._ (2023); Gorshkov  _et al._ (2009); Lis  _et al._
(2023); Omanakuttan  _et al._ (2021); Perlin  _et al._ (2022). The most
notable alkaline-earth species are Mg, Ca, Sr, but other species - in
particular Yb, and the group IIB transition metals (Hg McFerran  _et al._
(2012), Cd Yamaguchi  _et al._ (2019), Zn Wang and Ye (2007)) share these
properties Derevianko and Katori (2011); Dzuba and Derevianko (2019); Ludlow
_et al._ (2015). All of them but Zn have been laser-cooled and trapped in the
ground state. It is noteworthy that the number of nuclear spin states of their
fermionic isotopes can be quite large, such as $N=10$ for 87Sr, providing
beneficial features useful for quantum simulation and computation Omanakuttan
_et al._ (2021). The nuclear spin in these atoms thus represents a convenient
carrier of quantum information taking into account the long coherence times
and the ability of coherent spin control with magnetic and optical fields.

Manipulation of atomic spin states by off-resonant laser fields was studied
extensively in the context of both the alkali Chaudhury  _et al._ (2007);
Rosenbusch  _et al._ (2009); Deutsch and Jessen (2010); Rosenbusch  _et al._
(2009); Le Kien  _et al._ (2013); Goldman  _et al._ (2014) and alkaline earth
atoms Omanakuttan  _et al._ (2021); Perlin  _et al._ (2022). The light-induced
effective atomic spin Hamiltonian can be generally cast in terms of components
known as the scalar, vector and tensor light shifts Deutsch and Jessen (2010);
Le Kien  _et al._ (2013). Although the light shifts can be generally expressed
in terms of the Clebsch-Gordan coefficients describing virtual transitions to
various excited states, these expressions are complex to analyze and require
numerical calculation. On the other hand, the simple form of the effective
Hamiltonian and the scalar, vector and tensor components of the light shifts
are desired by practical experimental implementations and easier theoretical
understanding.

Here we develop a formalism allowing us to obtain compact analytical
expressions for the scalar, vector and tensor light shifts of the nuclear spin
states for the alkaline-earth-like atoms. The method involves the Green
operator (the resolvent) approach and the solution of the corresponding Dyson
equation, allowing to by-pass explicit analysis of transitions to the specific
hyperfine states of the excited state manifold. The technique is related to
the one used to analyze the vector light shifts of the alkali atoms presented
in Sec. 4 of Ref. Goldman  _et al._ (2014). Yet, there are significant
differences between the two treatments. In the present analysis of the
alkaline-earth-like atoms the electron spin is zero in the ground state, and
it is the hyperfine splitting of the excited electronic states which provides
the vector and tensor light shifts for the (nuclear) spin states of the ground
state atoms. On the other hand, for alkali atoms with a single unpaired outer
electron, it is the excited state fine structure splitting which is
responsible for the vector light shifts for the spin states in the ground
electronic manifold Goldman  _et al._ (2014). In our analytical calculations,
we have managed to include not only linear but also quadratic terms of the
hyperfine coupling operator. Therefore the theory provides a very accurate
description of the light shifts for the alkaline-earth-like atoms. To resolve
the hyperfine splitting without significant perturbation by spontaneous
emission we consider the light shifts due to the ${}^{3}P_{1}$
intercombination line characterized by a relatively large frequency of the
hyperfine splitting $A_{\rm HF}$ and a comparatively small spontaneous
emission rate $\Gamma$. As an example, for 87Sr atoms $A_{\rm HF}/2\pi$ is
around $2.6\times 10^{5}$ kHz, whereas $\Gamma/2\pi$ is of the order of 10 kHz
Miyake  _et al._ (2019) for the ${}^{3}P_{1}$ line. This allows an efficient
light-induced manipulation of the nuclear spin states with minimum
decoherence. In the subsequent text, we will consider 87Sr as a typical
example. Nevertheless, our analytical findings are broad-reaching and can be
applied to other alkaline earth (or alkaline earth-like) atoms with the same
structure.

The paper is organised as follows. In Section II we introduce the fine and
hyperfine structures of alkali earth or alkali earth-like atoms. Next, in
Section III, we present our theory for the derivation of effective Hamiltonian
describing the spin ground-state manifold when an atom is coupled with laser
radiation. In Section IV, we discuss typical examples of light configuration
and the resulting atom-light coupling Hamiltonian, as well as prospects for
the realization of spin-orbit coupling for fermionic atoms with nuclear spin.
Finally, the concluding Section V summarizes our findings.

##  II Manifolds of ground and excited states for alkaline-earth-like atoms

We consider alkali earth or alkali earth-like atoms with two $s$-shell valence
electrons characterised by a nuclear spin $i_{I}$. In the ground electronic
manifold denoted by ${}^{1}S_{0}$ the spins of outer electrons form a singlet,
so the quantum numbers of the electronic spin ${\bf S}$ and the orbital
angular momenta ${\bf L}$ are zero, as well as the total electronic momentum
${\bf J}$, namely $s=l=j=0$. It is a reason why the fine and hyperfine
couplings are not affecting the ground electronic state. The energy of this
state is taken to be zero $E_{g}=0$.

The excited electronic state of interest $|j,m_{j}\rangle$ with $j=0,1,2$
belongs to the manifold ${}^{3}P_{j}$ where spins of the outer electrons form
a triplet, and the state of the atom is characterized by the spin and orbital
angular momentum quantum numbers, $s=1$ and $l=1$, respectively
111Alternatively, the theory presented here is also applicable for virtual
transitions via the ${}^{1}P_{1}$ manifold. However, we will concentrate on
the ${}^{3}P_{1}$ manifold, as its higher ratio between hyperfine structure
and linewidth makes it more relevant experimentally due to much smaller rate
of spontaneous emission.. The state $|j,m_{j}\rangle$ is an eigenstate of
$\mathbf{S}^{2}$, $\mathbf{L}^{2}$, $\mathbf{J}^{2}$ and $\mathbf{J}_{z}$ with
the corresponding eigenvalues $\hbar^{2}s\left(s+1\right)=2\hbar^{2}$,
$\hbar^{2}l\left(l+1\right)=2\hbar^{2}$, $\hbar^{2}j\left(j+1\right)$ and
$\hbar m_{j}$. Without including the fine and hyperfine couplings the excited
state is degenerate concerning $j$. The fine structure coupling removes this
degeneracy.

The fine structure (FS) coupling between the total electron spin and orbital
angular momentum is described by the Hamiltonian $H_{\rm FS}=\frac{A_{\rm
FS}}{\hbar}\mathbf{L}\cdot\mathbf{S}=\frac{A_{\rm
FS}}{\hbar}\left(\mathbf{J}^{2}-\mathbf{L}^{2}-\mathbf{S}^{2}\right)/2$
providing $j$-dependent energies $E_{\rm FS}^{\left(j\right)}=\hbar A_{\rm
FS}\left[j\left(j+1\right)-4\right]/2$ for $l=s=1$. Therefore, one has

$E_{\rm FS}^{\left(j=0\right)}=-2\hbar A_{\rm FS}\,,\,\,E_{\rm
FS}^{\left(j=1\right)}=-\hbar A_{\rm FS}\,,\,\,E_{\rm
FS}^{\left(j=2\right)}=\hbar A_{\rm FS}\,.$ (1)

The state with $j=2$ is not accessible via the dipole transitions from the
ground state manifold ${}^{1}S_{0}$ with $j=0$, and thus will not be
considered by us in the further part of the text. The transition to the state
with $j=0$, also known as "clock transition", is double forbidden, both by
spin and orbital degrees of freedom. It is very weak, and only relevant if the
radiation frequency is very close to this state’s resonance. Therefore, in
what follows, we will consider only the transition to the excited state
manifold ${}^{3}P_{1}$ corresponding to $j=1$. This transition called the
intercombination line, is forbidden by spin flip, but opens due to $j-j$
coupling admixing with the ${}^{1}P_{1}$ state Foot (2004).

The hyperfine (HF) coupling between the total electronic angular momentum
$\bf{J}$ and the nuclear spin $\bf{I}$ is generally given by the Hamiltonian
(see e.g. ref. Le Kien  _et al._ (2013)):

$H_{\mathrm{HF}}=\frac{A_{{\rm
HF}}^{\prime}}{\hbar}\mathbf{I}\cdot\mathbf{J}+B_{{\rm
HF}}\frac{6\left(\mathbf{I}\cdot\mathbf{J}\right)^{2}+3\hbar^{2}\mathbf{I}\cdot\mathbf{J}-2\mathbf{I}^{2}\mathbf{J}^{2}}{4\hbar^{3}i_{I}(2i_{I}\text{\textminus}1)j(2j\text{\textminus}1)}\,,$
(2)

where in the present situation $j=1$. Here $A_{{\rm HF}}^{\prime}$ and
$B_{{\rm HF}}$ are HF constants. In the case of strontium atoms, their values
are $A_{{\rm HF}}^{\prime}/2\pi=-260085\pm 2\,$kHz and $B_{{\rm
HF}}/2\pi=-35667\pm 21\,$kHz Miyake  _et al._ (2019); zu Putlitz (1963). It is
convenient to represent $H_{\mathrm{HF}}$ as

$H_{\mathrm{HF}}=\frac{A_{{\rm
HF}}}{\hbar}\left[\mathbf{I}\cdot\mathbf{J}+\frac{\gamma}{\hbar^{2}}\left(\mathbf{I}\cdot\mathbf{J}\right)^{2}\right]\,,$
(3)

where $A_{\rm HF}=A_{{\rm HF}}^{\prime}+3B_{{\rm
HF}}/[4i_{I}(2i_{I}-1)j(2j-1)]$, and the dimensionless parameter
$\gamma=6B_{{\rm HF}}/[A_{{\rm HF}}4i_{I}(2i_{I}-1)j(2j-1)]$ describing a
relative strength of the quadratic term is generally much less than the unity.
For 87Sr atoms one has $\gamma\approx 0.006$, and we will use this value in
specific examples. In Eq. (3) a uniform shift of energy proportional to
$\mathbf{I}^{2}\mathbf{J}^{2}$ has been incorporated to the excitation energy
$E_{\rm FS}^{\left(j=1\right)}$, whereas the two terms of Eq. (2) proportional
to $\mathbf{I}\cdot\mathbf{J}$ have been merged to a single term characterised
by the modified frequency $A_{{\rm HF}}$.

The eigenstates of hyperfine Hamiltonian are also eigenstates of
$\mathbf{F}^{2}$, $\mathbf{I}^{2}$, $\mathbf{J}^{2}$ and $\mathbf{F}_{z}$, and
the corresponding eigenvalues are $\hbar^{2}f\left(f+1\right)$,
$\hbar^{2}i_{I}\left(i_{I}+1\right)$, $\hbar^{2}j\left(j+1\right)$, and $\hbar
m_{f}$. For the excited state ${}^{3}P_{1}$ manifold of interest one has
$j=1$, so $f$ can take the values $f=i_{I}$ and $f=i_{I}\pm 1$. The energies
of these excited states are:

$\displaystyle E_{i_{I}+1}$ $\displaystyle=\hbar A_{{\rm
HF}}\,i_{I}\left(1+\gamma i_{I}\right)\,,$ (4) $\displaystyle E_{i_{I}}$
$\displaystyle=-\hbar A_{{\rm HF}}\left(1-\gamma\right)\,,$ (5) $\displaystyle
E_{i_{I}-1}$ $\displaystyle=-\hbar A_{{\rm
HF}}\left(i_{I}+1\right)\left[1-\gamma\left(i_{I}+1\right)\right]\,.$ (6)

The hyperfine structure for the electronic line ${}^{3}P_{1}$ is schematically
represented in Fig. 1.

Figure 1:  A sketch of the atomic levels for alkaline-earth-like atoms
illustrating the hyperfine structure of the excited electronic line
${}^{3}P_{1}$. The order of the excited states depends on the sign of $A_{\rm
HF}$; here we take $A_{\rm HF}<0$ as in 87Sr. The probe frequency $\omega_{L}$
is tuned in a vicinity of the transition ${}^{1}S_{0}\to{}^{3}P_{1}$ manifold
with the detuning $\Delta$.

##  III Effective Hamiltonian

We consider the effects of light field $\mathbf{E}$ on atomic nuclear spin
states in the ground electronic manifold ${}^{1}S_{0}$ due to the virtual
transition to the excited state manifold ${}^{3}P_{1}$. In the dipole
approximation, the interaction between a monochromatic light field and an atom
can be described by the operator

$V=-\mathbf{d}\cdot\mathbf{E}=-\sum_{q}{d}_{q}\tilde{\mathcal{E}}_{q}\cos(\phi_{q}-\omega_{L}t).$
(7)

where the components of amplitude $\tilde{\mathcal{E}}_{q}$ and phase
$\phi_{q}$ of the radiation are generally position dependent. Here ${\bf
d}=-e\sum_{\alpha}{\bf r}_{\alpha}$ is the electric dipole operator, $e$ is an
electron charge and ${\bf r}_{\alpha}$ is the position of the $\alpha$th
electron within the atom. The interaction Hamiltonian, in a frame rotated at
the driving frequency $\omega_{L}$, reads

$V=-\frac{1}{2}\sum_{q}\left[\mathcal{E}^{*}_{q}P_{g}d_{q}P_{e}+\mathcal{E}_{q}P_{e}d_{q}P_{g}\right],$
(8)

where we introduced the complex electric field
$\mathcal{E}_{q}=\tilde{\mathcal{E}}_{q}e^{i\phi_{q}}$ and we applied the
rotated wave approximation (RWA) by removing terms oscillating at frequencies
$\omega_{L}$, $2\omega_{L}$, see e.g. Chapter 4 of Ref. Goldman  _et al._
(2014). In the above formula, $P_{g}$ and $P_{e}$ are projector operators onto
the ground and excited state manifolds ${}^{1}S_{0}$ and ${}^{3}P_{1}$,
respectively.

###  III.1 Analysis of effective Hamiltonian

If the incident light $\boldsymbol{\mathcal{E}}$ is sufficiently weak and/or
well detuned from atomic resonances, the saturation parameter is small, and
the atom-light interaction effectively provides coupling between the nuclear
spin states $|i_{I},m_{i_{I}}\rangle$ of the electronic ground state
$|g\rangle$ via virtual transitions to the excited states without populating
them. Specifically, the adiabatic elimination of the excited states can be
described in the second order via Schriffer-Wolf transformation, see e.g.
Supplementary Materials of Hernández Yanes  _et al._ (2022). This provides the
following effective atomic Hamiltonian for the ground state manifold:

$\hat{H}_{\rm eff}=P_{g}VGVP_{g}$ (9)

where

$G=(\Delta-H_{\rm HF})^{-1}$ (10)

is the Green operator and $\Delta=E_{g}+E_{\rm FS}^{(j=1)}-\hbar\omega_{L}$ is
detuning. Note that by adding an imaginary part to the detuning $\Delta$, one
can include effects of losses due to spontaneous emission. It is discussed in
Sec. III.2.

The effective atomic Hamiltonian can be expressed as

$\hat{H}_{\rm eff}=\frac{1}{4}\mathcal{E}_{s}^{*}D_{s,q}\mathcal{E}_{q}.$ (11)

Here and below the summation over repeated indices of vectors and tensors is
assumed. In Eq.(11) we introduced the tensor operator acting on the ground
state manifold:

$\displaystyle D_{s,q}$ $\displaystyle=P_{g}d_{s}P_{e}GP_{e}d_{q}P_{g}\,,$
(12)

where the Green operator $G$ obeys a Dyson-type equation (identity)

$\displaystyle G$
$\displaystyle=\frac{1}{\Delta}+\frac{1}{\Delta}H_{\mathrm{HF}}G,$ (13)

see Appendix A. Combining Eqs. (12) and (13), one has

$\displaystyle D_{s,q}$
$\displaystyle=\frac{\left|d_{ge}^{2}\right|}{\Delta}P_{g}\delta_{sq}+\frac{1}{\Delta}P_{g}d_{s}P_{e}H_{\rm
HF}GP_{e}d_{q}P_{g}\,,$ (14)

where

$\left|d_{ge}^{2}\right|=\frac{1}{3}\left\langle
g\right|d_{q}P_{e}d_{q}\left|g\right\rangle$ (15)

characterizes the strength of the dipole transition and equals to the third of
the reduced matrix element $\langle
g|d_{q}P_{e}d_{q}|g\rangle=\sum_{e}|\langle e|d_{q}|g\rangle|^{2}$. The dipole
transition is weak, but non-zero. This is because the excited state manifold
${}^{3}P_{1}$ contains a small admixture of the manifold ${}^{1}P_{1}$ with
the same $j=1$, but zero spin $s=0$, and the transitions to the latter
manifold ${}^{1}S_{0}\rightarrow{}^{1}P_{1}$ are dipole allowed.

According to the Landé projection theorem Deutsch and Jessen (2010), the
effective Hamiltonian (11) can be represented in the form

$\displaystyle\hat{H}_{\rm eff}$
$\displaystyle=\frac{b_{0}}{4}\left(\boldsymbol{\mathcal{E}}^{*}\cdot\boldsymbol{\mathcal{E}}\right)+i\frac{b_{1}}{4\hbar}\textbf{$\mathbf{I}$}\cdot\left(\boldsymbol{\mathcal{E}}^{*}\times\boldsymbol{\mathcal{E}}\right)$
$\displaystyle+\frac{b_{2}}{4\hbar^{2}}\left[\left(\boldsymbol{\mathcal{E}}^{*}\cdot\textbf{$\mathbf{I}$}\right)\left(\boldsymbol{\mathcal{E}}\cdot\textbf{$\mathbf{I}$}\right)-\frac{1}{3}\left|\boldsymbol{\mathcal{E}}\right|^{2}\mathbf{I}^{2}+{\rm
h.c.}\right].$ (16)

where the coefficients $b_{0},b_{1}$ and $b_{2}$ characterize the strengths of
the scalar, vector and tensor terms. Conventionally the effective Hamiltonian
is calculated using the explicit expressions for the transition matrix
elements involving the Clebsch-Gordon coefficients. Here, we perform
analytical calculations of the coefficients $b_{0},b_{1}$ and $b_{2}$ defining
the effective Hamiltonian to derive their simple forms in a self-consistent
way based on solution of Eq. (14) for $D_{s,q}$. To this end, it is convenient
to introduce a vector operator $\mathbf{K}$ with Cartesian components derived
from $D_{s,q}$:

$K_{s}=D_{s,q}\mathcal{E}_{q}\,.$ (17)

In terms of $\mathbf{K}$, the effective Hamiltonian of Eq (11) reads:

$\hat{H}_{{\rm
eff}}=\frac{1}{4}\boldsymbol{\mathcal{E}}^{*}\cdot\mathbf{K}\,.$ (18)

Multiplying Eq. (14) by $\mathcal{E}_{q}$ from the right-hand side, we write
down the relation for the vector-operator $\mathbf{K}$

$\displaystyle\mathbf{K}$
$\displaystyle=\frac{1}{\Delta}P_{g}\left|d_{ge}^{2}\right|\mathcal{\boldsymbol{\mathcal{E}}}$
$\displaystyle+\frac{A_{{\rm
HF}}}{\hbar\Delta}P_{g}\mathbf{d}\left(\mathbf{I}\cdot\mathbf{J}\right)\left(1+\frac{\gamma}{\hbar^{2}}\mathbf{I}\cdot\mathbf{J}\right)GP_{e}\left(\mathbf{d}\cdot\boldsymbol{\mathcal{E}}\right)P_{g}\,.$
(19)

As demonstrated in Appendix B.1, this provides the following equation for
$\mathbf{K}$:

$\displaystyle\eta\mathbf{K}$
$\displaystyle=\frac{1}{\Delta}P_{g}\left|d_{ge}^{2}\right|\boldsymbol{\mathcal{E}}+i\frac{A_{{\rm
HF}}}{\Delta}\left(1-\gamma\right)\left(\textbf{$\mathbf{I}$}\times\mathbf{K}\right)-$
$\displaystyle-\gamma\frac{A_{{\rm
HF}}}{\hbar\Delta}\mathbf{I}\left(\textbf{$\mathbf{I}$}\cdot\mathbf{K}\right)\,,$
(20)

where $\eta=1-\gamma A_{{\rm HF}}\hbar i_{I}\left(i_{I}+1\right)/\Delta$.

Anticipating the effective Hamiltonian of the form (16), we will look for the
following ${\bf K}$ ansatz

$\displaystyle\mathbf{K}$
$\displaystyle=\left[a_{0}\boldsymbol{\mathcal{E}}-i\frac{a_{1}}{\hbar}\textbf{$\mathbf{I}$}\times\boldsymbol{\mathcal{E}}+\frac{a_{2}}{\hbar^{2}}\textbf{$\mathbf{I}$}\left(\textbf{$\mathbf{I}$}\cdot\boldsymbol{\mathcal{E}}\right)\right]P_{g},$
(21)

where $a_{0}$, $a_{1}$ and $a_{2}$ are real coefficients. Inserting the ansatz
(21) into Eq. (20), one gets a system of equations for these coefficients,
giving

$\displaystyle a_{0}$ $\displaystyle=-\frac{|d_{ge}^{2}|}{\hbar A_{\rm
HF}}\frac{\bar{E}^{(+)}}{\left(\bar{\Delta}-\bar{E}_{i_{I}+1}\right)\left(\bar{\Delta}-\bar{E}_{i_{I}-1}\right)}\,,$
(22) $\displaystyle a_{1}$ $\displaystyle=\frac{|d_{ge}^{2}|}{\hbar A_{\rm
HF}}\frac{\bar{E}_{i_{I}}}{\left(\bar{\Delta}-\bar{E}_{i_{I}+1}\right)\left(\bar{\Delta}-\bar{E}_{i_{I}-1}\right)}\,,$
(23) $\displaystyle a_{2}$ $\displaystyle=\frac{|d_{ge}^{2}|}{\hbar A_{\rm
HF}}\frac{\gamma\bar{E}^{(-)}-\bar{E}_{i_{I}}^{2}}{\left(\bar{\Delta}-\bar{E}_{i_{I}+1}\right)\left(\bar{\Delta}-\bar{E}_{i_{I}}\right)\left(\bar{\Delta}-\bar{E}_{i_{I}-1}\right)}\,,$
(24)

with

$\bar{E}^{(\pm)}=(\bar{E}_{i_{I}+1}\pm\bar{E}_{i_{I}}+\bar{E}_{i_{I}-1})/2-\bar{\Delta}.$
(25)

where $\bar{\Delta}=\Delta/(\hbar A_{\rm HF})$ and $\bar{E}_{f}=E_{f}/(\hbar
A_{\rm HF})$ (with $f=i_{I}-1,i_{I},i_{I}+1$) are, respectively, the
dimensionless detuning and energies of hyperfine states, and $E_{f}$ is
defined by Eqs. (4)-(6). Details of derivation are in Appendix B.2.

Substituting the solution (21) for $\mathbf{K}$ into Eq. (18), one arrives at
the explicit result for the ground state effective Hamiltonian expressed in
terms of the coefficients $a_{0,1,2}$

$\displaystyle\hat{H}_{\rm eff}$
$\displaystyle=\frac{a_{0}}{4}\left(\boldsymbol{\mathcal{E}}^{*}\cdot\boldsymbol{\mathcal{E}}\right)+i\frac{a_{1}}{4\hbar}\textbf{$\mathbf{I}$}\cdot\left(\boldsymbol{\mathcal{E}}^{*}\times\boldsymbol{\mathcal{E}}\right)+$
$\displaystyle+\frac{a_{2}}{4\hbar^{2}}\left(\boldsymbol{\mathcal{E}}^{*}\cdot\textbf{$\mathbf{I}$}\right)\left(\boldsymbol{\mathcal{E}}\cdot\textbf{$\mathbf{I}$}\right)\,,$
(26)

where for brevity we have omitted the ground state projector operator $P_{g}$.

Note that the self-consistent method we used to find the vector ${\bf K}$ of
Eq. (21) can be applied on the level of the tensor $D_{s,q}$ featured in Eqs.
(12) and (14). In that case one finds

$D_{s,q}=a_{0}\delta_{sq}+i\frac{a_{1}}{\hbar}I_{k}\epsilon_{ksq}+\frac{a_{2}}{\hbar^{2}}I_{s}I_{q}\,.$
(27)

This leads to the same form of the effective Hamiltonian (26) with the same
coefficients $a_{0,1,2}$ defined by Eqs. (22)-(24).

Using the fact that
$\boldsymbol{\mathcal{E}}^{*}\cdot\left(\boldsymbol{I}\times\boldsymbol{\mathcal{E}}\right)=-\boldsymbol{\mathcal{E}}^{*}\cdot\left(\boldsymbol{\mathcal{E}}\times\boldsymbol{I}\right)=-\left(\boldsymbol{\mathcal{E}}^{*}\times\boldsymbol{\mathcal{E}}\right)\cdot\boldsymbol{I}$,
we express the effective Hamiltonian (26) in the symmetric and traceless form
as in Eq. (16). Consequently, we obtain the following relations for the
coefficients $b_{0}$, $b_{1}$ and $b_{2}$ defining the scalar, vector and
tensor light shifts

$b_{0}=a_{0}+\bar{i}_{I}^{2}a_{2}/3,\,\,b_{1}=a_{1}+a_{2}/2\,\,\mathrm{and}\,\,b_{2}=a_{2}/2.$
(28)

with

$\bar{i}_{I}^{2}\equiv\mathbf{I}^{2}/\hbar^{2}=i_{I}(i_{I}+1)\,,$ (29)

It is convenient to define the dimensionless $b$ coefficients

$\bar{b}_{u}=b_{u}\hbar A_{\rm HF}/|d_{ge}^{2}|\,\,,\quad\mathrm{with}\quad
u=0,1,2\,.$ (30)

Calling on Eqs. (22)-(24) and (28), one arrives at the explicit expressions
for the latter coefficients:

$\displaystyle\bar{b}_{0}$
$\displaystyle=\frac{-3\bar{E}^{(+)}\left(\bar{\Delta}-\bar{E}_{i_{I}}\right)+\bar{i}_{I}^{2}\left(\gamma\bar{E}^{(-)}-\bar{E}_{i_{I}}^{2}\right)}{3\left(\bar{\Delta}-\bar{E}_{i_{I}+1}\right)\left(\bar{\Delta}-\bar{E}_{i_{I}}\right)\left(\bar{\Delta}-\bar{E}_{i_{I}-1}\right)}\,,$
(31) $\displaystyle\bar{b}_{1}$
$\displaystyle=\frac{2\bar{E}_{i_{I}}\left(\bar{\Delta}-\bar{E}_{i_{I}}\right)+\left(\gamma\bar{E}^{(-)}-\bar{E}_{i_{I}}^{2}\right)}{2\left(\bar{\Delta}-\bar{E}_{i_{I}+1}\right)\left(\bar{\Delta}-\bar{E}_{i_{I}}\right)\left(\bar{\Delta}-\bar{E}_{i_{I}-1}\right)}\,,$
(32) $\displaystyle\bar{b}_{2}$
$\displaystyle=\frac{\left(\gamma\bar{E}^{(-)}-\bar{E}_{i_{I}}^{2}\right)}{2\left(\bar{\Delta}-\bar{E}_{i_{I}+1}\right)\left(\bar{\Delta}-\bar{E}_{i_{I}}\right)\left(\bar{\Delta}-\bar{E}_{i_{I}-1}\right)}\,.$
(33)

In Fig. 2 we plot the dimensionless coefficients $\bar{b}_{0,1,2}$ as a
function of the relative detuning $\bar{\Delta}$ using Eqs. (31)-(33) and the
conventional Clebsch-Gordan approach outlined in Appendix C. There is a
perfect agreement with a relative deviation due to finite computing precision
below $10^{-13}$.

Figure 2:  Dependence of the rescaled scalar $\bar{b}_{0}$ (a), vector
$\bar{b}_{1}$ and tensor $\bar{b}_{2}$ coefficients (b) on the rescaled
detuning $\bar{\Delta}$ for $i_{I}=9/2$ and $\gamma=0.006$, corresponding to
87Sr atoms. The analytical results given by (31), (32) and (33) are marked by
solid lines while squares, circles and triangles, respectively, represent
Clebsch-Gordan numerical calculations.

The analytical expressions (31)-(33) show explicitly the detuning dependence
of the coefficients $\bar{b}_{0,1,2}$ describing, respectively, the scalar,
vector and tensor light shifts. For example, by taking the 87Sr atom with
$i_{I}=9/2$ and choosing detuning comparable with the HF energies,
$\bar{\Delta}=\sum_{f=i_{I},i_{I}\pm 1}\bar{E}_{f}/3\approx-0.57$, the
dimensional coefficients are: $\bar{b}_{0}\approx 0.76$, $\bar{b}_{1}\approx
0.09$ and $\bar{b}_{2}\approx 0.05$. Thus the scalar term dominates for the
specific value of $\bar{\Delta}$ chosen. Yet by changing $\bar{\Delta}$ it is
possible to find a point where the scalar term $\bar{b}_{0}$ changes the sign
and thus goes through the zero point, as one sees in Fig. 2. Hence one can
eliminate the scalar term without cancelling the other terms.

Neglecting the quadratic hyperfine coupling ($\gamma=0$ and thus
$\bar{E}_{i_{I}}=-1$), Eqs. (32)-(33) shows that the ratio of the vector and
tensor light shifts goes linearly with the detuning:

$\bar{b}_{1}/\bar{b}_{2}=2\bar{\Delta}+3\,.$ (34)

Such a linear dependence holds well also after including the quadratic
hyperfine coupling which is typically small ($\gamma\ll 1$). From Eq (34) it
follows that the ratio $\bar{b}_{1}/\bar{b}_{2}$ goes to zero at
$\bar{\Delta}=-3/2$, as one can see in Fig. 2.

For far-off resonance radiation when detuning is much larger than the energy
of the hyperfine splitting $\bar{|\Delta}|\gg 1$, the dimensionless
coefficients read up to terms cubic in detuning

$\displaystyle\bar{b}_{0}$
$\displaystyle\approx\bar{\Delta}^{-1}+\frac{2}{3}\gamma\bar{i}_{I}^{2}\bar{\Delta}^{-2}+\frac{2}{3}\bar{i}_{I}^{2}\left[1+\gamma(\gamma\bar{i}_{I}^{2}-1)\right]\bar{\Delta}^{-3},$
(35) $\displaystyle\bar{b}_{1}$
$\displaystyle\approx-\left(1-\frac{\gamma}{2}\right)\bar{\Delta}^{-2}+\frac{1}{2}\left[1-4\gamma\bar{i}_{I}^{2}+\gamma^{2}\left(3\bar{i}_{I}^{2}-1\right)\right]\bar{\Delta}^{-3},$
(36) $\displaystyle\bar{b}_{2}$
$\displaystyle\approx-\frac{\gamma}{2}\bar{\Delta}^{-2}+\frac{1}{2}\left[-1+4\gamma-\gamma^{2}\left(\bar{i}_{I}^{2}+3\right)\right]\bar{\Delta}^{-3}.$
(37)

In this way, for large detuning, $\bar{\Delta}\gg 1$, the scalar light shift
$\propto\bar{\Delta}^{-1}$ decreases linearly with the inverse detuning,
whereas the vector light shift $\propto\bar{\Delta}^{-2}$ goes quadratically.
For intermediate detuning $1\ll\bar{\Delta}\ll\gamma^{-1}$ the tensor light
shift decreases as $\bar{\Delta}^{-3}$ and then cross over to an asymptotic
decay $\propto\gamma/\bar{\Delta}^{2}$ for $\bar{\Delta}\gg\gamma^{-1}$, with
$|\bar{b}_{2}/\bar{b}_{1}|\ll 1$ ensured by the weaker quadratic hyperfine
coupling (typically $\gamma\ll 1$).

###  III.2 Effects of spontaneous emission

Let us now include the effects of the finite radiative lifetime of the excited
state. This can be done by replacing the detuning $\Delta$ by
$\Delta-i\hbar\Gamma$ in our treatment including the expressions for the
coefficients $a_{0,1,2}$ and $b_{0,1,2}$, where $\Gamma$ is the excited state
linewidth due to spontaneous emission. Consequently, the effective Hamiltonian
acquires a non-Hermitian part representing the radiative losses.

To estimate the losses let us expand ${\rm Im}\,\bar{b}_{0,1,2}$ to first
order in the linewidth $\Gamma$ assuming that the detuning is not too close to
the hyperfine lines ($|\Delta-E_{f}|\gg\hbar\Gamma$) and neglecting the
quadratic contribution to the hyperfine coupling ($\gamma=0$):

${\rm
Im}\,\bar{b}_{0}\approx\frac{3(\bar{\Delta}+1)^{4}+2(\bar{\Delta}+1)\bar{i}_{I}^{2}+\bar{i}_{I}^{4}}{3\left(\bar{\Delta}-\bar{E}_{i_{I}+1}\right)^{2}\left(\bar{\Delta}-\bar{E}_{i_{I}}\right)^{2}\left(\bar{\Delta}-\bar{E}_{i_{I}-1}\right)^{2}}\bar{\Gamma}\,,$
(38) ${\rm
Im}\,\bar{b}_{1}\approx-\frac{4\bar{\Delta}^{3}+13\bar{\Delta}^{2}+12\bar{\Delta}-\bar{i}_{I}^{2}+3}{2\left(\bar{\Delta}-\bar{E}_{i_{I}+1}\right)^{2}\left(\bar{\Delta}-\bar{E}_{i_{I}}\right)^{2}\left(\bar{\Delta}-\bar{E}_{i_{I}-1}\right)^{2}}\bar{\Gamma}\,,$
(39) ${\rm
Im}\,\bar{b}_{2}\approx-\frac{3\bar{\Delta}^{2}+4\bar{\Delta}-\bar{i}_{I}^{2}+1}{2\left(\bar{\Delta}-\bar{E}_{i_{I}+1}\right)^{2}\left(\bar{\Delta}-\bar{E}_{i_{I}}\right)^{2}\left(\bar{\Delta}-\bar{E}_{i_{I}-1}\right)^{2}}\bar{\Gamma}\,,$
(40)

with $\bar{\Gamma}=\Gamma/A_{\rm HF}$.

Although these formulas do not hold when the detuning $|\Delta-E_{f}|$ becomes
comparable to $\hbar\Gamma$ or is smaller, such a situation is not interesting
because of a significant increase of losses compared to the light shifts. On
the other hand, the intermediate ($\hbar\Gamma\ll|\Delta-E_{f}|\ll\hbar A_{\rm
HF}$) and far detuned regimes are well described by Eqs. (38)-(40). For the
far detuned radiation, $|\Delta-E_{f}|\gg\hbar A_{\rm HF}$, the losses go as
${\rm Im}\,\bar{b}_{0}\propto\bar{\Delta}^{-2}$, ${\rm
Im}\,\bar{b}_{1}\propto\bar{\Delta}^{-3}$ and ${\rm
Im}\,\bar{b}_{2}\propto\bar{\Delta}^{-4}$, so the dominant loss term ${\rm
Im}\,\bar{b}_{0}$ and the vector light shift coefficient ${\rm Re}\
\bar{b}_{1}$ both go as $1/\Delta^{2}$. Hence, in the far detuned regime the
ratio ${\rm Im}\bar{b}_{0}/{\rm Re}\bar{b}_{1}$ becomes constant and is
determined by the ratio of the linewidth $\Gamma$ and the hyperfine structure
constant $A_{\rm HF}$:

${|{\rm Im}\bar{b}_{0}/{\rm Re}\ \bar{b}_{1}|\approx\Gamma/|A_{{\rm HF}}|\ll
1\,.}$ (41)

For example, for the ${}^{1}S_{0}\rightarrow^{3}P_{1}$ line of the species:
spin-7/2 43Ca; spin-9/2 87Sr; spin-1/2 171Yb; spin-5/2 173Yb, this ratio is,
respectively, $\Gamma/|A_{{\rm
HF}}|=2.\,10^{-6};3.\,10^{-5};5.\,10^{-5};2.\,10^{-4}$, showing that the
intercombination line of alkaline-earth-like atoms is very favorable to
engineer nuclear-spin-sensitive vector or tensor light shifts.

Yet in the far-detuned regime, the vector and tensor light shits decrease
considerably. To get their larger values, one needs to go closer to the
resonances or even between the resonance lines. In that case the relative
contribution of the losses can still have reasonably small (of the order of
$10^{-4}$) with a considerably increase of the vector light polarizability, as
one can see in Fig. 3.

Figure 3:  Dependence of the loss ratio ${\rm Im}\ \bar{b}_{0}/{\rm Re}\
\bar{b}_{1,2}$ (blue solid lines, left vertical axis) and the dimensionless
vector (a) and tensor (b) polarizabilities ${\rm Re}\,\bar{b}_{1,2}$ (dashed
yellow lines, the right vertical axis) on the dimensionless detuning
$\bar{\Delta}$ for $\Gamma/|A_{{\rm HF}}|=3\times 10^{-5}$. Red dash-dotted
horizontal lines represent the limiting value of the ratio ${\rm Im}\
\bar{b}_{0}/{\rm Re}\ \bar{b}_{1}$ for large detuning given by Eq. (41).

##  IV Specific configurations

The analytical form of the effective Hamiltonian (16) of the ground state
manifold, can be presented as

$\hat{H}_{{\rm eff}}=\hat{H}_{{\rm eff}(0)}+\hat{H}_{{\rm
eff}(1)}+\hat{H}_{{\rm eff}(2)}\,,$ (42)

where

$\hat{H}_{{\rm
eff}(0)}\coloneqq\frac{b_{0}}{4}\left(\boldsymbol{\mathcal{E}}^{*}\cdot\boldsymbol{\mathcal{E}}\right)=\frac{b_{0}}{4}\left|\boldsymbol{\mathcal{E}}\right|^{2}\,,$
(43)

is the scalar component,

$\hat{H}_{{\rm
eff}(1)}\coloneqq\frac{ib_{1}}{4\hbar}\textbf{$\mathbf{I}$}\cdot\left(\boldsymbol{\mathcal{E}}^{*}\times\boldsymbol{\mathcal{E}}\right)=-\frac{b_{1}}{2\hbar}\textbf{$\mathbf{I}$}\cdot\left(\boldsymbol{\mathcal{E}}^{\prime}\times\boldsymbol{\mathcal{E}}^{\prime\prime}\right)\,,$
(44)

is the vector component, and

$\displaystyle\hat{H}_{{\rm eff}(2)}\coloneqq$
$\displaystyle\frac{b_{2}}{4\hbar^{2}}\left[\left(\boldsymbol{\mathcal{E}}^{*}\cdot\textbf{$\mathbf{I}$}\right)\left(\boldsymbol{\mathcal{E}}\cdot\textbf{$\mathbf{I}$}\right)-\frac{1}{3}\left|\boldsymbol{\mathcal{E}}\right|^{2}{\mathbf{I}}^{2}+{\rm
h.c.}\right]$ $\displaystyle=$
$\displaystyle\frac{b_{2}}{2\hbar^{2}}\left[\left(\boldsymbol{\mathcal{E}}^{\prime}\cdot\textbf{$\mathbf{I}$}\right)^{2}+\left(\boldsymbol{\mathcal{E}}^{\prime\prime}\cdot\textbf{$\mathbf{I}$}\right)^{2}-\frac{1}{3}\left|\boldsymbol{\mathcal{E}}\right|^{2}{\mathbf{I}}^{2}\right]\,,$
(45)

is the tensor component. In the above expressions, we used the real and
imaginary parts of the complex light field
$\boldsymbol{\mathcal{E}}=\boldsymbol{\mathcal{E}}^{\prime}+i\boldsymbol{\mathcal{E}}^{\prime\prime}$.
This entirely analytical form allows one to explore its form to design desired
atom-light coupling. Below we discuss three example configurations including a
single beam and a pair of counter-propagating beams and a pair of
perpendicular beams. The last configuration is demonstrated to generate spin-
orbit coupling for the nuclear spin. Finally in Subsection IV.4 we discuss a
way of eliminating the tensor light shifts by combining light fields with
highly different frequencies.

###  IV.1 Single beam

Let us consider first a situation where the system is subjected to a single
linearly polarized laser beam

$\boldsymbol{\mathcal{E}}=\mathcal{E}e^{ikz}\boldsymbol{\mathrm{e}}_{z}\,.$
(46)

where the amplitude $\mathcal{E}$ is assumed to be real and no extra global
phase is added. This can be done by properly choosing the origin of spatial
coordinates in all the examples considered here in Secs. IV.1-IV.3.

For the single beam the vector light shift (44) is zero and effective
Hamiltonian contains only the position-independent scalar and tensor
components

$\hat{H}_{{\rm
eff}}=\frac{\mathcal{E}^{2}}{4}\left[b_{0}\mathbb{1}+\frac{2b_{2}}{\hbar^{2}}\left({\bf
I}_{z}^{2}-\frac{{\bf I}^{2}}{3}\right)\right]\,,$ (47)

where $\mathbb{1}$ is an identity operator. This shifts the ground-state
energies according to $\hat{H}_{{\rm
eff}}|g\rangle|i_{I},m_{I}\rangle=E^{(i_{I})}_{m_{I}}|g\rangle|i_{I},m_{I}\rangle$
with eigenvalues having a scalar term (independent of $m_{I}$) and tensor
shifts (quadratic in $m_{I}$), namely
$E^{(i_{I})}_{m_{I}}=\mathcal{E}^{2}/4\left[b_{0}+2b_{2}\left(m_{I}^{2}-\bar{i}_{I}^{2}/3\right)\right]$.

In the case of a circularly polarized laser beam,

$\boldsymbol{\mathcal{E}}=\mathcal{E}e^{ikz}\boldsymbol{e}_{\pm}\,,\quad\mathrm{with}\quad\boldsymbol{e}_{\pm}=\frac{1}{\sqrt{2}}\left(\boldsymbol{\mathrm{e}}_{x}\pm
i\boldsymbol{\mathrm{e}}_{y}\right)\,,$ (48)

the effective Hamiltonian contains additionally the (position-independent)
vector components $\propto{\bf I}_{z}$:

$\hat{H}_{{\rm
eff}}=\frac{\mathcal{E}^{2}}{4}\left[b_{0}\mathbb{1}\mp\frac{b_{1}}{\hbar}{\bf
I}_{z}+\frac{b_{2}}{\hbar^{2}}\left(\frac{1}{3}{\bf I}^{2}-{\bf
I}_{z}^{2}\right)\right]\,.$ (49)

In this configuration, eigenenergies of the ground states manifold
$|g\rangle|i_{I},m_{I}\rangle$ contain additional linear shifts,
$E^{(i_{I})}_{m_{I}}=\mathcal{E}^{2}/4\left[b_{0}\mp
b_{1}m_{I}-b_{2}\left(m_{I}^{2}-\bar{i}_{I}^{2}/3\right)\right]$.

###  IV.2 Two cross polarized counter-propagating laser beams

Consider next the effect of two cross-polarized counter-propagating Raman
beams where

$\boldsymbol{\mathcal{E}}=\frac{\mathcal{E}^{2}}{\sqrt{2}}\left(e^{ikz}\boldsymbol{\mathrm{e}}_{x}+e^{-ikz}\boldsymbol{\mathrm{e}}_{y}\right)\,.$
(50)

In that case, the scalar light shift is uniform

$\hat{H}_{{\rm eff}(0)}=\frac{\mathcal{E}^{2}}{4}b_{0}\mathbb{1}\,$ (51)

and there are spatially periodic vector and tensor components

$\hat{H}_{{\rm
eff}(1)}=\frac{\mathcal{E}^{2}}{4}\frac{b_{1}}{\hbar}\sin\left(2kz\right){\bf
I}_{z}\,.$ (52)

and

$\hat{H}_{{\rm
eff}(2)}=\frac{1}{2}\frac{b_{2}}{\hbar^{2}}\mathcal{E}^{2}\left(\cos^{2}\left(kz\right){\bf
I}_{\tilde{y}}^{2}+\sin^{2}\left(kz\right){\bf
I}_{\tilde{x}}^{2}-\frac{1}{3}{\bf I}^{2}\right)\,,$ (53)

with

${\bf I}_{\tilde{x}}\coloneqq\frac{{\bf I}_{x}-{\bf
I}_{y}}{\sqrt{2}}\quad\mathrm{and}\quad{\bf I}_{\tilde{y}}\coloneqq\frac{{\bf
I}_{x}+{\bf I}_{y}}{\sqrt{2}}\,.$ (54)

The latter $\hat{H}_{{\rm eff}(2)}$ can be also represented in terms of the
original nuclear spin operators ${\bf I}_{x,y,z}$:

$\hat{H}_{{\rm
eff}(2)}=\frac{1}{2}\frac{b_{2}}{\hbar^{2}}\mathcal{E}^{2}\left(\frac{1}{6}{\bf
I}^{2}-\frac{1}{2}{\bf I}_{z}^{2}+\frac{\cos\left(2kz\right)}{2}\left\\{{\bf
I}_{x},{\bf I}_{y}\right\\}\right)\,,$ (55)

where anti-commutators can be rephrased by using the commutation relations:
$\\{{\bf I}_{x},{\bf I}_{y}\\}=[{\bf I}_{y},{\bf I}_{x}]+2{\bf I}_{x}{\bf
I}_{y}=-i\hbar{\bf I}_{z}+2{\bf I}_{x}{\bf I}_{y}$.

###  IV.3 Two perpendicular light beams for spin orbit coupling

In the typical situations presented in the previous subsections, the effective
Hamiltonian contains both vector and tensor shifts. Here we will consider a
scheme in which these light shifts describe the spin-orbit coupling (SOC) of
the NIST type Goldman  _et al._ (2014); Lin  _et al._ (2011) supplied with an
extra quadratic term (tensor light shift). For additional flexibility, we will
allow a small frequency difference $\delta\omega$ with $|\delta\omega/A_{\rm
HF}|\ll 1$ between two beams enabling us to remove or modify the position-
independent linear shift in the effective Hamiltonian.

The scheme involves two light beams of equal amplitudes propagating at a right
angle. One of them is propagating along the $z$ axis and is circularly
$\boldsymbol{e}_{+}$ polarized. The second beam is propagating along the $y$
axis with a linear polarization $\boldsymbol{e}_{z}$ and a frequency
difference $\delta\omega$ with respect the the first beam, so that

$\boldsymbol{\mathcal{E}}=\frac{\mathcal{E}}{\sqrt{2}}\left(\boldsymbol{e}_{+}e^{ikz}+\boldsymbol{e}_{z}e^{i(ky-\delta\omega
t)}\right)\,.$ (56)

The scalar light shift is then uniform, $\hat{H}_{{\rm
eff}(0)}=\frac{1}{4}b_{0}\mathcal{E}^{2}$, and thus will be omitted by
redefining the detuning energy $\Delta$. The vector and tensor components are

$\displaystyle\hat{H}_{{\rm eff}(1)}$
$\displaystyle=-\frac{\mathcal{E}^{2}}{8}\frac{b_{1}}{\hbar}\left({\bf
I}_{z}-\sqrt{2}{\bf I}_{xy}\left(s\right)\right),$ (57)
$\displaystyle\hat{H}_{{\rm eff}(2)}$
$\displaystyle=-\frac{\mathcal{E}^{2}}{4}\frac{b_{2}}{\hbar^{2}}\left(\frac{{\bf
I}^{2}}{6}-\frac{{\bf I}_{z}^{2}}{2}-\frac{\\{{\bf I}_{xy}\left(s\right),{\bf
I}_{z}\\}}{\sqrt{2}}\right)\,,$ (58)

where

${\bf I}_{xy}\left(s\right)\coloneqq{\bf I}_{x}\cos\left(s\right)+{\bf
I}_{y}\sin\left(s\right)\,$ (59)

and $s=ky-kz-\delta\omega t$. Alternatively the operator ${\bf
I}_{xy}\left(s\right)$ can be cast in terms of the spin raising and lowering
operator ${\bf I}_{\pm}={\bf I}_{x}\pm i{\bf I}_{y}$ as

${\bf I}_{xy}\left(s\right)={\bf I}_{+}e^{is}+{\bf I}_{-}e^{-is}\,.$ (60)

The time dependence of linear and quadratic shifts can be eliminated by a
unitary transformation describing the time-dependent spin rotation around the
$z$ axis

$\hat{U}=\exp\left(-i{\bf I}_{z}\delta\omega t/\hbar\right)\,.$ (61)

Consequently ${\bf I}_{xy}(s)$ entering $\hat{H}_{{\rm eff}(1)}$ and
$\hat{H}_{{\rm eff}(2)}$ transforms into the time-independent oprator
$\hat{U}^{\dagger}{\bf I}_{xy}(s)\hat{U}={\bf I}_{xy}(ky-kz)$, and the light
shifts take the form:

$\hat{H}_{\mathrm{eff(1)}}=-\frac{\mathcal{E}^{2}}{8}\frac{b_{1}}{\hbar}\left[{\bf
I}_{z}-\sqrt{2}{\bf I}_{xy}(ky-kz)\right]+\delta\omega{\bf I}_{z}\,$ (62)

and

$\hat{H}_{\mathrm{eff(2)}}=\frac{\mathcal{E}^{2}}{4}\frac{b_{2}}{\hbar^{2}}\left[-\frac{{\bf
I}^{2}}{6}+\frac{{\bf I}_{z}^{2}}{2}+\frac{\\{{\bf I}_{xy}(ky-kz),{\bf
I}_{z}\\}}{\sqrt{2}}\right]\,,$ (63)

where an extra linear Zeeman shift term $\delta\omega{\bf I}_{z}$ appears in
$\hat{H}_{{\rm eff}(1)}$ due to the time-dependence of the unitary
transformation $\hat{U}$. In particular, by choosing
$\delta\omega=b_{1}\mathcal{E}^{2}/8\hbar$, the linear shift $\propto{\bf
I}_{z}$ is eliminated in Eq. (62), giving

$\hat{H}_{\mathrm{eff(1)}}=\frac{\mathcal{E}^{2}}{4\sqrt{2}}\frac{b_{1}}{\hbar}{\bf
I}_{xy}(ky-kz)\,.$ (64)

The operator ${\bf I}_{xy}(ky-kz)$ defined by Eq. (59) and featured in
$\hat{H}_{{\rm eff}(1)}$ and $\hat{H}_{{\rm eff}(2)}$ represents the spin
rotating in the $xy$ plane when the atomic position changes in the $y-z$
direction. This provides a linear SOC of the NIST type Goldman  _et al._
(2014); Lin  _et al._ (2011) in the vector component of the effective
Hamiltonian $\hat{H}_{{\rm eff}(1)}$ given by Eqs. (62) or (64), and a more
complex quadratic SOC in the tensor component $\hat{H}_{{\rm eff}(2)}$ given
by Eq. (63). Note that the position dependence of both linear and quadratic
terms can be eliminated by an additional unitary operation
$\hat{U}_{1}=\exp\left[-i(kz-ky){\bf I}_{z}/\hbar\right]$ transforming ${\bf
I}_{xy}(ky-kz)$ to the position-independent spin operator ${\bf I}_{x}$. As in
the case of spin-1/2 Goldman  _et al._ (2014); Lin  _et al._ (2011), the SOC
is then represented by the spin-dependent momentum shift when the atomic
kinetic energy is included.

###  IV.4 Combining significantly different frequencies for eliminating tensor
light shifts

As exemplified above, realising purely linear light shifts, in particular the
linear SOC, can be complicated for atoms with a large spin, as the tensor
terms are present in typical configurations. It has been proposed to tune
these out by stroboscopically alternating different polarisation
configurations at a rate exceeding the atomic response rate Perlin  _et al._
(2022). As an alternative relying on continuous illumination, we propose to
use a bichromatic light with two widely different frequencies, ensuring that
the resulting vector light shifts are of the same sign. Yet, perfectly
opposing tensor terms cancel each other.

In this way, let us consider two laser fields $\alpha$ and $\beta$ with the
same configurations, but significantly different detunings $\Delta_{\alpha}$
and $\Delta_{\beta}$. Specifically the frequency difference
$\Delta_{\alpha}-\Delta_{\beta}$ should be large compared to the
characteristic frequencies of the vector and tensor light shifts, so that the
fast oscillating cross terms of the effective Hamiltonian average out. The
resulting effective Hamiltonian is thus additive, namely is given by the sum
of the effective Hamiltonians due to the separate fields $\alpha$ and $\beta$,

$\hat{H}_{\rm eff}=\hat{H}_{{\rm eff},\alpha}+\hat{H}_{{\rm eff},\beta}.$ (65)

In all situations considered in Secs. IV.1-IV.3, the linear and quadratic
light shifts are characterised by the factors
$\mathcal{E}^{2}b_{1}\equiv\mathcal{E}^{2}b_{1}(\Delta)$ and
$\mathcal{E}^{2}b_{2}\equiv\mathcal{E}^{2}b_{2}(\Delta)$. For the bichromatic
field these factors are to be replaced by
$\mathcal{E}_{\alpha}^{2}b_{1}(\Delta_{\alpha})+\mathcal{E}_{\beta}^{2}b_{1}(\Delta_{\beta})$
and
$\mathcal{E}_{\alpha}^{2}b_{2}(\Delta_{\alpha})+\mathcal{E}_{\beta}^{2}b_{2}(\Delta_{\beta})$,
respectively, where the normalization condition
$\mathcal{E}_{\alpha}^{2}+\mathcal{E}_{\beta}^{2}=1$ applies to the amplitudes
of the electric field. The quadratic light shift can be eliminated by properly
choosing detunings and amplitudes of the electric fields such that
$\mathcal{E}_{\alpha}^{2}{\rm
Re}[b_{2}(\Delta_{\alpha})]+\mathcal{E}_{\beta}^{2}{\rm
Re}[b_{2}(\Delta_{\beta})]=0$.

Figure 4:  The ratio ${\rm Re}\,\bar{b}_{1,\Sigma}/{\rm
Im}\,\bar{b}_{0,\Sigma}$ (a, upper panel) and ${\rm Re}\,\bar{b}_{1,\Sigma}$
(b, lower panel) vs the dimensionless detuning imbalance
$\bar{\delta}=\delta/(\hbar A_{\rm HF})$, where
$\bar{b}_{i,{\rm\Sigma}}=\mathcal{E}_{\alpha}^{2}\bar{b}_{i}(\Delta_{\alpha})+\mathcal{E}_{\beta}^{2}\bar{b}_{i}(\Delta_{\beta})$
(with i=0,1,2) are the dimensionless bichromatic coefficients. Here the
normalization condition $\mathcal{E}_{\alpha}^{2}+\mathcal{E}_{\beta}^{2}=1$
is used and the calculations were carried out for $\Gamma/|A_{\rm HF}|\approx
3\times 10^{-5}$. The amplitudes of bichromatic light beams were chosen such
that the condition $\mathcal{E}_{\alpha}^{2}{\rm
Re}\,\bar{b}_{2}(\Delta_{\alpha})+\mathcal{E}_{\beta}^{2}{\rm
Re}\,\bar{b}_{2}(\Delta_{\beta})=0$ holds, and the tensor light shift is
eliminated. This requirement can be always satisfied for values of
$\bar{\delta}$ shown in this plot, but not for arbitrary $\bar{\delta}$, as
beyond this region signs of ${\rm Re}\,\bar{b}_{2}(\Delta_{\alpha})$ and ${\rm
Re}\,\bar{b}_{2}(\Delta_{\beta})$ may become the same.

The analytical expressions for the vector and tensor coefficients $b_{1}$ and
$b_{2}$ make it easy to determine suitable configurations. In particular, we
note that $b_{2}$ changes sign when crossing the central hyperfine line at
$\Delta=E_{i_{I}}$, while $b_{1}$ does not (apart from a tiny spectral
window), as one can see in Fig. 2 and 3. Thus one can eliminate the quadratic
light shift by taking the detunnings $\Delta_{\alpha}$ and $\Delta_{\beta}$ to
be on different sides of the central hyperfine line. For example, one can
choose detunings to be symmetric with respect to the central hyperfine line:
$\Delta_{\alpha}=E_{i_{I}}+\delta$ and $\Delta_{\beta}=E_{i_{I}}-\delta$
222The frequency difference $\delta$ in the bichromatic scheme should not be
confused with a much smaller frequency difference $\delta\omega$ included in
the scheme considered in Sec IV.3 to have an adjustable linear Zeeman shift..

As shown in Fig. 4, for $i_{I}=9/2$ corresponding to 87Sr, at $\delta_{\rm
optimal}\simeq 3|A_{\rm HF}|$ there is a local optimal ratio between vector
shift and loss rate due to spontaneous emission, ${\rm
Re}\,\bar{b}_{1,\Sigma}/{\rm Im}\,\bar{b}_{0,\Sigma}$, is $10^{4}$, where
$\bar{b}_{i,{\rm\Sigma}}=\mathcal{E}_{\alpha}^{2}\bar{b}_{i}(\Delta_{\alpha})+\mathcal{E}_{\beta}^{2}\bar{b}_{i}(\Delta_{\beta})$
(with i=0,1,2) are the dimensionless bichromatic coefficients.

Using two light fields even further detuned, away from the hyperfine
structure, one can also have the opposite signs of $b_{2}$ and same signs of
$b_{1}$ needed for the cancelation of the tensor light shifts without
cancelling the vector ones. There are however two drawbacks: i) vector shifts
smaller in absolute value, ii) a complication in the implementation of spin-
dependent lattices or vector spin-orbit coupling. Indeed, the schemes
considered in Sections IV.2 and IV.3 involve a local phase related to optical
path differences between two beams. When superimposing two spin-dependent
lattices created with light beams at differing wavelengths, the slight
difference $k=k_{\alpha}-k_{\beta}$ in wavevectors $k_{\alpha}$ and
$k_{\beta}$ might matter. These difficulties can be overcome by working at the
relatively small detuning within the hyperfine structure, as discussed above
($\delta_{\rm optimal}\simeq 3|A_{\rm HF}|$ for $i_{I}=9/2$). In that case,
one can ensure that the two polarization lattices are in phase over the
lengthscale of the atomic cloud. For example, in the counterpropagating beam
configuration, where the optical path difference $2k(z-z_{0})$ is typically
controlled by the position $z_{0}$ of a mirror, this rephasing occurs every
$z-z_{0}=n\times\pi c/2\delta$ with $n\in\mathbb{N}$. For 87Sr, one has $\pi
c/2\delta_{\rm optimal}\simeq 10\,$cm. This is a macroscopic distance, so it
is possible to tune the optical path length accurately enough and thus realize
a good cancellation of the tensor light shift.

##  V Summary and discussions

We have analyzed the effective Hamiltonian describing the scalar, vector, and
tensor light shifts of nuclear spin states in alkaline earth atoms. Our
approach, rooted in the Dyson identity, circumvents the need for explicit
scrutiny of transitions within the excited state manifold’s specific hyperfine
states. Instead, we have derived simple analytical expressions for the
coefficients defining the scalar, vector and tensor light shift in the
effective Hamiltonian. We have explored typical examples of light
configurations capable of inducing light shifts for the nuclear spin states in
the ground electronic state manifold. Additionally, we have discerned
configurations governing SOC that encompass both linear and quadratic shifts,
alongside those characterized by purely linear SOC effects.

The breadth of our analytical findings is broad-reaching and can be applied to
other alkaline earth (or alkaline earth-like) atoms with analogous structures.
We anticipate that the straightforward analytical expressions for scalar,
vector, and tensor terms we have provided will prove helpful in guiding future
experimental endeavours and fostering a more profound theoretical
comprehension of systems subjected to off-resonance light coupling.

## ACKNOWLEDGMENTS

We gratefully acknowledge discussions with Bruno Laburthe-Tolra. This work was
supported by DAINA project of the Polish National Science Center
UMO-2020/38/L/ST2/00375 and Lithuanian Research Council S-LL-21-3. E.W.
acknowledges support of the Polish National Science Centre project
UMO-2016/22/E/ST3/00045.

##  Appendix A Dyson equation/identity

Let us consider a Hamiltonian represented in terms of the zero-order
Hamiltonian $H_{0}$ and the interaction operator $V$ as: $H=H_{0}+V$. We
define the full Green operator $G=(\Delta-H)^{-1}$ and the zero-order Green
operator $G_{0}=(\Delta-H_{0})^{-1}$, where $\Delta$ is some parameter in
energy units; in the present situation it is the detuning. The Green operator
obeys the Dyson equation (identity)

$G=G_{0}+G_{0}VG.$ (66)

The validity of the Dyson identity can be proved by rearranging the Green
operator in the following way:

$G=G_{0}(\Delta-H_{0})G=G_{0}(\Delta-H+V)G=G_{0}+G_{0}VG.$ (67)

The Dyson equation (66) is exact and valid for any $H_{0}$ and $V$, no
approximation or assumption has been made. Note also that in this work we take
$H_{0}=0$ and $V=H_{\mathrm{HF}}$. In that case, Eq. (66) reduces to Dyson
equation presented in the main text in Eq. (13).

##  Appendix B Derivation of the the equation for $\mathbf{K}$ and its
solution

###  B.1 Equation for $\mathbf{K}$

The vector operator $\mathbf{K}$ with the Cartesian components $K_{i}$ is
defined as

$K_{s}=D_{s,q}\mathcal{E}_{q}\,,$ (68)

where $D_{s,q}$ is a tensor describing the second order light induced
interaction for the ground state atoms:

$D_{s,q}=P_{g}d_{s}P_{e}GP_{e}d_{q}P_{g}\,.$ (69)

Applying the Dyson identity for $G$ given by Eq. (13) or (66), one gets the
following relation for $D_{s,q}$:

$\displaystyle D_{s,q}$
$\displaystyle=\frac{\left|d_{ge}^{2}\right|}{\Delta}P_{g}\delta_{sq}+P_{g}d_{s}P_{e}H_{\rm
HF}GP_{e}d_{q}P_{g}\,.$ (70)

Combining it with the expression (68) for $\mathbf{K}$, one has:

$\displaystyle\mathbf{K}$
$\displaystyle=\frac{1}{\Delta}P_{g}\left|d_{ge}^{2}\right|\mathcal{\boldsymbol{\mathcal{E}}}+\frac{A_{{\rm
HF}}}{\hbar\Delta}P_{g}\mathbf{d}P_{e}\left(\mathbf{I}\cdot\mathbf{J}\right)\left(1+\frac{\gamma}{\hbar^{2}}\mathbf{I}\cdot\mathbf{J}\right)GP_{e}\left(\mathbf{d}\cdot\boldsymbol{\mathcal{E}}\right)P_{g}\,.$
(71)

Our aim is get a closed equation [Eq. (84)] for $\mathbf{K}$ containing only
the spin operators $\mathbf{I}$ and the electric field amplitude
$\mathcal{\boldsymbol{\mathcal{E}}}$. For this, we will carry out the
following steps. Since $\mathbf{J}$ commutes with $P_{e}$ and gives zero when
acting on the ground state manifold ($P_{g}{\bf J}=0$), one has

$P_{g}d_{i}P_{e}\left(\mathbf{J}\cdot\mathbf{I}\right)=P_{g}d_{i}\left(\mathbf{J}\cdot\mathbf{I}\right)P_{e}=P_{g}\left[d_{i},\left(\mathbf{J}\cdot\mathbf{I}\right)\right]P_{e}=-i\hbar\epsilon_{lik}P_{g}d_{k}P_{e}I_{l},$
(72)

where the use has been made of the commutator relations Landau and Lifshitz
(1981):

$[J_{i},d_{j}]=[L_{i},d_{j}]=i\hbar\epsilon_{ijk}d_{k},$ (73)

Equation (72) can be represented in a vector form as:

$P_{g}\mathbf{d}P_{e}\left(\mathbf{J}\cdot\mathbf{I}\right)=i\hbar
P_{g}\mathbf{I}\times\mathbf{d}P_{e}\,.$ (74)

With this Eq. (71) reduces to

$\displaystyle\mathbf{K}=\frac{1}{\Delta}P_{g}\left|d_{ge}^{2}\right|\mathcal{\boldsymbol{\mathcal{E}}}+i\frac{A_{{\rm
HF}}}{\Delta}P_{g}\left[\mathbf{I}\times\mathbf{d}P_{e}+\frac{\gamma}{\hbar^{2}}\left(\mathbf{I}\times\mathbf{d}\right)P_{e}\left(\mathbf{I}\cdot\mathbf{J}\right)\right]GP_{e}\left(\mathbf{d}\cdot\boldsymbol{\mathcal{E}}\right)P_{g}\,.$
(75)

Using the commutation relation (73), one has for the Cartesian components of
the operator
$P_{g}\left(\mathbf{I}\times\mathbf{d}\right)P_{e}\left(\mathbf{I}\cdot\mathbf{J}\right)$:

$P_{g}\left(\mathbf{I}\times\mathbf{d}\right)_{i}P_{e}\left(\mathbf{I}\cdot\mathbf{J}\right)=\epsilon_{ilk}I_{l}P_{g}d_{k}P_{e}\left(\mathbf{I}\cdot\mathbf{J}\right)=i\hbar\epsilon_{ilk}\epsilon_{kpq}I_{l}P_{g}d_{q}P_{e}I_{p}\,.$
(76)

Since
$\epsilon_{ilk}\epsilon_{kpq}=\epsilon_{ilk}\epsilon_{pqk}=\delta_{ip}\delta_{lq}-\delta_{iq}\delta_{lp}$,
then

$P_{g}\left(\mathbf{I}\times\mathbf{d}\right)_{i}P_{e}\left(\mathbf{I}\cdot\mathbf{J}\right)=i\hbar
P_{g}\left(I_{l}d_{l}I_{i}-I_{l}d_{i}I_{l}\right)P_{e}=i\hbar
P_{g}\left(I_{i}I_{l}d_{l}-d_{i}I_{l}I_{l}\right)P_{e}-\hbar^{2}\epsilon_{ikl}I_{k}P_{g}d_{l}P_{e}\,.$
(77)

In vector form this reads:

$\displaystyle
P_{g}\left(\mathbf{I}\times\mathbf{d}\right)P_{e}\left(\mathbf{I}\cdot\mathbf{J}\right)=i\hbar
P_{g}\left[\mathbf{I}\left(\mathbf{I}\cdot\mathbf{d}\right)-\mathbf{I}^{2}\mathbf{d}\right]P_{e}-\hbar^{2}P_{g}\left(\mathbf{I}\times\mathbf{d}\right)P_{e}\,.$
(78)

Therefore, the relation (75) transforms to

$\displaystyle\mathbf{K}=\frac{1}{\Delta}P_{g}\left|d_{ge}^{2}\right|\mathcal{\boldsymbol{\mathcal{E}}}+i\frac{A_{{\rm
HF}}}{\Delta}P_{g}\left[\left(1-\gamma\right)\left(\mathbf{I}\times\mathbf{d}\right)+i\frac{\gamma}{\hbar}\left(\mathbf{I}\left(\mathbf{I}\cdot\mathbf{d}\right)-\mathbf{I}^{2}\mathbf{d}\right)\right]P_{e}GP_{e}\left(\mathbf{d}\cdot\boldsymbol{\mathcal{E}}\right)P_{g}\,.$
(79)

The final step is to form the operator $\mathbf{K}$ from $\mathbf{d}$ and $G$
in Eq. (79). Using Eqs. (68) and (69), one can write

$P_{g}d_{i}GP_{e}\left(\mathbf{d}\cdot\boldsymbol{\mathcal{E}}\right)P_{g}=P_{g}d_{i}P_{e}GP_{e}d_{u}P_{g}\mathcal{E}_{u}=D_{iu}\mathcal{E}_{u}=K_{i}\,.$
(80)

Equivalently in the vector form one has

$P_{g}\mathbf{d}P_{e}GP_{e}\left(\mathbf{d}\cdot\boldsymbol{\mathcal{E}}\right)P_{g}=\mathbf{K}\,.$
(81)

Analogously, the following relations hold:

$P_{g}\left(\mathbf{I}\times\mathbf{d}\right)P_{e}GP_{e}\left(\mathbf{d}\cdot\boldsymbol{\mathcal{E}}\right)P_{g}=\textbf{$\mathbf{I}$}\times\mathbf{K},$
(82)
$P_{g}(\mathbf{I}\cdot\mathbf{d})P_{e}GP_{e}\left(\mathbf{d}\cdot\boldsymbol{\mathcal{E}}\right)P_{g}=\mathbf{I}\cdot\mathbf{K}\,.$
(83)

Hence the relation (79) provides the required equation for $\mathbf{K}$:

$\displaystyle\eta\mathbf{K}$
$\displaystyle=\frac{1}{\Delta}P_{g}\left|d_{ge}^{2}\right|\boldsymbol{\mathcal{E}}+i\frac{A_{{\rm
HF}}}{\Delta}\left(1-\gamma\right)\left(\textbf{$\mathbf{I}$}\times\mathbf{K}\right)-\gamma\frac{A_{{\rm
HF}}}{\hbar\Delta}\mathbf{I}\left(\textbf{$\mathbf{I}$}\cdot\mathbf{K}\right)\,,$
(84)

where

$\eta=1-\gamma A_{{\rm HF}}\hbar i_{I}\left(i_{I}+1\right)/\Delta\,.$ (85)

###  B.2 Solution of equation for $\mathbf{K}$

As explained in the main text, we are looking for the solution of Eq. (84) in
the form

$\displaystyle\mathbf{K}$
$\displaystyle=\left[a_{0}\boldsymbol{\mathcal{E}}-i\frac{a_{1}}{\hbar}\textbf{$\mathbf{I}$}\times\boldsymbol{\mathcal{E}}+\frac{a_{2}}{\hbar^{2}}\textbf{$\mathbf{I}$}\left(\textbf{$\mathbf{I}$}\cdot\boldsymbol{\mathcal{E}}\right)\right]P_{g}\,,$
(86)

Substituting this ansatz into Eq. (84), one finds

$\displaystyle\eta\mathbf{K}$
$\displaystyle=\frac{1}{\Delta}P_{g}\left|d_{ge}^{2}\right|\boldsymbol{\mathcal{E}}+i\frac{A_{{\rm
HF}}}{\Delta}\left(1-\gamma\right)P_{g}\left[a_{0}\mathbf{I}\times\boldsymbol{\mathcal{E}}-i\frac{a_{1}}{\hbar}\mathbf{I}\times\left(\mathbf{I}\times\boldsymbol{\boldsymbol{\mathcal{E}}}\right)+\frac{a_{2}}{\hbar^{2}}\mathbf{I}\times\mathbf{I}\left(\mathbf{I}\cdot\boldsymbol{\mathcal{E}}\right)\right]+\mathbf{K}_{\gamma}\,,$
(87)

where

$\mathbf{K}_{\gamma}=-\gamma\frac{A_{{\rm
HF}}}{\hbar\Delta}\mathbf{I}\left(\mathbf{I}\cdot\mathbf{K}\right)=-\gamma\frac{A_{{\rm
HF}}}{\hbar\Delta}P_{g}\mathbf{I}\left[a_{0}\left(\mathbf{I}\cdot\boldsymbol{\mathcal{E}}\right)-i\frac{a_{1}}{\hbar}\mathbf{I}\cdot\left(\mathbf{I}\times\boldsymbol{\mathcal{E}}\right)+\frac{a_{2}}{\hbar^{2}}\left(\boldsymbol{\mathbf{I}}\cdot\boldsymbol{\mathbf{I}}\right)\left(\mathbf{I}\cdot\boldsymbol{\mathcal{E}}\right)\right]\,.$
(88)

The above relation can be rewritten as

$\mathbf{K}_{\gamma}=-\gamma\frac{A_{{\rm
HF}}}{\hbar\Delta}P_{g}\left[a_{0}+a_{1}+a_{2}\frac{\mathbf{I}^{2}}{\hbar^{2}}\right]\mathbf{I}\left(\mathbf{I}\cdot\boldsymbol{\mathcal{E}}\right)\,,$
(89)

where we used cyclic commutation relations for the components of the nuclear
spin operator ${\bf I}$, that gives
$\textbf{$\mathbf{I}$}\times\textbf{$\mathbf{I}$}=i\hbar\textbf{$\mathbf{I}$}$
(a general property of quantum angular momenta operators), and also the
relation
$\mathbf{I}\cdot\left(\mathbf{I}\times\boldsymbol{\mathcal{E}}\right)=\boldsymbol{\mathcal{E}}\cdot\left(\boldsymbol{\mathbf{I}}\times\mathbf{I}\right)=i\hbar\boldsymbol{\mathbf{I}\cdot\boldsymbol{\mathcal{E}}}$.

Rearranging the vector operators entering the r.h.s. of Eq. (87), one has

$\displaystyle\eta\mathbf{K}$
$\displaystyle=\frac{1}{\Delta}P_{g}\left|d_{ge}^{2}\right|\boldsymbol{\mathcal{E}}+i\frac{A_{{\rm
HF}}}{\Delta}\left(1-\gamma\right)P_{g}\left[\left(a_{0}+a_{1}\right)\boldsymbol{\mathbf{I}}\times\boldsymbol{\mathcal{E}}+\frac{i}{\hbar}\left(-a_{1}+a_{2}\right)\boldsymbol{\mathbf{I}}\left(\boldsymbol{\mathbf{I}}\cdot\boldsymbol{\mathcal{E}}\right)+i\frac{a_{1}}{\hbar}\boldsymbol{\mathcal{E}}\mathbf{I}^{2}\right]+\mathbf{K}_{\gamma}\,,$
(90)

where in the term proportional to $a_{2}$ relation
$\textbf{$\mathbf{I}$}\times\textbf{$\mathbf{I}$}=i\hbar\textbf{$\mathbf{I}$}$
was applied, while in the term proportional to $a_{1}$ the calculation is as
follows:

$\big{(}\mathbf{I}\times(\mathbf{I}\times\boldsymbol{\mathcal{E}})\big{)}_{k}=\epsilon_{kij}\epsilon_{mnj}I_{i}I_{m}\mathcal{E}_{n}=I_{i}\mathcal{E}_{i}I_{k}-I_{i}I_{i}\mathcal{E}_{k}=I_{k}I_{i}\mathcal{E}_{i}-I_{i}I_{i}\mathcal{E}_{k}+i\hbar\epsilon_{ikl}\mathcal{E}_{i}I_{l}=I_{k}(\mathbf{I}\cdot\boldsymbol{\mathcal{E}})-\mathcal{E}_{k}\mathbf{I}^{2}+i\hbar(\mathbf{I}\times\boldsymbol{\mathcal{E}})_{k}\,.$
(91)

Putting to Eq. (90) the expression (89) for $\mathbf{K}_{\gamma}$ and grouping
terms proportional to $\boldsymbol{\mathcal{E}}$,
$\mathbf{I}\times\boldsymbol{\mathcal{E}}$ and
$\mathbf{I}(\mathbf{I}\cdot\boldsymbol{\mathcal{E}})$, the following equation
is obtained

$\displaystyle\eta\Delta\cdot\mathbf{K}$
$\displaystyle=\left[\left|d_{ge}^{2}\right|-A_{{\rm
HF}}\left(1-\gamma\right)\mathbf{I}^{2}\frac{a_{1}}{\hbar}\right]\boldsymbol{\mathcal{E}}\,P_{g}+\frac{i}{\hbar}\Big{[}\hbar
A_{{\rm
HF}}\left(1-\gamma\right)\left(a_{0}+a_{1}\right)\Big{]}\boldsymbol{\mathbf{I}}\times\boldsymbol{\mathcal{E}}\,P_{g}+$
$\displaystyle\phantom{=}\
+\frac{1}{\hbar^{2}}\left[A_{\mathrm{HF}}\hbar\,\bigg{(}\left(1-\gamma\right)\left(a_{1}-a_{2}\right)-\gamma\,\Big{(}a_{0}+a_{1}+\tfrac{\mathbf{I}^{2}}{\hbar^{2}}a_{2}\Big{)}\bigg{)}\right]\boldsymbol{\mathbf{I}}\left(\boldsymbol{\mathbf{I}}\cdot\boldsymbol{\mathcal{E}}\right)P_{g}\,.$
(92)

Substituting Eq. (86) into (92) and comparing the factors at terms
proportional to $\boldsymbol{\mathcal{E}}$,
$\mathbf{I}\times\boldsymbol{\mathcal{E}}$ and
$\mathbf{I}(\mathbf{I}\cdot\boldsymbol{\mathcal{E}})$, one arrives at the
following closed set of equations for the coefficients $a_{0}$, $a_{1}$ and
$a_{2}$:

$\displaystyle\eta\Delta a_{0}$ $\displaystyle=\left|d_{ge}^{2}\right|-A_{{\rm
HF}}\left(1-\gamma\right)\mathbf{I}^{2}\frac{a_{1}}{\hbar}\,,$ (93)
$\displaystyle\eta\Delta a_{1}$ $\displaystyle=-\hbar A_{{\rm
HF}}\left(1-\gamma\right)\left(a_{0}+a_{1}\right)\,,$ (94)
$\displaystyle\eta\Delta a_{2}$ $\displaystyle=\hbar
A_{\mathrm{HF}}\left[\left(1-\gamma\right)\left(a_{1}-a_{2}\right)-\gamma\left(a_{0}+a_{1}+\frac{\mathbf{I}^{2}}{\hbar^{2}}a_{2}\right)\right]\,.$
(95)

The solution of this system of equations reads

$\displaystyle a_{0}=$
$\displaystyle\frac{\left|d_{ge}^{2}\right|\left(\Delta+\hbar
A_{\mathrm{HF}}\left(1-\gamma\left(\mathbf{I}^{2}/\hbar^{2}+1\right)\right)\right)}{\left(\Delta-
E_{i_{I}+1}\right)\left(\Delta-E_{i_{I}-1}\right)}\,,$ (96) $\displaystyle
a_{1}=$ $\displaystyle-\frac{\hbar
A_{\mathrm{HF}}\left|d_{ge}^{2}\right|\left(1-\gamma\right)}{\left(\Delta-
E_{i_{I}+1}\right)\left(\Delta-E_{i_{I}-1}\right)}\,,$ (97) $\displaystyle
a_{2}=$ $\displaystyle-\frac{\hbar
A_{\mathrm{HF}}\left|d_{ge}^{2}\right|\left[\hbar
A_{\mathrm{HF}}+\gamma\left(\Delta-2\hbar
A_{\mathrm{HF}}\right)-\gamma^{2}\hbar
A_{\mathrm{HF}}\left(\mathbf{I}^{2}/\hbar^{2}-1\right)\right]}{\left(\Delta-
E_{i_{I}+1}\right)\left(\Delta-E_{i_{I}}\right)\left(\Delta-
E_{i_{I}-1}\right)}\,,$ (98)

where $E_{i_{I}}$ and $E_{i_{I}\pm 1}$ featured in denominators are
eigenenergies given by Eqs. (4)-(6) The numerators can also be represented in
terms of eigenenergies, as in Eqs. (22)-(24) of the main text.

##  Appendix C A Clebsch-Gordan approach for calculating effective Hamiltonian

Throughout this paper, we used the Dyson equation method for calculating the
Green operator and thus the effective Hamiltonian. The method has the
advantage of being analytically solvable and providing a physical intuition.
Here, we summarize another method, a commonly used Clebsch-Gordan approach
which works well for numerical analysis but is not convenient for analytics.
The calculations based on the Clebsch-Gordan approach will then be compared
with the ones based on the Dyson equation approach.

We are interested in the effective Hamiltonian $\hat{H}_{\mathrm{eff}}$
describing the light-induced coupling between the atomic nuclear spin states
in the ground electronic state manifold. It is given by Eqs. (10)-(12) in the
main text, i.e.

$\hat{H}_{\mathrm{eff}}=\frac{1}{4}\mathcal{E}_{s}^{*}D_{s,q}\mathcal{E}_{q}\,,$
(99)

where:

$D_{s,q}=P_{g}d_{s}P_{e}GP_{e}d_{q}P_{g}\,,\quad\mathrm{with}\quad
G=\left(\Delta-H_{{\rm HF}}\right)^{-1}\,.$ (100)

Here $P_{g}$ and $P_{e}$ are unit projector operators onto the ground state
(${}^{1}S_{0}$) and excited state (${}^{3}P_{1}$) manifolds, respectively. In
the latter ${}^{3}P_{1}$ manifold the eigenstates of the hyperfine Hamiltonian
$H_{{\rm HF}}$ are $|j_{e}\,i_{I}\,f\,m_{f}\rangle$, where $j_{e}=1$ is
electron total angluar momentum and $i_{I}-1\leq f\leq i_{I}+1$ is the total
angular momentum of the atom. The corresponding hyperfine energies $E_{f}$ are
given by Eqs. (4) and (6) in the main text. The Green operator projected onto
this excited state manifold is diagonal in the basis of eigenstates
$|j_{e}\,i_{I}\,f\,m_{f}\rangle$ and reads explicitly:

$P_{e}GP_{e}=\sum_{f=i_{I}-j_{e}}^{i_{I}+j_{e}}\sum_{m_{f}=-f}^{f}\frac{|j_{e}\,i_{I}\,f\,m_{f}\rangle\langle
j_{e}\,i_{I}\,f\,m_{f}|}{\Delta-E_{f}}\,.$ (101)

Expressing hyperfine eigenstates $|j_{e}\,i_{I}\,f\,m_{f}\rangle$ in terms of
the bare electron and spin states
$|j_{e}\,i_{I}\,m_{j_{e}}\,m_{i}\rangle\equiv|j_{e}\,m_{j_{e}}\rangle|i_{I}\,m_{i}\rangle$
via the Clebsch-Gordan coefficients $C_{j_{e}i_{I}m_{j_{e}}m_{i}}^{fm_{f}}$,
one has

$P_{e}GP_{e}=\sum_{m_{j},m_{i},m_{j}^{\prime},m_{i}^{\prime},f,m_{f}}|j_{e}\,m_{j_{e}}\rangle|i_{I}\,m_{i}\rangle\frac{C_{j_{e}i_{I}m_{j_{e}}m_{i}}^{fm_{f}}C_{j_{e}i_{I}m_{j_{e}}^{\prime}m_{i}^{\prime}}^{fm_{f}*}}{\Delta-
E_{f}}\langle j_{e}\,m_{j_{e}}^{\prime}|\langle i_{I}\,m_{i}^{\prime}|\,.$
(102)

Finally, since the dipole operator acts only on the electronic degrees of
freedom, Eq. (100) for the tensor $D_{s,q}$ takes the form:

$D_{s,q}=P_{g}\sum_{m_{i},m_{i}^{\prime}}|i_{I}\,m_{i}\rangle
D_{s,q}^{m_{i},m_{i}^{\prime}}\langle i_{I}\,m_{i}^{\prime}|\,,$ (103)

with

$D_{s,q}^{m_{i},m_{i}^{\prime}}=\sum_{m_{j},m_{j}^{\prime},f,m_{f}}\langle
j_{g}\,m_{j_{g}}|d_{s}|j_{e}\,m_{j}\rangle\frac{C_{j_{e}i_{I}m_{j_{e}}m_{i}}^{fm_{f}}C_{j_{e}i_{I}m_{j_{e}}^{\prime}m_{i}^{\prime}}^{fm_{f}*}}{\Delta-
E_{f}}\langle
j_{e}\,m_{j_{e}}^{\prime}|d_{q}|j_{g}\,m_{j_{g}}^{\prime}\rangle\,,$ (104)

where $j_{g}=0$ and $m_{j_{g}}=0$ are, respectively, the electron angular
momentum and its projection quantum numbers for the ground state manifold
${}^{1}S_{0}$, and the unit projection operator onto this manifold is
$P_{g}=|j_{g}\,m_{j_{g}}\rangle\langle j_{g}\,m_{j_{g}}|$.

We now define the dipole moment matrix elements $\langle
j_{g}\,m_{j_{g}}|d_{s}|j_{e}\,m_{j}\rangle\equiv\langle^{1}S_{0}|d_{s}|{}^{3}P_{1},m_{j}\rangle$
featured in Eq. (104). As discussed in ref. Le Bellac and Bers (2003) (chapt.
3 p. 50, chapt. 10 p. 350-351), if we identify the two-photon polarisation
states $\ket{\sigma^{\pm}}$ to photon angular momentum states $\ket{j=1,m=\pm
1}$, and defining the photon angular momentum operator $\vec{J}$, the states
must transform as $\exp{-i\pi J_{y}}\ket{\sigma^{+}}=+\ket{\sigma^{-}}$ to
satisfy the standard raising and lowering operator conventions. The commonly
used relations between circular and linear polarization states
$\ket{\sigma^{\pm}}=(\ket{x}\pm i\ket{y})/\sqrt{2}$ do not satisfy that, and
will lead to erroneous tensorial shift calculations. Instead one must use
$\ket{\sigma^{\pm}}=(\mp\ket{x}-i\ket{y})/\sqrt{2}$. Consequently, the matrix
elements of the atomic dipole operators
$d_{x,y,z}=\vec{d}\cdot\vec{e}_{x,y,z}$ are

$\langle^{1}S_{0}|d_{x}|{}^{3}P_{1},m_{j}\rangle=\frac{d_{ge}}{\sqrt{2}}\left(\delta_{m_{j},-1}-\delta_{m_{j},1}\right)\,,$
(105)
$\langle^{1}S_{0}|d_{y}|{}^{3}P_{1},m_{j}\rangle=-\frac{d_{ge}}{i\sqrt{2}}\left(\delta_{m_{j},1}+\delta_{m_{j},-1}\right)\,,$
(106)
$\langle^{1}S_{0}|d_{z}|{}^{3}P_{1},m_{j}\rangle=d_{ge}\delta_{m_{j},0}\,.$
(107)

Using Eqs. (99), (103)-(104) and (105)-(107), the matrix elements of the
effective Hamiltonian were numerically calculated and compared with analytical
results in Fig. 2 showing a perfect agreement between numerical and analytical
results.

## References

  * Daley  _et al._ (2008) A. J. Daley, M. M. Boyd, J. Ye, and P. Zoller, Phys. Rev. Lett. 101, 170504 (2008).
  * González-Cuadra  _et al._ (2023) D. González-Cuadra, D. Bluvstein, M. Kalinowski, R. Kaubruegger, N. Maskara, P. Naldesi, T. V. Zache, A. M. Kaufman, M. D. Lukin, H. Pichler, B. Vermersch, J. Ye, and P. Zoller, Proceedings of the National Academy of Sciences 120 (2023), 10.1073/pnas.2304294120.
  * Li  _et al._ (2023) Q. Li, C. Mukhopadhyay, and A. Bayat, Phys. Rev. Res. 5, 043175 (2023).
  * Gorshkov  _et al._ (2009) A. V. Gorshkov, A. M. Rey, A. J. Daley, M. M. Boyd, J. Ye, P. Zoller, and M. D. Lukin, Phys. Rev. Lett. 102, 110503 (2009).
  * Lis  _et al._ (2023) J. W. Lis, A. Senoo, W. F. McGrew, F. Rönchen, A. Jenkins, and A. M. Kaufman, Phys. Rev. X 13, 041035 (2023).
  * Omanakuttan  _et al._ (2021) S. Omanakuttan, A. Mitra, M. J. Martin, and I. H. Deutsch, Phys. Rev. A 104, L060401 (2021).
  * Perlin  _et al._ (2022) M. A. Perlin, D. Barberena, M. Mamaev, B. Sundar, R. J. Lewis-Swan, and A. M. Rey, Phys. Rev. A 105, 023326 (2022).
  * McFerran  _et al._ (2012) J. J. McFerran, D. V. M. aes, C. Mandache, J. Millo, W. Zhang, Y. L. Coq, G. Santarelli, and S. Bize, Opt. Lett. 37, 3477 (2012).
  * Yamaguchi  _et al._ (2019) A. Yamaguchi, M. S. Safronova, K. Gibble, and H. Katori, Phys. Rev. Lett. 123, 113201 (2019).
  * Wang and Ye (2007) G. Wang and A. Ye, Phys. Rev. A 76, 043409 (2007).
  * Derevianko and Katori (2011) A. Derevianko and H. Katori, Rev. Mod. Phys. 83, 331 (2011).
  * Dzuba and Derevianko (2019) V. A. Dzuba and A. Derevianko, Journal of Physics B: Atomic, Molecular and Optical Physics 52, 215005 (2019).
  * Ludlow  _et al._ (2015) A. D. Ludlow, M. M. Boyd, J. Ye, E. Peik, and P. O. Schmidt, Rev. Mod. Phys. 87, 637 (2015).
  * Chaudhury  _et al._ (2007) S. Chaudhury, S. Merkel, T. Herr, A. Silberfarb, I. H. Deutsch, and P. S. Jessen, Phys. Rev. Lett. 99, 163002 (2007).
  * Rosenbusch  _et al._ (2009) P. Rosenbusch, S. Ghezali, V. A. Dzuba, V. V. Flambaum, K. Beloy, and A. Derevianko, Phys. Rev. A 79, 013404 (2009).
  * Deutsch and Jessen (2010) I. H. Deutsch and P. S. Jessen, Optics Communications 283, 681 (2010).
  * Le Kien  _et al._ (2013) F. Le Kien, P. Schneeweiss, and A. Rauschenbeutel, The European Physical Journal D 67, 92 (2013).
  * Goldman  _et al._ (2014) N. Goldman, G. Juzeliūnas, P. Öhberg, and I. B. Spielman, Reports on Progress in Physics 77, 126401 (2014).
  * Miyake  _et al._ (2019) H. Miyake, N. C. Pisenti, P. K. Elgee, A. Sitaram, and G. K. Campbell, Phys. Rev. Res. 1, 033113 (2019).
  * Foot (2004) C. J. Foot,  _Atomic physics_ , Vol. 7 (OUP Oxford, 2004).
  * zu Putlitz (1963) G. zu Putlitz, Zeitschrift für Physik 175 (1963), https://doi.org/10.1007/BF01375346.
  * Hernández Yanes  _et al._ (2022) T. Hernández Yanes, M. Płodzień, M. Mackoit Sinkevičienė, G. Žlabys, G. Juzeliūnas, and E. Witkowska, Physical Review Letters 129 (2022), 10.1103/physrevlett.129.090403.
  * Lin  _et al._ (2011) Y.-J. Lin, K. Jiménez-García, and I. B. Spielman, Nature 471, 83 (2011).
  * Landau and Lifshitz (1981) L. D. Landau and L. M. Lifshitz,  _Quantum Mechanics Non-Relativistic Theory, Third Edition: Volume 3_, 3rd ed. (Butterworth-Heinemann, 1981).
  * Le Bellac and Bers (2003) M. Le Bellac and A. Bers,  _Physique quantique (1ere edition)_ (EDP sciences, 2003).