††thanks<EMAIL_ADDRESS>
# Remote Cooling of Spin-ensembles Coupled to a Spin-mechanical Hybrid
Interface

Yang Wang  3\. Physikalisches Institut, ZAQuant, University of Stuttgart,
Allmandring 13, 70569 Stuttgart, Germany  Durga Dasari  3\. Physikalisches
Institut, ZAQuant, University of Stuttgart, Allmandring 13, 70569 Stuttgart,
Germany

###### Abstract

We present here a protocol for the ground-state cooling of a hybrid quantum
system consisting of a spin ensemble, an oscillator, and a single probe spin.
In the weak dispersive coupling limit between the spins and the oscillator, a
back-action effect from frequent measurements of the probe spin allows
exponential cooling of the hybrid system. We identify the parameter regimes
necessary to cool the ensemble, the oscillator, or both to their thermal
ground states. Our results contribute to the development of a practical
solution for cooling/polarizing large spin ensembles. More importantly, we
show that even with weak coupling, significant changes in the dynamics of such
ensemble-oscillator hybrid systems can be achieved by manipulating a single
spin. This result provides a relatively simple means of tuning the dynamics of
a hybrid system, which might facilitate advances in many quantum technology
applications such as preparing entangled states in macroscopic objects.

††preprint: APS/123-QED

##  I Introduction

Quantum technology holds great potential to revolutionize areas such as
information processing [1, 2, 3], secure communication [4, 5], sensing [6, 7],
and fundamental testing [8]. However, translating this potential into real-
world applications remains a significant challenge due to the inherent
limitations of current experimental platforms [9, 10]. For example, nitrogen-
vacancy (NV) centers in diamond [11] have an extended coherence time [12, 13]
and an optical interface that enables long-range entanglement [14, 15], but
their limited control accuracy [16] hinders advanced applications such as
fault-tolerant quantum computation [17, 18]. Therefore, it is essential to
explore the hybrid system approach that fully exploits the unique capabilities
of different systems [19, 20, 21, 22]. Following this idea, large NV ensembles
have been used exclusively as memory for superconducting qubits [23, 24],
known for their superior control fidelity [1], reducing the need for high-
precision entangling gates over NV centers required for quantum computation.
The rich dynamics of hybrid quantum systems also opens up significant
prospects for sensing [25, 22], simulation [21, 26] and machine learning [27,
28, 29].

Establishing high-fidelity entanglement between remote spin qubits through a
bus is a well-studied paradigm [30, 22]. The idea is to transduce their common
coupling to macroscopic objects (bus) into indirect interactions that lead to
their entanglement [30, 31]. Nano-mechanical oscillators have been playing a
pivotal role in creating universal quantum transducers for spin–spin spin-
photon and phonon-photon interactions [32, 33]. Other transducers in the form
of magnons [34], superconducting qubits [35, 36], optical cavities [25],
optomechanical cavities [37], and superconducting resonator [38] have also
been considered. There is also a proposal of using a three-level spin as the
transducer to mediate entanglement between two oscillators [39].

Transducers have primarily been proposed to create indirect interactions
between the initially non-interacting systems, effectively implementing
entanglement gates. between remote qubits [37, 33, 34]. Nevertheless, there
are few transducer protocols that focus on heat extraction from quantum
systems to achieve ground state cooling or quantum state preparation. This
aspect is crucial, especially for large spin ensembles, since conventional
ensemble cooling or polarization techniques often rely on resonant energy
exchange methods, such as the Hartman-Hahn transfer [40, 41].

Polarization based on energy transfer typically exploit the interaction
between a single spin and all the ensemble spins. This single spin can be
repeatedly polarized, thus realizing a continuous removal of entropy from the
ensemble, leaving it finally polarized. However, when dealing with identical
spins (with similar g-factors), the risk of the ensemble being trapped in dark
states due to Hamiltonian symmetries is significantly high [41, 42], which
prevents full polarization. This dark state problem can be overcome by
performing projective measurements on a single qubit interacting with the
target ensemble and using its back action to extract heat. Such quantum state
projections of a single spin not only allow efficient heat extraction from an
ensemble, but can also work in a non-trivial situation where the original
system dynamics does not provide any energy exchange [43, 44].

In this work, we will consider a tripartite hybrid system consisting of a spin
ensemble, an oscillator, and a single spin (probe). Since we envision the
oscillator to be a macroscopic object such as a nano-mechanical oscillator,
the probe spin and the ensemble can be remotely entangled, as schematically
shown in Fig. 1. In the limit where the oscillator is dispersively coupled
(off-resonant) to the probe and ensemble spins, no energy exchange is possible
from the natural Hamiltonian dynamics. However, by exploiting the measurement
back action from the frequently measured probe spin, it is possible to
cool/polarize the spin ensemble from its initial fully mixed thermal states to
the ground states.

Our cooling protocol is based on indirect coupling between the probe spin and
the ensemble that is mediated by the oscillator. It has already been
demonstrated that multiple spins can be coupled to a nano-mechanical
oscillator [43, 44], which would be much easier than achieving strong coupling
between a single spin and all target spins in a large ensemble directly [24].
Moreover, since only weak dispersive coupling between the spins and the
oscillator is required, our protocol can be straightforwardly implemented on a
large variety of physical systems and is anticipated to be more feasible for
cooling down larger electron-spin ensembles. Furthermore, given the wide
operating temperature range of contemporary hybrid systems, from $4$K to
$300$K, the need for cooling down individual components is paramount [45].

The paper is organized as follows: In Sec. II, we outline the Hamiltonian of
the system and discuss possible experimental implementations of our proposal.
In Sec. III, we analyze the dynamics of the hybrid system and describe how
periodic dynamical decoupling (DD) pulses on the probe spin can modify the
dynamics of other parts in the system. In Sec. V, we study the parameters that
allow cooling of the spin ensemble, the oscillator, or both. Finally, Sec. VI
provides a summary and discusses possible directions for future research.

##  II SETUP

This work focuses on a tripartite quantum system consisting of an ancilla spin
(probe), a spin ensemble, and an oscillator. This hybrid system can be
realized as shown in Fig. 1 through a one-sided clamped cantilever embedding a
single NV center positioned such that it experiences maximal strain coupling.
To the edge of the cantilever is attached a magnetic tip that produces a field
gradient $\nabla B$. This field is then coupled to a spin ensemble on an
external substrate. The spin ensemble and the probe spin are both coupled to
the mechanical motion through the field gradients [33, 22] or longitudinal
strain [46, 47, 48].

For simplicity, we assume (though not necessary) that the ensemble spins
interact with the common mechanical mode uniformly with a coupling strength
denoted as $g$. In addition, the probe spin is assumed to have a distinct
coupling strength with the oscillator, denoted as $g_{0}$.

The hybrid system’s Hamiltonian is then formulated as follows:

$H=\omega
a^{\dagger}a+\sum_{k=0}^{N}\omega_{k}S_{k,z}+\left(g_{0}S_{0,z}+g\sum_{k=1}^{N}S_{k,z}\right)\left(a+a^{\dagger}\right),$
(1)

where $S_{i,z}$ signifies the z-component of the spin-1/2 operator for the ith
spin, oscillating at the Larmor frequency $\omega_{i}$. The oscillator’s
annihilation and creation operators are denoted by $a$ and $a^{\dagger}$, with
a frequency $\omega$.

Figure 1: (a) Schematic illustration of the control procedure detailed in the
paper. (a) A clamped nano/micro-mechanical cantilever hosting single defect
spin, that is coupled to its motion through strain and is manipulated by laser
and RF fields. A magnetic tip attached at the edge of the cantilever produces
a gradient field $\nabla\vec{B}$ on a nearby spin ensemble allowing for its
coupling to the mechanical motion of the cantilever. (b) Schematic
illustration of the coupling and control mechanism for the above described
spin-mechanical system. While the probe spin (green) is reset (to state
$\ket{0}$) after every projective measurement, the oscillator ($\rho_{o}$) and
ensemble states ($\rho_{S}$) continue to evolve conditioned on the measurement
of the probe.

The single probe spin is envisioned to be solid state spins in diamond such as
NV center electronic spins, due to good spin coherence properties, high degree
spin control and a well-resolved optical spectrum at low temperatures [49, 17,
16]. The optical interface allows optical excitation to various levels that
can either initialize or readout the spin state with fidelities exceeding
$98\%$ [50]. In addition, these spins can also (weakly) couple to external
spins, photons [43, 51] and phonons [44, 37].

Note that the NV electronic spin-1 can be operated as a qubit with its
$m_{s}=0,-1$ sub-levels by applying a 400 Gauss external magnetic field
aligned with its symmetry axis to separate the $m_{s}=\pm 1$ levels [17]. Such
an arrangement negates the need to consider transverse spin-strain coupling
due to its direct interaction between the $m_{s}=\pm 1$ levels [30, 47, 48],
so we can safely ignore the $m_{s}=+1$ sublevel in this work. What remains is
the longitudinal strain coupling acting on the qubit levels. This number can
be up to 10 kHz for NV centers [52], and higher for other defect centers [53].

The spin ensemble to be cooled can be chosen with more flexibility, as
stringent control over these spins is not a prerequisite for the current
study. Potential candidates for the spin ensemble include distant spins on
other coupled cantilevers, spins at the clamped edge that are strain-coupled
to the mechanical mode, or other surface electron spins. Even more complex
would be the situation where the spins coupled to the mechanical mode are not
the spins attached to the NV centers in diamond but are dark (electronic)
spins on the surface [54] of the NMO, or other ferromagnetic particles [55].
To harness those spins for spin-mechanics one should rely on indirect
polarization techniques as they cannot be directly polarized optically. Note
taht, as schematically shown in Fig. 1(a), the interaction between the
ensemble spins and the mechanical vibrational mode is generated by a magnetic
tip attached to the cantilever with high field gradients. Here the coupling
strength can be much higher [56], than the strain coupling between probe and
the oscillator.

##  III System dynamics

In the following, we will consider the system in the rotating frame of
precession frequencies of the oscillator and the spins. The rotating-frame
Hamiltonian reads

$\tilde{H}=\left(g_{0}S_{0,z}+g\sum_{k=1}^{N}S_{k,z}\right)\left(\hat{a}e^{-i\omega
t}+\hat{a}^{\dagger}e^{i\omega t}\right).$ (2)

Since the ensemble spins form a finite-dimensional Hilbert space, the system
evolution operators can be written in the form [57, 58]:

$\begin{split}U(t)&=\ket{0}\bra{0}\otimes\sum_{\scriptscriptstyle
k=0}^{\scriptscriptstyle N}\mathcal{D}_{k,+}(t)\otimes\mathrm{I}_{k}\\\
&\qquad\qquad+\ket{1}\bra{1}\otimes\sum_{\scriptscriptstyle
n=0}^{\scriptscriptstyle N}\mathcal{D}_{n,-}(t)\otimes\mathrm{I}_{k},\\\
&\mathcal{D}_{k,\pm}(t)=\mathcal{T}\exp\left(\pm
i\int_{0}^{t}\mathrm{d}t^{\prime}\ h_{\scriptscriptstyle
k,\pm}(t^{\prime})\right),\end{split}$ (3)

where $\mathcal{T}$ the time ordering operator, and
$\mathrm{I}_{\scriptscriptstyle k}$ projects the spin ensemble into the
subspace where $k$ spins are pointing up. For example, when $N=3$ and $k=2$,
the projector is written as

$\mathrm{I}_{\scriptscriptstyle
k=2}=\ket{\uparrow\uparrow\downarrow}\bra{\uparrow\uparrow\downarrow}+\ket{\uparrow\downarrow\uparrow}\bra{\uparrow\downarrow\uparrow}+\ket{\downarrow\uparrow\uparrow}\bra{\downarrow\uparrow\uparrow}.$
(4)

The term $h_{\scriptscriptstyle k,\pm}(t)$ is hence the oscillator Hamiltonian
conditioning on the states if the spins. It is expressed as

$h_{\scriptscriptstyle k,\pm}(t)=\frac{g(k-N/2)\pm
g_{0}}{2}\left(\hat{a}e^{-i\omega t}+\hat{a}^{\dagger}e^{i\omega t}\right).$
(5)

It has been shown that the Magnus expansion of the evolution operators can be
simplified as [57, 58]

$\begin{split}\mathcal{D}_{\scriptscriptstyle k,\pm}(t&)=\exp\left(\pm
i\int_{0}^{t}\mathrm{d}t^{\prime}h_{\scriptscriptstyle
k,\pm}(t^{\prime})\right)\\\
&\qquad\times\exp\left(\frac{1}{2}\int_{0}^{t}\mathrm{d}t^{\prime}\int_{0}^{t^{\prime}}\mathrm{d}t^{\prime\prime}\left[h_{\scriptscriptstyle
k,\pm}(t^{\prime}),h_{\scriptscriptstyle
k,\pm}(t^{\prime\prime})\right]\right)\\\ &=\exp\left[\frac{g(k-N/2)\pm
g_{0}}{2\omega}\left(\alpha(t)a^{\dagger}-\alpha^{*}(t)a\right)\right]\\\
&\qquad\quad\times\exp\big{(}it\theta_{\scriptscriptstyle
k}\big{)}.\end{split}$ (6)

where the functions $\alpha(t)$, $\beta(t)$ and $\theta_{\scriptscriptstyle
k}$ are given by

$\begin{split}\alpha(t)&=1-e^{i\omega t},\quad\beta(t)=\frac{\omega
t-sin(\omega t)}{\omega t},\\\ \theta_{\scriptscriptstyle
k}&=-\beta(t)\left(g(k-N/2)\pm g_{0}\right)^{2}/4\omega.\end{split}$ (7)

##  IV Probe spin control

In addition to the natural dynamics induced by the spin-oscillator interaction
given in Eq. (1), we also apply a periodic sequence of $\pi$ pulses on the
probe spin. These pulses significantly alter the dynamics of the system by
embedding a time-periodic function $f(t)$ into the Hamiltonian, described as
follows:

$f(t)=\begin{cases}+1&\text{if }2j\tau<t<(2j+1)\tau\\\ -1&\text{if
}(2j+1)\tau<t<(2j+2)\tau\end{cases}$ (8)

where $j=0,1,2,\cdots,n$ and $\tau$ denotes the interval between successive
pulses. By replacing $g_{0}$ with $f(t)g_{0}$ in $h_{\scriptscriptstyle
k,\pm}(t)$,the conditional evolution operators in Eq (6) are now written as

$\begin{split}\mathcal{D}^{f}_{\scriptscriptstyle
k,\pm}(t)=&\exp\Big{[}\frac{g(k-N/2)}{2\omega}\left(\alpha(t)a^{\dagger}-\alpha^{*}(t)a\right)\\\
&\;\pm\frac{f(t)g_{0}}{2\omega}\left(F(\epsilon,t)a^{\dagger}-F^{*}(\epsilon,t)a\right)\Big{]}\exp\left[i\theta_{\scriptscriptstyle
k}(t)\right],\end{split}$ (9)

where the filter function $F(t)$ is given by [57]

$\begin{split}F(t)&=i\omega\
\sum_{k=0}^{n}\int\limits_{k\tau}^{(k+1)\tau}\mathrm{d}t\;(-1)^{k}e^{-i\omega
t}\\\
&=-\sum_{k=0}^{n}(-1)^{k}\left[e^{-i(k+1)\omega\tau}-e^{-ik\omega\tau}\right]\\\
&=\left[1-(-1)^{n+1}e^{-i(n+1)\omega\tau}+2\sum_{k=1}^{n}(-1)^{k}e^{-ik\omega\tau}\right].\end{split}$
(10)

This filter function can be expressed in a well-known form, centered at
$\tau=\pi/\omega$ and with its width decreasing approximately as $1/n^{2}$
[58, 57, 59]:

$F^{2}(\tau,n,\omega)=4\tan^{2}\left(\frac{\omega\tau}{2n+2}\right)\cos^{2}\left(\frac{\omega(n+1)\tau}{2}\right),$
(11)

The asymptotic behavior of the filter function is given by

$\lim\limits_{n\rightarrow\infty}F(t)=2(n+1)\delta(\omega-\frac{\pi}{\tau}).$
(12)

In practice, to ensure only the mode at the particular frequency $\pi/\tau$
couples to the probe spin, the number of pulses on the probe should be
$n>\mathcal{O}(\omega/g_{0}$). Other phonon modes are effectively decoupled
[58], which supports our initial assumption of a single oscillator mode in Eq.
(1).

Furthermore, our analysis focuses on the scenario where the inversion pulses
resonate with the oscillator frequency. To achieve resonance, we set the
interval between successive inversion pulses as

$\tau=\frac{\pi}{\omega-\epsilon},\text{ where }\frac{\epsilon}{\omega}\ll 1,$
(13)

where the detuning parameter $\epsilon$ is intentionally kept significantly
smaller than the oscillator frequency $\omega$ to facilitate near-resonant
conditions. The filter function is now simplified to

$F(\epsilon,t=(n+1)\tau)=1-e^{i\epsilon t}+\frac{2(1-e^{i\epsilon
t})}{e^{-i\epsilon\tau}-1}.$ (14)

@C=.8em @R=1em Oscillator & Oscillator D^f_k, ±(t)

Ensemble  Ensemble  -1

Probe  —0⟩ Probe  —0⟩ H -1 H Z

Figure 2: Quantum circuit illustrating the implementation of a projector on
the oscillator, influenced by the spin ensemble states and the measurement
result of the probe spin. The process starts with the ancilla spin in the
$\ket{+}$ state, followed by the application of repeated DD pulses. This leads
to a $\mathcal{D}^{f}_{\scriptscriptstyle k,\pm}(t)$ controlled evolution
guided by the spin ensemble $P_{k}$. The procedure ends with a $\pi/2$ pulse
and measurement of the ancilla spin in the Pauli Z basis. The specific
projector realized is expressed mathematically in Eq (15).

In our investigation, we focus specifically on the weak dispersive coupling
between each spin and the oscillator, deliberately ignoring any direct
interactions between the spins. Nevertheless, by carefully choosing the
detuning $\epsilon$ and the number of pulses $n$, substantial changes in the
dynamics of the spin ensemble and the oscillator can be induced. The
modification of the dynamics is facilitated by a cyclic process involving the
projective measurement and reinitialization of the ancilla spin, as
schematically depicted in Fig. 1. Initially, the ancilla spin is prepared in
the state $\ket{+}$. Then the system undergoes a DD-modulated evolution for a
time $t=n\tau$, as described in Eq (3). Finally, the probe spin is measured in
the Pauli-Z basis, following the application of an additional $\pi/2$ pulse to
it.

The process described above is represented by the circuit in Fig. 2. It leads
to the implementation of a projector throughout the system depending on the
probe measurement outcome, as mathematically formulated:

$\mathcal{P}_{\scriptscriptstyle N,\pm}=\sum_{\scriptscriptstyle
k=0}^{\scriptscriptstyle N}V_{k,\pm}\otimes\mathrm{I}_{k},\quad
V_{k,\pm}=\frac{1}{2}\left(\mathcal{D}^{\scriptscriptstyle
f}_{\scriptscriptstyle k,+}(t)\pm\mathcal{D}^{\scriptscriptstyle
f}_{\scriptscriptstyle k,-}(t)\right),$ (15)

where $V_{k,\pm}$ is the projector on the oscillator, $k$ is the number of
ensemble spins pointing up, and $\pm$ is the result of the ancilla
measurement. The repeated application of such projectors is central to the
cooling of the oscillator in our previous work [58], which is now extended to
cooling/polarizing an entire spin ensemble in this work.

Figure 3: We plot the probability of projecting a spin ensemble into its
possible total magnetization states (from $-N/2$ to $N/2$) after repeated
conditional measurements of the probe spin in the state $\ket{0}$, as
described in Eq. (18). Here we consider an ensemble with $N=14$ spins, an
oscillator with frequency $\omega=1.2$MHz that is weakly coupled to the probe
and spin ensembles with strengths of $g_{0}=g=10$kHz. The probe is modulated
periodically with $n=128$ $\pi$-pulses with an interval
$\tau=\pi/{\omega-\epsilon}$. Here the detuning is set as $\epsilon=1.1$kHz.
The probability distribution reaches saturation after $M=100$ measurements.

##  V System cooling

In this section, we describe the protocol that imposes a thermal filter on the
hybrid system, allowing for controlled cooling of either the oscillator, the
spin ensemble, or both. We start with $N$ spins in the ensemble, which are
initially assumed to be in a fully mixed state, and the oscillator is assumed
to be in its thermal state. These initial states are written as

$\rho_{\text{ens}}=\sum_{k=0}^{N}\frac{\mathrm{I}_{k}}{2^{N}},\quad\rho_{\text{osc}}=\frac{1}{\sum_{n}e^{-\frac{n\omega}{k_{\scalebox{0.5}{{\hbox{B}}}}T}}}\sum_{n}e^{-\frac{n\omega}{k_{\scalebox{0.5}{{\hbox{B}}}}T}}\ket{n}\bra{n}.$
(16)

Here the projector $\mathrm{I}_{k}$ represents the states with $k$ spins
pointing up, as shown in Eq. (4); $k_{\scalebox{0.5}{$B$}}$ represents the
Boltzmann constant; and $T$ is the temperature. The ensemble spins are fully
unpolarized, i.e., $\langle\sum_{j=1}^{N}S_{z,j}\rangle=0$; and the initial
thermal occupancy of the oscillator is calculated as

$n_{0}=\operatorname{Tr}(\rho_{osc}\,a^{\dagger}a)=\frac{1}{e^{(n\omega/k_{\scalebox{0.5}{{\hbox{B}}}}T)}-1}.$
(17)

Following the circuit depicted in Fig. 2, we measure the probe spin, selecting
only outcomes of 0. This process effectively implements the projector
$\mathcal{P}_{\scriptscriptstyle N,+}$ as defined in Eq. (15). Upon performing
$M$ such measurements, the probability that the ensemble is in the
$\mathrm{I}_{k}$ sector is calculated as

$\mathrm{P}_{m}=\frac{1}{\sum_{m}\mathrm{P}_{m}}\mathrm{Tr}[\mathcal{P}_{\scriptscriptstyle
N,+}^{M}\,\rho_{B}(0)\,\mathcal{P}_{\scriptscriptstyle N,+}^{M}],$ (18)

where $m=(2k-N)/2$ is the total magnetization of the ensemble. This dependence
on $m$ shows that, with appropriate parameter choices, it is possible to
achieve $\mathrm{P}{\scriptscriptstyle m=N/2}\approx 1$, thereby polarizing
the ensemble. This ensemble polarization depends non-trivially on the filter
function imposed by the pulses on the probe spin, as can be seen from Eq.
(14). It gives us a new control parameter $\epsilon$ that would affect the
projection of the spin ensemble onto different $m$-states.

In this work, we will focus on the weak coupling regime, i.e., $g\ll\omega$,
and we choose the detuning $\epsilon\sim\mathcal{O}(g)$ [31] so that
$\epsilon\ll\omega$. Specifically, we consider the case where the spin
couplings $g,~{}g_{0}$ are $10$ kHz and the mechanical oscillation frequency
is chosen to be $1.2$ MHz. Consequently, for a spin ensemble of size $N=14$,
by fixing the detuning to be $\epsilon=1.16$ kHz, one can see that the sector
with all spins up ($m=7$) is predominantly left, effectively polarizing the
ensemble. To further illustrate the non-trivial dependence of the maximum
ensemble polarization on the detuning $\epsilon$, we vary its value from 0 to
2 kHz and also vary the number of ensemble spins. Note that we always assume
$M=100$ projections of the probe spin.

Figure 4: Non-trivial dependence of the final maximum polarization on the
detuning $\epsilon$ for varying spin-ensemble sizes ($N$). The y-axis
represents the probability of finding all ensemble spins pointing up after
$M=100$ probe projection measurements, denoted as
$\mathrm{P}_{\scriptscriptstyle m=N/2}(M)$. Note that the oscillator frequency
is set as $\omega=1.2$MHz and the spin-oscillator coupling strengths are set
as $g_{0}=g=10$kHz.

A critical factor in achieving ensemble polarization is to ensure that the
interval between successive pulses on the probe is close to resonance, i.e., a
negligible detuning, $\epsilon\approx 0$, as indicated in Eq. (13). The reason
for this is to prevent premature cooling of the oscillator, since its thermal
occupancy significantly affects the interaction between the probe and the
ensemble, as detailed in Eq. (1). To allow the oscillator to cool, the
detuning must be set to a higher value to facilitate its rapid cooling [58].

In Fig. 5 we show a pulse sequence that polarizes an ensemble of $N=4$ spins
and cools an oscillator starting from an initial thermal occupancy of
$n_{0}=45$. We first optimize the cooling of the spin-ensemble by varying the
detuning, to find an optimal point $\epsilon_{1}=1.1$kHz where the ensemble
polarizes perfectly, while the oscillator occupancy remains unchanged. This
behavior is similar to what one would expect for resonant case ($\epsilon=0$).
Now for another $50$ probe measurements and shifting the detuning to a higher
value of $\epsilon_{2}=5$kHz, would rapidly cool the oscillator as shown in
the right part of Fig. 5. Shifting $\epsilon_{2}$ to even higher values would
require more inversion pulses on the probe to achieve a similar effect, as can
be seen from Eq. (12). We would like to note that there would be no effect of
the dephasing caused by the surrounding spin-bath on the cooling dynamics, as
each measurement also partially collapses the spin-bath towards steady states
that eventually stabilize the cooling mechanism [60]. Other noise sources only
affect the cooling rates, but not the final steady states of the ensemble and
the oscillator.

Figure 5: Cooling the oscillator and spin ensemble. This figure shows the
polarization $\mathrm{P}_{z}$ and the thermal occupancy $n(\omega)$ for both
the spin ensemble and the oscillator, illustrating how the
cooling/polarization progresses with the number of measurements $M$. The
effective cooling of these quantum systems occurs in different parameter
regimes. Initially, the spin ensemble is cooled by adjusting the $\pi$ pulse
intervals to near-resonance conditions, specifically $\tau_{1}$ such that
$\epsilon_{1}=0.11$kHz, where the thermal occupancy of the oscillator remains
largely unaffected. After the spin ensemble has cooled, changing the pulse
interval to $\tau_{2}$, which produces a more significant off-resonance with
$\epsilon_{2}=1.5$kHz, leads to a rapid cooling of the oscillator to its
ground state. The coupling parameters are also $g_{0}=g=10$kHz, with an
oscillator frequency $\omega=1.2$MHz. Here, we have considered an initial
oscillator thermal occupation of $n_{0}\approx 45$, and a spin ensemble size
of $N=4$ so that $\mathrm{P}_{max}=2$.

A subtle but crucial aspect of cooling the oscillator is to choose the
detuning $\epsilon$ precisely so that the conditional displacement operator
$D(\theta)$ in Fig. 2 is not perfectly aligned with the position quadrature.
Note that when the detuning is zero, $D(\theta)$ is just perfectly aligned
with the momentum quadrature. In both cases the oscillator cannot be cooled
[58]. This observation is intriguingly related to Gottesman-Kitaev-Preskill
(GKP) state encoding by phase estimation, which uses the same control of
applying inversion pulses to the probe spin and its measurement. However, the
detuning $\epsilon$ must now be carefully chosen for each measurement to
ensure that the conditional oscillator displacement is exactly aligned with
either the momentum or position quadratures. In addition, we must carefully
tune the magnitude of the oscillator displacement. Starting with an oscillator
in its zero photon ground state, this process actually increases the number of
photons and allows the preparation of GKP states in the oscillator [61, 62].

In the Appendix A, we discuss this proposal in more detail, where it becomes
clear that encoding GKP states in a nanomechanical oscillator would be much
more challenging. This is because a much stronger spin-mechanical coupling
strength is needed, since the coherence time of mechanical vibration is rather
limited.

##  VI Discussion

In this research, we have emphasized the profound influence that the precise
control of a single spin can have on the dynamics of ensemble oscillator
systems, and demonstrated this effect in the regime of weak coupling. This
finding is central to a wide range of quantum technological applications.
Specifically, we have introduced a ground-state cooling protocol for a
tripartite hybrid quantum system that includes a spin ensemble, an oscillator,
and a single probe spin. By exploiting the weak dispersive coupling between
the spins and the oscillator, together with the feedback from frequent
measurements of the probe spin, our protocol can cool both the spin ensemble
and the oscillator to their ground states, either individually or in tandem.

Our numerical simulations over various parameter settings not only confirm the
practicality of our cooling protocol, but also provide a basis for its
experimental implementation. For example, our analysis shows that a spin
ensemble of $N=14$ spins could be cooled with a spin-mechanical coupling
strength as low as 1 kHz at a mechanical oscillation frequency of about 1 MHz.
Such coupling strengths, though are on the higher side of current state-of-
the-art experiments, can become realizable by further reducing dimensions of
the mechanical oscillator [63]. For other solid-state defects in silicon
carbide, the coupling strengths could reach a few MHz [53], allowing us to
attain higher cooling rates.

Another notable limitation of our cooling scheme is its reliance on post-
selection, which leads to an exponential decrease in success rate as the
number of probe measurements increases. Nevertheless, our protocol offers a
distinct advantage in its applicability to cooling large spin ensembles that
are otherwise difficult to control and cannot be cooled by conventional
optical or energy transfer methods. Furthermore, for spin ensembles that are
amenable to such traditional cooling methods, such as NV centers in diamond,
our approach adds to the toolbox and serves as a complementary polarization
method. For example, optical methods could initially polarize an NV ensemble
to a certain degree, after which our protocol could be used to further cool
the ensemble toward the ground state, while maintaining a reasonable success
rate. This will help us to alleviate the heating problem caused by intensive
optical initialization.

Finally, this work also explores the potential of extending the cooling scheme
to encode complex quantum states in oscillators, with a particular focus on
the GKP states. Despite the important role of these states in quantum error
correction [61, 62] and quantum sensing [64], their creation in macroscopic
objects, a such as a nano-mechanical oscillator considered in this work,
promises to deepen our understanding of the boundary between the classical and
quantum worlds [65, 66]. Despite the inherent challenges, in particular the
need for much stronger coupling strengths for experimental realization, our
results provide insights for future efforts in this direction.

## Acknowledgement

The authors would like to thank Ruoming Peng and Jixing Zhang for feedback on
the manuscript. This work is supported by German Federal Ministry of Education
and Research (BMBF) for the project QECHQS under Grant No. 16KIS1590K.

##  Appendix A Encoding a qubit into a macroscopic oscillator

In this section, we consider a hybrid system consisting of only the probe spin
and the oscillator. We show that we can prepare a macroscopic mechanical
oscillator in a complicated quantum state by controlling the probe. The
exploration of such quantum superposition states in macroscopic objects has
emerged as a significant area of interest, mainly because of its potential to
explore the frontier between classical and quantum physics [66, 65].

Our proposal is realized by carefully tuning the detuning ($\epsilon$) and
applying a certain number of inversion pulses to the auxiliary spin. The
quantum state in which we propose to place a macroscopic mechanical oscillator
is the so-called Gottesman-Kitaev-Preskill (GKP) state [61, 67], which are
defined as the simultaneous eigenstates of the two stabilizers, i.e.,

$S_{q}=e^{-i2\sqrt{\pi}\hat{p}},\,\,S_{p}=e^{i2\sqrt{\pi}\hat{q}}.$ (19)

The logical operators are correspondingly given by $X_{\scriptscriptstyle
L}=\sqrt{S_{p}}$ and $Z_{\scriptscriptstyle L}=\sqrt{S_{q}}$, respectively.

However, the ideal GKP states are unphysical, requiring infinite amount of
energy to prepare [61]. Instead, we will consider the approximate GKP states
that are defined as (up to normalization):

$\begin{split}&\ket{\tilde{0}}\propto\sum_{s=-\infty}^{\infty}\int_{-\infty}^{\infty}e^{-\frac{\Delta^{2}}{2}(2s)^{2}\pi}e^{-\frac{1}{2\Delta^{2}}(q-2s\sqrt{\pi})^{2}}\ket{q}\mathrm{d}q,\\\
&\ket{\tilde{1}}\propto\sum_{s=-\infty}^{\infty}\int_{-\infty}^{\infty}e^{-\frac{\Delta^{2}}{2}(2s+1)^{2}\pi}e^{-\frac{1}{2\Delta^{2}}(q-(2s+1)\sqrt{\pi})^{2}}\ket{q}\mathrm{d}q.\end{split}$
(20)

Note that the closer the squeezing parameter $\Delta$ is to 0, the higher is
the quality of the GKP states [61, 67].

The GKP code is a promising candidate for quantum error correction (QEC) in
continuous variable systems [68], such as photonic systems and cavity modes.
In addition, GKP states hold promise for enhanced sensing capabilities,
allowing simultaneous detection of an oscillator’s displacements in both its
position and momentum quadratures [64].

###  A Oscillator displacements along the momentum and position quadratures

In the following, we assume that the spin ensemble is effectively decoupled
from the rest of the system. The projectors given in Eq. (15) will then become
just a displacement operator of the oscillator conditioning on the measurement
result of the probe spin, as indicated by:

$\mathcal{P}_{\scriptscriptstyle N=0,\pm}=\frac{D(\theta)\pm
e^{i\phi}D(-\theta)}{2},$ (21)

where $D(\theta)=\exp{(\theta a^{\dagger}-\theta^{*}a)}$ represents the
displacement operator, with the amplitude $\theta=g_{0}F(\epsilon,t)/2\omega$.

The additional phase $\phi$ is introduced by rotating the ancilla spin by an
angle $\phi$ before the measurement, as shown in the circuit in Fig. 6. The
angle $\phi$ selects the effective measurement basis of the ancilla spin.

The proposed GKP encoding scheme involves repeated execution of this circuit,
adaptively setting the $\phi$ angle based on previous measurement results.
This is a phase estimation protocol that incrementally identifies the
eigenvalue of $D(2\theta)$ while simultaneously projecting the oscillator into
the corresponding eigenstate.

#### Displacement along position quadrature

For $\epsilon=0$, the filter function simplifies to $F(\epsilon,t)=2(n+2)$ so
that $\theta=g_{0}t/\pi$. This leads to conditional displacements exactly
along the position quadrature:

$D(\theta)=\exp\Big{[}\mp
i\frac{\sqrt{2}g_{0}t}{\pi}\hat{p}\Big{]},\,\,\text{with }\epsilon=0,$ (22)

where $\hat{p}=\frac{i}{\sqrt{2}}(\hat{a}-\hat{a}^{\dagger})$ denotes the
momentum operator of the oscillator. According to Ref. [33], this ability to
adjust the oscillator’s momentum based on the spin state facilitates a
mechanism for entangling two spins by the oscillator.

#### Displacement along momentum quadrature

Additionally, we also need to realize controlled oscillator-displacements
along the momentum quadrature. This is achieved by setting
$\epsilon=\omega/(n+1)$ so that $\epsilon t=\pi$, the filter function reads

$F(\epsilon,t)=2+\frac{4}{e^{i\epsilon\tau}-1}=\frac{4}{i\epsilon\tau}+\mathcal{O}(\epsilon\tau),$
(23)

where $\hat{q}=\frac{1}{\sqrt{2}}(\hat{a}+\hat{a}^{\dagger})$ is the position
operator of the oscillator. Correspondingly, we have
$\theta=ig_{0}t(n+1)/\pi^{2}(n+2)$. Furthermore, neglecting the higher order
terms of $\epsilon\tau\sim\mathcal{O}(1/n)$ leads to the following
displacement operator of the oscillator:

$D(\theta)=\exp\Big{[}\mp
i\frac{2\sqrt{2}g_{0}t}{\pi^{2}}\frac{n+1}{n+2}\hat{q}\Big{]},\,\,\text{with
}\epsilon=\frac{\omega}{n+1}.$ (24)

#### State preparation via phase estimation

@C=.8em @R=1em Oscillator& Oscillator D(±θ)

Ancilla spin  —0⟩ Ancilla spin  —0⟩ H -1 R_Z(ϕ) H Z

Figure 6: Single round of an adaptive phase estimation protocol that estimates
the eigenvalue of oscillator displacement operator, $D(2\theta)=\exp{(2\theta
a^{\dagger}-2\theta^{*}a)}$ with amplitude
$2\theta=g_{0}F(\epsilon,t)/\omega$. Initially, the ancillary spin is set in
the $\ket{+}$ state and subsequently measured in the x-basis after applying a
series of periodic inversion pulses, which facilitates the implementation of
controlled oscillator displacements. Before measuring the ancilla spin, it
undergoes a z-axis rotation by an angle $\phi$. This angle $\phi$ is
strategically determined based on the ancilla measurement results of previous
phase estimation cycles, with the goal of minimize the uncertainty in the
posterior probability distribution estimate of $D(2\theta)$’s eigenvalue.
Specifically, when $\phi=0$, this circuit simply executes the projectors given
in Eq. (21). Through repeated application of this phase estimation circuit,
the oscillator states outside the subspace spanned by the eigenstates of
$D(2\theta)$ are gradually suppressed.

Phase estimation has been shown to be effective for preparing such approximate
GKP logical states by determining the eigenvalues of stabilizers and logical
operators [62]. For the spin-oscillator hybrid system considered here, such
eigenvalue determination is realized by adjusting the detuning magnitude
$\epsilon$ and the total evolution duration $t$. In this way, we can choose
the magnitude by which the oscillator is shifted along the momentum or
position quadrature.

For example, to prepare the approximate GKP state $\ket{\tilde{0}}$, first
initialize the oscillator in its ground state with zero photons. It then sets
$D(\theta)=S_{q}$ and runs the phase estimation circuit several times to
estimate its eigenvalue with high confidence. A subsequent fit to
$D(\theta)=Z_{\scriptscriptstyle L}$ allows for similar eigenvalue estimation.
Determining the eigenvalues of $S_{q}$ and $Z_{\scriptscriptstyle L}$ projects
the oscillator into the desired approximate GKP state. This process can be
described as iterative and adaptive application of the circuit in Fig. 6.

#### GKP encoding and oscillator cooling

The GKP encoding scheme we discuss is similar to the oscillator cooling
technique described in Ref. [58], with both methods using the phase estimation
circuit shown in Fig. 6 many times. However, the encoding scheme starts from
the ground state, as opposed to the cooling method, which applies the circuit
to a thermal state and then cools it down to the ground state.

A crucial difference between the two lies in how the circuit is used in the
cooling process, where it is implemented with non-zero detuning, so that
neither Eq.(16) nor Eq.(18) hold. When these conditions are satisfied, instead
of cooling to a zero photon ground state, the oscillator, initially in a
thermal state, stabilizes in a state that maintains a finite number of photons
[58].

###  B Experimental feasibility

In contrast to the oscillator cooling scheme developed in Ref. [58], which
requires only a weak ancilla-oscillator coupling, the GKP encoding scheme is
much more demanding as it requires a strong coupling strength. Current
experimental configurations are clearly inadequate in this regard to implement
the preparation scheme. Furthermore, the GKP encoding scheme has been
extensively studied in Ref. [62].

Therefore, we choose not to perform detailed numerical simulations to evaluate
the practicality of the GKP encoding scheme in this study. Instead, we briefly
assess the challenges and discuss potential improvements to existing setups
that could facilitate the preparation of GKP states in macroscopic mechanical
oscillators. The time duration for each phase estimation round needs to be of
the order of $1/g_{0}$, so that the displacement magnitude can be of the order
of $|\theta|\sim O(\sqrt{\pi})$. Given the necessity of around 10 phase
estimation cycles to yield an approximate GKP state of acceptable quality
[62], the preparation scheme requires an overall time exceeding
$\mathcal{O}(10/g_{0})$.

From an experimental point of view, it’s crucial that the coherence time of
the oscillator is considerably long compared to the coupling strength $g_{0}$,
ideally at least on the order of $100/g_{0}$. Diamond NMO with a quality
factor $Q>10^{6}$ have been demonstrated at sub-Kelvin temperatures and
mechanical damping rates of about a few kHz [69]. This implies that the
oscillator coherence time needs to be significantly improved. Note that the
coherence time of the NV electronic spin in cryogenic conditions can exceed 1
second through the application of dynamical decoupling [12]. This is far from
sufficient, as the probe spin only needs to maintain coherence during each
phase estimation round.

## References

  * Acharya  _et al._ [2023] R. Acharya, I. Aleiner, R. Allen, T. I. Andersen, M. Ansmann, F. Arute, K. Arya, A. Asfaw, J. Atalaya, R. Babbush, D. Bacon, J. C. Bardin, J. Basso, A. Bengtsson, S. Boixo, G. Bortoli, A. Bourassa, J. Bovaird, L. Brill, M. Broughton, B. B. Buckley, D. A. Buell, T. Burger, B. Burkett, N. Bushnell, Y. Chen, Z. Chen, B. Chiaro, J. Cogan, R. Collins, P. Conner, W. Courtney, A. L. Crook, B. Curtin, D. M. Debroy, A. Del Toro Barba, S. Demura, A. Dunsworth, D. Eppens, C. Erickson, L. Faoro, E. Farhi, R. Fatemi, L. Flores Burgos, E. Forati, A. G. Fowler, B. Foxen, W. Giang, C. Gidney, D. Gilboa, M. Giustina, A. Grajales Dau, J. A. Gross, S. Habegger, M. C. Hamilton, M. P. Harrigan, S. D. Harrington, O. Higgott, J. Hilton, M. Hoffmann, S. Hong, T. Huang, A. Huff, W. J. Huggins, L. B. Ioffe, S. V. Isakov, J. Iveland, E. Jeffrey, Z. Jiang, C. Jones, P. Juhas, D. Kafri, K. Kechedzhi, J. Kelly, T. Khattar, M. Khezri, M. Kieferová, S. Kim, A. Kitaev, P. V. Klimov, A. R. Klots, A. N. Korotkov, F. Kostritsa, J. M. Kreikebaum, D. Landhuis, P. Laptev, K.-M. Lau, L. Laws, J. Lee, K. Lee, B. J. Lester, A. Lill, W. Liu, A. Locharla, E. Lucero, F. D. Malone, J. Marshall, O. Martin, J. R. McClean, T. McCourt, M. McEwen, A. Megrant, B. Meurer Costa, X. Mi, K. C. Miao, M. Mohseni, S. Montazeri, A. Morvan, E. Mount, W. Mruczkiewicz, O. Naaman, M. Neeley, C. Neill, A. Nersisyan, H. Neven, M. Newman, J. H. Ng, A. Nguyen, M. Nguyen, M. Y. Niu, T. E. O’Brien, A. Opremcak, J. Platt, A. Petukhov, R. Potter, L. P. Pryadko, C. Quintana, P. Roushan, N. C. Rubin, N. Saei, D. Sank, K. Sankaragomathi, K. J. Satzinger, H. F. Schurkus, C. Schuster, M. J. Shearn, A. Shorter, V. Shvarts, J. Skruzny, V. Smelyanskiy, W. C. Smith, G. Sterling, D. Strain, M. Szalay, A. Torres, G. Vidal, B. Villalonga, C. Vollgraff Heidweiller, T. White, C. Xing, Z. J. Yao, P. Yeh, J. Yoo, G. Young, A. Zalcman, Y. Zhang, and N. Zhu, Suppressing quantum errors by scaling a surface code logical qubit, Nature 614, 676 (2023).
  * Bluvstein  _et al._ [2023] D. Bluvstein, S. J. Evered, A. A. Geim, S. H. Li, H. Zhou, T. Manovitz, S. Ebadi, M. Cain, M. Kalinowski, D. Hangleiter, J. P. Bonilla Ataides, N. Maskara, I. Cong, X. Gao, P. Sales Rodriguez, T. Karolyshyn, G. Semeghini, M. J. Gullans, M. Greiner, V. Vuletić, and M. D. Lukin, Logical quantum processor based on reconfigurable atom arrays, Nature 626, 58 (2023).
  * Cai  _et al._ [2013] J. Cai, A. Retzker, F. Jelezko, and M. B. Plenio, A large-scale quantum simulator on a diamond surface at room temperature, Nat. Phys. 9, 168 (2013).
  * Wehner  _et al._ [2018] S. Wehner, D. Elkouss, and R. Hanson, Quantum internet: A vision for the road ahead, Science 362, 10.1126/science.aam9288 (2018).
  * Zhang  _et al._ [2022] W. Zhang, T. van Leent, K. Redeker, R. Garthoff, R. Schwonnek, F. Fertig, S. Eppelt, W. Rosenfeld, V. Scarani, C. C.-W. Lim, and H. Weinfurter, A device-independent quantum key distribution system for distant users, Nature 607, 687 (2022).
  * Aslam  _et al._ [2023] N. Aslam, H. Zhou, E. K. Urbach, M. J. Turner, R. L. Walsworth, M. D. Lukin, and H. Park, Quantum sensors for biomedical applications, Nat. Rev. Phys. 5, 157 (2023).
  * Shi  _et al._ [2015] F. Shi, Q. Zhang, P. Wang, H. Sun, J. Wang, X. Rong, M. Chen, C. Ju, F. Reinhard, H. Chen, J. Wrachtrup, J. Wang, and J. Du, Single-protein spin resonance spectroscopy under ambient conditions, Science 347, 1135 (2015).
  * Hensen  _et al._ [2015] B. Hensen, H. Bernien, A. E. Dréau, A. Reiserer, N. Kalb, M. S. Blok, J. Ruitenberg, R. F. L. Vermeulen, R. N. Schouten, C. Abellán, W. Amaya, V. Pruneri, M. W. Mitchell, M. Markham, D. J. Twitchen, D. Elkouss, S. Wehner, T. H. Taminiau, and R. Hanson, Loophole-free bell inequality violation using electron spins separated by 1.3 kilometres, Nature 526, 682 (2015).
  * Battistel  _et al._ [2023] F. Battistel, C. Chamberland, K. Johar, R. W. J. Overwater, F. Sebastiano, L. Skoric, Y. Ueno, and M. Usman, Real-time decoding for fault-tolerant quantum computing: progress, challenges and outlook, Nano Futures 7, 032003 (2023).
  * Lordi and Nichol [2021] V. Lordi and J. M. Nichol, Advances and opportunities in materials science for scalable quantum computing, MRS Bulletin 46, 589 (2021).
  * Wang [2023] Y. Wang,  _Using spins in diamond for quantum technologies_ , Ph.D. thesis, Delft University of Technology (2023).
  * Abobeih  _et al._ [2018] M. H. Abobeih, J. Cramer, M. A. Bakker, N. Kalb, M. Markham, D. J. Twitchen, and T. H. Taminiau, One-second coherence for a single electron spin coupled to a multi-qubit nuclear-spin environment, Nat. Commun. 9, 10.1038/s41467-018-04916-z (2018).
  * Wood  _et al._ [2022] B. D. Wood, G. A. Stimpson, J. E. March, Y. N. D. Lekhai, C. J. Stephen, B. L. Green, A. C. Frangeskou, L. Ginés, S. Mandal, O. A. Williams, and G. W. Morley, Long spin coherence times of nitrogen vacancy centers in milled nanodiamonds, Phys. Rev. B 105, 205401 (2022).
  * Hermans  _et al._ [2022] S. L. N. Hermans, M. Pompili, H. K. C. Beukers, S. Baier, J. Borregaard, and R. Hanson, Qubit teleportation between non-neighbouring nodes in a quantum network, Nature 605, 663 (2022).
  * Pompili  _et al._ [2021] M. Pompili, S. L. N. Hermans, S. Baier, H. K. C. Beukers, P. C. Humphreys, R. N. Schouten, R. F. L. Vermeulen, M. J. Tiggelman, L. dos Santos Martins, B. Dirkse, S. Wehner, and R. Hanson, Realization of a multinode quantum network of remote solid-state qubits, Science 372, 259 (2021).
  * Bradley  _et al._ [2019] C. E. Bradley, J. Randall, M. H. Abobeih, R. C. Berrevoets, M. J. Degen, M. A. Bakker, M. Markham, D. J. Twitchen, and T. H. Taminiau, A ten-qubit solid-state spin register with quantum memory up to one minute, Phys. Rev. X 9, 031045 (2019).
  * Abobeih  _et al._ [2022] M. H. Abobeih, Y. Wang, J. Randall, S. J. H. Loenen, C. E. Bradley, M. Markham, D. J. Twitchen, B. M. Terhal, and T. H. Taminiau, Fault-tolerant operation of a logical qubit in a diamond quantum processor, Nature 606, 884 (2022).
  * Waldherr  _et al._ [2014] G. Waldherr, Y. Wang, S. Zaiser, M. Jamali, T. Schulte-Herbrüggen, H. Abe, T. Ohshima, J. Isoya, J. F. Du, P. Neumann, and J. Wrachtrup, Quantum error correction in a solid-state hybrid spin register, Nature 506, 204 (2014).
  * Xiang  _et al._ [2013] Z.-L. Xiang, S. Ashhab, J. Q. You, and F. Nori, Hybrid quantum circuits: Superconducting circuits interacting with other quantum systems, Rev. Mod. Phys. 85, 623 (2013).
  * Kurizki  _et al._ [2015] G. Kurizki, P. Bertet, Y. Kubo, K. Mølmer, D. Petrosyan, P. Rabl, and J. Schmiedmayer, Quantum technologies with hybrid systems, Proc. Natl. Acad. Sci. 112, 3866 (2015).
  * Hei  _et al._ [2023] X.-L. Hei, P.-B. Li, X.-F. Pan, and F. Nori, Enhanced tripartite interactions in spin-magnon-mechanical hybrid systems, Phys. Rev. Lett. 130, 073602 (2023).
  * Kolkowitz  _et al._ [2012] S. Kolkowitz, A. C. Bleszynski Jayich, Q. P. Unterreithmeier, S. D. Bennett, P. Rabl, J. G. E. Harris, and M. D. Lukin, Coherent sensing of a mechanical resonator with a single-spin qubit, Science 335, 1603 (2012).
  * Marcos  _et al._ [2010] D. Marcos, M. Wubs, J. M. Taylor, R. Aguado, M. D. Lukin, and A. S. Sørensen, Coupling nitrogen-vacancy centers in diamond to superconducting flux qubits, Phys. Rev. Lett. 105, 210501 (2010).
  * Zhu  _et al._ [2011] X. Zhu, S. Saito, A. Kemp, K. Kakuyanagi, S.-i. Karimoto, H. Nakano, W. J. Munro, Y. Tokura, M. S. Everitt, K. Nemoto, M. Kasu, N. Mizuochi, and K. Semba, Coherent coupling of a superconducting flux qubit to an electron spin ensemble in diamond, Nature 478, 221 (2011).
  * Xie  _et al._ [2024] L. Xie, Z. Yan, L. Wang, D. Wang, J. Liu, Y. Song, W. Xiong, and M. Wang, Unitary and efficient spin squeezing in cavity optomechanics, arXiv:2401.15553 (2024).
  * Fedele  _et al._ [2024] F. Fedele, F. Cerisola, L. Bresque, F. Vigneau, J. Monsel, J. Tabanera, K. Aggarwal, J. Dexter, S. Sevitz, J. Dunlop, A. Auffèves, J. Parrondo, A. Pályi, J. Anders, and N. Ares, Coupling a single spin to high-frequency motion, arXiv:2402.19288 (2024).
  * Körber  _et al._ [2023] L. Körber, C. Heins, T. Hula, J.-V. Kim, S. Thlang, H. Schultheiss, J. Fassbender, and K. Schultheiss, Pattern recognition in reciprocal space with a magnon-scattering reservoir, Nat. Commun. 14, 10.1038/s41467-023-39452-y (2023).
  * Fujii and Nakajima [2017] K. Fujii and K. Nakajima, Harnessing disordered-ensemble quantum dynamics for machine learning, Phys. Rev. Appl. 8, 024030 (2017).
  * Negoro  _et al._ [2018] M. Negoro, K. Mitarai, K. Fujii, K. Nakajima, and M. Kitagawa, Machine learning with controllable quantum dynamics of a nuclear spin ensemble in a solid, arXiv:1806.10910 (2018).
  * Rabl  _et al._ [2010a] P. Rabl, S. J. Kolkowitz, F. H. L. Koppens, J. G. E. Harris, P. Zoller, and M. D. Lukin, A quantum spin transducer based on nanoelectromechanical resonator arrays, Nat. Phys. 6, 602 (2010a).
  * Bhaktavatsala Rao  _et al._ [2011] D. D. Bhaktavatsala Rao, N. Bar-Gill, and G. Kurizki, Generation of macroscopic superpositions of quantum states by linear coupling to a bath, Phys. Rev. Lett. 106, 010404 (2011).
  * Bochmann  _et al._ [2013] J. Bochmann, A. Vainsencher, D. D. Awschalom, and A. N. Cleland, Nanomechanical coupling between microwave and optical photons, Nat. Phys. 9, 712 (2013).
  * Rosenfeld  _et al._ [2021] E. Rosenfeld, R. Riedinger, J. Gieseler, M. Schuetz, and M. D. Lukin, Efficient entanglement of spin qubits mediated by a hot mechanical oscillator, Phys. Rev. Lett. 126, 250505 (2021).
  * Fukami  _et al._ [2024] M. Fukami, J. C. Marcks, D. R. Candido, L. R. Weiss, B. Soloway, S. E. Sullivan, N. Delegan, F. J. Heremans, M. E. Flatté, and D. D. Awschalom, Magnon-mediated qubit coupling determined via dissipation measurements, Proc. Natl. Acad. Sci. 121, 10.1073/pnas.2313754120 (2024).
  * Wang and Terhal [2021] Y. Wang and B. M. Terhal, Preparing dicke states in a spin ensemble using phase estimation, Phys. Rev. A 104, 032407 (2021).
  * Hakoshima and Matsuzaki [2020] H. Hakoshima and Y. Matsuzaki, Efficient detection of inhomogeneous magnetic fields from a single spin with dicke states, Phys. Rev. A 102, 042610 (2020).
  * Peng  _et al._ [2024] R. Peng, X. Wu, Y. Wang, J. Zhang, J. Geng, D. Dasari, A. N. Cleland, and J. Wrachtrup, Spin-phonon entanglement in sic optomechanical quantum oscillators (2024), manuscript in preparation.
  * Dijkema  _et al._ [2023] J. Dijkema, X. Xue, P. Harvey-Collard, M. Rimbach-Russ, S. L. de Snoo, G. Zheng, A. Sammak, G. Scappucci, and L. M. K. Vandersypen, Two-qubit logic between distant spins in silicon, arXiv:2310.16805 (2023).
  * Ansel  _et al._ [2023] Q. Ansel, A. D. Chepelianskii, and J. Lages, Enhancing quantum exchanges between two oscillators, Phys. Rev. A 107, 042609 (2023).
  * Hartmann and Hahn [1962] S. R. Hartmann and E. L. Hahn, Nuclear double resonance in the rotating frame, Phys. Rev. 128, 2042 (1962).
  * Rao  _et al._ [2020] D. D. B. Rao, A. Ghosh, D. Gelbwaser-Klimovsky, N. Bar-Gill, and G. Kurizki, Spin-bath polarization via disentanglement, New J. Phys. 22, 083035 (2020).
  * Sasaki and Abe [2024] K. Sasaki and E. Abe, Suppression of pulsed dynamic nuclear polarization by many-body spin dynamics, Phys. Rev. Lett. 132, 106904 (2024).
  * Yang  _et al._ [2016] S. Yang, Y. Wang, D. D. B. Rao, T. Hien Tran, A. S. Momenzadeh, M. Markham, D. J. Twitchen, P. Wang, W. Yang, R. Stöhr, P. Neumann, H. Kosaka, and J. Wrachtrup, High-fidelity transfer and storage of photon states in a single nuclear spin, Nat. Photonics 10, 507 (2016).
  * Browaeys and Lahaye [2020] A. Browaeys and T. Lahaye, Many-body physics with individually controlled rydberg atoms, Nat. Phys. 16, 132 (2020).
  * Aspelmeyer  _et al._ [2014] M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, Cavity optomechanics, Rev. Mod. Phys. 86, 1391 (2014).
  * Rabl  _et al._ [2010b] P. Rabl, S. J. Kolkowitz, F. H. L. Koppens, J. G. E. Harris, P. Zoller, and M. D. Lukin, A quantum spin transducer based on nanoelectromechanical resonator arrays, Nat. Phys. 6, 602 (2010b).
  * Bennett  _et al._ [2013a] S. D. Bennett, N. Y. Yao, J. Otterbach, P. Zoller, P. Rabl, and M. D. Lukin, Phonon-induced spin-spin interactions in diamond nanostructures: Application to spin squeezing, Phys. Rev. Lett. 110, 156402 (2013a).
  * Barfuss  _et al._ [2015] A. Barfuss, J. Teissier, E. Neu, A. Nunnenkamp, and P. Maletinsky, Strong mechanical driving of a single electron spin, Nat. Phys. 11, 820 (2015).
  * Robledo  _et al._ [2011] L. Robledo, L. Childress, H. Bernien, B. Hensen, P. F. A. Alkemade, and R. Hanson, High-fidelity projective read-out of a solid-state spin quantum register, Nature 477, 574 (2011).
  * Aslam  _et al._ [2017] N. Aslam, M. Pfender, P. Neumann, R. Reuter, A. Zappe, F. Fávaro de Oliveira, A. Denisenko, H. Sumiya, S. Onoda, J. Isoya, and J. Wrachtrup, Nanoscale nuclear magnetic resonance with chemical resolution, Science 357, 67 (2017).
  * Chikkaraddy  _et al._ [2017] R. Chikkaraddy, V. A. Turek, N. Kongsuwan, F. Benz, C. Carnegie, T. van de Goor, B. de Nijs, A. Demetriadou, O. Hess, U. F. Keyser, and J. J. Baumberg, Mapping nanoscale hotspots with single-molecule emitters assembled into plasmonic nanocavities using dna origami, Nano Lett. 18, 405 (2017).
  * Bennett  _et al._ [2013b] S. D. Bennett, N. Y. Yao, J. Otterbach, P. Zoller, P. Rabl, and M. D. Lukin, Phonon-induced spin-spin interactions in diamond nanostructures: Application to spin squeezing, Phys. Rev. Lett. 110, 156402 (2013b).
  * Whiteley  _et al._ [2019] S. J. Whiteley, G. Wolfowicz, C. P. Anderson, A. Bourassa, H. Ma, M. Ye, G. Koolstra, K. J. Satzinger, M. V. Holt, F. J. Heremans, A. N. Cleland, D. I. Schuster, G. Galli, and D. D. Awschalom, Spin–phonon interactions in silicon carbide addressed by gaussian acoustics, Nature Physics 15, 490–495 (2019).
  * Huillery  _et al._ [2020] P. Huillery, T. Delord, L. Nicolas, M. Van Den Bossche, M. Perdriat, and G. Hétet, Spin mechanics with levitating ferromagnetic particles, Phys. Rev. B 101, 134415 (2020).
  * Schomburg [2011] W. K. Schomburg,  _Introduction to Microsystem Design_ (Springer Berlin Heidelberg, 2011).
  * Fung  _et al._ [2023] F. Fung, E. Rosenfeld, J. D. Schaefer, A. Kabcenell, J. Gieseler, T. X. Zhou, T. Madhavan, N. Aslam, A. Yacoby, and M. D. Lukin, Programmable quantum processors based on spin qubits with mechanically-mediated interactions and transport (2023).
  * Agarwal [2010] G. S. Agarwal, Saving entanglement via a nonuniform sequence of $\pi$ pulses, Phys. Scr. 82, 038103 (2010).
  * Rao  _et al._ [2016] D. D. B. Rao, S. A. Momenzadeh, and J. Wrachtrup, Heralded control of mechanical motion by single spins, Phys. Rev. Lett. 117, 077203 (2016).
  * Uhrig [2007] G. S. Uhrig, Keeping a quantum bit alive by optimized pi-pulse sequences, Phys. Rev. Lett. 98, 100504 (2007).
  * Dasari  _et al._ [2022] D. B. R. Dasari, S. Yang, A. Chakrabarti, A. Finkler, G. Kurizki, and J. Wrachtrup, Anti-zeno purification of spin baths by quantum probe measurements, Nat. Commun 13, 10.1038/s41467-022-35045-3 (2022).
  * Gottesman  _et al._ [2001] D. Gottesman, A. Kitaev, and J. Preskill, Encoding a qubit in an oscillator, Phys. Rev. A 64, 012310 (2001).
  * Terhal and Weigand [2016] B. M. Terhal and D. Weigand, Encoding a qubit into a cavity mode in circuit QED using phase estimation, Phys. Rev. A 93, 012315 (2016).
  * Babinec  _et al._ [2010] T. M. Babinec, B. J. M. Hausmann, M. Khan, Y. Zhang, J. R. Maze, P. R. Hemmer, and M. Lončar, A diamond nanowire single-photon source, Nat. Nanotechnol 5, 195–199 (2010).
  * Duivenvoorden  _et al._ [2017] K. Duivenvoorden, B. M. Terhal, and D. Weigand, Single-mode displacement sensor, Phys. Rev. A 95, 012305 (2017).
  * Flühmann  _et al._ [2019] C. Flühmann, T. L. Nguyen, M. Marinelli, V. Negnevitsky, K. Mehta, and J. P. Home, Encoding a qubit in a trapped-ion mechanical oscillator, Nature 566, 513 (2019).
  * Bild  _et al._ [2023] M. Bild, M. Fadel, Y. Yang, U. von Lüpke, P. Martin, A. Bruno, and Y. Chu, Schrödinger cat states of a 16-microgram mechanical oscillator, Science 380, 274 (2023).
  * Wang [2017] Y. Wang,  _Quantum Error Correction with the GKP Code and Concatenation with Stabilizer Codes_ , Master’s thesis, RWTH-AACHEN UNIVERSITY (2017).
  * Vuillot  _et al._ [2019] C. Vuillot, H. Asasi, Y. Wang, L. P. Pryadko, and B. M. Terhal, Quantum error correction with the toric gottesman-kitaev-preskill code, Phys. Rev. A 99, 032344 (2019).
  * Tao  _et al._ [2014] Y. Tao, J. M. Boss, B. A. Moores, and C. L. Degen, Single-crystal diamond nanomechanical resonators with quality factors exceeding one million, Nat. Commun 5, 10.1038/ncomms4638 (2014).


# A measurement of the dark energy equation of state with 40 strong lenses

Natalie B. Hogg
Laboratoire Univers et Particules de Montpellier, Université de Montpellier,
CNRS, Montpellier, France E-mail<EMAIL_ADDRESS>
###### Abstract

Problems with the cosmological constant model of dark energy motivate the
investigation of alternative scenarios. I make the first measurement of the
dark energy equation of state using the hierarchical strong lensing time delay
likelihood provided by TDCOSMO. I find that the combination of seven TDCOSMO
lenses and 33 SLACS lenses is only able to provide an upper bound on the dark
energy equation of state, $w<-1.75$, which nevertheless implies the presence
of a phantom dark energy component. When the strong lensing time delay data is
combined with a collection of cosmic microwave background, baryon acoustic
oscillation and Type Ia supernova data, I find that the equation of state
remains consistent with a cosmological constant.

###### keywords:

dark energy – gravitational lensing: strong

††pubyear: 2023††pagerange: A measurement of the dark energy equation of state
with 40 strong lenses–References

##  1 Introduction

The better part of three decades has passed since one of the most profound
discoveries in modern cosmology was made: that the expansion rate of the
Universe is currently accelerating (Riess et al., 1998; Perlmutter et al.,
1999). In the standard model of cosmology, this acceleration is attributed to
the cosmological constant $\Lambda$ acting as a dark energy, with an equation
of state $w=P/\rho=-1$, where $P$ is the pressure and $\rho$ the energy
density of the dark energy. However, the true nature of dark energy, and in
particular whether its energy density is actually constant in time, remains a
subject of debate (Escamilla et al., 2023).

Various dark energy models have been proposed as alternatives to the
cosmological constant, most of which incorporate a dark energy whose density
changes over time (Copeland et al., 2006). A model-agnostic way of
investigating this dynamical kind of dark energy is to use a parameterisation
of the equation of state which allows for deviations from $w=-1$. A popular
choice is the Chevalier–Polarski–Linder (CPL) parameterisation (Chevallier &
Polarski, 2001; Linder, 2003),

$w(a)=w_{0}+(1-a)\,w_{a},$ (1)

a first-order Taylor expansion in the scale factor $a$. The parameters $w_{0}$
and $w_{a}$ can thus be constrained with data; deviations from $w_{0}=-1$ and
$w_{a}=0$ would indicate that the energy density of dark energy is evolving in
time.

Myriad observational data have been used to place constraints on the equation
of state of dark energy. For example, the Planck satellite’s measurements of
the cosmic microwave background (CMB) temperature anisotropies, E-mode
polarisation and CMB lensing, along with baryon acoustic oscillations (BAO)
measured by the 6dF Galaxy Survey and SDSS (Beutler et al., 2011; Ross et al.,
2015; Alam et al., 2017) and the Pantheon catalogue of Type Ia supernovae
(SNIa) (Scolnic et al., 2018), produced a measurement of $w=-1.028\pm 0.031$
(Aghanim et al., 2020b). A stated goal of Stage IV surveys such as Euclid is
to measure the dark energy equation of state with percent-level precision
(Amendola et al., 2018).

From a theoretical perspective, strong lensing time delays are a powerful
probe of cosmology due to their ability to provide cosmological constraints at
low redshift without reliance on the distance ladder. Strong lensing time
delays are the measurement of the different arrival times of the multiple
images of a gravitationally lensed source object (Refsdal, 1964). This time
delay depends on the angular diameter distances between the objects involved
and thus directly probes the expansion rate of the Universe, $H_{0}$.

Whilst there is no reliance on the distance ladder in strong lensing time
delay measurements of $H_{0}$, there can be a fairly significant dependence on
how the mass profile of the lens galaxy is modelled (Sonnenfeld, 2018). This
problem has been recognised since the studies of the very first strongly
lensed quasar for which time delays were measured, Q0957+561 (Walsh et al.,
1979; Falco et al., 1985).

As an example of the effect of different lens modelling on the measurement of
$H_{0}$, we can compare the H0LiCOW measurement of
$H_{0}=73.3^{+1.7}_{-1.8}$kms-1Mpc-1 (Wong et al., 2020), which comes from six
strong lensing time delays and is a 2.4% precision measurement, with the
TDCOSMO measurement, $H_{0}=74.5^{+5.6}_{-6.1}$kms-1Mpc-1 (Birrer et al.,
2020), which comes from the original six H0LiCOW lenses plus one new lens and
is a 9% precision measurement.

The reduction in precision is mainly the result of relaxing certain
assumptions in the H0LiCOW lens modelling. To combat this increase in
uncertainty, the TDCOSMO team added 33 additional lenses from the SLACS
catalogue in order to provide more information on the mass profiles of the
TDCOSMO lenses. This data and its nuisance parameters were combined with the
original seven lens dataset in a hierarchical Bayesian manner, leading to a
final measurement of $H_{0}=67.4^{+4.1}_{-3.2}$kms-1Mpc-1, which is in
statistical agreement with both the H0LiCOW measurement and the measurement
from the seven TDCOSMO lenses alone.

The TDCOSMO $H_{0}$ value has been used to constrain variations of the fine
structure constant (Colaço et al., 2021), axion–photon couplings (Buen-Abad et
al., 2022) and interacting dark energy models (Wang et al., 2022), but the
full likelihood has not yet been used to obtain cosmological constraints
beyond those on $H_{0}$ and the matter density parameter $\Omega_{\rm m}$.

In this work, I present the first measurement of the dark energy equation of
state using the full hierarchical TDCOSMO likelihood. Strong lensing time
delays are particularly useful for this task as the constraints in the
$w-\Omega_{\rm m}$ plane are typically orthogonal to those coming from
standard rulers i.e. the CMB and BAO (Motta et al., 2021). I begin by
reviewing strong lensing time delays and the construction of the hierarchical
TDCOSMO likelihood. I explain my analysis method and then present my results
and conclusions.

##  2 Strong lensing time delays

###  2.1 Theory

Gravitational lensing is the phenomenon which arises due to the deflection of
light by massive objects. The strong lensing regime may be defined as that in
which multiple images of a single source are produced. This typically occurs
on super-galactic scales, with lenses being galaxies and sources being distant
and bright objects such as quasars.

Light from the distant source is deflected by the lens galaxy. Different light
paths have different lengths, leading to measurable delays between the arrival
times of images. This strong lensing time delay is given by

$t(\bm{\theta},\bm{\beta})=\frac{(1+z_{\rm od})}{c}\frac{D_{\rm od}D_{\rm
os}}{D_{\rm
ds}}\left[\frac{(\bm{\theta}-\bm{\beta})^{2}}{2}-\psi(\bm{\theta})\right],$
(2)

where $z_{\rm od}$ is the redshift of the deflector; $D_{\rm od}$, $D_{\rm
os}$ and $D_{\rm ds}$ are the angular diameter distances between the observer
and deflector, observer and source and deflector and source; $\bm{\theta}$ is
the observed image position; $\bm{\beta}$ is the unknown source position; and
$\psi(\bm{\theta})$ is the lensing potential, which carries the information
about the mass density profile of the deflector. In a spatially flat Universe
($\Omega_{k}=0$), the time delay is inversely proportional to $H_{0}$ via the
angular diameter distances involved,

$D(z)=\frac{c}{H_{0}(1+z)}\int_{0}^{z}\frac{\mathrm{d}z^{\prime}}{E(z^{\prime})},$
(3)

where $E(z)\equiv H(z)/H_{0}$ is the dimensionless Hubble rate.

The difficulty faced by all strong lensing time delay inference is that under
any arbitrary linear re-scaling of the source position
$\bm{\beta}\rightarrow\lambda\bm{\beta}$, image positions $\bm{\theta}$ are
preserved. The lens model is also accordingly transformed,
$\bm{\alpha}\rightarrow\lambda\bm{\alpha}+(1-\lambda)\bm{\theta}$, where
$\bm{\alpha}$ is the deflection angle. This is known as the internal
mass–sheet transform or degeneracy (Falco et al., 1985). The only way that
such a degeneracy can be broken and a measurement of $H_{0}$ made is by
obtaining direct knowledge either of the absolute source size or of the
lensing potential itself. The former is quite obviously near-impossible, so in
order to constrain $H_{0}$ a choice must be made in how to model the deflector
mass density profile.

###  2.2 The TDCOSMO likelihood

Whilst the H0LiCOW team used analytic mass profiles for the deflectors to
break the mass-sheet degeneracy, the TDCOSMO work used stellar kinematics data
to provide information about the lensing potential. This led to the initial
reduction in precision of the $H_{0}$ measurement compared to the H0LiCOW
result. The precision was increased again by the addition of further stellar
kinematics data from a set of 33 SLACS lenses, specifically selected for their
similarity to the original seven TDCOSMO lenses. Note that these additional
lenses do not have time delay information. Furthermore, the combination of
data was made under the assumption that the TDCOSMO and SLACS lenses were
drawn from the same parent population.

The complete likelihood describing this dataset was constructed
hierarchically, meaning that a set of hyperparameters were defined which
allowed all constraints related to the mass-sheet degeneracy to be inferred on
a population level, whilst the lens and light model parameters, $\bm{\xi}_{\rm
mass}$ and $\bm{\xi}_{\rm light}$, could be inferred on a lens-by-lens basis.
Thus, all remaining uncertainty about the mass-sheet degeneracy is propagated
to the level of the $H_{0}$ inference.

Following Birrer et al. (2020), the posterior distribution of the cosmological
parameters of interest, $\bm{\pi}$, given $N$ sets of individual lens data
$\mathcal{D}_{i}$ and the model parameters $\bm{\xi}$, is given by

$\displaystyle P(\bm{\pi}|\\{\mathcal{D}_{i}\\}_{N})$
$\displaystyle\propto\mathcal{L}(\\{\mathcal{D}_{i}\\}_{N}|\bm{\pi})P(\bm{\pi}),$
$\displaystyle=\int\mathcal{L}(\\{\mathcal{D}_{i}\\}_{N}|\bm{\pi},\bm{\xi})P(\bm{\pi},\bm{\xi}),$
$\displaystyle=\int\prod_{i}^{N}\mathcal{L}(\mathcal{D}_{i}|\bm{\pi},\bm{\xi})P(\bm{\pi},\bm{\xi}).$
(4)

The nuisance parameter $\bm{\xi}$ is divided into the mass and light model
parameters which are constrained at the level of each individual lens, and the
set of mass-sheet degeneracy hyperparameters which are constrained at the
population level, $\bm{\xi}_{\rm pop}$. The hierarchical TDCOSMO likelihood is
thus given by

$\displaystyle\mathcal{L}(\mathcal{D}_{i}|D,\bm{\xi}_{\rm pop})=\int$
$\displaystyle\mathcal{L}(\mathcal{D}_{i}|D,\bm{\xi}_{\rm pop},\bm{\xi}_{\rm
mass},\bm{\xi}_{\rm light})$ $\displaystyle\times P(\bm{\xi}_{\rm
mass},\bm{\xi}_{\rm light})\;\mathrm{d}\bm{\xi}_{\rm
mass}\,\mathrm{d}\bm{\xi}_{\rm light},$ (5)

where $D$ is the set of angular diameter distances $\\{D_{\rm od},D_{\rm
os},D_{\rm ds}\\}$ from which all cosmological results are obtained. For a
complete discussion of the construction of the hierarchical likelihood,
including the details of the population hyperparameters, I refer the reader to
Section 3 of Birrer et al. (2020).

The final TDCOSMO measurement of $H_{0}$ is thus the most precise measurement
possible from strong lensing time delays which does not make any assumptions
about the deflector mass profiles in order to artificially break the mass-
sheet degeneracy. It is also the first strong lensing time delay measurement
of $H_{0}$ which used information from other datasets to improve the precision
of the measurement. I will now discuss how I used this hierarchical likelihood
to measure the dark energy equation of state.

##  3 Method

I wrote an external likelihood package for the cosmological modelling and
sampling software Cobaya (Torrado & Lewis, 2021), so that the hierarchical
TDCOSMO likelihood can be used to obtain constraints on cosmological model
parameters in combination with any other cosmological likelihood and with any
choice of Boltzmann code. This package is publicly
available.111https://github.com/nataliehogg/tdcosmo_ext.

Using the Markov chain Monte Carlo sampler provided by Cobaya, which is
adapted from CosmoMC (Lewis & Bridle, 2002; Lewis, 2013) and uses a fast-
dragging procedure to increase sampling speed (Neal, 2005), I obtained
constraints on three cosmological models: $\Lambda$CDM, $w$CDM and
$w_{0}w_{a}$CDM. The $w$CDM model allows for a dark energy with a constant
equation of state that may differ from $w=-1$; the $w_{0}w_{a}$CDM model
allows for a dynamical dark energy. In a $w$CDM cosmology (extendable to
$w_{0}w_{a}$CDM by the CPL parameterisation shown in Equation 1), and
recalling that I keep $\Omega_{k}=0$, the dimensionless Hubble rate is given
by

$E(z)=\left[\Omega_{\rm m}(1+z)^{3}+\Omega_{\rm
DE}(1+z)^{3(1+w)}\right]^{\frac{1}{2}},$ (6)

where $\Omega_{\rm m}=\Omega_{\rm b}+\Omega_{\rm c}$, the sum of the baryon
and cold dark matter densities, and $\Omega_{\rm DE}$ is the dimensionless
energy density of dark energy today.

For each cosmological model I considered, I sampled the posterior
distributions of $H_{0}$, $\Omega_{\rm b}h^{2}$ and $\Omega_{\rm c}h^{2}$,
along with the relevant dark energy equation of state parameter(s). In each
case, I also sampled the posterior distributions of the hyperparameters
associated with the likelihood, and marginalised over them to obtain the
posterior distributions on the cosmological parameters of interest. The
computation of the angular diameter distances required by the TDCOSMO
likelihood was done using the Boltzmann code CAMB (Lewis et al., 2000; Howlett
et al., 2012), but I emphasise that I designed the Cobaya interface so that
any theory code currently available in or added to Cobaya in the future can be
used for this task. Lastly, I used the Parameterised Post-Friedmann framework
in CAMB to allow $w$ to cross $-1$ (Fang et al., 2008).

I validated my implementation of the likelihood package in Cobaya by comparing
my $\Lambda$CDM results with the original TDCOSMO results, finding an
excellent agreement in terms of constraints on both the cosmological and the
nuisance parameters. This validation test, plus all of the code needed to
reproduce the results and figures in this paper is also publicly
available.222https://github.com/nataliehogg/slide.

Besides the TDCOSMO and SLACS datasets, I also obtained constraints on the
models with a dataset consisting of the Planck 2018 measurements of the CMB
temperature, polarisation and lensing (Aghanim et al., 2020a, c); the BAO
measurements from the 6dF Galaxy Survey (Beutler et al., 2011), the SDSS Main
Galaxy Sample (Ross et al., 2015) and the SDSS DR12 consensus catalogue (Alam
et al., 2017); the Pantheon catalogue of Type Ia supernovae (Scolnic et al.,
2018); and the TDCOSMO + SLACS lenses, in order to compare with the
constraints obtained just using the strong lensing time delay data and with
those of the Planck collaboration (Aghanim et al., 2020b). I will refer to
this combination of data as the “full combination” from now on.

The priors I used for the cosmological parameters are listed in Table 1.
Following Birrer et al. (2020), I also used the Pantheon prior on $\Omega_{\rm
m}=\mathcal{N}(0.298,0.022)$ when using the TDCOSMO and SLACS data alone,
though it is important to note that in my analysis $\Omega_{\rm m}$ was not
directly sampled since it is treated as a derived parameter in CAMB, the
posterior being obtained from those of $\Omega_{\rm b}$ and $\Omega_{\rm c}$.
In the dynamical dark energy cases, I ensured that acceleration occurs by
setting a prior on the dark energy equation of state such that
$w<-\frac{1}{3}$.

Parameter | Prior
---|---
$\Omega_{b}h^{2}$ | $[0.005,0.1]$
$\Omega_{c}h^{2}$ | $[0.001,0.99]$
$H_{0}$ | $[0,150]$
$w$ and $w_{0}$ | $\left[-3.0,-\frac{1}{3}\right]$
$w_{a}$ | $[-4.0,4.0]$
Table 1: Prior ranges of the parameters sampled in my analysis.

##  4 Results

In this section, I present the constraints obtained on $\Lambda$CDM, $w$CDM
and $w_{0}w_{a}$CDM using the hierarchical TDCOSMO likelihood and the full
combination of data listed above. Given the smaller uncertainty on $H_{0}$
that comes from combining the TDCOSMO + SLACS lenses, I will only present the
TDCOSMO alone results for $\Lambda$CDM, to demonstrate the replication of
Birrer et al. (2020). For the rest of the results, I will show constraints
obtained using the complete TDCOSMO + SLACS dataset of 40 lenses. I used
GetDist to make the figures and to compute the marginalised parameter values
quoted in the text (Lewis, 2019).

###  4.1 $\Lambda$CDM

Figure 1: The one and two-dimensional marginalised posterior distributions of
$H_{0}$ and $\Omega_{\rm m}$ in a $\Lambda$CDM cosmology from the seven
TDCOSMO lenses (blue) and the 40 TDCOSMO + SLACS lenses (red).

In Figure 1, I show the one and two-dimensional marginalised posterior
distributions of $H_{0}$ and $\Omega_{\rm m}$ in a $\Lambda$CDM cosmology. The
blue contours show the constraints from the seven TDCOSMO lenses alone, whilst
the red contours show the constraints from the full dataset of 40 lenses
(TDCOSMO + SLACS). As expected, this result replicates the findings of the
TDCOSMO paper, with $H_{0}=76.8^{+6.4}_{-5.6}$kms-1Mpc-1 from the TDCOSMO
lenses and $H_{0}=68.7^{+3.4}_{-3.9}$ from the TDCOSMO + SLACS lenses. These
values are fully consistent at $1\sigma$ with the TDCOSMO result.

###  4.2 $w$CDM

Figure 2: The one and two-dimensional marginalised posterior distributions of
$H_{0}$, $w$ and $\Omega_{\rm m}$ in a $w$CDM cosmology from the 40 TDCOSMO +
SLACS lenses.

In Figure 2, I show the one and two-dimensional marginalised posterior
distributions of $H_{0}$, $w$ and $\Omega_{\rm m}$ in a $w$CDM cosmology from
the 40 TDCOSMO + SLACS lenses. The likelihood is only able to provide an upper
bound on the equation of state of dark energy, $w<-1.75$. This bound implies
that the dark energy is ‘phantom’, the term applied when $w<-1$.

Some consideration must be given here to what is, at first glance, a
surprising result. It is known that in some datasets, such as Planck 2018,
large negative values of $w$ correlate with large positive values of $H_{0}$
(Escamilla et al., 2023). This is due to the so-called geometrical degeneracy,
where $H_{0}$, $\Omega_{\rm m}$ and $w$ can take various values which in
combination lead to the same value for the angular diameter distance to the
surface of last scattering and hence the same angular size of the sound
horizon, provided the physical sound horizon size is kept fixed (Efstathiou &
Bond, 1999).

Since strong lensing time delays also rely on angular diameter distances to
probe cosmology, I infer that a similar degeneracy exists here. This
conclusion is supported by the clear correlation between $H_{0}$ and $w$ in
Figure 2. Thus, the strongly negative dark energy equation of state is likely
driven by the high central value of $H_{0}=78.4^{+8.3}_{-6.3}$ obtained in
this case – which is nevertheless consistent with the $\Lambda$CDM value at
$1\sigma$. As expected, based on the results of Birrer et al. (2020), the
marginalised posterior distribution of the matter density $\Omega_{\rm m}$ is
largely informed by the Pantheon prior. This also acts to somewhat ameliorate
the aforementioned degeneracy.

###  4.3 $w_{0}w_{a}$CDM

Figure 3: The one and two-dimensional marginalised posterior distributions of
$H_{0}$, $w_{0}$, $w_{a}$ and $\Omega_{\rm m}$ in a $w_{0}w_{a}$CDM cosmology
from the 40 TDCOSMO + SLACS lenses.

In Figure 3, I show the one and two-dimensional marginalised posterior
distributions of $H_{0}$, $w_{0}$, $w_{a}$ and $\Omega_{\rm m}$ in a
$w_{0}w_{a}$CDM cosmology from the 40 TDCOSMO + SLACS lenses. The likelihood
is again only able to provide upper bounds on the dark energy equation of
state parameters, $w_{0}<-1.86$ and $w_{a}<0.102$. Again the value of $H_{0}$
is larger than but still consistent with the $\Lambda$CDM measurement:
$H_{0}=79.6^{+7.5}_{-6.0}$.

Figure 4: The two-dimensional marginalised posterior distributions of $H_{0}$
and $\Omega_{\rm m}$ in the three cosmologies studied from the 40 TDCOSMO +
SLACS lenses.

In Figure 4, I show the two-dimensional marginalised posterior distributions
of $H_{0}$ and $\Omega_{\rm m}$ in a $\Lambda$CDM cosmology (red), a $w$CDM
cosmology (orange) and a $w_{0}w_{a}$CDM cosmology (yellow) from the 40
TDCOSMO + SLACS lenses. From this plot, it is clear that whilst the $H_{0}$
values in the extended cosmologies are large, they are still consistent at
$1\sigma$ with the $\Lambda$CDM result.

###  4.4 Full combination of data

Figure 5: The two-dimensional marginalised posterior distributions of $H_{0}$
and $\Omega_{\rm m}$ in the three cosmologies studied from the full
combination of Planck 2018 + BAO + Pantheon SNIa + TDCOSMO + SLACS.

In Figure 5, I show the two-dimensional marginalised posterior distributions
of $H_{0}$ and $\Omega_{\rm m}$ in a $\Lambda$CDM cosmology (dark purple), a
$w$CDM cosmology (medium purple) and a $w_{0}w_{a}$CDM cosmology (light
purple) from the full combination of Planck 2018 + BAO + Pantheon SNIa +
TDCOSMO + SLACS data. In this case, I did not use the Pantheon prior on
$\Omega_{\rm m}$, since the inclusion of the Pantheon dataset provides the
same information as that prior. From this plot we can see that the values of
$H_{0}$ and $\Omega_{\rm m}$ in the extended cosmologies are completely
consistent with the $\Lambda$CDM values at $1\sigma$ when using the full
combination of data; the $w$CDM and $w_{0}w_{a}$CDM constraints are virtually
identical.

Furthermore, this combination of data inevitably provides a much stronger
constraint on $w$, $w_{0}$ and $w_{a}$ than the TDCOSMO + SLACS data alone,
and removes any hint of a phantom dark energy component. In a $w$CDM
cosmology, the dark energy equation of state is measured to be $w=-1.025\pm
0.029$, a marginal increase in precision compared to the Planck 2018 + BAO +
SNIa value of $w=-1.028\pm 0.031$. Both of these measurements are consistent
with a cosmological constant.

Figure 6: The two-dimensional marginalised posterior distributions of $w_{0}$
and $w_{a}$ from the TDCOSMO + SLACS data (red) and from the full combination
of Planck 2018 + BAO + Pantheon SNIa + TDCOSMO + SLACS (purple). The dashed
lines show the values of these parameters which correspond to a cosmological
constant, $w_{0}=-1$, $w_{a}=0$.

In a $w_{0}w_{a}$CDM cosmology, $w_{0}=-0.985^{+0.071}_{-0.091}$ and
$w_{a}=-0.18^{+0.33}_{-0.25}$, compared to $w_{0}=-0.957\pm 0.08$ and
$w_{a}=-0.29^{+0.32}_{-0.26}$ from Planck 2018 + BAO + SNIa. These are again
both almost equivalent in precision and consistent with a cosmological
constant. This can be clearly seen in Figure 6, where I show the two-
dimensional marginalised posterior distributions of $w_{0}$ and $w_{a}$ from
the TDCOSMO + SLACS data (red) and from the full combination of Planck 2018 +
BAO + Pantheon SNIa + TDCOSMO + SLACS (purple). The dashed lines in this plot
show the values of the dark energy equation of state parameters which
correspond to a cosmological constant, $w_{0}=-1$ and $w_{a}=0$.

Model | Data | $\Delta\chi^{2}$
---|---|---
$w$CDM | TDCOSMO + SLACS | $-0.5$
$w$CDM | Full combination | $-0.7$
$w_{0}w_{a}$CDM | TDCOSMO + SLACS |  $0.0$
$w_{0}w_{a}$CDM | Full combination | $-2.0$
Table 2: The $\Delta\chi^{2}$ for each case studied.

Lastly, I computed the
$\Delta\chi^{2}\equiv\chi^{2}-\chi^{2}_{\Lambda\rm{CDM}}$ for each case
studied. The computed values are shown in Table 2. Whilst the $\Delta\chi^{2}$
is negative for the $w$CDM cosmology obtained with both datasets, implying
that the extended cosmology is favoured over $\Lambda$CDM, the significance of
this decrease in $\chi^{2}$ must be evaluated using a difference table due to
the greater number of degrees of freedom in the $w$CDM cosmology with respect
to $\Lambda$CDM. With one additional degree of freedom, and for a 95% level of
significance, a $|\Delta\chi^{2}|>3.841$ is required for the improvement in
fit to be considered significant. This is clearly not the case here. From
these results it is also evident that $\Lambda$CDM is a better fit to the data
than the $w_{0}w_{a}$CDM cosmologies, since the $\Delta\chi^{2}$ is on parity
with $\Lambda$CDM for the TDCOSMO + SLACS data and still does not exceed
$3.841$ for the full combination of data.

##  5 Conclusions

In this work, I presented the first constraints on the equation of state of
dark energy from the seven TDCOSMO lenses plus 33 SLACS lenses using the
hierarchical likelihood provided by TDCOSMO. I wrote an external likelihood
package for the Cobaya software, making this likelihood readily available for
public use in combination with other cosmological likelihoods and theory
codes.

I replicated the original TDCOSMO results in a $\Lambda$CDM cosmology and then
explored two extended cosmologies, finding that the TDCOSMO + SLACS data was
not able to place strong constraints on the equation of state of dark energy,
obtaining only upper bounds on $w$, $w_{0}$ and $w_{a}$. The bounds I obtained
all implied the presence of a phantom dark energy component, $w<-1$.

The use of the TDCOSMO likelihood in combination with the Planck 2018
likelihood, BAO data and the Pantheon SNIa catalogue yielded more precise
constraints, with all the dark energy equation of state parameters consistent
with a cosmological constant, $w=-1$. I computed the $\Delta\chi^{2}$ to
evaluate the fit of each model, finding no preference for the extended
cosmologies over the $\Lambda$CDM case.

In conclusion, while strong lensing time delays are beginning to provide a
competitive (albeit lens-modelling-dependent) measurement of $H_{0}$, it is
clear that a larger dataset is needed before they can be useful when studying
dark energy; perhaps on the order of hundreds or thousands of lenses
(Shiralilou et al., 2020). Furthermore, an increase in precision on the
$H_{0}$ inference will naturally lead to a reduction in the effect of the
geometrical degeneracy which is contributing to the hints of phantom dark
energy in the TDCOSMO data. Fortunately a number of current or near-future
experiments, such as JWST, Roman, LSST and Euclid are likely to provide such
data in abundance (Collett, 2015).

## Acknowledgements

I am grateful to Simon Birrer, Martin Millon and Judit Prat for valuable
discussions about the TDCOSMO likelihood, and to Pierre Fleury for his
comments on the manuscript.

## Data Availability

All data associated with this article is publicly available.

## References

  * Aghanim et al. (2020a) Aghanim N., et al., 2020a, Astronomy & Astrophysics, 641, A5
  * Aghanim et al. (2020b) Aghanim N., et al., 2020b, Astronomy & Astrophysics, 641, A6
  * Aghanim et al. (2020c) Aghanim N., et al., 2020c, Astronomy & Astrophysics, 641, A8
  * Alam et al. (2017) Alam S., et al., 2017, Monthly Notices of the Royal Astronomical Society, 470, 2617
  * Amendola et al. (2018) Amendola L., et al., 2018, Living Reviews in Relativity, 21, 2
  * Beutler et al. (2011) Beutler F., et al., 2011, Monthly Notices of the Royal Astronomical Society, 416, 3017
  * Birrer et al. (2020) Birrer S., et al., 2020, Astronomy & Astrophysics, 643, A165
  * Buen-Abad et al. (2022) Buen-Abad M. A., Fan J., Sun C., 2022, Journal of High Energy Physics, 02, 103
  * Chevallier & Polarski (2001) Chevallier M., Polarski D., 2001, International Journal of Modern Physics D, 10, 213
  * Colaço et al. (2021) Colaço L. R., Gonzalez J. E., Holanda R. F. L., 2021, The European Physical Journal C, 81, 533
  * Collett (2015) Collett T. E., 2015, The Astrophysical Journal, 811, 20
  * Copeland et al. (2006) Copeland E. J., Sami M., Tsujikawa S., 2006, International Journal of Modern Physics D, 15, 1753
  * Efstathiou & Bond (1999) Efstathiou G., Bond J. R., 1999, Monthly Notices of the Royal Astronomical Society, 304, 75
  * Escamilla et al. (2023) Escamilla L. A., Giarè W., Di Valentino E., Nunes R. C., Vagnozzi S., 2023 (arXiv:2307.14802)
  * Falco et al. (1985) Falco E. E., Gorenstein M. V., Shapiro I. I., 1985, The Astrophysical Journal Letters, 289, L1
  * Fang et al. (2008) Fang W., Hu W., Lewis A., 2008, Physical Review D, 78, 087303
  * Howlett et al. (2012) Howlett C., Lewis A., Hall A., Challinor A., 2012, Journal of Cosmology and Astroparticle Physics, 1204, 027
  * Lewis (2013) Lewis A., 2013, Physical Review D, 87, 103529
  * Lewis (2019) Lewis A., 2019 (arXiv:1910.13970)
  * Lewis & Bridle (2002) Lewis A., Bridle S., 2002, Physical Review D, 66, 103511
  * Lewis et al. (2000) Lewis A., Challinor A., Lasenby A., 2000, The Astrophysical Journal, 538, 473
  * Linder (2003) Linder E. V., 2003, Physical Review Letters, 90, 091301
  * Motta et al. (2021) Motta V., García-Aspeitia M. A., Hernández-Almada A., Magaña J., Verdugo T., 2021, Universe, 7, 163
  * Neal (2005) Neal R. M., 2005 (arXiv:math/0502099)
  * Perlmutter et al. (1999) Perlmutter S., et al., 1999, The Astrophysical Journal, 517, 565
  * Refsdal (1964) Refsdal S., 1964, Monthly Notices of the Royal Astronomical Society, 128, 307
  * Riess et al. (1998) Riess A. G., et al., 1998, The Astronomical Journal, 116, 1009
  * Ross et al. (2015) Ross A. J., Samushia L., Howlett C., Percival W. J., Burden A., Manera M., 2015, Monthly Notices of the Royal Astronomical Society, 449, 835
  * Scolnic et al. (2018) Scolnic D. M., et al., 2018, The Astrophysical Journal, 859, 101
  * Shiralilou et al. (2020) Shiralilou B., Martinelli M., Papadomanolakis G., Peirone S., Renzi F., Silvestri A., 2020, Journal of Cosmology and Astroparticle Physics, 04, 057
  * Sonnenfeld (2018) Sonnenfeld A., 2018, Monthly Notices of the Royal Astronomical Society, 474, 4648
  * Torrado & Lewis (2021) Torrado J., Lewis A., 2021, Journal of Cosmology and Astroparticle Physics, 05, 057
  * Walsh et al. (1979) Walsh D., Carswell R. F., Weymann R. J., 1979, Nature, 279, 381
  * Wang et al. (2022) Wang L.-F., Zhang J.-H., He D.-Z., Zhang J.-F., Zhang X., 2022, Monthly Notices of the Royal Astronomical Society, 514, 1433
  * Wong et al. (2020) Wong K. C., et al., 2020, Monthly Notices of the Royal Astronomical Society, 498