Preparation of Metrological States in Dipolar-Interacting Spin Systems

Consider a two-body interaction Hamiltonian:

$\displaystyle H_{\text{int}}=\sum_{i<j}V_{ij}\left(J^{\text{I}}S_{zi}S_{zj}+J^{\text{S}}\bm{S}_{i}\cdot\bm{S}_{j}\right).$

(S1)

In this Hamiltonian, $\bm{S}=(S_{x},S_{y},S_{z})$ is the vector of spin-1/2 operators, $V_{ij}$ is the interaction strength between spin $i$ and $j$ which depends on their locations, and $J^{\text{I}}(\neq 0),J^{\text{S}}$ are the Ising and symmetric coupling constant respectively.

The elementary gates in each layer of the variational circuit for preparing metrological states (Fig.1(c) main text) include two free evolutions under the interaction Hamiltonian $D(\tau),D(\tau^{\prime})$, one global rotation along the $x$-axis $R_{x}(\vartheta)$ and two fixed $\pi/2$ rotations $R_{y}(-\frac{\pi}{2}),R_{y}(\frac{\pi}{2})$ along the $y$-axis. We define the interaction gate in the $z$-direction as

$\displaystyle D_{z}(\tau)\equiv\exp(-i\tau H_{\text{int}}/\hbar)=\exp[-i\tau\sum_{i<j}V_{ij}\left(J^{\text{I}}S_{zi}S_{zj}+J^{\text{S}}\bm{S}_{i}\cdot\bm{S}_{j}\right)/\hbar].$

(S2)

The interaction gates in other directions can be obtained by $\pi/2$ rotations:

	$\displaystyle D_{x,y}(\tau)$	$\displaystyle=R_{y,x}(\pi/2)D_{z}(\tau)R_{y,x}(-\pi/2)$
		$\displaystyle=\exp[-i\tau\sum_{i<j}V_{ij}\left(J^{\text{I}}S_{x,yi}S_{x,yj}+J^{\text{S}}\bm{S}_{i}\cdot\bm{S}_{j}\right)/\hbar].$		(S3)

In Eqs. (S2) (Entanglement generation gates from two-body interaction Hamiltonian and global rotations), the symmetric interaction term stay unchanged because inner product is conserved under global rotation and the ‘direction of interaction’ is only determined by the Ising term. Using these definitions, we simplify the gate set in each layer as

	$\displaystyle\mathcal{U}_{i}$	$\displaystyle=R_{y}\left(\frac{\pi}{2}\right)D\left(\tau^{\prime}_{i}\right)R_{y}\left(-\frac{\pi}{2}\right)R_{x}\left(\vartheta_{i}\right)D\left(\tau_{i}\right)$
		$\displaystyle=D_{x}\left(\tau^{\prime}_{i}\right)R_{x}\left(\vartheta_{i}\right)D_{z}\left(\tau_{i}\right).$		(S4)

In the next two subsections, it will be shown that the sequence in Eq. (Entanglement generation gates from two-body interaction Hamiltonian and global rotations) is the most general gate set that uses only global rotations and preserves the collective spin direction along $x$-direction.

Define the x-parity operator $P_{x}\equiv\Pi_{i}^{N}\sigma_{xi}=P_{x}^{\dagger}$, with $P_{x}^{2}=I$. This operator describes the parity of a state in $x$-direction and is related to the global $\pi$ rotation along $x$-axis up to a phase constant, $R_{x}(\pi)=\exp(-i\pi\sum_{i}S_{xi})=(-i)^{N}\Pi_{i}^{N}\sigma_{xi}$. Applying the x-parity operator onto individual spin’s angular momentum operator gives $P_{x}S_{\mu j}P_{x}=(\sigma_{x}S_{\mu}\sigma_{x})_{j}=\pm S_{\mu j}$. Thus the interaction gates along $x$- and $z$-direction conserve the x-parity, $P_{x}D_{x,z}P_{x}=D_{x,z}$. Similarly, the only global rotation that conserves x-parity for arbitrary angles is $R_{x}(\vartheta)$. Then, based on Eq.(1) in the main text, the unitary operator of the whole control sequence conserves the x-parity

	$\displaystyle P_{x}\mathcal{S}(\bm{\theta})P_{x}$	$\displaystyle=P_{x}\mathcal{U}_{m}...\mathcal{U}_{2}\mathcal{U}_{1}P_{x}$
		$\displaystyle=P_{x}\Pi_{i}[D_{x}\left(\tau^{\prime}_{i}\right)R_{x}\left(\vartheta_{i}\right)D_{z}\left(\tau_{i}\right)]P_{x}$
		$\displaystyle=\mathcal{S}(\bm{\theta}).$		(S5)

The initial spin state pointing to the $+x$-direction is an eigenstate of $P_{x}$: $P_{x}\ket{\uparrow_{x}}^{\otimes N}=\ket{\uparrow_{x}}^{\otimes N}$. Thus, any state produced by this variational circuit remains an eigenstate of $P_{x}$:

	$\displaystyle P_{x}\ket{\Psi(\bm{\theta})}$	$\displaystyle=P_{x}\mathcal{S}(\bm{\theta})\ket{\uparrow_{x}}^{\otimes N}$
		$\displaystyle=P_{x}\mathcal{S}(\bm{\theta})P_{x}P_{x}\ket{\uparrow_{x}}^{\otimes N}$
		$\displaystyle=\ket{\Psi(\bm{\theta})}.$		(S6)

Now consider the expectation value of the total spin angular momentum operator $J_{\mu}\equiv\sum_{i}S_{\mu i}$ $(\mu\in\{x,y,z\})$:

	$\displaystyle\langle J_{y,z}\rangle$	$\displaystyle=\bra{\Psi(\bm{\theta})}J_{y,z}\ket{\Psi(\bm{\theta})}$
		$\displaystyle=\bra{\Psi(\bm{\theta})}P_{x}P_{x}J_{y,z}P_{x}P_{x}\ket{\Psi(\bm{\theta})}$
		$\displaystyle=-\langle J_{y,z}\rangle=0.$		(S7)

Thus, the collective spin direction $\langle\bm{J}\rangle/|\langle\bm{J}\rangle|$ always points along the $x$-direction.

To preserve the collective spin direction along $x$-axis, the global rotation and interaction gates that can be chosen are $R_{x}$, $D_{x}$, $D_{\perp}$ where $D_{\perp}$ stands for the interaction gates along any direction perpendicular to the $x$-direction. Combining $R_{x}$ and $D_{z}$ can generate any $D_{\perp}$, thus the simplest gate set fulfilling all the requirements is $D_{x}R_{x}D_{z}$, as described by Eq.(1) in the main text.

The derivations and results in this section about selecting the variational sequence agree with ref.[Kaubruegger et al., 2019]. However, the interaction Hamiltonian we discuss here is more general. In Eq. (S1), when $J^{\text{I}}=1,J^{\text{S}}=0$, the interaction becomes Ising type interaction which is equivalent to the Rydberg interaction in ref.[Kaubruegger et al., 2019,Kaubruegger et al., 2021]. The Ising interaction can also describe spin systems with large local disorder. The optimization results are shown in the next section. When $J^{\text{I}}=3,J^{\text{S}}=-1$, Eq. (S1) becomes the dipolar interaction Hamiltonian between spin-1/2 particles. When $J^{\text{I}}=2,J^{\text{S}}=-1$, it becomes the dipolar interaction Hamiltonian between spin-1 particles (such as NV centers). The simulation results for this case are shown in the next section. When $J^{\text{I}}=1,J^{\text{S}}=-1$, the interaction can describe the dipolar interaction between cold molecules [Yan et al., 2013].

Here we consider NV ensemble and only $\ket{m_{s}=1}$ and $\ket{m_{s}=0}$ are used as a 2-level system. The spin-1 operators are

$\displaystyle S_{x}^{(1)}=\frac{1}{\sqrt{2}}\begin{pmatrix}0&1&0\\ 1&0&1\\ 0&1&0\\ \end{pmatrix},S_{y}^{(1)}=\frac{1}{\sqrt{2}}\begin{pmatrix}0&-i&0\\ i&0&-i\\ 0&i&0\\ \end{pmatrix},S_{z}^{(1)}=\begin{pmatrix}1&0&0\\ 0&0&0\\ 0&0&-1\\ \end{pmatrix}.$

(S10)

If we only take the $\ket{m_{s}=1}$, $\ket{m_{s}=0}$ subspace into consideration, the relations between the ‘truncated’ spin-1 operators and the spin-1/2 operators are:

$\displaystyle S_{y}^{(1)}=\sqrt{2}S_{y}^{(\frac{1}{2})},S_{y}^{(1)}=\sqrt{2}S_{y}^{(\frac{1}{2})},S_{z}^{(1)}=\frac{I}{2}+S_{x}^{(\frac{1}{2})}$

(S11)

Plugging Eq. (S11) into Eq. (S8), we get the effective dipole-dipole interaction Hamiltonian for NV ensemble $\ket{m_{s}=1}$, $\ket{m_{s}=0}$ subspace [Kucsko et al., 2018,Choi et al., 2017b]:

$\displaystyle H_{\text{DD,NV}}=\sum_{i<j}V_{ij}(S^{(\frac{1}{2})}_{zi}S^{(\frac{1}{2})}_{zj}-S^{(\frac{1}{2})}_{xi}S^{(\frac{1}{2})}_{xj}-S^{(\frac{1}{2})}_{yi}S^{(\frac{1}{2})}_{yj})$

(S12)

Fig. S1(a) shows the Classical Fisher Information (CFI) optimization results for 2D square lattice spin configuration. They are similar to the results we get in Fig.(2) of the main text for spin-1/2 systems.

When the system has large local disorder, the flip-flop terms in the dipolar interaction Hamiltonian Eq. (S8) are suppressed because of the large energy gap:

	$\displaystyle H_{\text{DD,Ising}}$	$\displaystyle=\sum_{i}\delta_{i}S^{(\frac{1}{2})}_{zi}+\sum_{i<j}V_{ij}(2S^{(\frac{1}{2})}_{zi}S^{(\frac{1}{2})}_{zj}-S^{(\frac{1}{2})}_{xi}S^{(\frac{1}{2})}_{xj}-S^{(\frac{1}{2})}_{yi}S^{(\frac{1}{2})}_{yj})$
		$\displaystyle=\sum_{i}\delta_{i}S^{(\frac{1}{2})}_{zi}+\sum_{i<j}2V_{ij}(S^{(\frac{1}{2})}_{zi}S^{(\frac{1}{2})}_{zj}-S^{(\frac{1}{2})}_{+i}S^{(\frac{1}{2})}_{-j}-S^{(\frac{1}{2})}_{-i}S^{(\frac{1}{2})}_{+j})$
		$\displaystyle\approx\sum_{i}\delta_{i}S^{(\frac{1}{2})}_{zi}+\sum_{i<j}2V_{ij}S^{(\frac{1}{2})}_{zi}S^{(\frac{1}{2})}_{zj}.$		(S13)

This location-dependent single-spin energy shift ($\delta_{i}$) can be canceled by spin-echo pulse sequence where the interaction gate $D(\tau)$ needs to be applied:

	$\displaystyle D(\tau)$	$\displaystyle=R_{x}(\pi)\exp[-i\tau H_{\text{DD,Ising}}]R_{x}(\pi)\exp[-i\tau H_{\text{DD,Ising}}]$
		$\displaystyle=\exp[-i\tau\sum_{i<j}2V_{ij}S^{(\frac{1}{2})}_{zi}S^{(\frac{1}{2})}_{zj}].$		(S14)

Eq. (Ising type interaction (large local disorder)) is also valid when the local disorder $\delta_{i}$ is comparable with the interaction strength $V_{ij}$. If there is local disorder in the dipolar-interacting spin ensemble, applying spin-echo will generate the interaction gate $D(\tau)$ where the local disorder terms are canceled.

The CFI optimization results by using the effective Ising type interaction Hamiltonian $H_{\text{DD,Ising}}=\sum_{i<j}2V_{ij}S^{(\frac{1}{2})}_{zi}S^{(\frac{1}{2})}_{zj}$ is shown in Fig. S1 (b).

From Fig. S1, the CFI results close to the Heisenberg Limit are observed, indicating that the variational method can be applied to different kinds of spins in solid state systems and generate highly entangled state for high-precision quantum metrology. We also observe that for shallow variational circuits, the CFI ‘oscillation’ between even and odd spin numbers only appears when there are flip-flop terms in the Hamiltonian. For Ising type interaction, the ‘oscillation’ disappears.

The complete data for dipolar-interacting spin systems’ CFI optimization is shown in this section. Fig. S3(a) shows the 50-cases averaged optimization results for 3D-random spin configurations, the variational circuit layer number $m$ is chosen from 1 to 10. The optimized CFI results are approaching to the Heisenberg Limit (HL) when more layers ($m$) are used. However, when $m>7$, the CFI results stop increasing. This CFI ‘saturation’ effect might be caused by two reasons. First, when $m$ is larger, the number of the local maximum points in the high dimensional parameter space increases. This could potentially cause the optimizer to stuck in the local maximum point. Sometimes, take $N=7,m=10$ data in Fig. S3(a) as an example, adding more variational layers even leads to a lower CFI optimization result. The ‘local maximum’ problem could be solved by more advanced and powerful optimization algorithms, such as reinforcement learning [Sutton and Barto, 2018; Silver et al., 2017; Peng et al., 2021], and more computational resources. Second, the ‘saturation’ effect reflects the global maximum CFI one can reach, no matter what kind of optimization algorithm is applied. It’s still an open question what is the highest CFI the spin ensemble could reach for a given configuration.

Fig. S3(b)-(d) show the CFI optimization result for 1D chain, 2D square lattice and 2D symmetric cycle spin configurations. The results of regular patterns are better than those of 3D-random pattern.

As shown by the schematic on the left, the distances between spin No.4 and spin No.5, 7, and 9 are the same, so the interaction strengths between each pair are the same. Similarly, the distance between spin No.4 and spin No.2, 3, 6, and 8 are the same (smaller). Therefore, the plateau features in Fig. S4 are likely due to this symmetry: adding one more spin to the lattice does not require an extra layer to reach a given percentage of the CFI.

Due to the decaying feature ($\frac{1}{r^{3}}$) of dipolar interaction strength, the resulting states might be mainly generated by nearest-neighbor interaction. For studying ‘how much’ interaction is essential for generating the resulting entangled states, we calculate the overlap (state fidelity [Nielsen and Chuang, 2010]) between the original state and the new state, which is generated by using the cutoff Hamiltonian and optimized parameters. A cutoff interaction strength $f_{\text{cutoff}}$ is chosen, and all the pairwise potential $V_{ij}$ smaller than $f_{\text{cutoff}}$ are set to zero in the cutoff Hamiltonian. Fig. S5 shows the relation between the state fidelity $F$ versus $f_{\text{cutoff}}$. A state fidelity value less than 1 is observed when $f_{\text{cutoff}}$ is set to be equal to the averaged nearest-neighbor interaction strength $f_{\text{dd}}$. This result reflects higher order interactions in the spins ensemble are utilized for the entangled state generation.

We run the optimization with imperfect readout for $N=4$ and $N=10$ 2D square lattice spin configurations. The optimized states resemble GHZ states (high RF), SSS (low RF), CSS (RF close to 50%). For $N=4$ case, the Gaussian state appears for RF lower than 92%, but for $N=10$ case, the Gaussian states appears when RF is about 96%. We expected that for large spin system with finite RF, Gaussian states (e.g. SSS) are advantageous for quantum-enhanced metrology.

Based on the simulation results shown in Fig.4(c) in main text, we need $\bar{f}_{dd}T_{2}\geq 5$ to generate metrological states that beat the SQL. It’s worth mentioning that the $T_{2}$ in this situation stands for the coherence time without the dipole-dipole interaction’s influence. During the state preparation step, the dipolar interactions between the spins are included in the system Hamiltonian for the entanglement generation ($D$ gate in Fig.1(c) in main text). Thus, the $T_{2}$^(DD) in Table S1 is a lower bound and $T_{2}$^(best) is a more precise estimation for the spin coherence time.

In general, the Classical Fisher Information (CFI) measures the sensitivity of a statistical model to small changes of a parameter $\theta$ [Braunstein and Caves, 1994, Strobel et al., 2014]. Let $Z$ be a random variable and $P_{z}(\theta)\equiv P(z|\theta)$ be its probability distribution which depends on $\theta$. Let $\Theta$ be an unbiased estimator of $\theta$, i.e.

$\displaystyle\theta=\langle\Theta\rangle=\sum_{z}\Theta\cdot P_{z}(\theta).$

(S17)

From Eq. (S17) and the fact that the sum of probabilities of all outcomes is $1$,

	$\displaystyle 1$	$\displaystyle=\frac{\partial\langle\Theta\rangle}{\partial\theta}=\frac{\partial}{\partial\theta}\sum_{z}\Theta P_{z}(\theta),$		(S18)
	$\displaystyle 0$	$\displaystyle=\frac{\partial}{\partial\theta}\sum_{z}P_{z}(\theta).$		(S19)

Subtracting Eq. (S19) multiplied by $\theta$ from Eq. (S18), we get

	$\displaystyle 1$	$\displaystyle=\sum_{z}(\Theta-\theta)\frac{\partial}{\partial\theta}P_{z}(\theta)$
		$\displaystyle=\sum_{z}P_{z}(\theta)(\Theta-\theta)\frac{1}{P_{z}(\theta)}\frac{\partial}{\partial\theta}P_{z}(\theta)$
		$\displaystyle=\langle(\Theta-\theta)\frac{1}{P_{z}(\theta)}\frac{\partial}{\partial\theta}P_{z}(\theta)\rangle$
		$\displaystyle=\langle(\Theta-\theta)\frac{\partial}{\partial\theta}\log P_{z}(\theta)\rangle.$		(S20)

Letting $X=\Theta-\theta$ and $Y=\frac{\partial}{\partial\theta}\log P_{z}(\theta)$, by the Cauchy-Schwartz inequality for random variables: $\langle XY\rangle^{2}\leq\langle X^{2}\rangle\langle Y^{2}\rangle$, we have

$\displaystyle\langle(\Theta-\theta)^{2}\rangle\bigg{\langle}\left(\frac{\partial}{\partial\theta}\log P_{z}(\theta)\right)^{2}\bigg{\rangle}\geq 1,$

(S21)

where

	$\displaystyle\langle(\Theta-\theta)^{2}\rangle$	$\displaystyle=\langle\Theta^{2}\rangle-(2\theta\langle\Theta\rangle-\langle\theta^{2}\rangle)$
		$\displaystyle=\langle\Theta^{2}\rangle-(2\theta^{2}-\theta^{2})$
		$\displaystyle=\langle\Theta^{2}\rangle-\langle\Theta\rangle^{2}$
		$\displaystyle=\Delta\Theta^{2}$		(S22)

is the variance of $\Theta$. Defining

$\displaystyle\mathrm{CFI}=\sum_{z}P_{z}(\theta)\left(\frac{\partial}{\partial\theta}\log P_{z}(\theta)\right)^{2},$

(S23)

we have

$\displaystyle\Delta\Theta^{2}\geq\frac{1}{\mathrm{CFI}}.$

(S24)

If the measurement is repeated $M$ times, then by the additive property of CFI, we obtain the Cramér-Rao bound:

$$\Delta\Theta^{2}\geq\frac{1}{M\cdot\mathrm{CFI}}.$$

(S25)

In our variational circuit, we use CFI with respect to an infinitesimal angle $\phi$ as the cost function to generate entangled states. In our program, we use a method similar to parameter shift to calculate the $\mathrm{CFI}_{\phi}$ of our optimized states [Cerezo et al., 2021,Meyer et al., 2021,Schuld et al., 2019]. In the following notation,

1. 1.

   $z$ represents a multi-qubit state in the z-basis;

2. 2.

   $\mathcal{U}(\phi)=e^{-i\phi J_{y}}$ is the rotation operator where $\phi$ is a small angle;

3. 3.

   $\psi$ is the state we create from the variational circuit;

4. 4.

   $P_{z}(\theta)$ is the probability of measuring the state $z$ with the state after rotation.

Then

	$\displaystyle\frac{\partial}{\partial\phi}P_{z}(\phi)\Bigr{\rvert}_{\phi\rightarrow 0}=$	$\displaystyle\frac{\partial}{\partial\phi}|\langle z|\mathcal{U}(\phi)|\psi\rangle|^{2}\Bigr{\rvert}_{\phi\rightarrow 0}$
	$\displaystyle=$	$\displaystyle\frac{\partial}{\partial\phi}\langle\psi|\mathcal{U}^{\dagger}(\phi)|z\rangle\langle z|\mathcal{U}(\phi)|\psi\rangle\Bigr{\rvert}_{\phi\rightarrow 0}$
	$\displaystyle=$	$\displaystyle\langle\psi|\mathcal{U}^{\dagger}(\phi)iJ_{y}|z\rangle\langle z|\mathcal{U}(\phi)|\psi\rangle\Bigr{\rvert}_{\phi\rightarrow 0}+\langle\psi|\mathcal{U}^{\dagger}(\phi)|z\rangle\langle z|\mathcal{U}(\phi)(-i)J_{y}|\psi\rangle\Bigr{\rvert}_{\phi\rightarrow 0}$
	$\displaystyle=$	$\displaystyle i\bra{\psi}\bigg{(}J_{y}\ket{z}\bra{z}-\ket{z}\bra{z}J_{y}\bigg{)}\ket{\psi}.$		(S26)

Note that assuming the rotation operator $\mathcal{U}(\phi)=e^{-i\phi J_{y}}\equiv\mathcal{U}_{y}(\phi)$ along $y$-axis is for calculation simplicity. In experiments, the signal (e.g. the external B-field) usually induces a rotation along $z$-axis, $\mathcal{U}_{z}(\phi)=e^{-i\phi J_{z}}$. It’s equivalent to assume that the prepared state is firstly rotated by a $R_{x}(\pi/2)$ pulse and then accumulates a signal $\phi$ along $y$-axis, or firstly accumulates a signal along $z$-axis and then rotated by $R_{x}(-\pi/2)$ pulse [Strobel et al., 2014]. In another word, $R_{x}(-\pi/2)\mathcal{U}_{z}(\phi)=\mathcal{U}_{y}(\phi)R_{x}(\pi/2)$, so the signal accumulation process we assumed in the calculation is able to simulate the experiments.

After creating the entangled states, we want to know how useful they are in a Ramsey spectroscopy, where the signal we want to detect is a frequency $\omega$. By the same calculation as above except the difference that we take derivative with respect to $\omega=\frac{\phi}{t_{\text{R}}}$ where $t_{\text{R}}$ is the Ramsey sensing time, we have

$\displaystyle\mathrm{CFI}_{\omega}=\mathrm{CFI}_{\phi}\cdot t_{\text{R}}^{2}.$

(S27)

We illustrate the Ramsey protocol for a single qubit.

1. 1.

   The qubit is initialized into the ground state $\ket{0}$.

2. 2.

   A $\frac{\pi}{2}$ pulse along the y-direction is applied to transform it into a superposition state $\frac{1}{\sqrt{2}}(\ket{0}+\ket{1})$. Its matrix form is

   $\displaystyle\rho(t)=\frac{1}{2}\begin{pmatrix}1&1\\ 1&1\end{pmatrix}.$

   (S28)

3. 3.

   After evolving under noise and a signal with frequency $\omega$ for time $t$, its state becomes

   $\displaystyle\rho(t)=\frac{1}{2}\begin{pmatrix}1&e^{-i\omega t}e^{-2\gamma t}\\ e^{i\omega t}e^{-2\gamma t}&1\end{pmatrix}$

   (S29)

   where $\gamma$ is the decoherence rate.

4. 4.

   A second $\frac{\pi}{2}$ pulse along the x-direction is applied for readout. The qubit is then in the state

   $\displaystyle R_{x}\left(\frac{\pi}{2}\right)\rho(t)R_{x}^{\dagger}\left(\frac{\pi}{2}\right)$

   (S30)

   .

5. 5.

   After the rotation, the probability of the qubit being in the ground state is

   $\displaystyle P_{0}=\frac{1}{2}+\frac{1}{2}e^{-2\gamma t}\sin{\omega t}.$

   (S31)

The CFI with respect to $\omega$ is

$\displaystyle\mathrm{CFI}_{\omega}$

$\displaystyle=\frac{1}{P_{0}}\left(\frac{\partial P_{0}}{\partial\omega}\right)^{2}+\frac{1}{P_{1}}\left(\frac{\partial P_{1}}{\partial\omega}\right)^{2}=\frac{t^{2}\cos^{2}{\omega t}}{e^{4\gamma t}-\sin^{2}{\omega t}}.$

(S32)

Assuming only quantum projection noise, the Signal-to-Noise Ratio (SNR) is $\frac{\delta P_{0}}{\sqrt{\frac{1}{M}P_{0}(1-P_{0})}}$ where $M$ is the total number of measurements. Then

$\displaystyle\mathrm{SNR}^{2}=\frac{Mt^{2}\cos^{2}{\omega t}\delta\omega^{2}}{e^{4\gamma t}-\sin^{2}{\omega t}}.$

(S33)

Assuming no time overhead, i.e., $M=\frac{T_{\text{tot}}}{t_{\text{R}}}$ where $T_{\text{tot}}$ is the total measurement time and $t_{\text{R}}$ is the time between Ramsey pulses, we obtain the relationship

$$\mathrm{CFI}_{\omega}\cdot\frac{T_{\text{tot}}}{t_{\text{R}}}\cdot\delta\omega^{2}=\mathrm{SNR}^{2}.$$

(S34)

In unit time ($T_{\text{tot}}=1$), when SNR$=1$, the smallest signal we can measure is

$\displaystyle\delta\omega=\frac{1}{\sqrt{M\cdot\mathrm{CFI}_{\omega}}},$

(S35)

leading to the saturated Cramér-Rao bound.

Since a measurement collapses a quantum state to an eigenstate of the observable, it’s impossible to directly measure $P(\theta)$. In experiments, we repeat the sequence to obtain the results for estimating the $P(\theta)$ and then get an estimate value of $\theta$. To understand the relation between the variance of the estimation and CFI, we introduce the Maximum Likelihood Estimator (MLE), which has asymptotic properties to saturate the Cramér-Rao bound. Below we summarize the proof given in [Panchenko, ].

Let $\bm{X}=\{X_{1},X_{2},...,X_{M}\}$ be a collection of independent and identically distributed (i.i.d.) random variables with a parametric family of probability distributions $\{P(X|\theta)|\theta\in\Theta\}$, where $\theta$ is an unknown parameter and $\Theta$ is the parameter space. Let $\bm{x}=\{x_{1},x_{2},...,x_{M}|x_{i}\in X_{i}\}$ be the experimental data set from $M$ repetitions. The goal is to estimate $\theta$ (the signal we want to measure) from $\bm{x}$, i.e., find $\theta$ that is most likely to produce the outcome $\bm{x}$. Thus, the normalized log-likelihood function is defined as

$\displaystyle L_{M}(\theta)$

$\displaystyle=\frac{1}{M}\log P(\bm{X}|\theta)=\frac{1}{M}\log\prod_{i=1}^{M}P(X_{i}|\theta)=\frac{1}{M}\sum_{i=1}^{M}\log P(X_{i}|\theta).$

(S36)

A MLE maximizes the log-likelihood function

$$\Theta_{\mathrm{MLE}}=\operatorname*{argmax}_{\theta\in\Theta}L_{M}(\theta).$$

(S37)

In the following, we first show that

1. 1.

   $\Theta_{\mathrm{MLE}}$ converges to the true parameter $\theta_{0}$;

2. 2.

   the distribution of $\sqrt{M}(\Theta_{\mathrm{MLE}}-\theta_{0})$ tends to a normal distribution $\mathcal{N}\left(0,\frac{1}{\mathrm{CFI}_{\theta_{0}}}\right)$ as $M$ increases.

In other words, not only does the MLE converge to the true parameter, it converges at a rate $\frac{1}{\sqrt{M}}$.

Define

$\displaystyle L(\theta)=\langle\log P(\bm{X}|\theta)\rangle_{\theta_{0}}$

(S38)

which denotes the expected log-likelihood function with respect to $\theta_{0}$, then by the Weak Law of Large Numbers (WLLN), the average outcomes from a large number of trials should approach the expected value:

$$\forall\theta,L_{M}(\theta)\xrightarrow{\mathit{M\rightarrow\infty}}L(\theta).$$

(S39)

In fact, $\theta_{0}$ maximizes $L(\theta)$:

	$\displaystyle\forall\theta,L(\theta)-L(\theta_{0})$	$\displaystyle=\langle\log P(\bm{X}|\theta)\rangle_{\theta_{0}}-\langle\log P(\bm{X}|\theta_{0})\rangle_{\theta_{0}}$
		$\displaystyle=\bigg{\langle}\log\frac{P(\bm{X}|\theta)}{P(\bm{X}|\theta_{0})}\bigg{\rangle}_{\theta_{0}}$
		$\displaystyle\leq\bigg{\langle}\frac{P(\bm{X}|\theta)}{P(\bm{X}|\theta_{0})}-1\bigg{\rangle}_{\theta_{0}}$
		$\displaystyle=\sum_{\bm{x}\in\bm{X}}\left(\frac{P(\bm{x}|\theta)}{P(\bm{x}|\theta_{0})}-1\right)P(\bm{x}|\theta_{0})$
		$\displaystyle=1-1=0.$		(S40)

Moreover, we show that $\theta_{0}$ is the unique maximizer. Jensen’s inequality states that for a strictly convex function $f$ and a random variable $Y$,

$\displaystyle\langle f(Y)\rangle>f(\langle Y\rangle).$

(S41)

Taking $f(y)=-\log y$ and $P(\bm{X}|\theta)\neq P(\bm{X}|\theta_{0})$, we have

$$\bigg{\langle}-\log\frac{P(\bm{X}|\theta)}{P(\bm{X}|\theta_{0})}\bigg{\rangle}_{\theta_{0}}>-\log\langle\frac{P(\bm{X}|\theta)}{P(\bm{X}|\theta_{0})}\bigg{\rangle}_{\theta_{0}}=0,$$

(S42)

or

$$L(\theta_{0})>L(\theta).$$

(S43)

Therefore, since

1. 1.

   $\Theta_{\mathrm{MLE}}$ maximizes $L_{M}(\theta)$,

2. 2.

   $\theta_{0}$ maximizes $L(\theta)$, and

3. 3.

   $L_{M}(\theta)\xrightarrow{\mathit{M\rightarrow\infty}}L(\theta)$,

$\Theta_{\mathrm{MLE}}$ converges to $\theta_{0}$.