#  Analysis of a Cesium lattice optical clock

A. Sharma, S. Kolkowitz, and M. Saffman  Department of Physics, University of
Wisconsin-Madison, Madison, Wisconsin 53706

(August 29, 2024)

###### Abstract

We propose and analyze a Cesium lattice optical clock (CLOC) which has the
potential for high performance and simple operation in a compact form factor
using a forbidden optical transition in Cs atoms at 685 nm. Cs atoms are
trapped in a 3D optical lattice using a magic trap wavelength of
$\lambda_{\textrm{m}}=803\leavevmode\nobreak\ \textrm{n}m$. To reduce
sensitivity to magnetic fields the atoms are probed on two cycling transitions
with equal magnitude, but opposite magnetic shifts. Operation of the clock
requires only three diode lasers at 685, 803, 852 nm and is simplified
compared to other higher-performance optical clocks that rely on alkaline
earth atoms or single trapped ions. Analysis shows a quantum noise limited
stability of $8.4\times 10^{-16}/\sqrt{\tau}$ and the potential for reaching 1
ns timing uncertainty at 1 month with realistic system parameters.

##  I Introduction

Optical atomic clocks are the most precise timekeepers ever developedLudlow
_et al._ (2015). These clocks have now reached fractional uncertainties
roughly two orders of magnitude smallerBeloy  _et al._ (2021) than the current
SI time-keeping standard based on the microwave transition in the Cs atomic
fountain clockJefferts  _et al._ (2007). High performance state-of-the-art
optical atomic clocks remain complex and laboratory-scale devices that
typically operate with low uptimes, in large part due to the complexity of the
laser systems involved in cooling and repumping the atoms. For many well
established and emerging applications of precision time keeping, reducing the
size, weight, and power (SWaP) of the clock is critical. In this work we
propose and analyze a cesium lattice optical clock (CLOC) based on a forbidden
optical transition in Cs atoms. Our proposed clock therefore takes advantage
of the relative simplicity of the level structure of alkali atoms and the
well-established tool-box that has been developed to control them, while
sacrificing some of the stability and accuracy offered by ultra-narrow
forbidden transitions in more complex atoms and ions. We show that as a
result, the proposed CLOC has the potential for performance exceeding current
microwave standards in a low SWaP form factor and with simpler operation than
current optical clocks.

Several recent reviews have surveyed progress in transportable high
performance atomic clocksDelehaye and Lacroûte (2018); Gellesch  _et al._
(2020). Possible approaches include ion clocks Cao  _et al._ (2017); Hannig
_et al._ (2019); Burt  _et al._ (2021), microwave frequency cold atom clocks
Farkas  _et al._ (2010); Esnault  _et al._ (2013); Liu  _et al._ (2017); Elvin
_et al._ (2019), and the optical frequency atomic clock based on a two-photon
transition in Rb atomsPerrella  _et al._ (2019); Maurice  _et al._ (2020);
Gerginov and Beloy (2018); Martin  _et al._ (2018, 2019); Newman  _et al._
(2019). Of these approaches, our CLOC architecture is most similar to the two-
photon Rb clock, which relies on lock-in detection of the fluorescence signal
scattered from the upper level of the clock transition as opposed to Ramsey or
Rabi interrogation. Recent demonstrations of this type of Rb clock have
reached fractional uncertainties at integration time $\tau$ of $4\times
10^{-13}/\sqrt{\tau}$ Martin  _et al._ (2018) and $4.4\times
10^{-12}/\sqrt{\tau}$ in a compact form factor with photonic integrationNewman
_et al._ (2019).

Figure 1: CLOC architecture. a) The clock transition at 685 nm connects the
$\ket{6s_{1/2},f=4}$ level to $\ket{5d_{5/2},f=6}$. The upper state decays
with near 100% branching ratio to $\ket{6p_{3/2},f=5}$ while emitting a 3490
nm photon, followed by decay back to $\ket{6s_{1/2},f=4}$ accompanied by
emission of a 852 nm photon. Atoms are cooled at 852 and 685 nm and trapped in
a 3D magic lattice with 803 nm light. The signal derived from detection of the
852 nm photons is used to stabilize the clock laser. b) Detailed atomic level
diagram showing the participating states. The atomic quantization axis is
defined by a bias magnetic field $B$. Optical pumping and repumping light
polarized along $z$ prepares the atoms in $\ket{6s_{1/2},f=4,m=\pm 4}$. The
clock probe light propagates along $x$ and is $y$ polarized. The probe drives
the closed cycle $6s_{1/2}\ket{4,\pm 4}\rightarrow 5d_{5/2}\ket{6,\pm
6}\rightarrow 6p_{3/2}\ket{5,\pm 5}\rightarrow 6s_{1/2}\ket{4,\pm 4}$. The
$\pm$ transitions at frequencies $\nu_{\pm}=\nu_{\textrm{c}}\pm\nu_{B}$ are
alternately sampled at the switching rate $f_{\textrm{s}}$.

##  II Clock operation

The proposed optical clock operates on the Cs $\ket{g}=\ket{6s_{1/2},f=4,m=\pm
4}\rightarrow\ket{e}=\ket{5d_{5/2},f=6,m=\pm 6}$ quadrupole (E2) transition
near 685 nm, as shown in Fig. 1. As we describe below, this cold atom optical
clock can be operated continuously, and requires only three laser wavelengths,
685, 803, and 852 nm, together with optical detection of emitted photons at
852 nm.

Atoms are cooled in a magneto-optical trap (MOT) with 852 nm light followed by
second stage cooling using the clock laser on the clock transition which has a
Doppler temperature of $3.0\leavevmode\nobreak\ \mu\textrm{K}$Carr (2014).
They are then loaded into a projected 3D optical latticeHuft  _et al._ (2022)
which is at a “magic wavelength” Takamoto  _et al._ (2005) for the clock
transition at a wavelength of $\lambda_{\textrm{m}}=803\leavevmode\nobreak\
\textrm{n}m$. The 3D lattice only requires optical access from a single side
which facilitates a compact form factor for the vacuum assembly. The atoms are
then pumped into the $6s_{1/2}\ket{4,\pm 4}$ stretched states using 852 nm $x$
polarized light resonant with the $6s_{1/2},f=4\rightarrow 6p_{3/2},f=4$
transition. The stretched states with $m=\pm 4$ are dark states which
minimizes heating during optical pumping.

To drive the E2 transition a 685 nm clock laser propagates along $x$ and is
linearly polarized along $y$. This configuration also couples $\ket{g}$ to
$f=6,m=\pm 2$ and $f=6,m=\pm 4$ states but these transitions are strongly
suppressed by Zeeman shifts from a bias magnetic field $B$. Sensitivity to
magnetic noise is suppressed by probing both the
$\ket{4,4}\rightarrow\ket{6,6}$ and $\ket{4,-4}\rightarrow\ket{6,-6}$
transitions at frequencies $\nu_{\pm}=\nu_{\textrm{c}}\pm\nu_{B}$, where
$\nu_{\textrm{c}}$ is the unshifted clock frequency and $\nu_{B}$ is the
Zeeman shift from the bias field. The sum $\nu_{+}+\nu_{-}=2\nu_{\textrm{c}}$
is insensitive to magnetic noise, and the two transitions are alternately
probed at a switching rate $f_{\textrm{s}}\sim 10\leavevmode\nobreak\
\textrm{kHz}$.

The excited states cycle back to the $6s_{1/2}\ket{4,\pm 4}$ states emitting
3490 and 852 nm photons. The 852 nm photons are detected and used to feedback
on the clock laser frequency and stabilize it to the atomic transition
frequency $\nu_{\textrm{c}}.$ Selection rules ensure that the decay is cycling
so that all atoms return to the $6s_{1/2}\ket{4,\pm 4}$ states. To maintain a
constant signal, atoms lost due to background collisions are periodically
replenished by loading from the MOT. By interleaving loading with alternate
probing of two samples continuous operation of the clock can be maintained
Biedermann  _et al._ (2013); Schioppo  _et al._ (2017).

###  II.1 Clock stability analysis

The CLOC architecture uses detection of the scattered 852 nm signal for
stabilization of the 685 nm laser to the quadrupole transition. The clock
uncertainty is fundamentally limited by detector shot noise resulting in a
stability in fractional frequency units of

$\sigma(\tau)=\frac{\Delta\nu}{\nu_{\textrm{c}}}\sqrt{\frac{1}{\dot{N}\tau}}.$
(1)

Here $\Delta\nu$ is the full width at half maximum (FWHM) transition
linewidth, $\nu_{\textrm{c}}$ is the clock frequency, $\dot{N}$ is the number
of detection events per second, and $\tau$ is the integration time.

To estimate the shot noise limited performance we use
$\nu_{\textrm{c}}=4.376\times 10^{14}\leavevmode\nobreak\ \textrm{Hz}$ for the
clock transition frequency. At unit saturation
$\Delta\nu=\sqrt{2}/(2\pi\tau_{\textrm{a}})=1.75\times
10^{5}\leavevmode\nobreak\ \textrm{Hz}$, with
$\tau_{\textrm{a}}=1.28\leavevmode\nobreak\ \mu\textrm{s}$ the $5d_{5/2}$
radiative lifetime. The detection rate is

$\dot{N}=\eta_{\textrm{c}ol}\eta_{\textrm{d}et}\frac{1}{4\tau_{\textrm{a}}}N_{\textrm{a}}$
(2)

where $N_{\textrm{a}}$ is the number of atoms, the optical efficiency is
$\eta_{\textrm{c}ol}=0.2$ corresponding to collection optics with numerical
aperture NA$=0.8$, the detector quantum efficiency is
$\eta_{\textrm{d}et}=0.9$ for a Si detector at 852 nm, and
$1/4\tau_{\textrm{a}}=1.95\times 10^{5}\leavevmode\nobreak\ \textrm{s}^{-1}$
is the scattering rate per atom at unit saturation.

With $N_{\textrm{a}}=6.6\times 10^{6}$ atoms (see Sec. III.2 for details of
the trap design) and $\dot{N}=2.3\times 10^{11}\leavevmode\nobreak\
\textrm{s}^{-1}$ eq. (1) evaluates to

$\sigma(\tau)=\frac{8.4\times 10^{-16}}{\sqrt{\tau}}.$ (3)

This is the standard quantum limit for the clock stability. The corresponding
frequency uncertainty at 1 s is
$\sigma\nu_{\textrm{c}}=0.37\leavevmode\nobreak\ \textrm{Hz}$.

In addition to the quantum limit the short term stability is impacted by the
intermodulation effectAudoin  _et al._ (1991) at frequencies of
$2f_{\textrm{s}}$ corresponding to switching between the $\nu_{\pm}$
transitions and $2f_{\textrm{m}}$ corresponding to the modulation frequency
used for lock in detection of the photodetector current. We anticipate using
$f_{\textrm{s}}=10\leavevmode\nobreak\ \textrm{kHz}$ and
$f_{\textrm{m}}=100\leavevmode\nobreak\ \textrm{kHz}$. The intermodulation
frequency contribution to the frequency deviation is

$\sigma^{(\textrm{I}M)}=\frac{\left[S(2f)\right]^{1/2}}{2\nu_{c}\sqrt{\tau}}$
(4)

with $S(2f)$ the power spectral density of frequency noise at $2f$ in units of
$\textrm{Hz}^{2}/Hz$. There are two primary contributions to $S_{y}$: the
frequency noise of the probe laser and the shot noise of the detector signal.
For the probe laser we assume that at frequencies above 20 kHz there is white
frequency noise of $1\leavevmode\nobreak\ \textrm{Hz}^{2}/Hz$. This level of
frequency noise has been achieved with optically injected slave laser
diodesLiu  _et al._ (2021). The correspondsing local oscillator contribution
to the stability is

$\sigma^{(\textrm{I}M,LO)}=1.1\times 10^{-15}/\sqrt{\tau}.$ (5)

The detector photocurrent is $I=e\dot{N}=37\leavevmode\nobreak\ \textrm{n}A$.
This implies a shot noise limited signal to noise ratio of SNR$=3.4\times
10^{5}\leavevmode\nobreak\ \textrm{(}Hz)^{-1/2}$. We estimate the
corresponding frequency deviation as

$\sigma^{(\textrm{I}M,shot)}=\frac{\left[\Delta\nu/{\textrm{S}NR}\right]^{1/2}}{2\nu_{c}\sqrt{\tau}}=8.2\times
10^{-16}/\sqrt{\tau}.$ (6)

Adding eqs. (3, 5, 6) in quadrature gives a short term stability limit of
$\sigma=1.6\times 10^{-15}/\sqrt{\tau}$. The achievable stability at long
times depends on the level of control over systematic shifts, and is discussed
in Sec. IV.

##  III Experimental Design

Figure 2: Continuous clock operation using two cells that are alternately
probed at frequencies $\nu_{\pm}$. A periodic MOT and optical pumping phase
replenishises lattice atoms.

We envision a clock based on alternating interrogation between two separate
vacuum cells, or time alternating interrogation of a single atomic sample. The
timing sequence is shown in Fig. 2. A pre-cooled, and occupied atomic lattice
will be prepared in each cell. As discussed in the following the lattice will
be addressed with alternating probing and cooling cycles. Using two cells and
out of phase probing and cooling a continuous clock signal will be available.
The Zeeman shifted $\pm$ transitions are alternately probed at a switching
rate of $f_{\textrm{s}}$. Each transition is detected by modulating the probe
light at frequency $f_{\textrm{m}}\gg f_{\textrm{s}}$ followed by lock in
detection of the scattered 852 nm light to provide an error signal for the
laser frequency. On a time scale of a few seconds one of the lattices will be
refreshed from a MOT followed by second stage cooling and optical pumping. In
this way a continuous clock measurement can be maintained indefinitely. The
continuous and high bandwidth photon signal relaxes the requirements on the
clock laser, potentially obviating the need for a ultra-high-stability bulky
and vibration sensitive cavity.

###  III.1 Laser Cooling

We propose a two-step cooling process starting with a standard magneto-optical
trap (MOT) at 852 nm driving the $6s_{1/2},f=4\rightarrow 6p_{3/2},f=5$
transition. This can be performed with a single diode laser by adding a 8.94
GHz sideband for repumping from $f=3$. Cold clouds of $>10^{6}$ atoms can be
prepared in $\ll 1\leavevmode\nobreak\ \textrm{s}$ using this approach.

Second stage cooling uses frequency shifted clock laser light tuned to the red
of the $6s_{1/2},f=4\rightarrow 5d_{5/2},f=6$ transition which has a Doppler
temperature of $3.0\leavevmode\nobreak\ \mu\textrm{K}$. Preliminary
demonstrations of cooling on this transition as well as analysis of sub-
Doppler cooling can be found in Carr (2014). Decay from the excited state
proceeds via the closed cascade $5d_{5/2},f=6\rightarrow
6p_{3/2},f=5\rightarrow 6s_{1/2},f=4$ with essentially 100% branching ratio.
The depumping fraction is $\sim 10^{-7}$ (see Sec. III.4) so a repumper is not
needed. We conservatively predict cooling to $T\sim 1\leavevmode\nobreak\
\mu\textrm{K}$ which corresponds to about $10$ times the recoil temperature.

###  III.2 Trap design

We propose to use a three dimensional array of dark optical traps that confine
atoms near zeros of the optical intensity. Such arrays are traditionally
prepared using six counterpropagating beams. The implementation can be
simplified using a a recently demonstrated projected “hole” arrayHuft  _et
al._ (2021). The array is prepared using a partially transmitting optical mask
followed by $4f$ Fourier filtering to create a two-dimensional pattern of near
Gaussian holes in a uniform background. The magnification of the optics used
to project the light into the vacuum cell is chosen to give an array period of
$d$. The Talbot effect leads to repetition of the planar array at axial
spacings of the Talbot length
$L_{\textrm{T}}=\lambda/[1-\sqrt{1-\lambda^{2}/d^{2}}]$, where $\lambda$ is
the wavelength of the trapping light. Taking $d=0.9\leavevmode\nobreak\
\mu\textrm{m}$ and $\lambda=0.803\leavevmode\nobreak\ \mu\textrm{m}$ gives
$L_{\textrm{T}}=1.46\leavevmode\nobreak\ \mu\textrm{m}$. We assume a cube of
size $w=250\leavevmode\nobreak\ \mu\textrm{m}$ on a side. This gives a total
of $w^{3}/(d^{2}L_{\textrm{T}})=1.3\times 10^{7}$ trapping sites. Assuming a
single atom filling fraction of 0.5 we obtain $N_{\textrm{a}}=6.6\times
10^{6}$ atoms. In this regime broadening due to atomic collisions is strongly
suppressed.

At the magic wavelength of $\lambda_{\textrm{m}}=803\textrm{ nm}$ the
$6s_{1/2}$ ground state has a polarizability of
$\alpha_{0}=-374.\leavevmode\nobreak\ \textrm{\AA}^{3}$. With the trap region
size of $250\leavevmode\nobreak\ \mu\textrm{m}$ on a side an optical power of
$P=2\leavevmode\nobreak\ \textrm{W}$ gives a trap depth of
$18\leavevmode\nobreak\ \mu\textrm{K}$ which is about 145 times the recoil
energy. This power level is available from semiconductor tapered amplifier
devices which provide power up to 3 W at 800 nm. In order to keep the optical
power requirement in the range of readily available tapered amplifier
components we propose a compact lattice and a relatively high atomic density.
This in turn implies a large optical density which will inhibit full optical
pumping and also lead to depumping from reabsorption of scattered 852 nm
photons during clock operation. Cycling photons emitted during the
$6p_{3/2}\ket{5,5}\rightarrow 6s_{1/2}\ket{4,4}$ decay are reabsorbed on the
same transition with cross section $\sigma_{0}=3\lambda_{852}^{2}/(2\pi).$
Since the emitted photons are distributed over the full solid angle they can
be reabsorbed by driving $\pi$ and $\sigma_{-}$ transitions to
$6p_{3/2}\ket{5,4}$ and $\ket{5,3}.$ The absorption cross section for these
transitions is $\sigma_{0}\times 1/5$ and $\sigma_{0}\times 1/45$. The
unwanted transitions are further suppressed by Zeeman shifts relative to
$\ket{5,5}.$ Assuming a bias magnetic field of $1$ mT We estimate the
effective optical depth of the medium for a photon emitted at the center of
the lattice as 6.5 with the optical depth due to depolarizing $\pi$ and
$\sigma_{-}$ transitions only 1.3. The depolarizing rate can in principle be
suppressed further at the cost of higher optical power to support a more
dilute optical ensemble. We defer a more detailed analysis of the impact of
photon reabsorption on clock performance to future work.

In order to suppress Doppler broadening and recoil heating the clock should be
operated in the resolved sideband limit with a small Lamb-Dicke parameter
$\eta$. The excited state linewidth is $\Delta\nu/2\pi=124\leavevmode\nobreak\
\textrm{kHz}$ which we require to be small compared to the trap vibrational
frequency along the axis of the probe beam. To achieve a high vibrational
frequency we propose to use an additional standing wave 1D lattice with cavity
build up, as shown in Fig. 3. With a 1D lattice beam power of 2.2 W, cross
sectional area of $w^{2}$, and modest cavity build up factor of 50 we obtain
an axial vibrational frequency of 1.05 MHz. When the clock transition is
probed at unit saturation the relative absorption at the first vibrational
sideband is a negligible 0.007. The Lamb-Dicke parameter is 0.078 which adds
further suppression to the sideband excitation.

The projected lattice and 1D lattice beams propagate in orthogonal directions
with the same electric field polarization. The beams can be detuned by a few
MHz to remove interference effects, without compromising the long time clock
stability (see Sec. IV.3).

Figure 3: Absorption profile of clock transition at unit saturation. Red and
blue vibrational sidebands are shown for the parameters given in the text.

###  III.3 Quadrupole clock transition

The quadrupole clock transition will be probed with 685 nm light that
propagates along $x$ and is $y$ polarized. With this geometry transitions
$m\rightarrow m^{\prime}=m-2,m+2$ are allowed. The matrix element for the
transitions $\ket{4,\pm 4}\rightarrow\ket{6,\pm 6}$ is larger than that for
the transitions $\ket{4,\pm 4}\rightarrow\ket{6,\pm 2}$ by a factor of
$3\sqrt{55}=22.2.$ Thus the transitions between stretched states are strongly
favored and the saturation intensity is $570.\leavevmode\nobreak\
\textrm{m}W/cm^{2}$. The unwanted transitions are further suppressed by tensor
shifts in the $5d_{5/2}$ level and by Zeeman shifts. For a beam of cross
sectional area $w^{2}$ matching the size of the atomic lattice this
corresponds to a modest probe power of 0.36 mW.

Since the clock transition is cycling the depumping rate into the ground $f=3$
level is dominated by scattering from $5d_{5/2},f=5$ at an approximate rate of

$\gamma_{\textrm{d}p}=\frac{\gamma}{2}\frac{1}{2+4\Delta_{65}^{2}/\gamma^{2}},$

where we have assumed unit saturation of the clock transition and the excited
state hyperfine splitting between $f=5$ and $6$ is $\Delta_{65}=2\pi\times
127\leavevmode\nobreak\ \textrm{MHz}$. We find
$\gamma_{\textrm{d}p}/\gamma=1.2\times 10^{-7}$ which implies a depumping rate
per atom of $0.004\leavevmode\nobreak\ \textrm{s}^{-1}$ or a depumping time of
$270\leavevmode\nobreak\ \textrm{s}.$ This is significantly longer than the
expected vacuum lifetime in current generation small form factor UHV cells, so
the cycle time for replenishing the atomic lattice will be governed by vacuum
quality, not depumping.

###  III.4 Magic optical trapping

The dynamic Stark shifts of the clock states can be readily calculated using a
standard sum over states formalism accounting for all electric dipole allowed
transitions. Matrix elements and transition frequencies were taken from
Barakhshan  _et al._ ; Sansonetti (2009). A core correction to the scalar
polarizabilities was included using the value of
$\alpha_{\textrm{c}}=15.84\leavevmode\nobreak\ a_{0}^{3}$ from Safronova  _et
al._ (2016), and $a_{0}$ is the Bohr radius. The hyperfine resolved
polarizabilities of the $5d_{5/2}$ states were calculated using

$\displaystyle\alpha(f,m)$ $\displaystyle=$
$\displaystyle\alpha_{0}^{(j)}+\frac{3m_{f}^{2}-f(f+1)}{f(2f-1)}\alpha_{2}^{(f)}$
(7) $\displaystyle\alpha_{2}^{(f)}$ $\displaystyle=$
$\displaystyle\alpha_{2}^{(j)}\frac{3K(K-1)-4f(f+1)j(j+1)}{(2f+3)(2f+2)j(2j-1)}$
(8)

where $\alpha_{0}^{(j)},\alpha_{2}^{(j)}$ are calculated in the fine structure
basis as outlined above and $K=f(f+1)+j(j+1)-I(I+1)$, with $I=7/2$ the nuclear
spin.

Figure 4: Polarizability of ground $6s_{1/2}$ and excited $5d_{5/2}$ states.
The solid blue, dashed yellow, and dotted green lines show the
$5d_{5/2}\alpha_{0}^{(j)}$, $5d_{5/2}\alpha_{2}^{(j)}$, and
$6s_{1/2}\alpha_{0}^{(j)}$ polarizabilities respectively.

Figure 4 shows that there are a large number of magic wavelengths for trapping
the $6s_{1/2},5d_{5/2}$ states between 0.6 and 0.9 $\mu\textrm{m}$ at which
the states have a negative polarizability suitable for dark optical traps.
Magic conditions occur for wavelengths near where $\alpha_{0}$ is the same for
the ground and excited state. In order to minimize the optical power needed
for trapping we choose the magic condition with the largest polarizability
near $0.8\leavevmode\nobreak\ \mu\textrm{m}$. The Zeeman state resolved
polarizabilities for the $5d_{5/2},f=6$ level are shown in Fig. 5. The
resonance behavior of the $5d_{5/2}$ polarizability is due to the transition
to the $5f$ levels at 808 nm.

The tensor polarizability of the $5d_{5/2}$ level results in differential
shifts between the Zeeman substates of the $f=6$ level. To estimate the
severity of these shifts we assume that the atoms are cooled to
$T=1\leavevmode\nobreak\ \mu\textrm{K}$, which is $3\times$ lower than the
$3\leavevmode\nobreak\ \mu\textrm{K}$ Doppler temperature of the
$6s_{1/2}-5d_{5/2}$ cooling transition. Calculations have confirmed the
presence of sub-Doppler cooling on this transitionCarr (2014). The fractional
spread of the polarizability at the magic wavelength for the $m_{f}=0$ excited
state is $\eta=|[\alpha(6,6)-\alpha(6,0)]/\alpha(6,0)|=0.39$. This implies
that an atom in a dark trap sees a spread of transition frequencies of
$\eta\frac{k_{B}T}{2\pi\hbar}=8.1\leavevmode\nobreak\ \textrm{kHz}$. This is a
factor of 15 smaller than the excited state linewidth and has a negligible
impact on the short term clock stability. However, the tensor shifts result in
a strong sensitivity to trap light intensity that impacts long term
systematics. For this reason we will only use transitions between the
stretched states $\ket{4,\pm 4}\rightarrow\ket{6,\pm 6}$ which have a magic
condition at 803 nm as shown in Fig. 5.

###  III.5 Magic magnetic field

Magic magnetic field conditions can be found from a Breit-Rabi type analysis
applied to the $6s_{1/2}$ and $5d_{5/2}$ levels. The $6s_{1/2},f=4$ level has
$g_{f}=1/4$ and $5d_{5/2},f=6$ has $g_{f}=1/2.$ We therefore expect magic
conditions to occur at low fields for transitions from $f=4$ to $f^{\prime}=6$
for Zeeman substates

$\displaystyle m_{f}=\pm 2\rightarrow m_{f}^{\prime}=\pm 1,$ $\displaystyle
m_{f}=\pm 4\rightarrow m_{f}^{\prime}=\pm 2.$

There are also magic conditions at higher field values that rely on
substantial mixing of the $f^{\prime}=5,6$ levels.

We have identified multiple magic magnetic field values for a variety of pairs
of ground and excited Zeeman states. However, while it is possible to
selectively tune the probe laser to a magnetic field dependent pair of states
at a magic field condition, the atomic decay back to the ground state is not
cycling with respect to the ground Zeeman substate which would lead to rapid
depumping of the probed transition. Repumping of the clock states is an
option, but this would result in additional dead time and increased
complexity.

We therefore propose to instead probe the $6s_{1/2}\ket{4,\pm 4}\rightarrow
5d_{5/2}\ket{6,\pm 6}$ transitions with a beam that propagates along $x$ and
is polarized along $y$. This geometry couples the ground states $\ket{4,\pm
4}$ to the excited states $\ket{6,\pm 2},\ket{6,\pm 6}.$ Using a moderate
magnetic field we can ensure that only the $\ket{6,\pm 6}$ states are excited
which decay in a cycling fashion back to $\ket{4,\pm 4}$. With a magnetic
field oriented along $z$ these states do not suffer Zeeman mixing so the
sensitivity of the clock transitions to magnetic noise is

$\delta\nu_{\pm}(m)=\pm(6g_{5d_{5/2},f=6}-4g_{6s_{1/2},f=4})\frac{\mu_{B}}{2\pi\hbar}\delta
B=\pm\frac{\mu_{B}}{\pi\hbar}\delta B$

with $\mu_{B}$ the Bohr magneton and $\delta B$ the magnetic field
fluctuation. A moderate magnetic field stability of $\delta
B=10^{-8}\leavevmode\nobreak\ \textrm{T}$ results in $\delta\nu_{\pm}=\pm
280\leavevmode\nobreak\ \textrm{Hz}$. This is much less than the excited state
linewidth and gives a negligible contribution to the short term clock
stability, but is a limiting factor for long term stability. To circumvent
this issue we will probe both stretched states and use the average frequency
which is insensitive to the magnetic field strength.

Figure 5: Polarizability of ground and excited states near $0.8\leavevmode\nobreak\ \mu\textrm{m}$. The dashed line shows the magic wavelength at $\lambda_{m}=0.8033\leavevmode\nobreak\ \mu\textrm{m}$ for the transition $6s_{1/2}\rightarrow 5d_{5/2},f=6,m_{f}=6$. The inset shows the differential polarizability as a function of detuning from the magic point.  Table 1: Systematic uncertainties affecting the clock transition. The right hand column gives the stability requirement for 1 ns timing uncertainty at 30 days. The requirement on magnetic field variation corresponds to probing a single transition on two separate clock ensembles. Alternating probing of $\nu_{\pm}$ on a single or two ensembles removes this requirement.  shift | fractional coefficient $\beta/\nu_{\textrm{c}}$ | 1 ns stability at 30 days
---|---|---
probe beam power | $2.4\times 10^{-16}\leavevmode\nobreak\ /\%$ | 1.6 %
lattice laser frequency shift from magic value | $1.8\times 10^{-17}/\textrm{MHz}$ | 21.6 MHz
magnetic field spatial variation | $3.2\times 10^{-12}\leavevmode\nobreak\ \textrm{/}(10^{-7}\leavevmode\nobreak\ \textrm{T})$ | $12\leavevmode\nobreak\ \textrm{pT}$
dc electric field | $-1.1\times 10^{-14}\leavevmode\nobreak\ /\textrm{(V/m)}$ | $0.036\leavevmode\nobreak\ \textrm{V/m}$
blackbody radiation | $2.7\times 10^{-16}\leavevmode\nobreak\ /\textrm{K}$ | $1.4\leavevmode\nobreak\ \textrm{K}$

##  IV Systematic uncertainties

The clock frequency is sensitive to parameters $a_{i}$ according to
$\delta\nu=\sum_{i}\beta_{i}\delta a_{i}$ where
$\beta_{i}=d\nu_{\textrm{c}}/da_{i}$ and $a_{i}$ are environmental parameters.
The required stability of the $a_{i}$ to ensure a timing uncertainty of
$\delta t$ after integration time $\tau$ is found from

$\delta t=\frac{\delta\nu}{\nu_{\textrm{c}}}\tau=\frac{\sum_{i}\beta_{i}\delta
a_{i}}{\nu_{\textrm{c}}}\tau$

so

$\delta a_{i}=\frac{\nu_{\textrm{c}}}{\beta_{i}\tau}\delta t_{i}$

and $\delta t=\sum_{i}\delta t_{i}$. Table 1 lists sources of systematic
uncertainty, and the required stability for each parameter to reach a timing
uncertainty of 1 ns at 30 days, corresponding to a long-term precision and
accuracy of $3.9\times 10^{-16}$. The fractional coefficients listed in Table
1 are $\beta_{i}/\nu_{\textrm{c}}$. Precise, but not unrealistic, control of
environmental parameters will be needed to reach this level of performance.

###  IV.1 Magnetic field

As described above the clock transition will be probed on
$6s_{1/2}\ket{4,4}\rightarrow 5d_{5/2}\ket{6,6}$ and
$6s_{1/2}\ket{4,-4}\rightarrow 5d_{5/2}\ket{6,-6}$. With a bias magnetic field
$B$ along the $z$ axis these states have a purely linear Zeeman shift. The
differential shifts of the $\ket{4,\pm 4}\rightarrow\ket{6,\pm 6}$ transitions
are

$\delta\nu_{\pm}=\pm\frac{\mu_{B}}{\pi\hbar}B.$

At fixed magnetic field the average of the shifts is zero, which can be used
to determine the unshifted clock frequency.

We envision two different approaches to probing both transitions. Applying $z$
polarized pumping light on the $6s_{1/2},f=4\rightarrow 6p_{3/2},f=3$
transition, together with a repump sideband driving $6s_{1/2},f=3\rightarrow
6p_{3/2},f=4$, we pump the atoms into a 50/50 mix of $\ket{4,4}$ and
$\ket{4,-4}$. Using a modulator to rapidly switch the probe frequency seen by
the atoms between $\nu_{\textrm{c}}+\delta\nu_{+}$ and
$\nu_{\textrm{c}}+\delta\nu_{-}$ we can acquire data from both transitions,
and average to find the unshifted clock frequency. In this approach the atom
number for each transition is reduced by a factor of two so $\sigma$ in eq.
(1) is increased by $\sqrt{2}$. If the alternation is implemented at frequency
$f_{\textrm{s}}\sim 10\leavevmode\nobreak\ \textrm{kHz}$ this will introduce
intermodulation noise at $2f_{\textrm{s}}$Audoin  _et al._ (1991). This has
has been accounted for in Sec. II.1.

Alternatively we may probe two different atomic samples, either in nearby
regions of one vacuum cell, or in separate vacuum cells, each on one of the
two stretched state transitions. In this case there is a potential sensitivity
to the difference in the magnetic field at each location. If the magnetic
field variation is $\pm\delta B/2$ the average of the measured frequencies
will be

$\nu_{\textrm{c}}+\delta\nu=\nu_{\textrm{c}}+\frac{\mu_{B}}{2\pi\hbar}\delta
B.$

The sensitivity coefficient is $\beta=1400\leavevmode\nobreak\
\textrm{Hz}/(10^{-7}\leavevmode\nobreak\ \textrm{T})$. This effect can be
suppressed by alternating back and forth between which region probes which
transition and again taking the average. However, this will not cancel the
effect of time varying magnetic field gradients.

###  IV.2 DC electric field

The static scalar polarizability of the Cs 6s ground state has been accurately
measured to beAmini and Gould (2003)
$\alpha_{0}^{(6s)}=59.4\leavevmode\nobreak\ \textrm{\AA}^{3}.$ The tensor
polarizability of the ground state is negligible for the following analysis.
Accounting for scalar and tensor contributions we calculate the static
polarizability of the $5d_{5/2}\ket{6,6}$ state to be
$\alpha_{5d}=30.8\leavevmode\nobreak\ \textrm{\AA}^{3}$ The sensitivity
coefficient for the clock transition is $\beta=\leavevmode\nobreak\
4.7\leavevmode\nobreak\ \textrm{Hz/(V/m)}.$

###  IV.3 Lattice wavelength

As is shown in Fig. 5 the differential polarizability of the excited and
ground states is approximately $1.4\times 10^{-4}\leavevmode\nobreak\
\textrm{\AA}^{3}$ per MHz of lattice laser shift from the magic point. This
translates into a sensitivity coefficient $\beta=0.0078\leavevmode\nobreak\
\textrm{Hz}/\textrm{MHz}$ assuming the atoms have a temperature of
$1\leavevmode\nobreak\ \mu\textrm{K}$.

###  IV.4 Probe beam power

Figure 6: Polarizability of ground and excited states near the clock
wavelength.

The clock frequency is sensitive to the probe beam power due to differential
Stark shifts of the ground and excited states. This is a relatively small
effect since the clock wavelength is close to a magic wavelength as can be
seen in Fig. 6. The differential polarizability at the clock wavelength is
$\delta\alpha=-57.5\leavevmode\nobreak\ \textrm{\AA}^{3}$. For a probe beam
intensity that provides unit saturation we find a sensitivity coefficient per
% change of the probe power of $\beta=0.11\leavevmode\nobreak\
\textrm{Hz}/\%$.

###  IV.5 Blackbody radiation

Accurate calculations of dynamic Stark shifts due to blackbody radiation can
be found in Farley and Wing (1981). At $T=300\leavevmode\nobreak\ \textrm{K}$
the shifts of the $6s_{1/2}$ and $5d_{5/2}$ states are -3.589 and 5.315 Hz,
respectively. The differential shift is thus
$\delta\nu=8.904\leavevmode\nobreak\ \textrm{Hz}$. Since the blackbody shift
scales as $T^{4}$ the sensitivity at 300 K is $\beta=0.12\leavevmode\nobreak\
\textrm{Hz}/\textrm{K}$.

##  V Discussion

We have introduced a cesium lattice optical clock using a forbidden quadrupole
transition. System complexity is considerably reduced compared to a state of
the art atomic optical clock using alkaline earth atoms. Predicted performance
is intermediate between microwave clocks and existing atomic optical clocks.
The use of a high atom number, combined with magic optical trapping, provides
a path to high stability that exceeds that of Cs fountain clocks at short
timesHeavner  _et al._ (2014), as well as competitive long term stability, all
in a more compact package.

The proposed cesium clock requires three lasers at wavelengths of 685, 803,
and 852 nm, all of which are available from semiconductor sources. Acousto-
optic and electro-optic modulators can be used to provide the frequency shifts
needed for cooling, optical pumping, and repumping operations. Cesium can be
conveniently laser cooled directly from a room temperature source in a compact
vacuum package. We envision implementing the lattice build-up resonator
described in Sec. III.2 in a compact monolithic form factor using a fiber
Fabry-Perot resonator. To complete the system, an optical frequency comb is
referenced to the 685 nm probe laser. Small form factor, yet high performance
architectures for frequency combsSinclair  _et al._ (2015) together with
progress in high stability optical oscillators without the need for bulky and
vibration sensitive cavity stabilizationLoh  _et al._ (2020) point towards the
potential for a complete clock system in a compact form factor. Progress in
engineering of state of the art Sr clocks has also resulted in high
performance in a reduced form factorOhmae  _et al._ (2021); Guo  _et al._
(2021). Compared to the CLOC, Sr requires five laser wavelengths (and often
multiple lasers per wavelength), as well as a more complex atomic source based
on an oven and a Zeeman slower. The narrow width of the Sr clock transition,
while advantageous for reaching the highest possible performance, also implies
the need for a more stable probe laser, which in turn requires relatively
bulky and vibration sensitive cavity stabilization.

Rubidium has a similar atomic structure to Cs and may in principle be used in
an analogous fashion. The corresponding clock transition in 87Rb is
$5s_{1/2},f=2\rightarrow 4d_{5/2},f=4$ at 517 nm with cycling decay via
$5p_{3/2},f=3$ back to the ground state. Analysis shows that there are several
magic wavelengths for both dark and bright optical traps. However, the atomic
parameters are not as favorable as they are in Cs. The lifetime of the
$4d_{5/2}$ level is 89 ns, corresponding to a linewidth of 1.8 MHz, which is
more than an order of magnitude broader than the upper clock state in Cs. This
is partially compensated for by the increase in $\dot{N}$ in eq. (1). Assuming
the same number of atoms and detection efficiency as for Cs we find
$\sigma(\tau)=2.4\times 10^{-15}/\sqrt{\tau}$ which is a factor of 2.9 larger
than for Cs. While still attractive for a clock implementation the requirement
of resolved sideband excitation to suppress Doppler broadening would imply
more than an order of magnitude higher trap vibrational frequencies. Since the
optical trap frequencies scale with the square root of the optical intensity
times the polarizability, a trap intensity more than two orders of magnitude
higher would be needed, which would be difficult to implement.

Finally we note that the atomic sample with individually trapped Cs atoms
arranged at $\mu\textrm{m}$ scale spacings is closely related to systems for
quantum simulation and computation that use Rydberg interactions to generate
entanglement. Recent advances have demonstrated multi-particle, highly
entangled states with both analogEbert  _et al._ (2015); Zeiher  _et al._
(2015); Omran  _et al._ (2019) and digital control modalitiesGraham  _et al._
(2021). Combining the architecture described here with Rydberg mediated
entanglement may provide a feasible route towards quantum enhanced clock
performance.

##  VI Acknowledgements

We gratefully acknowledge Andrew Jayich for a critical reading of the
manuscript. This material is based upon work supported by the Air Force Office
of Scientific Research under award number FA9550-21-1-0034.

## References

  * Ludlow  _et al._ (2015) A. D. Ludlow, M. M. Boyd, J. Ye, E. Peik, and P. O. Schmidt, “Optical atomic clocks,” Rev. Mod. Phys. 87, 637 (2015).
  * Beloy  _et al._ (2021) K. Beloy, M. I. Bodine, T. Bothwell, S. M. Brewer, S. L. Bromley, J.-S. Chen, J.-D. Deschênes, S. A. Diddams, R. J. Fasano, T. M. Fortier, Y. S. Hassan, D. B. Hume, D. Kedar, C. J. Kennedy, I. Khader, A. Koepke, D. R. Leibrandt, H. Leopardi, A. D. Ludlow, W. F. McGrew, W. R. Milner, N. R. Newbury, D. Nicolodi, E. Oelker, T. E. Parker, J. M. Robinson, S. Romisch, S. A. Schäffer, J. A. Sherman, L. C. Sinclair, L. Sonderhouse, W. C. Swann, J. Yao, J. Ye, and X. Zhang, “Frequency ratio measurements at 18-digit accuracy using an optical clock network,” Nature 591, 564 (2021).
  * Jefferts  _et al._ (2007) S. R. Jefferts, T. P. Heavner, T. E. Parker, and J. H. Shirley, “NIST Cesium fountains - current status and future prospects,” Act. Phys. Pol. A 112, 759 (2007).
  * Delehaye and Lacroûte (2018) M. Delehaye and C. Lacroûte, “Single-ion, transportable optical atomic clocks,” J. Mod. Opt. 65, 622 (2018).
  * Gellesch  _et al._ (2020) M. Gellesch, J. Jones, R. Barron, A. Singh, Q. Sun, K. Bongs, and Y. Singh, “Transportable optical atomic clocks for use in out-of-the-lab environments,” Adv. Opt. Techn. 9, 313 (2020).
  * Cao  _et al._ (2017) J. Cao, P. Zhang, J. Shang, K. Cui, J. Yuan, S. Chao, S. Wang, H. Shu, and X. Huang, “A compact, transportable single-ion optical clock with $7.8\times 10^{-17}$ systematic uncertainty,” Appl. Phys. B 123, 12 (2017).
  * Hannig  _et al._ (2019) S. Hannig, L. Pelzer, N. Scharnhorst, J. Kramer, M. Stepanova, Z. T. Xu, N. Spethmann, I. D. Leroux, T. E. Mehlstäubler, and P. O. Schmidt, “Towards a transportable aluminium ion quantum logic optical clock,” Rev. Sci. Instr. 90, 053204 (2019).
  * Burt  _et al._ (2021) E. A. Burt, J. D. Prestage, R. L. Tjoelker, D. G. Enzer, D. Kuang, D. W. Murphy, D. E. Robison, J. M. Seubert, R. T. Wang, and T. A. Ely, “Demonstration of a trapped-ion atomic clock in space,” Nature 595, 43 (2021).
  * Farkas  _et al._ (2010) D. M. Farkas, A. Zozulya, and D. Z. Anderson, “A compact microchip atomic clock based on all-optical interrogation of ultra-cold trapped Rb atoms,” Appl. Phys. B 101, 705 (2010).
  * Esnault  _et al._ (2013) F.-X. Esnault, E. Blanshan, E. N. Ivanov, R. E. Scholten, J. Kitching, and E. A. Donley, “Cold-atom double-$\Lambda$ coherent population trapping clock,” Phys. Rev. A 88, 042120 (2013).
  * Liu  _et al._ (2017) X. Liu, E. Ivanov, V. I. Yudin, J. Kitching, and E. A. Donley, “Low-drift coherent population trapping clock based on laser-cooled atoms and high-coherence excitation fields,” Phys. Rev. Applied 8, 054001 (2017).
  * Elvin  _et al._ (2019) R. Elvin, G. W. Hoth, M. Wright, B. Lewsi, J. P. McGilligan, A. S. Arnold, P. F. Griffin, and E. Riis, “Cold-atom clock based on a diffractive optic,” Opt. Expr. 27, 38359 (2019).
  * Perrella  _et al._ (2019) C. Perrella, P.S. Light, J.D. Anstie, F.N. Baynes, R.T. White, and A.N. Luiten, “Dichroic two-photon rubidium frequency standard,” Phys. Rev. Applied 12, 054063 (2019).
  * Maurice  _et al._ (2020) V. Maurice, Z. L. Newman, S. Dickerson, M. Rivers, J. Hsiao, P. Greene, M. Mescher, J. Kitching, M. T. Hummon, and C. Johnson, “Miniaturized optical frequency reference for next-generation portable optical clocks,” Opt. Expr. 28, 24708 (2020).
  * Gerginov and Beloy (2018) V. Gerginov and K. Beloy, “Two-photon optical frequency reference with active ac Stark shift cancellation,” Phys. Rev. Applied 10, 014031 (2018).
  * Martin  _et al._ (2018) K. W. Martin, G. Phelps, N. D. Lemke, M. S. Bigelow, B. Stuhl, M. Wojcik, M. Holt, I. Coddington, M. W. Bishop, and J. H. Burke, “Compact optical atomic clock based on a two-photon transition in Rubidium,” Phys. Rev. Applied 9, 014019 (2018).
  * Martin  _et al._ (2019) K. W. Martin, B. Stuhl, J. Eugenio, M. S. Safronova, G. Phelps, J. H. Burke, and N. D. Lemke, “Frequency shifts due to Stark effects on a rubidium two-photon transition,” Phys. Rev. A 100, 023417 (2019).
  * Newman  _et al._ (2019) Z. L. Newman, V. Maurice, T. Drake, J. R. Stone, T. C. Briles, D. T. Spencer, C. Fredrick, Q. Li, D. Westly, B. R. Ilic, B. Shen, M.-G. Suh, K. Y. Yang, C. Johnson, D. M. S. Johnson, L. Hollberg, K. J. Vahala, K. Srinivasan, S. A. Diddams, J. Kitching, S. B. Papp, and M. T. Hummon, “Architecture for the photonic integration of an optical atomic clock,” Optica 6, 680–685 (2019).
  * Carr (2014) A. W. Carr, “Improving quantum computation with neutral Cesium: readout and cooling on a quadrupole line, conditions for double magic traps and novel dissipative entanglement scheme,” PhD thesis Univ. Wisconsin-Madison (2014).
  * Huft  _et al._ (2022) P. Huft, Y. Song, K. Jooya, T. Graham, S. Deshpande, C. Fang, M. Kats, and M. Saffman, “A projected optical trap array,” to be published (2022).
  * Takamoto  _et al._ (2005) M. Takamoto, F.-L. Hong, R. Higashi, and H. Katori, “An optical lattice clock,” Nature 435, 321–324 (2005).
  * Biedermann  _et al._ (2013) G. W. Biedermann, K. Takase, X. Wu, L. Deslauriers, S. Roy, and M. A. Kasevich, “Zero-dead-time operation of interleaved atomic clocks,” Phys. Rev. Lett. 111, 170802 (2013).
  * Schioppo  _et al._ (2017) M. Schioppo, R. C. Brown, W. F. McGrew, N. Hinkley, R. J. Fasano, K. Beloy, T. H. Yoon, G. Milani, D. Nicolodi, J. A. Sherman, N. B. Phillips, C. W. Oates, and A. D. Ludlow, “Ultrastable optical clock with two cold-atom ensembles,” Nat. Photon. 11, 48–52 (2017).
  * Audoin  _et al._ (1991) C. Audoin, V. Candelier, and N. Dimarcq, “A limit to the frequency stability of passive frequency standards due to an intermodulation effect,” IEEE Tr. Instr. Meas. 40, 121 (1991).
  * Liu  _et al._ (2021) Y.-Y. Liu, Z. F, P. Xu, X.-D. He, J. Wang, and M.-S. Zhan, “Spectral filtering of dual lasers with a high-finesse length-tunable cavity for rubidium atom Rydberg excitation,” Chin. Phys. B 30, 074203 (2021).
  * Huft  _et al._ (2021) P. Huft, Y. Song, K. Jooya, T. Graham, S. Deshpande, C. Fang, M. Kats, and M. Saffman, “A simplified projected optical trap array,” arXiv:2112.xxxxx (to be submitted Dec. 2021) (2021).
  * (27) P. Barakhshan, A. Marrs, B. Arora, R. Eigenmann, and M. S. Safronova, Portal for High-Precision Atomic Data and Computation (version 1.0). University of Delaware, Newark, DE, USA. URL: https://www.udel.edu/atom [January 2021].
  * Sansonetti (2009) J. E. Sansonetti, “Wavelengths, transition probabilities, and energy levels for the spectra of Cesium CsI \- CsLV,” J. Phys. Chem. Ref. Data 38, 761 (2009).
  * Safronova  _et al._ (2016) M. S. Safronova, U. I. Safronova, and Charles W. Clark, “Magic wavelengths, matrix elements, polarizabilities, and lifetimes of Cs,” Phys. Rev. A 94, 012505 (2016).
  * Amini and Gould (2003) J. M. Amini and H. Gould, “High precision measurement of the static dipole polarizability of cesium,” Phys. Rev. Lett. 91, 153001 (2003).
  * Farley and Wing (1981) J. W. Farley and W. H. Wing, “Accurate calculation of dynamic Stark shifts and depopulation rates of Rydberg energy levels induced by blackbody radiation. Hydrogen, helium, and alkali-metal atoms,” Phys. Rev. A 23, 2397 (1981).
  * Heavner  _et al._ (2014) T. P. Heavner, E. A. Donley, F. Levi, G. Costanzo, T. E. Parker, J. H. Shirley, N. Ashb, S. Barlow, and S. R. Jefferts, “First accuracy evaluation of NIST-F2,” Metrologia 51, 174 (2014).
  * Sinclair  _et al._ (2015) L. C. Sinclair, J.-D. Deschênes, L. Sonderhouse, W. C. Swann, I. H. Khader, E. Baumann, N. R. Newbury, and I. Coddington, “Invited article: A compact optically coherent fiber frequency comb,” Rev. Sci. Instr. 86, 081301 (2015).
  * Loh  _et al._ (2020) W. Loh, J. Stuart, D. Reens, C. D. Bruzewicz, D. Braje, John Chiaverini, P. W. Juodawlkis, J. M. Sage, and R. McConnell, “Operation of an optical atomic clock with a Brillouin laser subsystem,” Nature 588, 244 (2020).
  * Ohmae  _et al._ (2021) N. Ohmae, M. Takamoto, Y. Takahashi, M. Kokubun, K. Araki, A. Hinton, I. Ushijima, T. Muramatsu, T. Furumiya, Y. Sakai, N. Moriya, N. Kamiya, K. Fujii, R. Muramatsu, T. Shiimado, and H. Katori, “Transportable strontium optical lattice clocks operated outside laboratory at the level of $10^{-18}$ uncertainty,” Adv. Quantum Technol. 4, 2100015 (2021).
  * Guo  _et al._ (2021) F. Guo, W. Tan, C. h. Zhou, J. Xia, Y. x. Chen, T. Liang, Q. Liu, Y. Liu, D. j. He, Y. z. Zhou, W. h. Wang, Y. Shen, H. x. Zou, and H. Chang, “A proof-of-concept model of compact and high-performance 87Sr optical lattice clock for space,” AIP Adv. 11, 125116 (2021).
  * Ebert  _et al._ (2015) M. Ebert, M. Kwon, T. G. Walker, and M. Saffman, “Coherence and Rydberg blockade of atomic ensemble qubits,” Phys. Rev. Lett. 115, 093601 (2015).
  * Zeiher  _et al._ (2015) J. Zeiher, P. Schauß, S. Hild, T. Macrí, I. Bloch, and C. Gross, “Microscopic characterization of scalable coherent Rydberg superatoms,” Phys. Rev. X 5, 031015 (2015).
  * Omran  _et al._ (2019) A. Omran, H. Levine, A. Keesling, G. Semeghini, T. T. Wang, S. Ebadi, H. Bernien, A. S. Zibrov, H. Pichler, S. Choi, J. Cui, M. Rossignolo, P. Rembold, S. Montangero, T. Calarco, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, “Generation and manipulation of Schrödinger cat states in Rydberg atom arrays,” Science 365, 570 (2019).
  * Graham  _et al._ (2021) T. M. Graham, Y. Song, J. Scott, C. Poole, L. Phuttitarn, K. Jooya, P. Eichler, X. Jiang, A. Marra, B. Grinkemeyer, M. Kwon, M. Ebert, J. Cherek, M. T. Lichtman, M. Gillette, J. Gilbert, D. Bowman, T. Ballance, C. Campbell, E. D. Dahl, O. Crawford, N. S. Blunt, B. Rogers, T. Noel, and M. Saffman, “Demonstration of multi-qubit entanglement and algorithms on a programmable neutral atom quantum computer,” arXiv:2112.14589 (2021).


#  Density functionals with spin-density accuracy for open shells

Timothy Callow1,2,3,4 Benjamin Pearce1 Nikitas Gidopoulos1
1Department of Physics Durham University  South Road  Durham  DH1 3LE  UK
2 Max-Planck-Institut für Mikrostrukturphysik  Weinberg 2  D-06120 Halle
Germany
3 Center for Advanced Systems Understanding (CASUS)  D-02826 Görlitz  Germany
4 Helmholtz-Zentrum Dresden-Rossendorf  D-01328 Dresden  Germany

###### Abstract

Electrons in zero external magnetic field can be studied with density
functional theory (DFT) or with spin-DFT (SDFT). The latter is normally used
for open shell systems because its approximations appear to model better the
exchange and correlation (xc) functional, but also because so far the
application of DFT implied a closed-shell-like approximation. Correcting this
error for open shells allows the approximate DFT xc functionals to become as
accurate as those in SDFT. In the limit of zero magnetic field, the Kohn-Sham
equations of SDFT emerge as the generalised KS equations of DFT.

Keywords— local density approximation, implicit density functionals, semi-
local density approximation

##  1 Introduction

Density functional theory (DFT) [12, 13] and spin-density functional theory
(SDFT) [22, 9] are two reformulations of the many-electron problem. Their
computational advantage is replacing the solution of Schrödinger’s equation
for the multi-dimensional ground state of the physical electronic system of
interest by the solution of a much less demanding equation that yields the
ground-state (total) density $\rho$ (DFT) or the spin-density
$(\rho^{\uparrow},\rho^{\downarrow})$ (SDFT) of the real system. For
simplicity, we discuss the special case of SDFT for a collinear magnetisation
density $m=-\mu_{0}(\rho^{\uparrow}-\rho^{\downarrow})$; the generalisation to
the noncollinear case is straightforward.

Both theories can be applied equally well to study electronic systems when the
external magnetic field, $B$, vanishes. Then, obviously, the predicted values
for any observable quantity must be the same between the two, formally exact,
theories. Specifically, the ground-state (total) density and energy of the
electronic system must be the same in both theories.

In DFT/SDFT the value of any observable quantity becomes a functional of the
density/spin-density, but so far little is known about the functional
dependence of general observable quantities, other than the energy and
identically the total density (DFT) or the spin-density (SDFT). For example,
the spin-density of the exact SDFT-KS system is equal to the true spin-density
of the real system. However, this does not hold for the spin-density of the
exact DFT-KS system when $B=0$. In order to obtain the real (observable) spin-
density from the density, we need the so far unknown density functional for
the spin-density.

In the absence of a magnetic field, $B=0$, by definition, the virtual KS
systems in DFT and SDFT of non-interacting electrons share with the real
(interacting) system the same total density (DFT and SDFT) and the same spin-
density (SDFT). For open shell systems, although the two KS systems have equal
total densities, they are not the same and the single-particle KS equations of
exact SDFT do not reduce to those of exact DFT.

On reflection it seems paradoxical that, for $B=0$ the virtual, non-
interacting KS system of SDFT does not reduce to the virtual, non-interacting
KS system of DFT, and that the two have different spin-densities when both
share the same total density. The apparent paradox is related to the limit of
the SDFT single-particle equations when $B=0$ and is addressed in the second
part of this paper.

Even though for $B=0$ both exact theories can equally well describe an
electronic system, this changes when density functional approximations (DFAs)
are introduced. The flexibility of two basic variables in SDFT is thought to
offer a better modelling of the exchange-correlation energy, $E_{\textrm{\rm
xc}}[\rho^{\uparrow},\rho^{\downarrow}]$, compared with just a single variable
in DFT, $E_{\textrm{\rm xc}}[\rho]$. As Parr and Yang write [15], in local and
semi-local density functional approximations “the exchange-correlation energy
of the electrons is approximated locally by results for the homogeneous spin-
compensated electron gas. Such a procedure is not appropriate for systems with
unpaired electrons, like open-shell molecules. A better description for such
systems will be obtained through the use of the exchange-correlation energy of
the homogeneous spin-polarised electron gas”. As a result, open shell systems
are in practice usually treated in SDFT rather than DFT.

For $B=0$, the main error in modelling open shell electronic systems as
unpolarised arises in the exchange energy functional, $E_{\textrm{x}}[\rho]$
(DFT) vs $E_{\textrm{x}}[\rho^{\uparrow},\rho^{\downarrow}]$ (SDFT). When a
local or semi-local approximate exchange energy expression is employed,
depending explicitly only on the density $\rho$ and its derivatives, it is
essentially assumed that the spin-up and spin-down densities are equal to each
other and to one half of the total density. For open shell systems, this
amounts to mixing partly the spin-up $\rho^{\uparrow}$ with the spin-down
$\rho^{\downarrow}$ densities in the exchange energy expression:

$E_{\textrm{x}}[\rho]\simeq E_{\textrm{x}}[\rho/2,\rho/2],\ \
\rho=\rho^{\uparrow}+\rho^{\downarrow}.$ (1)

This mixing of the spin-densities in the exchange energy leads to a spurious
error that we call the ghost exchange energy error, in analogy to the ‘ghost
interaction’ error of Ref. [8]. We define it by the difference of the exchange
energies with mixed spin-densities from the reference exchange energy where
the spin-densities are separate, as they must be in the KS Slater determinant:

$G_{\rm x}=E_{\rm x}[\rho/2,\rho/2]-E_{\rm
x}[\rho^{\uparrow},\rho^{\downarrow}].$ (2)

In DFT, the ghost exchange error is not caused by the approximate or exact
expression for the exchange energy functional but by the mixing of two spin-
densities when an open shell system is treated as if it were closed-shell,
considering half its electrons as spin-up and half as spin-down. It follows
that the ghost exchange error would even be present for the exact exchange
energy density functional, were the relation (1) to be (incorrectly) assumed.

In the following first part of this paper, we show that open shell systems in
zero magnetic field can already be modelled within DFT, without any mixing of
spin-up and spin-down KS orbitals and avoiding the ghost-exchange error. In
the literature there are alternative approaches for open-shell systems based
on a spin-restricted formalism, for example we mention the restricted open-
shell KS (ROKS) method [6, 7]. This has no ghost exchange error, but uses a
linear combination of KS determinants to satisfy symmetry constraints in a
similar spirit to the ROHF method [18]. Our approach, on the other hand,
requires only a single Slater determinant and is thus formally defined in KS-
DFT.

##  2 The xc energy as an implicit density functional

To continue, we consider the exact exchange energy in KS theory. It is given
by the Fock expression in terms of the spin-orbitals of the KS Slater
determinant,

$E_{\textrm{x}}[\rho^{\uparrow},\rho^{\downarrow}]=-\frac{1}{2}\iint\\!\\!\differential{\mathbf{r}}\differential{\mathbf{x}}\left\\{\frac{|\rho^{\uparrow}(\mathbf{r},\mathbf{x})|^{2}}{|\mathbf{r}-\mathbf{x}|}+\frac{|\rho^{\downarrow}(\mathbf{r},\mathbf{x})|^{2}}{|\mathbf{r}-\mathbf{x}|}\right\\}.$
(3)

where $\rho^{\sigma}(\mathbf{r},\mathbf{x})$, $\sigma=\uparrow,\downarrow$, is
the spin-$\sigma$, one-body reduced density matrix of KS orbitals.

The exact exchange energy functional separates in two disjoint terms,

$E_{\textrm{x}}[\rho^{\uparrow},\rho^{\downarrow}]=E_{\textrm{x}}[\rho^{\uparrow},0]+E_{\textrm{x}}[0,\rho^{\downarrow}],$
(4)

in which the subsets of spin-up and spin-down KS orbitals do not mix. Hence,
in order to avoid cross exchange effects between opposite spin-electrons also
in approximations, the approximate exchange energy density functionals must
also satisfy equality (4).

Since modelling the approximate exchange energy in terms of the total density
and its derivatives, as in LDA and in semi-local DFAs, always violates
equality (4), we have to model the exchange energy using the KS spin-density
111In DFT the KS spin-density need not coincide with the true spin-density..
This modelling is still within DFT (not SDFT) since in DFT the KS spin-density
(the spin-density of the KS determinant) is an implicit functional of the
total density.

Equality (4) is satisfied by the local spin-density approximation (LSDA), and
as far as we know all spin-dependent density functionals (DFs) for the
exchange energy, approximate or the exact one. Therefore, we shall use
expressions for the exchange energy which depend on the total density
indirectly, or implicitly, via the KS spin-density (spin-density functionals,
SDFs),

$E_{\textrm{x}}^{\textrm{iDF}}[\rho]=E_{\textrm{x}}^{\textrm{SDF}}\big{[}\rho^{\uparrow}[\rho],\rho^{\downarrow}[\rho]\big{]},$
(5)

where $E_{\rm x}^{\rm SDF}$ must be the sum of two disjoint terms (4).

The acronyms iDF and SDF denote that the exchange energy on the lhs is an
_implicit_ density-functional and that the exchange energy expression on the
rhs depends on the KS spin-density.

Given that the exchange-functional is now an implicit functional of the
density, and that the correlation energy is modelled better in terms of the
spin-density, we choose to write also the correlation energy functional as an
implicit functional, depending on the total density via the KS spin-density,
like exchange 222The correlation energy is not a sum of two disjoint terms and
does not satisfy (4). Indeed in approximations there is a small gain in
accuracy in treating correlation using spin-polarised DFs, though we have
observed this is typically a minor contribution compared to the ghost exchange
error. The whole exchange-correlation (xc) density functional is thus written
as

$E_{\textrm{\rm xc}}^{\textrm{iDF}}[\rho]=E_{\textrm{\rm
xc}}^{\textrm{SDF}}\big{[}\rho^{\uparrow}[\rho],\rho^{\downarrow}[\rho]\big{]}.$
(6)

Hence, the total energy DF is now given by (for the exact or approximate DFs)

$E_{v_{\rm
en}}^{\textrm{iDF}}[\rho]=T_{\textrm{s}}[\rho]+\int\differential{\mathbf{r}}v_{\textrm{en}}(\mathbf{r})\rho(\mathbf{r})+U[\rho]+E_{\textrm{\rm
xc}}^{\textrm{iDF}}[\rho],$ (7)

where $U[\rho]$ is the Hartree energy and $v_{\rm en}(\mathbf{r})$ is the
external (electron-nuclear) potential of the interacting system.

As the xc energy functional is an implicit DF, the xc-potential, given by the
functional derivative

$v_{\textrm{\rm xc}}[\rho](\mathbf{r})=\functionalderivative{E_{\textrm{\rm
xc}}^{\textrm{iDF}}[\rho]}{\rho(\mathbf{r})},$ (8)

must be determined using the optimised effective potential (OEP) method [20,
21],

$\int{\differential{\mathbf{r}^{\prime}}}\sum_{\sigma}\chi^{\sigma}(\mathbf{r},\mathbf{r}^{\prime})v_{\textrm{\rm
xc}}[\rho](\mathbf{r}^{\prime})=\\\
\int{\differential{\mathbf{r}^{\prime}}}\sum_{\sigma}\chi^{\sigma}(\mathbf{r},\mathbf{r}^{\prime})v_{\rm
xc}^{\sigma}\big{[}\rho^{\uparrow}[\rho],\rho^{\downarrow}[\rho]\big{]}(\mathbf{r}^{\prime}),$
(9)

with

$\displaystyle v_{\rm
xc}^{\sigma}\big{[}\rho^{\uparrow}[\rho],\rho^{\downarrow}[\rho]\big{]}(\mathbf{r}^{\prime})=\frac{\partial
E_{\rm xc}^{\rm
SDF}[\rho^{\uparrow},\rho^{\downarrow}]}{\partial\rho^{\sigma}({\bf
r}^{\prime})}\Bigg{|}_{\rho^{\uparrow}=\rho^{\uparrow}[\rho]\atop\rho^{\downarrow}=\rho^{\downarrow}[\rho]},$
(10)
$\displaystyle\chi^{\sigma}(\mathbf{r},\mathbf{r}^{\prime})=-2\sum_{i=1}^{N^{\sigma}}\sum_{a=N^{\sigma}+1}^{\infty}\frac{\phi_{i}(\mathbf{r})\phi_{a}(\mathbf{r})\phi_{i}(\mathbf{r}^{\prime})\phi_{a}(\mathbf{r}^{\prime})}{\epsilon_{a}-\epsilon_{i}}.$
(11)

This equation is our first important result. The KS-DFT xc potential $v_{\rm
xc}$ is given as the weighted sum of the spin-dependent xc potentials $v_{\rm
xc}^{\sigma}$, with weighting factors the spin-dependent response functions
$\chi^{\sigma}$. We argue that this is the only consistent way (that avoids
the ghost-exchange error) to evaluate the KS-DFT xc potential for open shell
systems. For closed shell systems, the solution reduces to the familiar KS xc
functional derivative. For fully spin-polarised systems (eg spin up,
$\rho^{\uparrow}=\rho,\rho^{\downarrow}=0$), the spin-down response function
vanishes, and the solution of (9) reduces to the spin-up xc potential, $v_{\rm
xc}[\rho]=v_{\rm xc}^{\uparrow}[\rho,0]$.

##  3 Results

We have implemented these equations in the Gaussian basis set code HIPPO
333For information, contact NL at<EMAIL_ADDRESS>For more detail on the
computational implementation of the OEP equation please see Refs. [1] and [2].
Unless otherwise stated, all results use cc-pVTZ orbital bases and
uncontracted cc-pVDZ auxiliary bases [5, 23].

Ground-state energies calculated with our iDF KS method compare favourably
with those from SDFT-KS. This is demonstrated in Table 1 for some atoms and
molecules at their equilibrium geometries, using the L(S)DA functional. We
emphasise again that we seek only the ground-state energy in our approach, and
thus where a state is referred to as ‘doublet’ this simply means

$\displaystyle\Delta N=N^{\uparrow}-N^{\downarrow}=1,$ (12)

and likewise for triplet and so on. We also note the effect of the ghost
exchange (and correlation) error, which is particularly pronounced for triplet
states; we shall now explore some examples in which this error significantly
affects results.

| $E_{\textrm{LDA}}\ (\textrm{Har})$ | $E_{\textrm{iLDA}}\ (\textrm{Har})$ | $E_{\textrm{LSDA}}\ (\textrm{Har})$
---|---|---|---
$\ch{Li}$ | -7.388721 | -7.398145 | -7.398155
$\ch{B}$ | -24.43315 | -24.44669 | -24.44747
$\ch{N}$ | -54.12891 | -54.14996 | -54.15110
$\ch{Na}$ | -161.6491 | -161.6571 | -161.6572
$\ch{Si}$111Triplet state | -288.4640 | -288.4905 | -288.4910
$\ch{LiH^{+}}$ | -7.652062 | -7.685603 | -7.685608
$\ch{O2}$111Triplet state | -149.6038 | -149.6383 | -149.6403
$\ch{OH}$ | -75.34813 | -75.37077 | -75.37208
$\ch{NH4}$ | -56.79800 | -56.80389 | -56.80404
Avg % diff | 0.0827 | 0.00101 | -
Table 1: L(S)DA ground-state total energies calculated with: (i) standard LDA;
(ii) implicit LDA (iLDA); and (iii) spin-LDA (LSDA). All states are doublets
unless specified otherwise; auxiliary bases are uncontracted cc-pVTZ. Figure
1: Energy dissociation curves for the H2 molecule. Top, iLDA energies for
different values of $\Delta N$; bottom, comparison of LDA, spin LDA (LSDA),
and iLDA minimum energies. Figure 2: Energy dissociation curves for the OH
radical. Top, implicit PBE (iPBE) energies for different values of $\Delta N$;
bottom, comparison of PBE, spin PBE (SPBE), and iPBE minimum energies.

In Fig 1, we have calculated the energy dissociation curve for the H2 molecule
with the L(S)DA functional. As is well-known in the literature [3] and can be
seen in this figure, the standard restricted solution yields a qualitatively
incorrect dissociation curve; however, in our new method, once the bond is
stretched enough the triplet state becomes lower in energy than the singlet,
and the energy tends to the correct limit. This transition does not occur when
the ghost exchange energy is present: in this case the triplet energy is
higher than the singlet at all bond distances.

We see a similar picture emerge for the stretched OH radical with the Perdew-
Burke-Ernzehof (PBE) [16] functional in Fig. 2. In this case, the energy
dissociation curve is again qualitatively inaccurate for the solution
contaminated with ghost-exchange errors; in our method which removes the
ghost-exchange error, the energy correctly becomes equal to the energy of the
two atoms in the infinitely stretched limit. In both these examples, the
transition region is interesting: the unrestricted solution yields a smooth
dissociation curve, whereas in our single-determinant method the transition is
abrupt. Of course, our method does not correct the ubiquitous inability of
typical semi-local functionals to capture static correlation effects [4].

##  4 Limit of SDFT-KS equations
for $B=0$.

We now turn our attention to the limit of SDFT when $B=0$ and compare with
DFT. We start with the universal internal energy density functional,

$F[\rho]=\min_{\Psi\rightarrow\rho}\langle\Psi|{\hat{T}}+{\hat{V}}_{\rm
ee}|\Psi\rangle\,.$ (13)

The minimising state $\Psi_{\rho}$ depends on the total density $\rho$. When
$\rho$ is the ground state density of an interacting electronic system bound
by a local potential $v$ (i.e. when $\rho$ is interacting $v$-representable),
$\Psi_{\rho}$ is the corresponding ground state.

To proceed, we follow Levy [14] and separate the minimisation in (13) into two
separate minimisations with the same minimum,

$F[\rho]=\min_{(\rho^{\uparrow},\rho^{\downarrow})\rightarrow\rho}\bigg{[}\min_{\Psi\rightarrow(\rho^{\uparrow},\rho^{\downarrow})}\langle\Psi|{\hat{T}}+{\hat{V}}_{\rm
ee}|\Psi\rangle\bigg{]}\,.$ (14)

The inner minimisation, performed first, is over all states $\Psi$ with spin-
density $(\rho^{\uparrow},\rho^{\downarrow})$. The outer minimisation is over
all spin-densities that add up to the total density $\rho$. The first
minimisation inside the brackets on the rhs defines SDFT’s internal energy
functional; we conclude that the DFT internal energy functional can be
obtained from the SDFT functional with an extra optimisation over spin-
densities with the same $\rho$ (see also Appendix B in Ref. [17]),

$F[\rho]=\min_{(\rho^{\uparrow},\rho^{\downarrow})\rightarrow\rho}F[\rho^{\uparrow},\rho^{\downarrow}]\,.$
(15)

The two minimisations in (14) have the same overall minimum as the
minimisation in (13). Therefore, the minimising spin-density,
$(\rho^{\uparrow}_{\rho},\rho^{\downarrow}_{\rho})$, is equal to the spin-
density of $\Psi_{\rho}$ in (13) and hence equal to the spin-density of the
real interacting system, even though the minimisation (15) is at fixed total
density (not fixed spin-density), i.e. within DFT, not SDFT.

The minimisation in (15) is worth investigating further. We invoke SDFT’s KS
system with spin-density $(\rho^{\uparrow},\rho^{\downarrow})$, to expand the
internal energy functional in the usual way,
$F[\rho^{\uparrow},\rho^{\downarrow}]=T_{\rm
s}[\rho^{\uparrow},\rho^{\downarrow}]+E_{\rm
xc}[\rho^{\uparrow},\rho^{\downarrow}]+U[\rho^{\uparrow}+\rho^{\downarrow}]$.
Dropping the superscript SDF from the functionals of SDFT, we obtain

$F[\rho]=\min_{(\rho^{\uparrow},\rho^{\downarrow})\rightarrow\rho}\Big{\\{}T_{\rm
s}[\rho^{\uparrow},\rho^{\downarrow}]+E_{\rm
xc}[\rho^{\uparrow},\rho^{\downarrow}]\Big{\\}}+U[\rho]\,.$ (16)

Since every SDFT-KS system has a different spin-density, the minimisation
effectively searches over all SDFT-KS systems with common total density $\rho$
and returns that with the correct spin-density
$(\rho^{\uparrow}_{\rho},\rho^{\downarrow}_{\rho})$. The KS Slater determinant
state is $\Phi_{\rho^{\uparrow}_{\rho},\rho^{\downarrow}_{\rho}}$. The
minimising SDFT-KS system depends only on the total density $\rho$. The
corresponding SDFT-KS potential (functional of the spin-density) at the
specific spin-density $(\rho^{\uparrow}_{\rho},\rho^{\downarrow}_{\rho})$ is
also an implicit functional of $\rho$. In the Supplementary Material we show
how this unrestricted KS potential emerges directly from the minimisation in
(16).

The SDFT-KS system that minimises (16) defines a new non-interacting system in
DFT (when $B=0$) represented by the state
$\Phi_{\rho^{\uparrow}_{\rho},\rho^{\downarrow}_{\rho}}$. We call it the
generalised KS (GKS) system, as its derivation is analogous to the well-known
GKS scheme in the literature [19].

We conclude that the SDFT-KS equations reduce to the DFT-GKS equations for
$B=0$. The DFT-GKS Slater determinant state
$\Phi_{\rho^{\uparrow}_{\rho},\rho^{\downarrow}_{\rho}}$ gives not only the
true total density but also the true spin-density (or magnetisation density).
The elusive exact density functional for the spin-density is the spin-density
of the DFT-GKS system.

##  5 The spin-symmetry dilemma

One of the oft-cited issues with SDFT in the absence of an external magnetic
field is the spin-symmetry dilemma, so-called because SDFT approximate results
yield either accurate total energies but with a poor (broken symmetry)
prediction for the spin-density, or an accurate prediction for the spin-
density with poor total energies [17, 11, 10]. The dilemma lies in that the
spin-density is the key quantity and hence SDFT (exact) results are supposed
to yield both the exact total energy and the exact spin-density. For example,
$\ch{H}_{2}$ should dissociate into two Hydrogen atoms with zero net
magnetisation (difference of spin densities): indeed the lowest energy
solution of local and semi-local approximations in SDFT gets the energy right,
but it wrongly dissociates the molecule into two atoms of opposite spin
(broken symmetry solution). There is also a self-consistent solution where the
spin-densities are correctly equal to each other but the total energy is too
high.

Perdew, Savin and Burke [17] do not consider that the spin-symmetry dilemma
results from the approximation and to address the issue they reinterpret exact
SDFT for $B=0$. They write that the KS spin-densities “are not physical spin
densities, but are only intermediate objects (like the Kohn-Sham orbitals…”.
Instead, they argue that the on-top electron pair density
$P(\mathbf{r},\mathbf{r})$, alongside the total density, are the two
fundamental variables of the theory. Since the KS spin-density is no longer a
basic variable of SDFT [17], its value is not a prediction for the real spin-
density of the interacting system. We comment that considering the SDFT-KS
spin-density to give the observable spin-density for any value of $B$ except
at $B=0$, where it needs re-interpretation according to Ref. [17], suggests
that SDFT may have a singular limit at $B=0$.

We share the principle of Ref. [17], namely that the spin-density is no longer
the basic quantity in the absence of a magnetic field; but in our approach,
there is no need to consider the auxiliary on-top pair density, with just the
total density being the fundamental variable: our (ghost-exchange error free)
theory yields the correct ground-state energy and the correct total density in
the dissociation limit, and thus the spin-symmetry dilemma is weakened as a
fundamental problem, since the KS spin-density is not expected equal to the
observable spin-density.

However, we also propose that the DFT-GKS spin-density is the exact density
functional of the observable spin-density; knowledge of the (exact) SDFT xc
functional $E_{\rm xc}[\rho^{\uparrow},\rho^{\downarrow}]$ is necessary to
obtain from the density the (exact) DFT-GKS spin-density. The value of the
latter is the same as the SDFT-KS spin-density 444Hence SDFT does not have a
singular limit at $B=0$. which, in local and semi-local approximations,
demonstrates the spin-symmetry dilemma. We believe the explanation is not a
singular limit of SDFT when $B=0$ but that the approximate DFT-GKS spin-
density is not a sufficiently accurate density functional for the real spin-
density. According to Cohen, Sanchez and Yang [4, 3] on the fractional spin
error, the crux of the problem is the degeneracy of different Slater
determinants with the same energy and total density but different spin-
densities. This kind of strong correlation cannot be captured with local or
semi-local functionals.

##  6 Conclusion

We have addressed two related conceptual or fundamental theoretical questions
in DFT and spin-DFT. The first is the challenge how to apply the KS equations
of DFT for open shells, avoiding a serious qualitative error, which we called
the "ghost-exchange error". Associated with this problem is the widespread
belief that spin-DFT approximations are inherently more flexible and hence can
model more accurately the xc energy for open shells than the corresponding
approximations of DFT. In the literature this common view seemed plausible,
almost obvious. We demonstrate that actually density functionals can be as
accurate as spin-density functionals, provided the ghost-exchange error is
corrected.

For zero external magnetic field (B=0), both DFT and spin-DFT can obviously be
applied to study open shell electronic systems and both exact theories should
give identical results for any observable quantity. Then, intuitively, one
would expect in the limit $B=0$ the single particle KS equations of spin-DFT
to reduce to a set of single-particle equations of DFT. We demonstrate this
intuition is correct and that the SDFT equations reduce to the generalised KS
equations of DFT.

## Acknowledgments

We thank N. Lathiotakis for help with the code HIPPO. We acknowledge
enlightening discussions with K. Burke, E.K.U. Gross and N. Lathiotakis.
N.I.G. acknowledges financial support from The Leverhulme Trust, through a
Research Project Grant with number RPG-2016-005. B.P. acknowledges financial
support from EPSRC with grant number EP/R513039/1.
During the final stage of the project T.C. was funded by the Center for
Advanced Systems Understanding (CASUS) which is financed by Germany’s Federal
Ministry of Education and Research (BMBF) and by the Saxon Ministry for
Science, Culture and Tourism (SMWK) with tax funds on the basis of the budget
approved by the Saxon State Parliament.

## Data Availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * [1] Timothy J Callow.  Systematic routes to improved approximations in kohn-sham theory.  Available at Durham E-Theses Online: http://etheses.dur.ac.uk/13817/, 2020.
  * [2] Timothy J. Callow, Benjamin J. Pearce, Tom Pitts, Nektarios N. Lathiotakis, Matthew J. P. Hodgson, and Nikitas I. Gidopoulos.  Improving the exchange and correlation potential in density-functional approximations through constraints.  Faraday Discuss., 224:126–144, 2020.
  * [3] Aron J. Cohen, Paula Mori-Sánchez, and Weitao Yang.  Insights into current limitations of density functional theory.  Science, 321(5890):792–794, 2008.
  * [4] Aron J. Cohen, Paula Mori-Sánchez, and Weitao Yang.  Fractional spins and static correlation error in density functional theory.  J. Chem. Phys., 129(12):121104, 2008.
  * [5] Thom H. Dunning.  Gaussian basis sets for use in correlated molecular calculations. i. the atoms boron through neon and hydrogen.  J. Chem. Phys., 90(2):1007–1023, 1989.
  * [6] Michael Filatov and Sason Shaik.  Spin-restricted density functional approach to the open-shell problem.  Chemical Physics Letters, 288(5):689 – 697, 1998.
  * [7] Michael Filatov and Sason Shaik.  Application of spin-restricted open-shell kohn–sham method to atomic and molecular multiplet states.  J. Chem. Phys., 110(1):116–125, 1999.
  * [8] N. I. Gidopoulos, P. G. Papaconstantinou, and E. K. U. Gross.  Spurious interactions, and their correction, in the ensemble-kohn-sham scheme for excited states.  Phys. Rev. Lett., 88:033003, Jan 2002.
  * [9] Nikitas I Gidopoulos.  Potential in spin-density-functional theory of noncollinear magnetism determined by the many-electron ground state.  Phys. Rev. B, 75(13):134408, 2007.
  * [10] Andreas Görling.  Symmetry in density-functional theory.  Phys. Rev. A, 47:2783–2799, Apr 1993.
  * [11] O. Gunnarsson and B. I. Lundqvist.  Exchange and correlation in atoms, molecules, and solids by the spin-density-functional formalism.  Phys. Rev. B, 13:4274–4298, May 1976.
  * [12] P. Hohenberg and W. Kohn.  Inhomogeneous electron gas.  Phys. Rev., 136:B864–B871, Nov 1964.
  * [13] W. Kohn and L. J. Sham.  Self-consistent equations including exchange and correlation effects.  Phys. Rev., 140:A1133–A1138, Nov 1965.
  * [14] Mel Levy.  Universal variational functionals of electron densities, first-order density matrices, and natural spin-orbitals and solution of the v-representability problem.  Proceedings of the National Academy of Sciences, 76(12):6062–6065, 1979.
  * [15] Robert G Parr and Weitao Yang.  Density Functional Theory of Atoms and Molecules, chapter 8, pages 173–174.  Oxford University Press, 1989.
  * [16] John P. Perdew, Kieron Burke, and Matthias Ernzerhof.  Generalized gradient approximation made simple.  Phys. Rev. Lett., 77:3865–3868, Oct 1996.
  * [17] John P. Perdew, Andreas Savin, and Kieron Burke.  Escaping the symmetry dilemma through a pair-density interpretation of spin-density functional theory.  Phys. Rev. A, 51:4531–4541, Jun 1995.
  * [18] C. C. J. Roothaan.  Self-consistent field theory for open shells of electronic systems.  Rev. Mod. Phys., 32:179–185, Apr 1960.
  * [19] A. Seidl, A. Görling, P. Vogl, J. A. Majewski, and M. Levy.  Generalized kohn-sham schemes and the band-gap problem.  Phys. Rev. B, 53:3764–3774, Feb 1996.
  * [20] R. T. Sharp and G. K. Horton.  A variational approach to the unipotential many-electron problem.  Phys. Rev., 90:317–317, Apr 1953.
  * [21] James D. Talman and William F. Shadwick.  Optimized effective atomic central potential.  Phys. Rev. A, 14:36–40, Jul 1976.
  * [22] Ulf von Barth and Lars Hedin.  A local exchange-correlation potential for the spin polarized case. i.  Journal of Physics C: Solid State Physics, 5(13):1629, 1972.
  * [23] David E. Woon and Thom H. Dunning.  Gaussian basis sets for use in correlated molecular calculations. iii. the atoms aluminum through argon.  J. Chem. Phys., 98(2):1358–1371, 1993.