Supplementary Material for "Theory of Metastability in Discrete-Time Open Quantum Dynamics"

Abstract

In this supplementary material, we first provide more details about fundamentals of quantum channels, and properties of fixed points and extreme metastable states (EMSs) in Ramsey interferometry measurements (RIMs). Then we present Monte Carlo simulation results for a target system composed of multiple qubits, with the ancilla qubit under both RIMs and dynamical decoupling (DD) sequences. Moreover, we numerically demonstrate that quantum metastability is robust even when the target system suffers additional dissipation.

I Fundamentals of quantum channels

I.1 Definition of quantum channels

Denote the linear operators acting on a Hilbert space ${\mathcal{H}}$ as ${\mathcal{B(H)}}$, then a map $\Phi:{\mathcal{B(H)}}\to\mathcal{B(H)}$ is a quantum channel if it satisfies the following conditions [1]

• Linear map: For any $A,B\in{\mathcal{B(H)}}$ and complex number $c$, $\Phi(A+cB)=\Phi(A)+cT(B)$;

• Trace preserving: For any $A\in{\mathcal{B(H)}}$, $\Tr[\Phi(A)]=\Tr(A)$. This implies unitality of $\Phi^{\dagger}$, i.e., $\Phi^{\dagger}(\mathbb{I})=\mathbb{I}$, where $\Phi^{\dagger}$ is defined by $\Tr[B\Phi(A)]=\Tr[\Phi^{\dagger}(B)A]$;

• Completely positive: For any positive operator $\rho$, $(\Phi\otimes\mathbb{I})(\rho)$ is still a positive operator, with $\mathbb{I}$ is the identity operator on an additional system with dimension ${\rm dim}({\mathcal{H}})$.

So a quantum channel is a completely positive and trace-preserving (CPTP) linear map, which maps a state to another one by $\Phi(\rho)=\rho^{\prime}$.

I.2 Representations of quantum channels

Every quantum channel has four different representations: the Kraus representation, the Stinespring representation, the natural representation, and the Choi representation. In this paper, we use the first three representations of quantum channels.

I.2.1 Kraus representation

Of the four representations, the Kraus representation is the most commonly used one. In this representation, a quantum channel is fully characterized by a collection of Kraus operators $\{M_{\alpha}\}_{\alpha=1}^{r}$ satisfying $\sum_{\alpha}M^{\dagger}_{\alpha}M_{\alpha}=\mathbb{I}$ so that

$$\Phi(\cdot)=\sum_{\alpha=1}^{r}M_{\alpha}(\cdot)M_{\alpha}^{\dagger}=\sum_{\alpha=1}^{r}\mathcal{M}_{\alpha}(\cdot),$$ (S1)

where $\mathcal{M_{\alpha}}(\cdot)=M_{\alpha}(\cdot)M_{\alpha}^{\dagger}$ is a superoperator.

I.2.2 Stinespring representation

The stinespring representation is a dilation of a quantum channel. The dilation can be realized by coupling the target system to an ancilla system, and letting the composite system undergoing a unitary evolution and then tracing over the ancilla system,

$$\Phi(\rho)=\Tr_{a}[U(\rho_{a}\otimes\rho)U^{\dagger}],$$ (S2)

where $\rho_{a}$ is the initial state of the ancilla system and $U$ is a unitary of the composite system and ${\rm Tr}_{a}$ denotes the trace over the ancilla system. For a $r$-dimensional ancilla system with an orthonormal basis $\{|\alpha\}_{\alpha=1}^{r}$ and the initial state $\rho_{a}=|1\rangle_{a}\langle 1|$, the Kraus operator can be easily obtained as $M_{\alpha}=\langle r|U|1\rangle_{a}$.

I.2.3 Natural representation

For sequential applications of the same channel, it is convenient to use the natural representation of a quantum channel. To understand this representation, we need to introduce Hilbert-Schmidt (HS) space. In the HS space, an operator on a $d$-dimensional Hilbert space (represented by a $d\times d$ matrix) is transformed to a $d^{2}\times 1$ vector,

$$A=\matrixquantity[a_{11}&\dots&a_{1d}\\ \vdots&\ddots&\vdots\\ a_{d1}&\dots&a_{dd}]=\matrixquantity[{\bf\it a}_{1}\\ \vdots\\ {\bf\it a}_{d}]\Rightarrow\left.\ket{A}\right\rangle=\matrixquantity[\,{\bf\it a}^{T}_{1}\,\\ \vdots\\ \,{\bf\it a}^{T}_{d}\,],$$ (S3)

where the superscript $T$ denotes the matrix transposition. Such a transformation can also be represented by $A=\sum_{ij}a_{ij}\ket{i}\bra{j}$$\Rightarrow$$\left.\ket{A}\right\rangle=\sum_{ij}a_{ij}\left.\ket{ij}\right\rangle$ with $\left.\ket{ij}\right\rangle=\ket{i}\otimes\ket{j}$. The inner product in the HS space is defined as

$$\left\langle\innerproduct{A}{B}\right\rangle =\sum_{ijpq}a^{*}_{ij}b_{pq}\left\langle\innerproduct{ij}{pq}\right\rangle$$ (S4)

$$=\sum_{ij}a^{*}_{ij}b_{ij}$$

$$=\sum_{j}(A^{\dagger}B)_{jj}$$

$$=\Tr(A^{\dagger}B).$$

Then a superoperator (assuming $\mathcal{O}(\rho)=A\rho B$) becomes a single matrix $A\otimes B^{T}$ acting on a vector $\left.\ket{\rho}\right\rangle$ in HS space,

$$A\rho B=\sum_{ijpq}a_{ij}\rho_{jk}b_{pq}\ket{i}\bra{q}\,\,\Rightarrow\,\,A\otimes B^{T}\left.\ket{\rho}\right\rangle,$$ (S5)

so the natural representation of the channel in Eq. (S1) in HS space is

$$\hat{\Phi}=\sum_{\alpha=1}^{r}\hat{\mathcal{M}}_{\alpha}$$ (S6)

where $\hat{\mathcal{M}}_{\alpha}=M_{\alpha}\otimes M^{*}_{\alpha}$. Note that we add hats on operators acting on HS space.

I.3 Jordan decomposition of a quantum channel

The natural representation of a quantum channel on the HS space is a $d^{2}\times d^{2}$ square operator. This operator may not be diagonalizable, but can always be converted to a Jordan normal form as

$$\hat{\Phi} =S\left(\bigoplus_{k=1}^{\kappa}\mathcal{J}_{d_{k}}(\lambda_{k})\right)S^{-1}$$ (S7)

$$=S\left(\sum_{\left|{\lambda}_{j}\right|=1}{\lambda}_{j}{\mathcal{P}}_{j}+\sum_{\left|{\lambda}_{k}\right|<1}({\mathcal{P}}_{k}+{\mathcal{N}}_{k})\right)S^{-1},$$

where $S$ is an invertible $d^{2}\times d^{2}$ matrix, $\kappa$ is the total number of Jordan blocks, and $\mathcal{J}_{d_{i}}(\lambda_{i})$ represents a $d_{i}$-dimensional Jordan block corresponding to the eigenvalue $\lambda_{i}$, ${\mathcal{P}}_{j}$ is a projection operator and ${\mathcal{N}}_{k}$ is a nilpotent operator satisfying ${\mathcal{N}}_{k}^{d_{k}}=0$. Considering that $\sum_{k}^{\kappa}d_{k}=d$, then the channel is diagonalizable iff $\kappa=d$. Note that the Jordan blocks corresponding to the fixed points or rotating points (with eigenvalues $\absolutevalue{\lambda_{i}}=1$) are all rank-one projectors .

If the channel is diagonalizable, we have

$$\hat{\Phi} =\sum_{i}\lambda_{i}|R_{i}\rangle\rangle\langle\langle L_{i}|,$$ (S8)

where $\{|R_{i}\rangle\rangle,|L_{i}\rangle\rangle\}$ is a complete biorthogonal basis satisfying $\left\langle\innerproduct{L_{i}}{R_{j}}\right\rangle=\delta_{ij}$. The trace-preserving property of $\hat{\Phi}$ implies the unitality of $\hat{\Phi}^{\dagger}$, i.e., $\hat{\Phi}^{\dagger}\left.\ket{\mathbb{I}}\right\rangle=\left.\ket{\mathbb{I}}\right\rangle$ or there exists a left eigenvector $\langle\langle L_{i}|=\langle\langle\mathbb{I}|$ for eigenvalue 1, then for right eigenvectors corresponding to $\absolutevalue{\lambda_{j}}<1$, we have

$$0=\left\langle\innerproduct{L_{i}}{R_{j}}\right\rangle={\Tr}(L_{i}^{\dagger}R_{j})=\Tr(R_{j}),$$ (S9)

I.4 Sequential quantum channels with rotating points

If the channel has $n$ eigenvectors corresponding to eigenvalue $\absolutevalue{\lambda}=1$, in which there are $n_{0}$ fixed points and $n-n_{0}$ rotating points with $\lambda_{j}={e}^{{i}\varphi_{j}}$, and $l-n$ eigenvectors with $\lambda_{k}=\absolutevalue{\lambda_{k}}{e}^{{i}\varphi_{k}}\approx 1$, then after sequentially applying the quantum channel for $m$ times, we have

$$\hat{\Phi}^{m}\left.\ket{\rho}\right\rangle =\sum_{i=1}^{n_{0}}c_{i}\left.\ket{\rho_{\rm fix}^{i}}\right\rangle+\sum_{j=n_{0}+1}^{n}c_{j}\lambda_{j}^{m}|{R_{\rm rot}^{j}}\rangle\rangle+\sum_{k=n+1}^{l}c_{k}\lambda_{k}^{m}\left.\ket{R_{k}}\right\rangle+\cdots$$ (S10)

$$\simeq\sum_{i=1}^{n}c_{i}\left.\ket{\rho_{\rm fix}^{i}}\right\rangle+\sum_{j=n_{0}+1}^{n}c_{j}{e}^{{i}m\varphi_{j}}|{R_{\rm rot}^{j}}\rangle\rangle\sum_{k=n+1}^{l}c_{k}{e}^{m(\ln{\absolutevalue{\lambda_{k}}+{i}\varphi_{k})}}\left.\ket{R_{k}}\right\rangle,$$

where we truncate the equation as we did in the main text.

Then, in the metastable region, ${e}^{m\ln{\absolutevalue{\lambda_{k}}}}\approx 1$, by considering conjugate pairs of complex eigenvalues and absorbing the phase in $c$, we have

$$\hat{\Phi}^{m}\left.\ket{\rho}\right\rangle\simeq\sum_{i=1}^{n_{0}}c_{i}\left.\ket{\rho_{\rm fix}^{i}}\right\rangle+\sum_{j=n_{0}+1}^{l}c_{j}^{\prime}(m)\left.\ket{R_{j}^{\prime}}\right\rangle,$$ (S11)

which is approximately the same as Eq.(4) in the main text, except that the summation of metastable points also include the rotating points.

II Fixed points of the channel induced by RIMs

Now we consider the fixed points of quantum channel on a target system induced by a RIM of an ancilla qubit. The anciila is coupled to the target system with a pure-dephasing Hamiltonian,

$$H=\sigma_{q}^{z}\otimes B+\gamma\,\mathbb{I}_{q}\otimes C,$$ (S12)

where $\sigma_{q}^{i}$ is the Pauli-$i$ operator of the ancilla, $B$ and $C$ are both operators on the target system and $\gamma$ controls the magnitude of the second term. For a RIM, the Kraus representation of the channel on the target system is

$$\hat{\Phi}=\hat{\mathcal{M}}_{0}+\hat{\mathcal{M}}_{1}=(\hat{\mathcal{U}}_{0}+\hat{\mathcal{U}}_{1})/2,$$ (S13)

where $\mathcal{\hat{U}_{\alpha}}=U_{\alpha}\otimes U_{\alpha}^{*}$ with $U_{\alpha}=e^{-i[(-1)^{\alpha}B+\gamma C]}$, and $\mathcal{\hat{M}_{\alpha}}=M_{\alpha}\otimes M_{\alpha}^{*}$ with the Kraus operator $M_{\alpha}=[U_{0}-(-1)^{\alpha}e^{i\Delta\phi}U_{1}]/2$ and $\Delta\phi=\phi_{1}-\phi_{2}$. Then the fixed points induced by such a unital channel is given by the following proposition.

Proposition 1. The fixed points of the channel in Eq. (S13) depends on the commutativity of $B$ and $C$. If $[B,C]=0$, the fixed points are spanned by a set of rank-one projections $\{|j\rangle\langle j|\}_{j=1}^{d}$; if $[B,C]\neq 0$, the fixed points are spanned by a set of projection operators $\{\Pi_{j}\}_{j=1}^{r}$ ($r\leq d$), satisfying $\sum_{j=1}^{r}\Pi_{j}=\mathbb{I}$.

Proof.

It has been proven that $\rho$ is a fixed point of a unital channel if and only if it commutes with every Kraus operator [3], i.e., $[\rho,M_{\alpha}]=0$ for any $\alpha$. This implies that $[\rho,U_{0}]=[\rho,U_{1}]=0$. If the above condition is always satisfied for any $\alpha$, then $[\rho,B]=[\rho,C]=0$.

If $[B,C]=0$, then $B$ and $C$ can be diagonalized simultaneously, $B=\sum_{j=1}^{d}b_{j}\ket{j}\bra{j}$ and $C=\sum_{j=1}^{d}c_{j}\ket{j}\bra{j}$. So the fixed points must include the rank-one projections $\{|j\rangle\langle j|\}_{j=1}^{d}$ and their linear combinations.

If $[B,C]\neq 0$, we can block diagonalize them simultaneously by unitary transformation,

$$B=W\left(\bigoplus_{j=1}^{r}B_{j}\right)W^{\dagger},\quad C=W\left(\bigoplus_{j=1}^{r}C_{j}\right)W^{\dagger}$$ (S14)

where $r\leq d$ is the number of blocks (with equality occurring only when $[B,C]=0$ and all of blocks are one-dimensional), $W$ is unitary matrix and should be chosen so that $B_{j}$ and $C_{j}$ for any $j$ cannot be reduced further to have more blocks. There must be at least one subspace ${\mathcal{H}}_{j}$ in which $[B_{j},C_{j}]\neq 0$ to make $[B,C]\neq 0$. Such a block diagonalization partitions the Hilbert space of the target system into the direct sum of $r$ subspaces $\mathcal{H}=\bigoplus_{j=1}^{r}\mathcal{H}_{j}$, and $[B_{j},{C}_{j}]\neq 0$ for at least one subspace $\mathcal{H}_{j}$ with ${\rm dim}(\mathcal{H}_{j})\geq 2$. Thus the Kraus operator is also transformed to a block-diagonal form as $M_{\alpha}=\oplus_{j=1}^{r}M_{\alpha}^{j}$. Then the fixed points must include the set of projections $\{\Pi_{j}\}_{j=1}^{r}$ ($r\leq d$) and their linear combinations, where $\Pi_{j}=WP_{j}W^{-1}$ with $P_{j}$ being the projector to $\mathcal{H}_{j}$. Note that the case $[B,C]=0$ can be regarded as a special case of $[B,C]\neq 0$.

Now we prove that there are no other fixed points for the case $[B,C]\neq 0$, where there is at least one block with $[B_{j},C_{j}]\neq 0$ and $[M_{0}^{j},M_{1}^{j}]\neq 0$. Suppose there is another density matrix satisfying $[\rho^{\prime},M_{\alpha}^{j}]=0$. If $\rank{\rho^{\prime}}=d_{j}$, then $[\rho^{\prime},M_{0}^{j}]=[\rho^{\prime},M_{1}^{j}]=0$. Since the positive operator $\rho^{\prime}$ can be diagonalized, this implies that $[M_{0}^{j},M_{1}^{j}]=0$. If $\rank(\rho^{\prime})<d_{j}$, then formulate another fixed point $\rho^{\prime\prime}=\rho^{\prime}+\eta\mathbb{I}$ with $\eta$ being a positive number such that $\rank(\rho^{\prime\prime})=d_{j}$, then the proof is similar to the former case.

∎

III Quantum metastability in sequential RIMs for a target qubit

III.1 Construction of EMSs

If target system is a single qubit with $B=\sigma_{z}$ and $C=0$, the fixed points are spanned by $\left.\ket{00}\right\rangle$ and $\left.\ket{11}\right\rangle$, where $\{|0\rangle,|1\rangle\}$ are eigenstates of $\sigma_{z}$. When there exists a small perturbing Hamiltonian $\gamma C=\gamma\sigma_{x}$, the channel has a single fixed point $\left.\ket{\rho_{\rm fix}}\right\rangle=\left.\ket{\mathbb{I}}\right\rangle/2$, and another metastable point $\left.\ket{R_{2}}\right\rangle$. So the dimension of MM is $2-1=1$. In the metastable region, $\lambda_{2}^{m}\approx 1$ (assuming $\lambda_{2}$ is real), then a metastable state is

$$\left.\ket{\rho_{\rm MS}}\right\rangle=\left.\ket{\rho_{\rm fix}}\right\rangle+c_{2}\left.\ket{R_{2}}\right\rangle,$$ (S15)

with $c_{2}={\Tr}(L_{2}^{\dagger}\rho)$. Thus the metastable state is determined only by $c_{2}$, and the extremal points in the one-dimensional MM are

$$\left.\ket{\rho_{1}}\right\rangle=\left.\ket{\rho_{\rm fix}}\right\rangle+c_{2}^{M}/h\left.\ket{R_{2}}\right\rangle,\quad\left.\ket{\rho_{2}}\right\rangle=\left.\ket{\rho_{\rm fix}}\right\rangle+c_{2}^{m}/h\left.\ket{R_{2}}\right\rangle,$$ (S16)

where $c_{2}^{M}$ ($c_{2}^{m}$) is the maximal (minimal) value of $c_{2}$, and $h=\sqrt{\left\langle\innerproduct{L_{2}}{L_{2}}\right\rangle\left\langle\innerproduct{R_{2}}{R_{2}}\right\rangle}$ is a normalization coefficient that will be derived below. For real $\lambda_{2}$, $R_{2}$ and $L_{2}$ can be chosen to be Hermitian [4], then $c_{2}=\Tr(\rho L_{2})$ and $c_{2}^{M}$ and $c_{2}^{m}$ are two eigenvalues of $L_{2}$. So we have $L_{2}=c_{2}^{M}\ket{M}\bra{M}+c_{2}^{m}\ket{m}\bra{m}$ with $\{\ket{M},\ket{m}\}$ forming a complete orthonormal basis.

Let $\rho=k_{1}\ket{M}\bra{M}+k_{2}\ket{m}\bra{m}$ with $k_{1}+k_{2}=1$, then $\rho L_{2}=c_{2}^{M}k_{1}\ket{M}\bra{M}+c_{2}^{m}k_{2}\ket{m}\bra{m}$. This means that only when $k_{1}=1,\,k_{2}=0$ (or $\rho=\ket{M}\bra{M}$), we have $c_{2}=c_{2}^{M}$. Since $\rho_{1,2}$ should be pure states, this means $\Tr(\rho_{1,2}^{2})=1$,

$$1=\Tr(\rho_{1}^{2}) =\Tr(\mathbb{I}^{2}/4+c_{2}^{M}/hR_{2}+(c_{2}^{M}/h)^{2}R_{2}^{2})$$ (S17)

$$=1/2+(c_{2}^{M}/h)^{2}\Tr(R_{2}^{2})$$

$$\Downarrow$$

$$1/2 =(c_{2}^{M}/h)^{2}\Tr(R_{2}^{2}).$$

We can derive a similar expression for $\rho_{2}$, i.e., $1/2=(c_{2}^{m}/h)^{2}\Tr(R_{2}^{2})$, so

$$h^{2}=[(c_{2}^{M})^{2}+(c_{2}^{m})^{2}]\Tr(R_{2}^{2})=\Tr(L_{2}^{2})\Tr(R_{2}^{2})=\left\langle\innerproduct{L_{2}}{L_{2}}\right\rangle\left\langle\innerproduct{R_{2}}{R_{2}}\right\rangle.$$ (S18)

We can illustrate the above analysis with a simple example. Consider a fixed point $\left.\ket{\rho_{\rm fix}}\right\rangle=\left.\ket{\mathbb{I}}\right\rangle/2=(\left.\ket{00}\right\rangle+\left.\ket{11}\right\rangle)/2$, and a metastable point with $\left.\ket{R_{2}}\right\rangle=(\left.\ket{00}\right\rangle-\left.\ket{11}\right\rangle)/2$, then $\left.\ket{L_{2}}\right\rangle=\left.\ket{00}\right\rangle-\left.\ket{11}\right\rangle$, where $c_{2}^{M}=1,\,c_{2}^{m}=-1$, and $h=\sqrt{\left\langle\innerproduct{L_{2}}{L_{2}}\right\rangle\left\langle\innerproduct{R_{2}}{R_{2}}\right\rangle}=1$. Thus the EMSs are

$$|\rho_{1}\rangle\rangle=|\rho_{\rm fix}\rangle\rangle+|R_{2}\rangle\rangle=|00\rangle\rangle\langle\langle 00|,\quad|\rho_{2}\rangle\rangle=|\rho_{\rm fix}\rangle\rangle-|R_{2}\rangle\rangle=|11\rangle\rangle\langle\langle 11|,$$ (S19)

For general cases, the EMSs are the eigenstates of $B$ up to some corrections.

III.2 Construction of metastable states

Any metastable state can be represented by a mixture of EMSs,

$$\left.\ket{\rho_{\rm MS}}\right\rangle=p_{1}\left.\ket{\rho_{1}}\right\rangle+p_{2}\left.\ket{\rho_{2}}\right\rangle$$ (S20)

where $p_{1,2}=\left\langle\innerproduct{P_{1}}{\rho}\right\rangle=\Tr(\rho P_{1,2})$ with $P_{1,2}$ being observables satisfying $\langle\langle{P_{v}}|{\rho_{u}}\rangle\rangle=\delta_{v,u}$, $P_{1}+P_{2}=\mathbb{I}$ and $P_{1,2}\geq 0$. Below we show how to construct such observables.

Since $\left.\ket{\rho_{\rm MS}}\right\rangle=\Phi^{m}\left.\ket{\rho}\right\rangle$ in the metastable region $\mu^{\prime\prime}\ll m\ll\mu^{\prime}$, so we have

$$\Phi^{m}\simeq|\rho_{1}\rangle\rangle\langle\langle P_{1}|+|\rho_{2}\rangle\rangle\langle\langle P_{2}|,$$ (S21)

Comparing the above equation with $\Phi^{m}\simeq\frac{1}{2}|\mathbb{I}\rangle\rangle\langle\langle\mathbb{I}|+|R_{2}\rangle\rangle\langle\langle L_{2}|$ in metastable region, we have

$$\frac{1}{2}|\mathbb{I}\rangle\rangle\langle\langle\mathbb{I}|+|R_{2}\rangle\rangle\langle\langle L_{2}| =\left(\frac{1}{2}|\mathbb{I}\rangle\rangle+\frac{c_{2}^{M}}{h}|R_{2}\rangle\rangle\right)\langle\langle P_{1}|+\left(\frac{1}{2}|\mathbb{I}\rangle\rangle+\frac{c_{2}^{m}}{h}|R_{2}\rangle\rangle\right)\langle\langle P_{2}|$$ (S22)

$$\Downarrow$$

$$|R_{2}\rangle\rangle\langle\langle L_{2}| =\frac{c_{2}^{M}}{h}|R_{2}\rangle\rangle\langle\langle P_{1}|+\frac{c_{2}^{m}}{h}|R_{2}\rangle\rangle(\langle\langle\mathbb{I}|-\langle\langle P_{1}|)$$

$$=|R_{2}\rangle\rangle\left[\frac{c_{2}^{M}}{h}\langle\langle P_{1}|+\frac{c_{2}^{m}}{h}(\langle\langle\mathbb{I}|-\langle\langle P_{1}|)\right],$$

so we have

$$P_{1}=\frac{hL_{2}-c_{2}^{m}\mathbb{I}}{\Delta c_{2}},\quad P_{2}=\frac{-hL_{2}+c_{2}^{M}\mathbb{I}}{\Delta c_{2}},$$ (S23)

with $\Delta c_{2}=c_{2}^{M}-c_{2}^{m}$.

IV Quantum metastability in sequential quantum channels for multiple target qubits

In this section, we consider a practical example to show quantum metastability in sequential quantum channels, that is, an NV center electron spin (ancilla qubit) in a ${}^{13}$C nuclear spins bath (target system). The coupling Hamiltonian has a form similar to Eq. (S12),

$$H=\sigma_{q}^{z}\otimes B+C,$$ (S24)

with

$$B=f(t)\sum_{k=1}^{K}{\bf\it A}_{k}\cdot{\bf\it I}_{k},\quad C=\omega_{L}\sum_{k=1}^{K}I_{k}^{z}+\sum_{k<j}D_{jk}\left[\bm{I}_{k}\cdot\bm{I}_{j}-\frac{3\left(\bm{I}_{k}\cdot\bm{r}_{kj}\right)\left(\bm{I}_{j}\cdot\bm{r}_{kj}\right)}{r_{kj}^{2}}\right],$$ (S25)

where $f(t)$ is a modulation function accounting for possible DD control of the ancilla qubit, ${\bf\it A}_{k}=(A_{k}^{x},A_{k}^{y},A_{k}^{z})$ is the hyperfine interaction vector of the $k$${\rm th}$ nuclear spin, ${\bf\it I}_{k}$ is the nuclear spin operator vector, $d=2^{K}$ is the dimension of Hilbert space of target system, $\omega_{L}=\gamma_{n}B_{z}$ is Larmor precession frequency of nuclear spins, $D_{jk}=\frac{\mu_{0}\gamma_{\mathrm{n}}^{2}}{4\pi r_{kj}^{3}}$ with $\gamma_{n}$ is the gyromagnetic ratios of the target spins, ${\bf\it r}_{kj}$ is the displacement from the $i$${\rm th}$ target spin to the $j$${\rm th}$ target spin and $\mu_{0}$ is vacuum permeability. Note that we let $\gamma=1$ and use practical parameters below.

IV.1 RIM sequences

For RIM sequences, we assume that the external magnetic field is zero or very weak, so that

$$B=\sum_{k=1}^{K}{\bf\it A}_{k}\cdot{\bf\it I}_{k},\quad C=\sum_{k<j}D_{jk}\left[\bm{I}_{k}\cdot\bm{I}_{j}-\frac{3\left(\bm{I}_{k}\cdot\bm{r}_{kj}\right)\left(\bm{I}_{j}\cdot\bm{r}_{kj}\right)}{r_{kj}^{2}}\right],$$ (S26)

To show metastable polarization of nuclear spins, we also measure the evolution of fidelity, which is defined as

$$F_{i}(\rho)=\Tr\sqrt{\sqrt{\rho_{i}}\rho\sqrt{\rho_{i}}},$$ (S27)

where $\rho_{i}=\ket{i}\bra{i}$ with $\ket{i}$ being the eigenstate of $B$ with $i\in\{1,2,3,4\}$. Fidelity here can show the distance between the state of system and unperturbed eigenstates. We can see the target system are polarized to nearby with the unperturbed eigenstates, which are our EMSs, and the fidelity plateaus can also confirm the metastability behaviors [Fig. S1]. When there is an external magnetic field involving, metastability also remains if the field is weak (Fig. S2), i.e., $\gamma_{n}B\ll\absolutevalue{A_{k}}$.

Figure S1: Metastability in measurement statistics of sequence RIMs. (a-d)The Monte Carlo simulation for two target qubits, in which we present four stages of the system evolution. (a) When the time of measurements is small, the width of peaks are relatively large and the peaks are not absolutely distinguishable. (b) For larger $m$, the four peaks corresponding to four EMSs are clear. (c) The EMSs relevant to $\left.\ket{22}\right\rangle$ and $\left.\ket{33}\right\rangle$ vanish faster, with that relevant to $\left.\ket{11}\right\rangle$ and $\left.\ket{44}\right\rangle$ are still existing. (d) All EMSs are collapsed. (e) The trajectory of evolution of fidelity, in which the stages in (a-d) are labeled by dashed line. In all Monte Carlo simulations, we use $10^{4}$ samples, parameters are $\Delta\phi=\frac{\pi}{2}$, $\absolutevalue{A_{1}}=0.585\mathrm{kHz},\,\absolutevalue{A_{2}}=0.890\mathrm{kHz},\absolutevalue{D_{12}}=11.6\mathrm{Hz}$, then $||H_{e}||/||V||=0.0316$. (f) Spectra of the quantum channel.

Figure S2: Monte Carlo simulation for two target qubits. The dynamics of polarization and depolarization can be seen from (a) to (d), which proves that our metastability theory also applies to the system involving Zeeman term. In all Monte Carlo simulation, we use $10^{4}$ samples, $\Delta\phi=\frac{\pi}{2}$ is chosen. $\absolutevalue{A_{1}}=0.138\mathrm{MHz},\,\absolutevalue{A_{2}}=0.517\mathrm{MHz},\absolutevalue{D_{12}}=60.5\mathrm{Hz}$, $B=20$G, then $\omega_{L}=13.5$ KHz, and $||H_{e}||/||V||=0.0412$.

IV.2 DD sequences

DD sequences are a generalization of RIM sequences, with additional $N$ $\pi$-flips of the ancilla spin during each cycle [Fig. S3]. Then we have

$$B=f(t)\sum_{k=1}^{K}{\bf\it A}_{k}\cdot{\bf\it I}_{k},\quad C=\omega_{L}\sum_{k=1}^{K}I_{k}^{z},$$ (S28)

where $f(t)$ is the DD modulation function jumping between $+1$ and $-1$ every time the sensor is flipped by a DD pulse. Here for simplicity, we neglect the dipolar interaction term in $C$, which means these target spins are spatially far away from each other.

For $N$-pulse Carr-Purcell-Meiboom-Gill (CPMG) control, $f(t)=f(t+T)$ with $T=4\tau$. Specifically, $f(t)=-1$ when $\tau+\beta T\leq t<3\tau+\beta T$ with $\beta=0,1,\dots,\frac{N}{2}-1$ and $f(t)=1$ otherwise. Then $f(t)$ can be expanded into a Fourier series as $f(t)=\sum_{n=1}^{\infty}c_{n}\cos(n\omega_{T}t)$, where $c_{n}$ is $n$th Fourier expansion coefficient and $\omega_{T}=\frac{2\pi}{4\tau}$ is the angular frequency. For the weak coupling condition, i.e., $\absolutevalue{{\bf\it A}_{k}}\ll\omega_{L}$, then we can approximate $f(t)$ by keeping only the first term in the Fourier series, $f(t)\approx c_{1}\cos(\omega_{T}t)=\frac{4}{\pi}\cos(\omega_{T}t)$. Then $B$ becomes,

$$B=\frac{4}{\pi}\cos(\omega_{T}t)\sum_{k=1}^{K}{\bf\it A}_{k}\cdot{\bf\it I}_{k}\approx\frac{4}{\pi}\cos(\omega_{T}t)\sum_{k=1}^{K}A^{\perp}_{k}I^{\perp}_{k},$$ (S29)

where $A^{\perp}_{k}=\sqrt{{A_{k}^{x}}^{2}+{A_{k}^{y}}^{2}}$ and $I_{k}^{\perp}=\cos\alpha I_{k}^{x}+\sin\alpha I_{k}^{y}$, and the longitude component $A_{z}I_{z}$ can be neglected in the weak coupling region [5].

Now we move to the rotating frame with respect to $\omega_{T}I_{z}$, by using rotating wave approximation, the effective Hamiltonian becomes time-independent,

$$B=\frac{2}{\pi}\sum_{k=1}^{K}A^{\perp}_{k}I^{\perp}_{k},\quad C=\Delta_{\omega}\sum_{k=1}^{K}I_{k}^{z},$$ (S30)

where $\Delta_{\omega}=\omega_{L}-\omega_{T}$ is the detuning away from Larmor precession frequency. Then we can see from Eq.(S30) that, if we choose $\omega_{T}=\omega_{L}$, also called resonance condition, then $C=0$ and the target qubits can be polarized. However, practically the adjustment of frequency may not be precise. With a small detuning, quantum metastability theory shows that the target can still be polarized with appropriate rounds of DD, as numerically verified in Fig. S4.

Figure S3: Schematic of an $N$-pulse CPMG sequence with $N$ $\pi-$flips at $t=\tau,3\tau,\cdots,(2N-1)\tau$.

Figure S4: Monte Carlo simulation for two target spins with DD sequences and growing time of measurements $m$. As $m$ grows, the peaks emerge and become apparent (a-b) since they are corresponding to the four EMSs, then they vanish gradually as $m$ approaching the boundary of metastable region (c), which finally collapse to the maximally mixed state (d). In all Monte Carlo simulation, we use $10^{4}$ samples, parameters are $\Delta\phi=\frac{\pi}{2}$, $\absolutevalue{A_{1}}=5\mathrm{kHz},\,\absolutevalue{A_{2}}=6\mathrm{kHz}$, $B=200$G, $\Delta_{\omega}=10^{-3}\omega_{L}$, $N=32$.

V Quantum metastability for dissipative target system

In this section, we show that quantum metastability is robust even when the target system suffers additional dissipations. Suppose that the target system (a single qubit) suffers dephasing and relaxation noise, the evolution of the composite systems can be described by the following Lindblad master equation,

$$\derivative{\rho_{\rm tot}}{t}=-i[H,\rho_{\rm tot}]+\sum_{k}\Gamma_{k}\left(L_{k}\rho_{\rm tot}L_{k}^{\dagger}-\frac{1}{2}\left\{L_{k}^{\dagger}L_{k},\rho_{\rm tot}\right\}\right),$$ (S31)

where $\rho_{\rm tot}$ is the density matrix of the composite system, $H=\sigma_{z}\otimes B+\gamma\,\mathbb{I}_{q}\otimes C$ , $L_{1}=\sigma_{z}$ denotes the target dephasing, $L_{2}=\sigma^{-}=|1\rangle_{q}\langle 0|$ denotes the target relaxation, and $\Gamma_{k}$ is the dissipation rate.

We perform Monte Carlo simulations for each kind of dissipation separately, in order to examine their influence. The results show that the dephasing noise of target system does not influence the measurement statistics. However, the relaxation of target system maps $\ket{0}$ to $\ket{1}$, which makes the peak corresponding to $\left.\ket{00}\right\rangle$ (when $\gamma=0$) or that corresponding to EMS perturbed from $\left.\ket{00}\right\rangle$ (when $\gamma\neq 0$) transfer to the other peak. This does not influence our discussion about metastability, if the dissipate rate is very small or when $m$ is relatively small (see upper panel of Fig. S5(b)).

Figure S5: Monte Carlo simulations for a target qubit with (a) dephasing noise and (b) relaxation noise. (a) is like Fig. 2 in our main text, showing that there is no influence in our discussion. In (b) it can be seen that the dissipation of the target system remains independent of the conclusion of the metastable, both for the case of $\gamma=0$ (left panel) and for $\gamma=0.05$ (right panel), the peak on the left side is gradually vanished as the round of the measurements increasing due to $\sigma^{-}$ makes $\ket{0}$ into $\ket{1}$. In all Monte Carlo simulation, we use $10^{4}$ samples, $\Delta\phi=\frac{\pi}{2}$ is chosen.