Pulsed Quantum-State Reconstruction of Dark Systems

TL;DR

This paper introduces a new method for reconstructing the quantum states of systems that cannot be directly measured, using a controlled probe and pulse sequences to extract detailed quantum information.

Contribution

It presents a novel pulse-based scheme for quantum-state reconstruction of dark systems, applicable to spins and oscillators, with robustness against noise.

Findings

Successfully reconstructs dark spin states and harmonic oscillator Wigner functions

Scheme is robust against slow noise in the probe

Feasible in solid-state spins and trapped ions

Abstract

We propose a novel strategy to reconstruct the quantum state of dark systems, i.e., degrees of freedom that are not directly accessible for measurement or control. Our scheme relies on the quantum control of a two-level probe that exerts a state-dependent potential on the dark system. Using a sequence of control pulses applied to the probe makes it possible to tailor the information one can obtain and, for example, allows us to reconstruct the density operator of a dark spin as well as the Wigner characteristic function of a harmonic oscillator. Because of the symmetry of the applied pulse sequence, this scheme is robust against slow noise on the probe. The proof-of-principle experiments are readily feasible in solid-state spins and trapped ions.

Equations

$$H=|0{⟩}_p⟨0|V_0+|1{⟩}_p⟨1|V_1,$$

$$\displaystyle|\Psi(t)\rangle=\frac{1}{\sqrt{2}}\big{(}|0\rangle_{p}U_{0}|\psi\rangle+|1\rangle_{p}U_{1}|\psi\rangle\big{)},$$

$$u_0=$$

$$u_1=$$

$$⟨{\sigma }_x^p⟩=$$

$$⟨{\sigma }_y^p⟩=$$

$$⟨{\sigma }_y^p⟩=-\sin \varphi \mathbf{n}\cdot \mathbf{r}.$$

$$V_0=\frac{{\omega }_0}{2}{\sigma }_z,V_1=\frac{{\omega }_1}{2}\left(v_x{\sigma }_x+v_z{\sigma }_z\right),$$

$$\cos \varphi =1-{\sin }^2(N\theta )(1-{\mathbf{n}}_0\cdot {\mathbf{n}}_1),$$

$$\mathbf{n}=\frac{\sqrt{2(1-{\mathbf{n}}_0\cdot {\mathbf{n}}_1)}{\sin }^2(N\theta )}{\sin \varphi \sin \theta }\left[\begin{array}{c}\sin \alpha \cos \beta +v_z\cos \alpha \sin \beta \\ \sin \theta \cot (N\theta )\\ -v_x\cos \alpha \sin \beta \end{array}\right]$$

$$V_0=\nu a^{†}a-\frac{g}{2}(a+a^{†}),V_1=\nu a^{†}a+\frac{g}{2}(a+a^{†}).$$

$$\xi (\tau ,N)=-2\frac{g}{\nu }\sin (N\nu \tau )\tan \left(\frac{\nu \tau }{2}\right)e^{iN\nu \tau },$$

$$⟨{\sigma }_x^p⟩+i⟨{\sigma }_y^p⟩=\chi \left(\xi \right),$$

$$\begin{array}{rl}{u}_0=& e^{-i\beta (v_z{\sigma }_z+v_x{\sigma }_x)}{e}^{-i\alpha {\sigma }_z},\\ u_1=& e^{-i\alpha {\sigma }_z}{e}^{-i\beta (v_z{\sigma }_z+v_x{\sigma }_x)},\end{array}$$

$${\text{n}}_0=\frac{1}{\sin \theta }\left[\begin{array}{c}{v}_x\cos \alpha \sin \beta \\ v_x\sin \alpha \sin \beta \\ \sin \alpha \cos \beta +v_z\cos \alpha \sin \beta \end{array}\right]$$

$$\cos \varphi ={\cos }^2(N\theta )+{\sin }^2(N\theta )({\text{n}}_0\cdot {\text{n}}_1),$$

$$\text{n}=\frac{\sin (N\theta )}{\sin \varphi }\left[\cos (N\theta )({\text{n}}_1-{\text{n}}_0)-\sin (N\theta )({\text{n}}_0\times {\text{n}}_1)\right].$$

$${\tau }_1=\frac{2\pi }{{\omega }_1+{\omega }_0},N=\frac{\pi }{4v_x}.$$

$$b(t+\Delta t)=b(t)e^{-\Delta t/t_b}+b_0\sqrt{1-e^{-2\Delta t/t_b}}\mathcal{N},$$

$$H_d=\sum _{k=1,2}\frac{{\omega }_0}{2}{\sigma }_z^{(k)}+\frac{A_x}{2}{\sigma }_x^{(1)}{\sigma }_x^{(2)},$$

$$H_{1}=\sum_{k=1,2}\bigg{[}\frac{a_{z}^{(k)}}{2}\sigma_{z}^{(k)}+\frac{a_{x}^{(k)}}{2}\sigma_{x}^{(k)}\bigg{]}.$$

$$V_0=\frac{A_x}{2}{\tilde{\sigma }}_x,V_1=\frac{A_x}{2}{\tilde{\sigma }}_x+\frac{A_z}{2}{\tilde{\sigma }}_z,$$

$$\begin{array}{r}{\tilde{r}}_x=\frac{1}{2}⟨{\sigma }_x^{(1)}{\sigma }_x^{(2)}+{\sigma }_y^{(1)}{\sigma }_y^{(2)}⟩,\\ {\tilde{r}}_y=\frac{1}{2}⟨{\sigma }_y^{(1)}{\sigma }_x^{(2)}-{\sigma }_x^{(1)}{\sigma }_y^{(2)}⟩,\end{array}$$

$$D^{†}(ϵ)V_{0}D(ϵ)=D(ϵ)V_1{D}^{†}(ϵ)=\nu a^{†}a-{ϵ}^2\nu .$$

$$u_0=e^{i{\phi }_1}{D}^{†}(ϵ)e^{-2i\nu \tau a^{†}a}D\left(2ϵe^{i\nu \tau }\right)D^{†}(ϵ),$$

$$U_0^{†}=e^{i{\phi }_2}{D}^{†}(ϵ){\left[D\left(2ϵ\left(1-e^{i\nu \tau }\right)\right)e^{2i\nu \tau a^{†}a}\right]}^{N}D(ϵ)$$

$$\begin{array}{rl}{\left[D(x)e^{iy{a}^{†}a}\right]}^N& =\left[\prod _{n=0}^{N-1}D\left(x{e}^{iny}\right)\right]e^{iNy{a}^{†}a}\\ & =e^{i\vartheta }D\left(x\frac{1-e^{iNy}}{1-e^{iy}}\right)e^{iNy{a}^{†}a},\end{array}$$

$$\begin{array}{rl}{U}_0^{†}=e^{i{\phi }_3}{D}^{†}(ϵ)D(& \zeta )e^{2iN\nu \tau a^{†}a}D(ϵ),\\ U_1=e^{-i{\phi }_3}{D}^{†}(ϵ)D(& {\zeta }^{\ast })e^{-2iN\nu \tau a^{†}a}D(ϵ),\end{array}$$

$$\zeta =\frac{g}{\nu }\left(1-e^{i\nu \tau }\right)\frac{1-e^{2iN\nu \tau }}{1-e^{2i\nu \tau }}$$

$$U_0^{†}{U}_1=D\left(-2\frac{g}{\nu }\sin (N\nu \tau )\tan \left(\frac{\nu \tau }{2}\right)e^{iN\nu \tau }\right),$$

$$⟨n|\rho |m⟩=\frac{1}{\pi }\int d^2\xi ⟨n|D(-\xi )|m⟩\chi (\xi ).$$

$$⟨n|D(\eta )|m⟩=\sqrt{\frac{m!}{n!}}{\eta }^{n-m}{L}_m^{(n-m)}\left(|\eta {|}^2\right)e^{-\frac{|\eta {|}^2}{2}},$$

Full text

*Introduction.—*The measurement of the quantum state of a system is a prerequisite ingredient in most modern quantum experiments, ranging from fundamental tests of quantum mechanics Aspelmeyer et al. (2014); Neukirch and Vamivakas (2015) to various quantum-information-processing tasks Nielsen and Chuang (2003); Briegel et al. (2009); Ladd et al. (2010). However, even with the rapid progress in the coherent manipulation and quantum-state tomography of several quantum systems, such as photons Altepeter et al. (2005); Deléglise et al. (2008), electron spins Elzerman et al. (2004); Barthel et al. (2009); Morello et al. (2010), atomic qubits Bochmann et al. (2010), superconducting circuits Steffen et al. (2006); Fedorov et al. (2014), and mechanical resonators Vanner et al. (2013); Rashid et al. (2017), many quantum systems still remain difficult to access for a direct observation of their state, systems we will refer to as dark. In order to circumvent the requirement of such a direct access, a promising technique is to employ an auxiliary quantum system as a measurement probe, on which measurements as well as coherent manipulations can be performed Leibfried et al. (1996); Bertet et al. (2002); Wallraff et al. (2005); Burgarth et al. (2011, 2015); Sone and Cappellaro (2017a, b); Zhang et al. . Interferometry Cronin et al. (2009) based on such a measurement probe allows us to extract information on a target system Recati et al. (2005); Quan et al. (2006); Dorner et al. (2013); Peng et al. (2015); Correa et al. (2015); Asadian et al. (2014). Nevertheless, it still remains a key challenge to achieve a full quantum-state tomography of dark systems without requiring any direct control.

In this Letter, we propose a general scheme for a probe-measurement based quantum-state reconstruction of dark systems, where the obtainable dark-system quantities can be tailored by a pulsed control of the two-level probe we employ. The scheme does not require any manipulation of the dark systems or the controllability of the coupling to the probe. This is a requirement on which, for example, many previous reconstruction methods for continuous-variable systems depend Lutterbach and Davidovich (1997); Kim et al. (1998); Singh and Meystre (2010); Casanova et al. (2012); Mazzola et al. (2013); Taketani et al. (2014). Additionally, it inherits the feature of robustness against noise on the probe from pulsed dynamical decoupling Uhrig (2007); de Lange et al. (2010); Naydenov et al. (2011), making it suitable also for noisy environments such as solid-state platforms. The proposed strategy is exemplified at the quantum-state tomography of a dark spin-1/2 and a dark harmonic oscillator by reconstructing their density operator and Wigner characteristic function, respectively. We discuss the feasibility of the proof-of-principle experiments to testify the distinct advantages of the present proposal in solid-state spin Taminiau et al. (2012); Zhao et al. (2011); Kolkowitz et al. (2012); Liu et al. (2017) and trapped-ion systems Johanning et al. (2009); Khromova et al. (2012); Sriarunothai et al. (2018) within state-of-the-art experimental capabilities. The scheme is applicable to other dark systems as in a variety of physical settings, making it a versatile tool for quantum measurements.

*Pulsed state-reconstruction scheme.—*The probe we consider is a generic two-level system described by the Hamiltonian $H_{p}=(\omega_{p}/2)\sigma_{z}^{p}$, with the Pauli operator $\sigma_{z}^{p}=|1\rangle_{p}\langle 1|-|0\rangle_{p}\langle 0|$. We denote the Hamiltonian of the dark system by $H_{d}$. The underlying idea of the proposed scheme is that the interaction between the probe and the dark system is given by probe-state-dependent potentials $H_{0}$ and $H_{1}$ for the dark system, i.e., an interaction of the form $H_{\rm int}=|0\rangle_{p}\langle 0|H_{0}+|1\rangle_{p}\langle 1|H_{1}$. The dynamics of the combined system is then governed by $H_{p}+H_{d}+H_{\rm int}$ and in the interaction picture with respect to $H_{p}$ this Hamiltonian has the form

[TABLE]

with $V_{0}=H_{d}+H_{0}$ and $V_{1}=H_{d}+H_{1}$ acting on the dark system. The dynamics generated by this Hamiltonian is used to obtain information about the state of the dark system.

In order to do so, the probe is initialized in the superposition state $|+\rangle_{p}=(|0\rangle_{p}+|1\rangle_{p})/\sqrt{2}$, such that the initial state of the full system has the form $|\Psi(0)\rangle=|+\rangle_{p}|\psi\rangle$, with the dark-system state $|\psi\rangle$. The free evolution of this state under the Hamiltonian (1) generates entanglement between the probe and the dark system and thereby allows to connect measurements on the probe with quantities of the dark system. However, as we will show, appreciably more information can be obtained by the application of a series of pulses that manipulate the probe Li et al. (2011); Zhao et al. (2014). Explicitly, we apply a series of $2N$ $\pi$ pulses, all separated by the free-evolution time $\tau$, and thereby modulate the effective potential acting on the dark system. After such an evolution of total time $t=2N\tau$ the state of the full system has evolved into

[TABLE]

where the state-dependent dark-system time-evolution operators are given by $U_{0}=u_{0}^{N}$ and $U_{1}=u_{1}^{N}$, with the single pulse-segment evolution operators

[TABLE]

In order to obtain any information on the dark-system state from this dynamics, a necessary condition is that the operators $V_{0}$ and $V_{1}$ do not commute, since otherwise the above evolution operators coincide. Following this time evolution, a measurement of the probe Pauli vector $\bm{\sigma}_{p}=(\sigma_{x}^{p},\sigma_{y}^{p},\sigma_{z}^{p})$ can be performed resulting in $\langle\bm{\sigma}_{p}\rangle={\rm Tr}\{\bm{\sigma}_{p}\varrho(t)\}$, with $\varrho(t)=|\Psi(t)\rangle\langle\Psi(t)|$. The generalization to initially separable density operators of the form $\varrho(0)=|+\rangle_{p}\langle+|\,\rho$, with the possibly mixed initial dark-system density operator $\rho$, is straightforward and yields

[TABLE]

and $\langle\sigma_{z}^{p}\rangle=0$. As we see, the probe-measurement outcomes $\langle\sigma_{x}^{p}\rangle$ and $\langle\sigma_{y}^{p}\rangle$ are, respectively, equal to the expectation value of the real and imaginary part of the operator $U_{0}^{\dagger}U_{1}$ in the initial dark-system state $\rho$. These expectation values are the information we can extract through Pauli measurements on the probe and by changing the pulse-sequence parameters $\tau$ and $N$ we can control to which dark-system quantity they correspond. The information is extracted by measuring the coherence of the probe and its dephasing thus affects the reconstruction fidelity. The scheme is feasible as long as the extended probe coherence time by the pulsed strategy is longer than the total measurement time. Up to this point we make no assumption on the nature of the dark system. In the following we give two explicit examples, one discrete and one continuous variable system, and demonstrate in both cases how the unitary $U_{0}^{\dagger}U_{1}$ can be engineered in order to perform a state reconstruction of these dark systems.

*State reconstruction of a spin-1/2 system.—*As a first case we consider a dark spin-1/2 system. Its density operator can be written as $\rho=(\mathbb{1}+\mathbf{r}\cdot\bm{\sigma})/2$, with the unity operator $\mathbb{1}$, the dark-spin Pauli vector $\bm{\sigma}$, and the Bloch vector $\mathbf{r}$ fulfilling $\mathbf{r}={\rm Tr}\{\bm{\sigma}\rho\}$. On the other hand, the unitary $U_{0}^{\dagger}U_{1}$ takes the form $U_{0}^{\dagger}U_{1}=\cos\phi\,\mathbb{1}-i\sin\phi\,\mathbf{n}\cdot\bm{\sigma}$, with a unit vector $\mathbf{n}$. Comparing this expression with Eq. (4) reveals $\langle\sigma_{x}^{p}\rangle=\cos\phi$ and allows us to connect the probe-measurement outcome $\langle\sigma_{y}^{p}\rangle$ with the dark-spin Bloch vector $\mathbf{r}$ according to

[TABLE]

Three independent measurements for different pulse-sequence parameters $\tau$ and $N$ are thereby sufficient for a full state reconstruction of the dark spin. It can also be seen that for a faithful measurement of any component of the dark-spin Bloch vector a crucial condition is $\cos\phi=0$. This additionally makes it possible to engineer the pulse sequences such as to obtain the three components $r_{\kappa}$ separately by ensuring $\sin\phi\,n_{\kappa}=-1$, for $\kappa=x,y,z$.

For a general dark-spin Hamiltonian $H_{d}=(\omega_{0}/2)\sigma_{z}$ this reconstruction can be achieved by the probe-state-dependent potentials $H_{0}=0$ and $H_{1}=(a_{z}/2)\sigma_{z}+(a_{x}/2)\sigma_{x}$, where $a_{z}$ and $a_{x}$ arise from the coupling between the probe and the dark spin. This results in the Hamiltonians

[TABLE]

with $v_{x}=a_{x}/\omega_{1}$, $v_{z}=(\omega_{0}+a_{z})/\omega_{1}$, and $\omega_{1}^{2}=(\omega_{0}+a_{z})^{2}+a_{x}^{2}$. The above Hamiltonians represent one effective longitudinal field of strength $\omega_{0}$ associated with the probe ground state and the other one of strength $\omega_{1}$ associated with the probe excited state, which is tilted from the $z$ direction by the angle $\arctan(v_{x}/v_{z})$. From these state-dependent effective fields and the pulse sequence applied to the probe one can obtain the explicit form of the operator $U_{0}^{\dagger}U_{1}$.

Every pulse-sequence segment of length $2\tau$, i.e., first an evolution under $V_{0}$ and then under $V_{1}$, or vice versa, produces a state-dependent rotation given by the unitaries $u_{k}=\exp(-i\theta\mathbf{n}_{k}\cdot\bm{\sigma})$, for $k=0,1$, respectively. Here, the angle $\theta$ satisfies $\cos\theta=\cos\alpha\cos\beta-v_{z}\sin\alpha\sin\beta$ and the two rotation axes fulfill $\mathbf{n}_{0}\cdot\mathbf{n}_{1}=1-2v_{x}^{2}\sin^{2}\alpha\sin^{2}\beta/\sin^{2}\theta$, with $\alpha=\omega_{0}\tau/2$ and $\beta=\omega_{1}\tau/2$. The operators $U_{k}=u_{k}^{N}$ are then rotations around the same axis, but by the angle $N\theta$ and one obtains the expressions Kolkowitz et al. (2012); Taminiau et al. (2012); Sup

[TABLE]

for the quantities $\phi$ and $\mathbf{n}$. This is their functional dependence on the pulse-sequence parameters $\tau$ and $N$, which can be used to fully determine the Bloch vector $\mathbf{r}$, according to Eq. (5), from three independent probe measurements.

Among the possible choices for the parameters $\tau$ and $N$, which ensure $\cos\phi=0$, we can choose to measure the observables corresponding to the three components $r_{\kappa}$ of the Bloch vector, for which the additional condition $\sin\phi\,n_{\kappa}=-1$ has to be fulfilled, for $\kappa=x,y,z$. As an example, in the $y$ case, these two conditions are $\mathbf{n}_{0}\cdot\mathbf{n}_{1}=-1$ and $\sin^{2}(N\theta)=1/2$. Here, the first one is fulfilled for the evolution time $\tau_{1}=2\pi/(\omega_{1}+\omega_{0})$ and the second one for $N=\pi/4v_{x}$, yielding $\sin\phi\,n_{y}=-1$ Sup . The results for the measurement of $r_{y}$ are illustrated in Fig. 1, where we show $|\sin\phi\,n_{y}|$ as a function of $\tau$ and $N$. Our further simulations demonstrate that the measurement protocol is very robust against noise acting on the probe Sup . As a brief note, we mention that without the application of pulses, one would have the unitary $U_{0}^{\dagger}U_{1}=\exp(iV_{0}\tau)\exp(-iV_{1}\tau)$ and the reachable points within the Bloch sphere are confined to a cylinder of radius $v_{x}$ around the $z$ axis, making a measurement of $r_{x}$ and $r_{y}$ impossible. The further generalization to multispin dark systems is feasible by employing sufficient controllability conditions and the technique of Cartan decomposition Albertini and D’Alessandro (2002). We also remark that the measurement of some specific observable may already be of significant interest, e.g., for entanglement and quantum-criticality detection Sup ; Cai et al. (2013); Chen et al. (2013); Xu and del Campo .

*State reconstruction of a harmonic oscillator.—*As a second case we consider a continuous-variable dark system which is formed by a harmonic oscillator of frequency $\nu$ with the annihilation operator $a$. The interaction between the probe and the harmonic oscillator is assumed to be of the form $H_{\rm int}=(g/2)\sigma_{z}^{p}(a+a^{\dagger})$, leading to the state-dependent Hamiltonians

[TABLE]

Using the multiplication properties of the displacement operator $D(\eta)=\exp(\eta a^{\dagger}-\eta^{\ast}a)$, in a picture displaced by $g/2\nu$, the operator $U_{0}^{\dagger}$ can be brought into the form $D(\zeta)\exp(2iN\nu\tau a^{\dagger}a)$, while $U_{1}$ similarly is of the form $D(\zeta^{\ast})\exp(-2iN\nu\tau a^{\dagger}a)$, where $\zeta$ is a function of $\tau$ and $N$ Sup . In this way, we can obtain the unitary $U_{0}^{\dagger}U_{1}=D(\xi)$ in the simple form of a single displacement operator with the quantity

[TABLE]

which depends on the pulse-sequence parameters $\tau$ and $N$ Sup . Equation (4) then yields

[TABLE]

with the Wigner characteristic function $\chi$, which is defined as $\chi(\xi)={\rm Tr}\{D(\xi)\rho\}$ Cahill and Glauber (1969). This function over reciprocal phase space is the complex Fourier transform of the Wigner function Wigner (1932) and contains all information on the initial density operator $\rho$ of the harmonic oscillator. For a full knowledge of the characteristic function, the completeness of the displacement operators Cahill and Glauber (1969) allows an exact reconstruction of the density operator itself according to $\rho=\int{\rm d}^{2}\xi\,\chi(\xi)D^{\dagger}(\xi)/\pi$. For example, in the Fock basis, the matrix elements $\langle n|\rho|m\rangle$ can easily be obtained using this expression and the matrix elements of the displacement operator $\langle n|D^{\dagger}(\xi)|m\rangle$ Cahill and Glauber (1969).

In our scheme, every fixed pulse-sequence parameter $N$ corresponds to a closed curve $\xi(\tau,N)$ in reciprocal phase space, shown in Fig. 2(a) for several pulse numbers. Their periodicity in $\tau$ is determined by the harmonic oscillator frequency $\nu$, requiring a maximal necessary evolution time of $\tau=2\pi/\nu$. The maximal distance from the origin is reached for $\tau=\pi/\nu$ and has the value $4Ng/\nu$; i.e., it scales linearly with the pulse number. By varying $N$ we can sample the characteristic function along this manifold of curves, as shown by blue lines in Fig. 2(b) for the example of a squeezed vacuum state $S(\lambda)|0\rangle$, with $S(\lambda)=\exp[(\lambda^{\ast}a^{2}-\lambda a^{\dagger 2})/2]$. For a real squeezing parameter $\lambda=\log(1/2)$, as chosen here, and $\xi=\xi_{\rm r}+i\xi_{\rm i}$ the corresponding characteristic function is $\chi(\xi)=\exp(-\xi_{\rm r}^{2}/8-2\xi_{\rm i}^{2})$. The property $\chi(-\xi)=\chi(\xi)^{\ast}$ allows us to obtain the values of $\chi$ along these curves mirrored around the origin by complex conjugation, as represented by red curves in Fig. 2(b).

Figure 2(c) shows the density matrices reconstructed from the characteristic function for the squeezed vacuum state from Fig. 2(b) and a coherent state by an interpolation of $\chi$ using $N=(1,...,20)$. The results are in good agreement with the exact density matrices, showing trace distances of the order $10^{-3}$ Sup . As a comparison, for the case of no pulses applied to the probe only points on the circle $\xi=g[\exp(i\nu\tau)-1]/\nu$ in reciprocal phase space can be sampled, which would be insufficient for a satisfactory state reconstruction.

The fact that the characteristic function has its maximum $\chi(0)=1$ at the origin, and is mostly centered in this region, is favorable in experiments since the density of reachable points $\xi$ is high close to the origin. In contrast, the Wigner function, which contains the same information, can have its maximum at any point in phase space, making it necessary to scan over larger areas with schemes for its measurement. We also stress that contrary to other methods for the measurement of phase-space distributions, this scheme does not require any manipulation of the harmonic oscillator, such as a displacement operation prior to the measurement procedure, or a control of the coupling strength. This advantage would become particularly important for systems in which direct manipulation on the harmonic oscillator is hard to achieve.

*Potential experimental implementations.—*As an example for the state tomography of a spin-1/2 system, we use a nitrogen-vacancy (NV) center as a probe and a dark spin of a weakly coupled 13C nucleus in diamond Doherty et al. (2013). Under an external magnetic field of strength $B$ along the NV axis, i.e., the $z$ direction, the Hamiltonians then read $H_{\rm NV}=DS_{z}^{2}+\gamma_{e}BS_{z}$ and $H_{\rm C}=\gamma_{\rm C}BI_{z}$, where $\gamma_{e}$ and $\gamma_{\rm C}$ are the gyromagnetic ratios of the NV-center spin and the 13C nuclear spin, respectively, and $D/2\pi=2.87$ GHz is the electron-spin zero-field splitting. Furthermore, the components of $S_{\kappa}$ and $I_{\kappa}$ denote their respective spin-$1$ and spins-$1/2$ operators, for $\kappa=x,y,z$. One can choose the $x$ direction such that the hyperfine interaction between the NV center and the nuclear spin is of the form $H_{\rm hf}=A_{\parallel}S_{z}I_{z}+A_{\perp}S_{z}I_{x}$ Taminiau et al. (2012). The external magnetic field allows us to address specific transitions of the NV-center electronic states and thereby, for example, to use the two states $|0\rangle_{p}=|m_{s}\!=\!0\rangle$ and $|1\rangle_{p}=|m_{s}\!=\!1\rangle$ as our probe. This results in the state-dependent effective fields acting on the 13C nuclear spin given by Eq. (6) with $\omega_{0}=\gamma_{\rm C}B$, $a_{z}=A_{\parallel}$, $a_{x}=A_{\perp}$, and $\sigma_{\kappa}=2I_{\kappa}$, for $\kappa=x,y,z$. As an example, we consider a weakly coupled 13C with $A_{\parallel}/2\pi=2.54$ and $A_{\perp}/2\pi=13.22$ kHz under a magnetic field $B=15.4$ mT, which are the parameters used in Fig. 1. The assumption of instantaneous $\pi$ pulses is well justified, since pulse durations $t_{\pi}$ of tens of nanoseconds have been realized, and the free-evolution time for the measurement of Bloch vector components thereby fulfill $t_{\pi}\ll\tau$ de Lange et al. (2010); Robledo et al. (2010); de Lange et al. (2011). The condition $2N\tau<T_{2p}$ for the total evolution time can also be satisfied for achievable long probe coherence times $T_{2p}$ Ryan et al. (2010); Naydenov et al. (2011); Bar-Gill et al. (2013).

To show the feasibility of an experimental demonstration for a continuous-variable dark system, we consider the motional-state reconstruction of a single trapped ion in a magnetic field gradient Johanning et al. (2009); Khromova et al. (2012). We orient ourselves at parameters from Ref. Sriarunothai et al. (2018) with single 171Yb+ ions trapped in a linear Paul trap with an axial frequency $\nu/2\pi=117$ kHz. In this setup, a magnetic field of the form $B(z)=B_{0}+B_{1}z$ is applied along the trap axis $z$. As a probe two-level system one can choose the two sublevels $|0\rangle_{p}=|m_{s}\!=\!-1/2\rangle$ and $|1\rangle_{p}=|m_{s}\!=\!1/2\rangle$ with $m_{F}=0$ of the ${}^{2}S_{1/2}$ state, whose coherence time can be longer than 1000 s Fisk et al. (1997). The linear magnetic field gradient $B_{1}$ induces a coupling between the ion motion and the probe of the form $H_{\rm int}=(g/2)\sigma_{z}^{p}(a+a^{\dagger})$, yielding the Hamiltonians $V_{0}$ and $V_{1}$ given in Eq. (8). The coupling strength is determined by $g=2\mu_{\rm B}B_{1}/\sqrt{2M\hbar\nu}$, with the Bohr magneton $\mu_{\rm B}$ and the ion mass $M$. For the reported magnetic field gradient $B_{1}=19$ T/m, this results in a coupling ratio of $g/\nu=0.072$, which is roughly the one we used above in Fig. 2. In this system, high-fidelity $\pi$ pulses with durations on the order of tens of picoseconds have also been demonstrated for these ${}^{2}S_{1/2}$ states Campbell et al. (2010).

*Conclusions and outlook.—*We present a general scheme for the quantum-state reconstruction of a dark system, which is inaccessible for direct control and measurements. The scheme only relies on the pulsed control and readout of a probe-two-level system, while requiring no manipulation of the target system. We illustrate our idea at the state tomography of a spin and a harmonic oscillator. For both examples, we show the feasibility to implement proof-of-principle demonstrations in currently available experimental setups. Moreover, the measurement scheme is intrinsically robust against slow noise acting on the probe due to the incorporated dynamic decoupling. The present idea provides a versatile tool for quantum-state measurement and can be extended to more general scenarios, such as dark systems formed by higher spins, many-body systems, and novel mechanical systems. A further generalization to a sequence of nonequidistant pulses and continuous processes is possible and may increase the information that can be obtained of the dark system.

The authors thank C. Arenz, Y. Chu, S. B. Jäger, G. Morigi, M. B. Plenio, and R. Said for helpful discussions and comments. This work was supported by the National Natural Science Foundation of China (Grants No. 11574103, No. 11874024, No. 11690030, and No. 11690032), the National Key R&D Program of China (Grant No. 2018YFA0306600), the China Postdoctoral Science Foundation (Grant No. 2017M622398), and the National Young 1000 Talents Plan.

S.1 S.1. Spin-$\frac{1}{2}$ system

S.1.1 S.1.1. Basic principle

For the pulse sequence applied to the probe in the presented state-reconstruction scheme a single pulse-sequence segment, which is applied $N$ times, has the form $[\pi\!-\!\tau\!-\!\pi\!-\!\tau]$, i.e., an evolution time $\tau$, a $\pi$ pulse, a second evolution time $\tau$ followed by a final $\pi$ pulse. For such a segment the probe-state-dependent evolution operators for a dark spin 1/2 read

[TABLE]

where $v_{x}=a_{x}/\omega_{1}$, $v_{z}=(\omega_{0}+a_{z})/\omega_{1}$, $\omega_{1}^{2}=(\omega_{0}+a_{z})^{2}+a_{x}^{2}$, $\alpha=\omega_{0}\tau/2$ and $\beta=\omega_{1}\tau/2$. Since these operators are composed of two rotations, they can both be combined into a single rotation $u_{k}=\exp(-i\theta\textbf{n}_{k}\cdot\bm{\sigma})$, with $k=0,1$, where the angle fulfills $\cos\theta=\cos\alpha\cos\beta-v_{z}\sin\alpha\sin\beta$, and the two rotation axes are

[TABLE]

and $\mathbf{n}_{1}=(n_{0}^{x},-n_{0}^{y},n_{0}^{z})$. For the full pulse sequence we can thereby write $U_{k}=\exp(-iN\theta\,\textbf{n}_{k}\cdot\bm{\sigma})$ Kolkowitz et al. (2012); Taminiau et al. (2012), for $k=0,1$, and the unitary $U_{0}^{\dagger}U_{1}$ can likewise be written as a single rotation $U_{0}^{\dagger}U_{1}=\exp(-i\phi\,\textbf{n}\cdot\bm{\sigma})$, with

[TABLE]

Using $\textbf{n}_{0}\cdot\textbf{n}_{1}=1-2v_{x}^{2}\sin^{2}\alpha\sin^{2}\beta/\sin^{2}\theta$ this corresponds to the expressions given in the main text. On the other hand, for a time evolution of time $\tau$ without the application of pulses the unitary corresponding to $U_{0}^{\dagger}U_{1}$ merely has the form of $u_{1}$ with $\alpha$ replaced by $-\alpha$.

As an example to illustrate our idea, we now consider measurements of a single components of the Bloch vector $\mathbf{r}$, under the typical situation of weak coupling, viz. $a_{x},a_{z}\ll\omega_{0}$. In this case we find $v_{z}=1-v_{x}^{2}/2$. And in first order of $v_{x}$ the $z$ component $\sin\phi\,n_{z}$ can be neglected, which is valid for relatively short evolution times $\tau$. This leads to the measurement condition for $x$ and $y$, namely the proper parameters $\tau$ and $N$, given in the main text: (i) $n_{\kappa}=-1$ with $\kappa=x,y$ and (ii) $\cos\phi=0$. First we focus on the $y$ component, where condition (i) implies $n_{x}=0$, which means $\sin\alpha\cos\beta+v_{z}\cos\alpha\sin\beta=0$. In first order of $v_{x}$ this leads to the condition $\tau_{k}=2k\pi/(\omega_{1}+\omega_{0})$, with odd integers $k$. One also has $\theta=v_{x}$ and $\sin\alpha\sin\beta=1$, resulting in two anti-parallel rotation axes, viz. $\textbf{n}_{0}\cdot\textbf{n}_{1}=-1$. For condition (ii) this yields $\cos\phi=\cos(2N\theta)\approx\cos(2Nv_{x})=0$ and thereby $N=(l+\pi/2)/2v_{x}$. Here we choose $k=1$ and $l=0$ to minimize the total measurement time. In summary, the $r_{y}$-measurement conditions are

[TABLE]

Similarly, in the $x$ case, conditions (i) and (ii) read $n_{y}=0$ and $\cos(N\theta)=0$, which is fulfilled for perpendicular rotation axes, $\textbf{n}_{1}\cdot\textbf{n}_{0}=0$. Figure S1 shows three examples of the points $\sin\phi\,\mathbf{n}$ for different pulse-segment numbers (red, blue, and gray).

S.1.2 S.1.2. Measurement of unknown coupling constants

A measurement of the probe-expectation value $\langle\sigma_{x}^{p}\rangle$, which is equal to $\cos\phi$, includes the possibility to extract the coupling constants $a_{z}$ and $a_{x}$ in case they are unknown Grinolds et al. (2014); Zopes et al. (2018). For every value of $N$ the function $\cos\phi$ shows a minimum at the free-evolution time $\tau_{1}=2\pi/(\omega_{1}+\omega_{0})$, which means in principle a single swipe over $\tau$, for a fixed $N$, is sufficient to determine $\tau_{1}$. In a second step the pulse number can be varied for the fixed evolution time $\tau_{1}$ in order to find the $N$ that minimizes $|\cos\phi|$. As in the case for a measurement of the $y$ component this is fulfilled for the integer closest to $\pi/4|v_{x}|$.

This procedure is illustrated in Fig S2 for the values from the main text, viz. $a_{z}=0.015\,\omega_{0}$ and $a_{x}=0.08\,\omega_{0}$. An experimental determination of $\tau_{1}$ and the optimal $N$ allows to infer the coupling constants $a_{z}$ and $a_{x}$ in dependence of the known $\omega_{0}$. These may then be used to locate the dark spin if, for example, the coupling is given by a magnetic dipole-dipole interaction.

S.1.3 S.1.3. Robustness against noise

In this section, we demonstrate the robustness of the proposed scheme against noise acting on the probe. We exemplify its robustness at the experimental implementation for the state tomography of a 13C nuclear spin weakly coupled to an NV-center probe. The main source of noise in this solid-state system comes from the spin bath surrounding the NV center. An established description of the collective effect of this complex environment on the probe is that of a fluctuating magnetic field Dobrovitski et al. (2008, 2009); de Lange et al. (2010); Maze et al. (2012); Wang and Takahashi (2013), resulting in the noise Hamiltonian $H_{\rm noise}=(b(t)/2)\sigma_{z}^{p}$. Here, the noise $b(t)$ is a random variable obeying a zero-mean Gaussian distribution with the autocorrelation $\langle b(t)b(0)\rangle=b_{0}^{2}\exp(-t/t_{b})$, where $b_{0}^{2}$ denotes the variance and $t_{b}$ the correlation time of the noise. The time evolution of this random variable can be modelled as an Ornstein-Uhlenbeck process Van Kampen (2007) with the update formula as follows Gillespie (1991, 1996a, 1996b)

[TABLE]

with a normally distributed random variable $\mathcal{N}$ and the time discretization $\Delta t$. Duo to the symmetrical structure of the pulse sequence, our scheme will be robust against slow noise, i.e., noise with a correlation time $t_{b}$ which is long compared to the free evolution time $\tau$ between pulses. This is evident in Fig. S3, where we show simulations for a $\langle\sigma_{y}\rangle$ measurement of the dark spin under the influence of the noise on the probe. Here, we use the parameters given in the main text and perform simulations for the three typical noise correlation times $t_{b}=(0.2,0.5,1)$ ms Dobrovitski et al. (2008, 2009); de Lange et al. (2010); Maze et al. (2012); Wang and Takahashi (2013) and different noise variances $b_{0}/2\pi\in[9,112]$ kHz, which are related to the probe coherence time $T_{2p}^{\ast}$, as measured in free-induction decay experiments, via $T_{2p}^{\ast}=\sqrt{2}/b_{0}$ Dobrovitski et al. (2008, 2009); de Lange et al. (2010); Maze et al. (2012); Wang and Takahashi (2013). In our simulations, we assume the 13C spin to be initially in the state $\rho=(\mathbb{1}+\mathbf{r}\cdot\bm{\sigma})/2$ with $\mathbf{r}=(0,0.4,0)$ and average the random process over 1000 realizations. The parameters are the ones from the main text, $N=10$ and $\tau=3\ \mu$s, leading to a total measurement time of $60\ \mu$s, which lies well below probe coherence times $T_{2p}$ achieved with pulsed dynamical decoupling de Lange et al. (2010); Naydenov et al. (2011).

S.1.4 S.1.4. Two interacting spin-$\frac{1}{2}$

As a further example for the applicability of the proposed scheme, we consider a dark system formed by two spin-1/2, whose Hamiltonian is given by

[TABLE]

with the Pauli operators $\sigma_{\lambda}^{(k)}$, for $k=1,2$ and $\lambda=x,y,z$, and the inter-dark-spin coupling strength $A_{x}$. Analogous to the single-spin case in the main text, the probe-state-dependent potentials are assumed to be $H_{0}=0$ and

[TABLE]

If we now consider the weak-coupling regime, where $a_{z}^{(k)},a_{x}^{(k)},A_{x}\ll\omega_{0}$ is fulfilled, the transitions in the manifold $\{|0\rangle_{1}|0\rangle_{2},|1\rangle_{1}|1\rangle_{2}\}$ of the bipartite dark system can be neglected due to their large energy separation. We therefore focus on the subspace spanned by the two states $|\tilde{0}\rangle=|0\rangle_{1}|1\rangle_{2}$ and $|\tilde{1}\rangle=|1\rangle_{1}|0\rangle_{2}$, which can be considered as a pseudo spin with the Pauli vector $\tilde{\bm{\sigma}}$. In this subspace, the resulting state-dependent potentials, in the interaction picture with respect to the probe Hamiltonian, are then simply given by

[TABLE]

with $A_{z}=a_{z}^{(1)}-a_{z}^{(2)}$. For the probe-measurement outcome we find $\langle\sigma_{y}^{p}\rangle=-\sin\phi\,\tilde{\mathbf{n}}\cdot\langle\bm{\tilde{\sigma}}\rangle$, where we can apply Eq. (5) of the main text, with $\tilde{\mathbf{n}}=(n_{z},n_{y},n_{x})$ and the substitutions $\alpha=A_{x}\tau/2$, $\beta=A\tau/2$, $A^{2}=A_{x}^{2}+A_{z}^{2}$, $v_{x}=A_{z}/A$, and $v_{z}=A_{x}/A$. The Bloch vector $\tilde{\mathbf{r}}=\langle\bm{\tilde{\sigma}}\rangle$ of the pseudo spin in this subspace can be determined by varying the pulse-sequence parameters $\tau$ and $N$, as described in the main text. Our strategy thereby yields a feasible method to measure correlations between the two dark spins, since the Bloch vector components $\tilde{r}_{x}$ and $\tilde{r}_{y}$ correspond to

[TABLE]

which may be used as entanglement witnesses.

S.2 S.2. Harmonic oscillator

S.2.1 S.2.1. Basic principle

For a dark harmonic oscillator we show the derivation of the explicit form of the the operator $U_{0}^{\dagger}$, the operator $U_{1}$ can be derived along the same lines. Defining $\epsilon=g/2\nu$ the probe-state-dependent potentials $V_{0}$ and $V_{1}$ can be displaced according to

[TABLE]

The multiplication property of displacement operators, $D(x)D(y)=\exp(i\,{\rm Im}\,xy^{\ast})D(x+y)$, and the identity $\exp(ixa^{\dagger}a)D(y)\exp(-ixa^{\dagger}a)=D(y\exp(ix))$ allows to write the single pulse-segment evolution operator $u_{0}$ in the form

[TABLE]

with $\varphi_{1}=2\epsilon^{2}\nu\tau$. By Hermitian conjugating and taking the operator to the $N$th power we can write $U_{0}^{\dagger}=(u_{0}^{\dagger})^{N}$ as

[TABLE]

where $\varphi_{2}$ incorporates $-N\varphi_{1}$ and an additional phase that we abstain from writing explicitly, since the operator $U_{1}$ carries exactly the opposite one. In a next step we rewrite the powers using the identities

[TABLE]

where we also do not need to determine the phase $\vartheta$. In these identities, the first equality can be easily proven by induction and the second one follows from the multiplication property of the displacement operator and the geometric sum. This leaves us with

[TABLE]

where $\varphi_{3}$ includes $\varphi_{2}$ and the phase corresponding to $\vartheta$ from identity S13. The quantity $\zeta$ mentioned in the main text thereby has the form

[TABLE]

and the relevant unitary $U_{0}^{\dagger}U_{1}$ finally reads

[TABLE]

where the arguments of the displacement operator are the curves $\xi(\tau,N)$.

S.2.2 S.2.2. Reconstruction example

In the Fock basis, the matrix elements $\langle n|\rho|m\rangle$ of a density operator $\rho$, with the associated characteristic function $\chi(\xi)={\rm Tr}\{D(\xi)\rho\}$, can be written as

[TABLE]

In order to calculate their explicit value, one can then use the form of the displacement-operator matrix elements Cahill and Glauber (1969), which is given by

[TABLE]

with the generalized Laguerre polynomials $L$. This expression is valid for $n\geq m$, while for $n<m$ one has to replace $\eta$ by $-\eta^{\ast}$ as well as exchange $n$ and $m$ on the right-hand side of Eq. (S18).

Here, we demonstrate the state reconstruction using the matrix elements from Eq. (S17) obtained from an interpolated characteristic function. As an example, we assume a dark harmonic oscillator initially prepared in the Fock state $|n\rangle$. The characteristic function of this state has the form $\chi_{n}(\xi)=L_{n}^{(0)}\left(|\xi|^{2}\right)\exp(-|\xi|^{2}/2)$, as can be see from Eq. (S18) for $n=m$. In Fig. S4(a) we show $\chi_{1}(\xi)$ [gray] with the contour lines $\xi(\tau,N)$ [blue] for $N=(1,...,20)$ and $g/\nu=3/40$. In our case $\chi(\xi)$ can be interpolated from samples of the exact characteristic functions along these curves, and Eqs. (S17) and (S18) can then be used to calculate the matrix elements in the Fock basis. To show the performance of such a reconstruction, we increase the maximum number $\tilde{N}$ of pulse-sequence segments and reconstruct a density operator $\rho_{n,\tilde{N}}$ from interpolating over all curves $\xi(\tau,N)$ with $N\leq\tilde{N}$. As a measure for the accuracy of $\rho_{n,\tilde{N}}$ we employ its trace distance $d_{n,\tilde{N}}$ to the exact density matrix $\rho_{n}=|n\rangle\langle n|$ of the Fock state, which is given by $d_{n,\tilde{N}}={\rm Tr}\{\sqrt{(\rho_{n,\tilde{N}}-\rho_{n})^{2}}\}/2$. Figure S4(b) shows this trace distance for $\tilde{N}=(1,...,20)$, $g/\nu=3/40$, and three values of $n$.