Symmetry-protected dissipative preparation of matrix product states

TL;DR

This paper introduces a symmetry-protected dissipative protocol for efficiently preparing matrix product states with topological order, robust against certain errors, and scalable with system size.

Contribution

It presents a novel dissipative method leveraging symmetry for preparing complex quantum states, with improved scalability and a concrete experimental scheme.

Findings

Preparation time scales polynomially with system size

Parallel preparation reduces scaling to logarithmic squared

Protocol is robust against symmetry-respecting errors

Abstract

We propose and analyze a method for efficient dissipative preparation of matrix product states that exploits their symmetry properties. Specifically, we construct an explicit protocol that makes use of driven-dissipative dynamics to prepare a many-body quantum state that features symmetry-protected topological order and non-trivial edge excitations. The preparation protocol is protected from errors that respect the symmetry, allowing for robust experimental implementation without fine-tuned control. Numerical simulations show that the preparation time scales polynomially in system size $n$. Furthermore, we demonstrate that this scaling can be improved to $\mathcal{O}(\log^2 n)$ by using parallel preparation of individual segments and fusing them via quantum feedback. A concrete scheme using excitation of trapped neutral atoms into Rydberg states via electromagnetically induced transparency is proposed, and generalizations to a broader class of matrix product states are discussed.

Equations

$$H_{\text{AKLT}}=\sum _i[{\vec{S}}_i\cdot {\vec{S}}_{i+1}+\frac{1}{3}({\vec{S}}_i\cdot {\vec{S}}_{i+1}{)}^2],$$

$$|G_{ab}^n⟩=\sum _{\{s_i\}}⟨a⟩A^{(s_1)}{A}^{(s_2)}\cdots A^{(s_n)}b|s_1{s}_2\cdots s_n⟩,$$

$$A^{(1)}=\sqrt{\frac{2}{3}}{\sigma }^{+},A^{(0)}=-\sqrt{\frac{1}{3}}{\sigma }^z,A^{(-1)}=-\sqrt{\frac{2}{3}}{\sigma }^{-}.$$

$$⟨G_{ab}^n⟩G_{a^{\prime }{b}^{\prime }}^n=⟨a{a}^{\prime }⟩{\mathbf{T}}^{n}b{b}^{\prime }=$$

$$\frac{1}{2}{\delta }_{a{a}^{\prime }}{\delta }_{b{b}^{\prime }}(1-(-1{)}^{{\delta }_{ab}}{ϵ}^n)+{\delta }_{ab}{\delta }_{a^{\prime }{b}^{\prime }}(1-{\delta }_{a{a}^{\prime }}){ϵ}^n,$$

$$\dot{\rho }=\mathcal{L}\rho \equiv -i[H,\rho ]+\sum _{\mu }{\Gamma }_{\mu }\mathcal{D}[c_{\mu }]\rho ,$$

$${\mathcal{L}}_{\text{MP}}=(\Gamma /5)\sum _i\sum _{\nu =0}^4\mathcal{D}[{\bar{c}}_{\nu }^{(i)}].$$

$${\mathcal{L}}_{\text{CW}}=\frac{\omega }{2\pi }{\int }_0^{2\pi /\omega }dt\Gamma \sum _i\mathcal{D}[e^{i{H}_{\text{CW}}t}{c}_2^{(i)}{e}^{-i{H}_{\text{CW}}t}].$$

$$|\varphi ⟩=\sum _{a,b=\downarrow }^{\uparrow }\sum _{z=-2}^2{\varphi }_{abz}|G_{ab}^{n-1}⟩|z⟩.$$

$$|{\varphi }_1^{\perp }⟩$$

$$|{\varphi }_2^{\perp }⟩$$

$$0$$

$$=\sum _{abz}{\varphi }_{abz}⟨G_{pq}^{n+1}⟩{\bar{c}}_{\nu }^{(n)}{G}_{ab}^{n-1}|z⟩\forall p,q,\nu .$$

$$p={\left|⟨{\psi }_f^{2m}⟩{\psi }_0\right|}^2=\frac{|\vec{\alpha }\cdot \vec{\beta }{|}^2}{2|\vec{\alpha }{|}^2|\vec{\beta }{|}^2}+\mathcal{O}(3^{-m}).$$

$$\mathrm{Pr}\{\text{completion}\}=(1-(1-p{)}^K{)}^{N_{\mathrm{chain}}-1}.$$

$$T(n)\le T_0+K_{\ast }{\tau }_c+(K_{\ast }-1){\tau }_r,$$

$$T(n)=\mathcal{O}(\log (n)\log (n/\mathcal{E}))=\mathcal{O}({\log }^{2}n).$$

$${\Gamma }_{DD}\approx -2\frac{g^4}{{\Omega }^4}\text{Im}\left[\frac{U}{1+i\chi U}\right],$$

$$|A_{ab}^n⟩=\sum _{\{s_i\}}⟨a⟩A^{(s_1)}{A}^{(s_2)}\cdots A^{(s_n)}b|s_1{s}_2\cdots s_n⟩,$$

$$V_g^{\otimes n}|A_{\circ }^n⟩=e^{i{\theta }_g}|A_{\circ }^n⟩.$$

$$\mathfrak{B}=\mathrm{range}\left(\sum _{\mu =1}^{k_{\min }}\sum _{g\in \mathcal{G}}{V}_g^{†}{c}_{\mu }^{†}{c}_{\mu }{V}_g\right).$$

$$k_{\min }=\underset{\begin{array}{c}\text{irrep }r\text{ of }\mathcal{G}\text{ in }\mathfrak{B}\end{array}}{\max }\lceil \frac{\text{\# of copies of irrep }r\text{ in }\mathfrak{B}}{\text{dimension of irrep }r}\rceil .$$

$$(\mathbf{M}{)}_{pq\nu }^{abz}=\sum _{r=1}^D⟨pa⟩{\mathbf{T}}^{n-2}rb\times ⟨G_{rq}^2⟩c_{\nu }z,$$

$$p$$

$$=\frac{|\vec{\alpha }\cdot \vec{\beta }{|}^2}{D|\vec{\alpha }{|}^2|\vec{\beta }{|}^2}+\mathcal{O}({ϵ}_2^m)\text{if}{C}_{bc}={\alpha }_b{\beta }_c,$$

$${\tilde{C}}_{b^{\prime }{c}^{\prime }}=\sum _{bc}⟨\psi |(|A_{b^{\prime }b}^1⟩\otimes |A_{c{c}^{\prime }}^1⟩)C_{bc}.$$

$$\mathrm{tr}(\tilde{C})$$

$$=\sum _{abc}⟨\psi |\mathcal{P}(|A_{ca}^1⟩\otimes |A_{ab}^1⟩)C_{bc}$$

$$=\sum _{bc}\pm ⟨\psi ⟩A_{cb}^2{C}_{bc}=0.$$

$${\left|\mathrm{tr}(u_g^{†}\tilde{C})\right|}^2$$

$$\displaystyle\mathds{E}_{g\in\mathcal{G}}\left[\big{|}\operatorname{tr}(u_{g}^{\dagger}\tilde{C})\big{|}^{2}\right]$$

$$=\frac{1}{D}\sum _i{\lambda }_i=\frac{1}{D}\mathrm{tr}({\tilde{C}}^{†}\tilde{C}).$$

$$\begin{array}{rl}{A}^{(0)}=|\uparrow ⟩⟨\uparrow |& =\left(\begin{array}{cc}1& 0\\ 0& 0\end{array}\right)\\ \text{and}{A}^{(1)}=|\downarrow ⟩⟨\downarrow |& =\left(\begin{array}{cc}0& 0\\ 0& 1\end{array}\right).\end{array}$$

$${\mathcal{L}}_{\text{MP}}(t)=(V_{\theta }^{†}{)}^k{\mathcal{L}}_0{V}_{\theta }^k\text{for}k\tau \le t<(k+1)\tau .$$

Full text

I Introduction

Entangled many-body states play a central role in understanding strongly correlated quantum matter and constitute the key resource for quantum information science. Matrix product states (MPSs) Perez-Garcia et al. (2007) form an important class of many-body entangled states that can describe a variety of one-dimensional quantum systems. Notably, MPSs include states featuring symmetry-protected topological (SPT) order Schuch et al. (2011); Gu and Wen (2009); Chen et al. (2011); Senthil (2015), corresponding to exotic quantum phases beyond the standard Landau paradigm of symmetry breaking. Such states can be prepared either through a sequence of unitary quantum gate operations, or by first engineering the parent Hamiltonian and subsequently preparing its ground state via adiabatic evolution or cooling Bohnet et al. (2016); Cohen et al. (2015); van Bijnen and Pohl (2015); Mazurenko et al. (2017); Cai et al. (2013); Yan et al. (2013); Barends et al. (2016). However, generating entanglement among many particles using these approaches is challenging, as it typically requires high-fidelity control of individual interactions while maintaining low entropy for intrinsically out-of-equilibrium systems. In particular, unavoidable coupling to environment limits the lifetime of these states and hinders their potential applications.

In this paper, we propose and analyze an alternative method to efficiently prepare an MPS by engineering couplings between a system and its environment such that the desired quantum state is obtained as the unique steady state of time evolution. Such approaches to prepare entangled states have been described previously Kraus et al. (2008); Verstraete et al. (2009); Ticozzi and Viola (2012); Cho et al. (2011); Albert and Jiang (2014); Iemini et al. (2016); Barreiro et al. (2010); Reiter et al. (2013); Rao and Mølmer (2014); Roghani and Weimer (2018). It has also been shown that under certain conditions, the dissipative method outperforms corresponding unitary schemes Kastoryano et al. (2011). In practice, however, the implementation of these suggested schemes in many-body systems is challenging as it requires engineering of complex interactions and decay channels with environment. Here, we show how symmetries can be used to design a simple, translation-invariant dissipative process that only requires a single decay channel and global manipulations to create a desired MPS. Remarkably, similar to symmetry-protected equilibrium phases, this dissipative dynamics is protected from imperfections that respect the symmetry, allowing for robust experimental realizations in large systems with current technologies. The symmetry protection implies that our method does not require high fidelity in controls such as global spin-rotations or interaction strength, which are otherwise essential in conventional unitary schemes.

The paper is organized as follows. In Sec. II, we describe our method for a well-known example of an MPS that exhibits SPT order. We elaborate on how to construct a driven-dissipative dynamics from given symmetry properties of the MPS and rigorously show that the engineered dynamics deterministically prepares the desired states. Using numerical simulations, we find the state-preparation time scales polynomially with system size $n$. In Sec. III, we show that this scaling of state-preparation time can be further improved to $\mathcal{O}(\log^{2}n)$ by first preparing multiple chains in parallel and then connecting them via repeated measurements with feedback. This corresponds to an exponential improvement over previously known scaling $\mathcal{O}(n^{\log n})$ for generic MPSs Verstraete et al. (2009). We also provide a detailed analysis of the effect of imperfections in quantum feedback. We emphasize that our example scheme utilizes the most natural, generic types of environmental couplings as a resource, and hence it can be directly implemented in cold atom systems with existing technologies. In Sec. IV, we propose a concrete scheme involving excitation of trapped neutral atoms Endres et al. (2016); Barredo et al. (2016); Bernien et al. (2017); Labuhn et al. (2016); Lester et al. (2015) into Rydberg states via electromagnetically induced transparency (EIT) Petrosyan et al. (2011), in which spontaneous emissions of photons from atomic excited states are harnessed as resources. In Sec. V, we generalize our protocol to a broader class of MPSs, including the ground states of all one-dimensional SPT phases, and derive a lower bound on the number of required decay channels that may be saturated by an explicit construction. We discuss our results and outlook in Sec. VI.

II Dissipative Preparation of AKLT states

II.1 Ground states of AKLT Hamiltonian

We illustrate our scheme by starting with an example to deterministically prepare a chain of spin-1 particles into the ground states of a gapped, frustration-free Hamiltonian

[TABLE]

where $\vec{S}_{i}$ is the spin-1 vector operator acting on a particle at site $i$. First studied by Affleck, Kennedy, Lieb, and Tasaki (AKLT), the ground states of $H_{\textrm{AKLT}}$ are paradigmatic examples of MPSs and model states for the Haldane phase Affleck et al. (1987); Haldane (1983); Pollmann et al. (2010). While under periodic boundary condition, $H_{\textrm{AKLT}}$ has a unique ground state; under open boundary condition, the ground states are fourfold degenerate due to two fractionalized degrees of freedom on the edges. These constitute a signature of symmetry-protected topological order, which can be experimentally verified by measuring a non-local string order parameter den Nijs and Rommelse (1989); Pérez-García et al. (2008).

The ground states of AKLT Hamiltonian have exact MPS representations Perez-Garcia et al. (2007). For a system of $n$ spin-1 particles, the unnormalized AKLT ground states can be written as

[TABLE]

where $s_{i}\in\{\pm 1,0\}$ runs over three possible spin projections along the $\hat{z}$-axis for a particle at site $i$. The quantum amplitude for each many-body basis state $\ket{s_{1}s_{2}\dots s_{n}}$ is obtained from the products of $D\times D$ matrices $A^{(s)}$, and the boundary conditions for the matrix products are specified by a row (column) vector $\bra{a}$ ($\ket{b}$) of dimension $D$. Specifically, for AKLT ground states, $A^{(s)}$ can be concisely written using ($D$ = 2) Pauli matrices $\sigma^{x}$, $\sigma^{y}$, $\sigma^{z}$, and $\sigma^{\pm}=(\sigma^{x}\pm i\sigma^{y})/2$:

[TABLE]

The choices of vectors $a,b\in\{\uparrow,\downarrow\}$ distinguish the four degenerate ground states with different fractionalized edge states under open boundary conditions 111We note that our convention here for labeling the edge states $a,b$ is different from the convention where the AKLT state is understood as projecting pairs of virtual spin-1/2 particles in singlet bonds back into the spin-1 particles, up to a basis change by $\mathds{1}\otimes i\sigma^{y}$.. Under periodic boundary condition, the unique ground state is $\ket{G^{n}_{\circ}}\equiv\sum_{a}\ket{G_{aa}^{n}}$.

We note a few properties of these AKLT states. With the notation $\ket{G^{n}_{ab}}$, one can conveniently rewrite the quantum state of an $n$-particle system as a linear superposition of composite systems, each with $m$ and $n-m$ particles, i.e. $\ket{G_{ab}^{n}}=\sum_{c}\ket{G_{ac}^{m}}\ket{G_{cb}^{n-m}}$. Moreover, the overlap between two AKLT states with different edge states can be analytically evaluated; using the transfer matrix $\boldsymbol{T}=\sum_{s}A^{(s)*}\otimes A^{(s)}$, we obtain exponentially small overlaps between distinct states with a normalization factor $1/2$:

[TABLE]

where $\epsilon=-1/3$. Finally, the AKLT ground state exhibits topological order that is protected by the symmetry groups $D_{2}$ (dihedral group corresponding to permutation of spin axes), $\mathcal{T}$ (time-reversal symmetry), and $\mathcal{P}$ (bond-inversion symmetry). While the presence of any of these symmetries can protect the non-trivial order reflected in the double degeneracy in the entanglement spectrum, the $D_{2}$ symmetry is necessary and sufficient to protect the string order parameter Pollmann et al. (2010). In addition to $D_{2}$, the parent Hamiltonian $H_{\textrm{AKLT}}$ respects a larger symmetry group of $\mathrm{SO}(3)$, corresponding to global rotation of spins.

II.2 Constructing driven-dissipative dynamics

Our key idea is to use $\mathrm{SO}(3)$ symmetry of the parent Hamiltonian $H_{\textrm{AKLT}}$ for preparation of an exact AKLT ground state. By converting energy penalties imposed by $H_{\textrm{AKLT}}$ into dissipative penalties in the form of decay channels, we can engineer a process that effectively cools to the ground states. More specifically, we start with a dissipative dynamics that eliminates one type of excitation in $H_{\textrm{AKLT}}$. Then, all other types of excitations can be eliminated using global spin-rotations in $\mathrm{SO}(3)$. Since spin-rotations are symmetries of $H_{\textrm{AKLT}}$, their implementations are robust against imperfections in control parameters such as durations, phases, or strengths of electromagnetic driving.

We consider a Markovian driven-dissipative dynamics described by a quantum master equation:

[TABLE]

where $\rho$ is the density operator of a system, $H$ is a Hamiltonian governing coherent dynamics, and $\mathcal{D}[c_{\mu}]\rho\equiv c_{\mu}\rho c_{\mu}^{\dagger}-\{c_{\mu}^{\dagger}c_{\mu},\rho\}/2$ characterizes incoherent dynamics by jump operators (i.e., decay channel) $c_{\mu}$ at rate $\Gamma_{\mu}$. We can interpret the dynamics of $\mathcal{L}$ as the system evolving with a non-Hermitian Hamiltonian $H_{\text{eff}}$ = $H-i\sum_{\mu}\Gamma_{\mu}c_{\mu}^{\dagger}c_{\mu}/2$, while stochastically undergoing quantum jumps $\rho$ $\mapsto$ $c_{\mu}\rho c_{\mu}^{\dagger}$ at rates $\operatorname{tr}(\Gamma_{\mu}c_{\mu}^{\dagger}c_{\mu}\rho)$ for each jump operator Gardiner and Zoller (2004).

In order to construct a simplest possible $\mathcal{L}$ that prepares an AKLT ground state, we exploit the $\mathrm{SO}(3)$ symmetry that conserves total angular momentum. In particular, each term in $H_{\textrm{AKLT}}$ can be written as $2P_{i}-2/3$, where $P_{i}$ is the projection operator onto the subspace of total angular momentum $J_{i}$ = $S_{i}+S_{i+1}$ $=$ 2 for the pair of particles $(i,i+1)$. Hence, a state $\ket{G}$ minimizes the energy if it has no population in the $J_{i}$ = 2 manifold, i.e. $P_{i}\ket{G}$ = [math] for every nearest-neighboring pair. Under open boundary condition, there are four such states $\ket{G_{ab}}$, labelled by two spin-1/2 edge degrees of freedom $a,b$ $\in$ $\{\uparrow,\downarrow\}$. Under periodic boundary condition, only a unique state $\ket{G_{\circ}}$ $\propto$ $\ket{G_{\uparrow\uparrow}}+\ket{G_{\downarrow\downarrow}}$ satisfies the constraints. Below, we use $\ket{G}$ to denote the ground state(s) when boundary conditions are not specified.

To prepare $\ket{G}$, we use jump operators to depopulate $J_{i}$ = 2 manifold of every neighboring pair. For example, we can set $H$ = 0 and use five types of jump operators, $c_{m}^{(i)}$ = $|{\phi_{m}}\rangle\langle{J=2,J^{z}=m}|_{i,i+1}$, where $\{|J$ = $2,J^{z}$ = $m\rangle_{i,i+1}$ :  $m$ = $-2,\dots,2\}$ is an orthonormal basis spanning the $J_{i}$ $=$ $2$ manifold for the pair of spins $(i,i+1)$, and $\ket{\phi_{m}}$ is any other quantum state with nonzero population in $J_{i}$ = $0,1$ manifolds ($\bra{\phi_{m}}P_{i}\ket{\phi_{m}}$ $<$ 1). With these jump operators, quantum jumps occur at rate $\Gamma_{\textrm{total}}$ = $\Gamma\sum_{i,m}\operatorname{tr}{(\rho c_{m}^{(i){\dagger}}c_{m}^{(i)})}$ = $\Gamma\sum_{i}\operatorname{tr}{(\rho P_{i})}$, which vanishes only for the ground state $\ket{G}$. This implies that $\ket{G}$ is a steady state of $\mathcal{L}$, and any other quantum state will undergo a series of quantum jumps.

Using $\mathrm{SO}(3)$ symmetry, this construction can be effectively realized with only one type of jump operator via global coherent manipulations $H$. More specifically, let us consider a dynamics with only one jump operator, $c_{2}$ = $\ket{00}\langle J$=$2,J^{z}$=$2|$ = $|{00}\rangle\langle{++}|$, written in the $S^{z}$ basis $\{\ket{+},\ket{0},\ket{-}\}$. After time evolution over duration $\tau/5$, we apply a fast global pulse $V$ $=$ $\exp[i(2\pi/5)\sum_{i}S_{i}^{y}]$, rotating the entire spin ensemble by an angle $2\pi/5$ about the $y$-axis. In a rotating frame, this operation implements the jump operator $V^{\dagger}c_{2}V$. Repeated multiple times, we obtain five distinct jump operators $\bar{c}_{\nu}$ $\equiv$ $(V^{\dagger})^{\nu}c_{2}(V)^{\nu}$ for $\nu$ $\in$ $\{0,\dots,4\}$ after the $\nu$-th (modulo 5) pulses. For a sufficiently short $\tau$ $\ll$ $1/\Gamma$, the effective Liouvillian of the five-pulse cycle can be well approximated using leading-order Magnus expansion by

[TABLE]

Note that the purpose of global rotations is to use a single jump operator for depopulating different states; different choices of angles and axes are equally effective as long as states rotated from $\ket{++}$ span the entire $J$ = 2 manifold. We may also employ a time-independent Hamiltonian $H_{\textrm{CW}}$ = $\omega\sum_{i}S_{i}^{y}$ to continuously rotate the ensemble, leading to an effective Liouvillian

[TABLE]

In both cases, the corresponding quantum jump rates vanish if and only if the system is in $\ket{G}$.

II.3 Proving uniqueness of steady states

While our construction of $\mathcal{L}_{\textrm{MP}}$ and $\mathcal{L}_{\textrm{CW}}$ ensures that $\ket{G}$ is a steady state, one can imagine an undesired mixed steady state that forms in dynamical equilibrium from the combination of coherent evolution and incoherent quantum jumps (Fig. 1). Such mixed steady states may arise only if there exists a subspace $\mathcal{S}$ orthogonal to $\ket{G}$ and closed under jump operators, $c_{\mu}\mathcal{S}$ $\subseteq$ $\mathcal{S}$ Kraus et al. (2008). Physically, this means that states in $\mathcal{S}$ cannot reach $\ket{G}$ even with arbitrarily many applications of jump operators $c_{\mu}$, allowing an equilibrium to form by their mixtures. In our scheme, we prove the following lemma that guarantees that the desired state $\ket{G}$ is the unique steady state:

Lemma**.**

For any finite system with size $n$ $\geq$ $2$ under open boundary condition, all states can reach $\ket{G}$ with some application of jump operators in $\mathcal{L}_{\textnormal{MP}}$ or $\mathcal{L}_{\textnormal{CW}}$, implying $\ket{G}$ is the unique steady state.

We now sketch the proof of this Lemma for open boundary condition, where the four $\ket{G_{ab}}$ states are unique steady states; more details on the proof can be found in Appendix C. The proof uses induction on system size $n$. Note that under an open boundary condition, the four $\ket{G_{ab}}$ states are the desired steady states. Let us define $D_{n}=\text{span}\{\ket{G_{ab}^{n}}\}$ as the subspace of AKLT ground states of $n$ spins. For $n$ = $2$ and 3, the uniqueness of steady states can be checked by exact diagonalization. Our induction hypothesis is that any $n$-spin input state can reach $D_{n}$ with some applications of jump operators, via some polynomial function $f_{[n]}(\{\bar{c}_{\nu}\})$ (which may depend on the input state). For the sake of contradiction, let us assume that a state $\ket{\psi^{n+1}}$ cannot reach $D_{n+1}$ with any sequence of jump operators $\bar{c}_{\nu}^{(i)}$ in $\mathcal{L}_{\textrm{MP}}$. (The same argument holds for $\mathcal{L}_{\textrm{CW}}$.) We will then construct a sequence of jump operators, involving some $f_{[n]}$ on the first $n$ spins followed by some $\bar{c}_{\nu}^{(n)}$ acting on spin $n$ and $n+1$, so that $\bar{c}_{\nu}^{(n)}f_{[n]}\ket{\psi^{n+1}}$ reaches $D_{n+1}$, leading to a contradiction.

To begin, we know that by our induction hypothesis, there exists $f_{[n]}(\{\bar{c}_{\nu}\})$ so that $f_{[n]}\ket{\psi^{n+1}}$ has nonzero population in $D_{n}$. Since the AKLT Hamiltonian is frustration-free, $f_{[n]}\ket{\psi^{n+1}}$ must also have nonzero population in $D_{n-1}$, i.e. the AKLT ground states on the first $n-1$ spins. We then do a general decomposition of $f_{[n]}\ket{\psi^{n+1}}=\ket{\phi}+\ket{\phi^{\perp}_{1}}+\ket{\phi^{\perp}_{2}}$, where

[TABLE]

Here $\ket{z}$ runs through the 5 states in the $J=2$ manifold on the last two spins. The remaining parts of the wavefunction $f_{[n]}\ket{\psi^{n+1}}$ are

[TABLE]

where $\ket{E_{\mu}^{n-1}}$ runs through all the excited eigenstates of the $H_{\textrm{AKLT}}$ on the first $n-1$ spins, and $s$ runs through all 9 possible 2-spin states. We now consider two cases:

Case (i): Suppose $\phi_{abz}=0$. Then $f_{[n]}\ket{\psi^{n+1}}=\ket{\phi^{\perp}_{1}}+\ket{\phi^{\perp}_{2}}$. From our inductive hypothesis, $f_{[n]}\ket{\psi^{n+1}}$ must have nonzero population in $D_{n}$ and $D_{n-1}$; this can only come from the $\ket{\phi^{\perp}_{1}}$ part, and thus there must be nonzero coefficients $\phi_{abst}^{\perp}\neq 0$. On the other hand, since we have assumed that $\ket{\psi^{n+1}}$ cannot reach $D_{n+1}$, we must have $0=\braket{G_{pq}^{n+1}}{\phi^{\perp}_{1}}$, which is only possible if $\ket{\phi^{\perp}_{1}}$ has the state of spins at $n-1$ and $n$-th sites in the $J$ = $2$ manifold. However, that would then imply that $\ket{\phi^{\perp}_{1}}$ has zero population in $D_{n}$, contradicting our inductive hypothesis.

Case (ii): Now suppose $\phi_{abz}\neq 0$ for some $a,b,z$. Note that for any jump operator $\bar{c}_{\nu}^{(n)}$ acting on spins $n$ and $n+1$, we have $\braket{G_{pq}^{n+1}}{\bar{c}_{\nu}^{(n)}}{\phi^{\perp}_{j}}=0$ for $j=1,2$, since $\bar{c}_{\nu}\ket{G_{st}^{2}}=0$ and $\braket{G_{ab}^{n-1}}{E_{\mu}^{n-1}}=0$. Then our assumption that $\ket{\psi^{n+1}}$ cannot reach $D_{n+1}$ gives the following set of linear equations for $\phi_{abz}$:

[TABLE]

Since $\ket{G^{n}_{ab}}$ and $\ket{G^{n+1}_{pq}}$ have explicit MPS formulas, one can analytically compute these expressions and find that only the trivial solution $\phi_{abz}$ = [math] are allowed for $n\geq 3$. This yields a contradiction and implies that all states $\ket{\psi^{n+1}}$ can reach at least one of the states $\ket{G^{n+1}_{pq}}$ in $D_{n+1}$ with some application of jump operators.

II.4 Numerical simulations and scaling

We numerically study the efficiency of our protocol via stochastic wavefunction method for systems of up to $n=25$ particles. We use both exact diagonalization (for $n$ $\leq$ $8$) and time-evolving block decimation (TEBD) algorithm Vidal (2003); *Vidal2 in MPS representations (for $n$ $\leq$ $25$); more details are discussed in Appendix D. We initialize the system in a random product state (representing a maximally mixed state), and evolve under $\mathcal{L}_{\mathrm{MP}}$ with open boundary condition. We then monitor the energy density with respect to $H_{\text{AKLT}}$, as well as the fidelity of state preparation $\mathcal{F}$ $=$ $\braket{P_{G}}$, where $P_{G}$ is the projector onto the ground states. The results in Figs. 2(a) and (b) demonstrate that both observables exponentially converge to their corresponding values for AKLT states in all system sizes. We extract the state-preparation time $T$ by first fitting $1-\mathcal{F}$ to an exponential in the long-time regime and extrapolating $\mathcal{F}(T)$ $=$ $0.9$. We find that $T$ generally increases with system size $n$. Plotted as a function of $n$ [Fig. 2(c)], we find a polynomial scaling $T$ $\sim$ $\mathcal{O}(n^{2.97})$. This scaling is consistent with the more complicated protocol in Ref. Verstraete et al. (2009) that requires up to $O(n^{\log_{2}n})$ time, up to system sizes simulated in this work.

III Improving Scaling via Parallelization and Quantum Feedback

While the time our protocol needs to prepare AKLT states is already shown to have an efficient polynomial scaling from our numerical simulations, we now provide a method to exponentially improve this scaling to $\mathcal{O}(\log^{2}n)$. Similar to approaches used in quantum repeaters Briegel et al. (1998), this exponential speedup is possible by preparing multiple chains in parallel, which are subsequently connected into a single long chain [Fig. 3(a)]. The key ingredient is the ability to efficiently connect or fuse two AKLT chains into a single entangled state [Fig. 3(b)], which we now describe.

III.1 Connecting two AKLT chains

Suppose we have independently prepared two chains of $m$ spins in AKLT states. As they are initially unentangled, their state can be written as $\ket{\psi_{0}}=\mathcal{N}_{0}(\sum_{b}\alpha_{b}\ket{G^{m}_{ab}})\otimes(\sum_{c}\beta_{c}\ket{G^{m}_{cd}})$, where $\vec{\alpha},\vec{\beta}\in\mathds{C}^{2}$ characterizes the edge states at the interface, and $\mathcal{N}_{0}$ is some normalization constant. The edge states represented by the indices $a$ and $b$ are unimportant for the connection. To connect the two chains into one, we can turn on the jump operators $\{\bar{c}_{\alpha}^{(m)}\}$ acting on spins $m$ and $m+1$ at the interface, and evolve for some time $\tau_{c}$. Then, we monitor quantum jump events to determine if we have succeeded in creating an AKLT state $\ket{\psi_{f}^{2m}}$ $=$ $\mathcal{N}_{f}\ket{G_{ad}^{2m}}$ with doubled length ($\mathcal{N}_{f}$ is some normalization constant).

A successful connection is heralded by the absence of quantum jumps, in which case evolution under the non-Hermitian Hamiltonian guides the system into an AKLT state of the combined chain. For sufficiently long $\tau_{c}$, the success probability is given by the overlap between initial and desired states (see Appx. B). This can be computed using Eq. (II.1) to be

[TABLE]

Note that if the edge states $\vec{\alpha}$ and $\vec{\beta}$ are random vectors in $\mathds{C}^{2}$, we have on average $p\simeq 1/4$. When the edge states are aligned, i.e. $\vec{\alpha}\propto\vec{\beta}^{*}$, we obtain the maximum success probability of $p_{\text{max}}=1/2$. The resultant state has an exponentially small error $\epsilon$ $\leq$ $e^{-\mathcal{O}(\tau_{c})}$.

The failure of the connection is signaled by detection of a quantum jump $\bar{c}_{\alpha}$, after which the state of the system changes according to $\ket{\psi_{0}}\mapsto\bar{c}_{\alpha}^{(m)}\ket{\psi_{0}}$. In this case, one can discard the pair of spins $(m,m+1)$ and then attempt the connection procedure again with two chains of length $m-1$. However, it turns out that quantum jumps affect the success probability of subsequent connection attempts, which in fact vanishes for this protocol without additional intervention, a phenomenon that we will explain in the following paragraph. Nevertheless, we can restore the success probability to $p_{\max}=1/2$ by applying a global spin rotations $U$ = $(e^{i\pi S^{y}})^{\otimes{m-1}}$ to one of the chains. This quantum feedback makes the procedure very efficient, since multiple repeated failures are exponentially unlikely and only one success is sufficient to fuse two chains. The number of attempts necessary follows the geometric distribution, and on average we need $(1-p)/p$ attempts with the loss of $2(1-p)/p$ particles per connection. By performing these connection procedures in parallel, we can quickly prepare an AKLT state of $n$ spins in $\mathcal{O}(\log^{2}n)$ time, as we show in Sec. III.2.

We now explain why a single failed connection attempt will cause subsequent attempts to fail, unless appropriate quantum feedback is applied. It turns out that this problem occurs whenever the matrix product state we want to prepare respects bond-inversion symmetry, but here we first focus on the example of AKLT states which has an intuitive explanation (see Sec. V.2 for the general case). Note that we can interpret the dynamics under $\mathcal{L}=\mathcal{D}[\bar{c}_{\theta}]$ as a continuous measurement of whether the pair of spins has total angular momentum $J_{\theta}=+2$, where $J_{\theta}=e^{-i\theta J_{y}}J_{z}e^{i\theta J_{y}}$. To be more specific, consider four spin-1 particles $\vec{S}_{1},\vec{S}_{2},\vec{S}_{3},\vec{S}_{4}$, and imagine that we are performing a connection between spins 2 and 3 by continuously measuring $\vec{J}=\vec{S}_{2}+\vec{S}_{3}$. Suppose we decompose each spin-1 into two virtual spin-$\frac{1}{2}$ particles: $\vec{S}_{i}=\vec{s}_{i,L}+\vec{s}_{i,R}$. It is known Affleck et al. (1987) that an AKLT state can be constructed by starting with singlet bonds of virtual spin-$\frac{1}{2}$ particles where $s_{i,R}+s_{i+1,L}=0$ for all $i$, and then projecting back into the triplet subspace of the original pairs of virtual spin-$\frac{1}{2}$ particles where $s_{i,L}+s_{i,R}=1$. The detection of a quantum jump $c_{\theta}$ in a failed connection attempt implies that $J_{\theta}=2$, which is only possible if $s_{2L}^{\theta}=s_{2R}^{\theta}=s_{3L}^{\theta}=s_{3R}^{\theta}=+\frac{1}{2}$. Due to the singlet bond conditions, this automatically implies that $s_{1R}^{\theta}=s_{4L}^{\theta}=-\frac{1}{2}$. Subsequently, when we retry the connection with spins 1 and 4 (after discarding 2 and 3), the two virtual spin-$\frac{1}{2}$ particles at the interface are in the state $\ket{s_{1R}^{\theta}=-\frac{1}{2}}\ket{s_{4L}^{\theta}=-\frac{1}{2}}$, which has no overlap with the desired singlet bond state $\ket{+\frac{1}{2}}\ket{-\frac{1}{2}}-\ket{-\frac{1}{2}}\ket{+\frac{1}{2}}$. Hence the overlap with the AKLT state is zero, and the connection is certain to fail. Now, applying $U$ = $(e^{i\pi S^{y}})^{\otimes{m-1}}$ to the first chain flips $s_{1R}^{\theta}$ so that the resultant state is $\ket{+\frac{1}{2}}\ket{-\frac{1}{2}}$. This has an overlap of 1/2 with the singlet state, restoring our success probability to roughly 1/2.

III.2 Scaling of preparation time of the parallelized protocol

We are now ready to describe and analyze the full parallelized protocol to prepare AKLT states on large system sizes $n$ in $\mathcal{O}(\log^{2}n)$ time. To prepare a length-$n$ chain, we first prepare $\mathcal{O}(n/n_{0})$ AKLT chains of length $n_{0}$, and then apply the above connection procedure in parallel for adjacent chains, as illustrated in Fig. 3(a). Since on average we lose $n_{c}=2(1-p)/p$ particles per connection, we should choose $n_{0}\gg n_{c}$, and prepare $N_{\rm chains}=(n-n_{c})/(n_{0}-n_{c})$ such chains. We can imagine attempting the connection simultaneously for all junctions between chains, and we will typically have some successes and some failures. By noting the locations of the failures and applying the quantum feedback $e^{i\pi S^{y}}$ to every other segment, we can then re-attempt connections on the junctions that have failed, which will then succeed with $p\simeq 1/2$ as shown above. This process is repeated for multiple rounds until all junctions become connected, and the system becomes one connected chain. The successful connection of each junction occurs independently and probabilistically; hence their order may be arbitrary. The probability of completing all $N_{\rm chains}-1$ connections after $K$ rounds is

[TABLE]

To achieve a completion probability of $p_{\rm comp}$, we need $K_{*}=\log[1-p_{\rm comp}^{1/(N_{\rm chain}-1)}]/\log(1-p)$ = $\mathcal{O}(\log n)$ rounds. Hence, the time required to successfully complete all connection with constant probability and obtain a length-$n$ AKLT chain is

[TABLE]

where $T_{0}$ is the preparation time of the length-$n_{0}$ chains, $\tau_{c}$ is the time for each connection attempt, and $\tau_{r}$ is the time required for feedback after each failed attempt. Recall that each successful connection induces an error of $\epsilon$ $\leq$ $e^{-\mathcal{O}(\tau_{c})}$ in the quantum state. Thus, a total of $\mathcal{O}(n/n_{0})$ connections yield a final error of $\mathcal{E}\leq\mathcal{O}(n/n_{0})e^{-\mathcal{O}(\tau_{c})}$, which means we should choose $\tau_{c}$ $=$ $\mathcal{O}(\ln(n/n_{0}\mathcal{E}))$ to achieve a final error of $\mathcal{E}$. Assuming arbitrarily fast classical communication and control, the quantum feedback of applying homogeneous spin-rotation $e^{i\pi S_{y}}$ to a subset of the chains can be done in a system-size independent time $\tau_{r}=\mathcal{O}(1)$. Thus, the average time necessary to prepare an AKLT state of length $n$ with bounded error $\mathcal{E}$ in this parallelized protocol is

[TABLE]

The required number of spins that we need initially is $n_{0}(n-n_{c})/(n_{0}-n_{c})$, indicating that an $\mathcal{O}(n)$ (i.e., constant factor) spatial overhead is sufficient.

While so far we have assumed that we can detect the occurrence of quantum jumps perfectly, this parallelized protocol remains very efficient even when the detection is imperfect. We describe two possible methods to address detection inefficiency. The first method is to use jump operators of the form $c=|{++}\rangle\langle{++}|$ for connection, which will cause an indefinite number of quantum jumps to occur once the first quantum jump occurs, enhancing the quantum jump signal. Alternatively, we can also slowly turn on additional jump operators near the interface if we do not detect any quantum jump initially; this serves to confirm that the two original chains have been successfully connected, as any failed connection attempt that evaded detection would cause more quantum jumps. In both cases, the scaling of our parallelized protocol is not significantly altered compared to the ideal case, and can be largely accounted for by modifying the effective success probability $p$. We provide a more detailed analysis of these two methods in Appendix E.

IV Experimental realization

The key task in implementing our protocol is to engineer the nearest-neighbor jump operators. Such engineering has been previously demonstrated in systems of trapped ions Barreiro et al. (2010). Here, we provide an explicit method to realize our scheme in systems of trapped atoms Endres et al. (2016); Barredo et al. (2016); Bernien et al. (2017); Labuhn et al. (2016); Lester et al. (2015) based on the Rydberg-EIT mechanism Petrosyan et al. (2011). We consider a five-level system consisting of a metastable Rydberg state $\ket{r}$, a short-lived excited state $\ket{e}$, and three long-lived hyperfine ground states $\ket{+},\ket{0}$, and $\ket{-}$ as shown in Fig. LABEL:fig:(a). Using lasers, we coherently couple the ground state $\ket{+}$ to the excited state with a time-dependent Rabi frequency $g(t)$. The excited state is further coupled to the Rydberg state with Rabi frequency $\Omega$. Owing to large dipole moments, simultaneous excitations of two Rydberg states within distance $R$ are suppressed by an interaction energy shift that decays as $1/R^{6}$ Lukin et al. (2001).

In the absence of interactions, our coherent driving ensures that every atom supports three stable states $\ket{-},\ket{0}$, and $\ket{D(t)}$ $\propto$ $\Omega\ket{+}-g(t)\ket{r}$ for arbitrary choices of $g$ and $\Omega$. We use these three states to encode the spin-1 degree of freedom. When $g(t)$ slowly increases starting from zero, $\ket{+}$ = $|D(t$=$0)\rangle$ adiabatically follows $\ket{D(t)}$ without populating any excited states. In the presence of strong interactions, however, population in the Rydberg state of one atom prevents another Rydberg excitation in its vicinity. Thus, as one gradually turns on $g(t)$, any neighboring atoms initially in $\ket{++}$ necessarily populate the excited states, followed by their decay into one of the three ground states. When [math] $<$ $g$ $\ll$ $\Omega$, this dissipative dynamics produces effective jump operators of the form $c_{\phi}$ = $\ket{\phi}\bra{DD}$ with a total rate

[TABLE]

where $\chi$ $\approx$ $1/\gamma+\gamma/(4\Omega^{2})$, and $\phi$ is one of 9 different combinations of two-particle ground states (see Ref. Reiter and Sørensen (2012) and Appendix F). To engineer the full Liouvillian, we can apply microwave pulses to the three ground states and globally rotate the spin-1 particles by $\theta$ $\approx$ $2\pi/5$ [Fig. 4(b)]. When dissipative interactions and global rotations are alternated, this protocol effectively realizes a dynamics similar to $\mathcal{L}_{\text{MP}}$ and deterministically prepares AKLT states. The experimental platform with a rearrangeable atom array in Refs. Endres et al. (2016); Barredo et al. (2016); Bernien et al. (2017); Labuhn et al. (2016); Lester et al. (2015) is particularly well-suited to parallelize the implementation and exponentially shorten preparation times for large systems.

In practice, unwanted dissipations or interactions can affect the fidelity of our protocol by perturbing the steady state of dissipative dynamics. There are two main imperfections in our proposed implementation: i) atomic states have finite dephasing time $T_{2}$, and ii) long-range Rydberg interaction can lead to dissipative coupling with particles beyond nearest neighbors. For the latter, we find that pairs of particles separated by distance $R$ with interaction $U$ $\sim$ $1/R^{6}$ acquire decay rates $\Gamma_{DD}(R)$ $\sim$ $1/R^{12}$. We study the effects of these imperfections by numerical simulations of long-range effective Hamiltonians and stochastic quantum jumps that now include dephasing operators $|{s}\rangle\langle{s}|$ for $s$ $=$ $+,0,-$. We introduce an effective temperature $T_{\textrm{eff}}$ defined by $\operatorname{tr}[P_{G}\rho(T_{\textrm{eff}})]$ $=$ $\mathcal{F}_{\textrm{SS}}$, where $\mathcal{F}_{\textrm{SS}}$ is the steady-state fidelity and $\rho(T_{\textrm{eff}})$ = $\exp(-H_{\textrm{AKLT}}/T_{\textrm{eff}})/Z$ is the Gibbs ensemble with $Z$ = $\textrm{tr}[\exp(-H_{\textrm{AKLT}}/T_{\textrm{eff}})]$ Deutsch (1991); Eisert et al. (2015). When the steady state is near the gapped ground state, $T_{\textrm{eff}}$ characterizes the quality of prepared state in the thermodynamic limit. The results in Fig. 4(c) show that the temperature decreases with increasing dephasing time $T_{2}$ and eventually saturates due to long-range interactions. While $T_{\textrm{eff}}$ also depends on system size $n$, we find that it stays below the gap of $H_{\textrm{AKLT}}$ for all $n$ studied in the present work (up to $n=8$). We note that the effect of long-range interaction is mitigated in our parallelized protocol, where jump operators are turned on only for a few spin pairs well separated by the length of connected chains. Also, throughout the parallelized protocol, the effective temperature does not increase, since our connection procedure ensures that $1-\mathcal{F}$ scales linearly in system sizes while the density of excited states grows at least as fast (see Appendix G).

V Generalization to MPS with symmetry

Our symmetry-based approach can be generalized to efficiently prepare a broader class of matrix product states. In general, any translation-invariant MPS of $n$ spins can be written as

[TABLE]

where $s_{i}\in\{1,2,\ldots,d\}$ runs over the physical spin basis for the $i$-th particle, and $\ket{a},\ket{b}\in\mathds{C}^{D}$ indicate the “boundary conditions” in the virtual bond space. For the case of AKLT states, we have $d=3$ and $D=2$. Under periodic boundary condition, the unique MPS is given by $\ket{A_{\circ}^{n}}=\sum_{a}\ket{A_{aa}^{n}}$. We say that such an MPS respects an internal symmetry group $\mathcal{G}$ if for every $g\in\mathcal{G}$ and some unitary representation $V_{g}:\mathcal{G}\to U(d)$ we have

[TABLE]

Our protocol can be generalized to prepare any such translation-invariant MPS with internal symmetry, which is a class of states that includes ground states of all one-dimensional SPT phases Sanz et al. (2009); Schuch et al. (2011).

V.1 Generalized protocol with a minimal set of decay channels

We now show how to generalize our protocol to any translation-invariant MPS with internal symmetry $\mathcal{G}$. Specifically, we design a dissipative dynamics that deterministically prepares the ground state(s) of the MPS parent Hamiltonian. This protocol uses a set of global coherent manipulations corresponding to symmetry operations on the MPS, as well as a minimal number $k_{\textrm{min}}$ of decay channels $\{c_{1},c_{2},\dots,c_{k_{\textrm{min}}}\}$ acting on pairs of neighboring particles. We are able to derive a lower bound for $k_{\textrm{min}}$ based on irreducible representations of the symmetry group $\mathcal{G}$, and provide an explicit construction of a minimal set of jump operators saturating the bound. Given such a set $\{c_{\mu}\}$, the uniqueness of steady states can be efficiently verified via the same inductive proof technique in Sec. II.3. For the purpose of preparing ground states of SPT phases, we note that the symmetry $\mathcal{G}$ of a parent Hamiltonian of the MPS may be larger than the minimal symmetry $\mathcal{G}_{p}\subset\mathcal{G}$ that protects the topological order. For example, $\mathcal{G}=\mathrm{SO}(3)$ while $\mathcal{G}_{p}=\mathds{Z}_{2}\times\mathds{Z}_{2}$ in the case of AKLT states  Pollmann et al. (2010).

We now describe our results on the minimum complexity on the decay channels (i.e. jump operators) necessary to prepare a general translation-invariant MPS with symmetry. Without loss of generality, we may assume that the desired states are ground states of a gapped, frustration-free parent Hamiltonian $H_{p}=\sum_{i}h^{(i)}$, where $h^{(i)}$ is a translation-invariant, nearest-neighbor projector that respects the internal symmetry $\mathcal{G}$ Perez-Garcia et al. (2007); Sanz et al. (2009). Each term $h$ can be written in a block diagonal form, corresponding to different irreducible representations of $\mathcal{G}$. We refer to the two-particle subspace that $h$ projects onto as “bright manifold” $\mathfrak{B}\equiv\operatorname{range}(h)\subset\mathds{C}^{d^{2}}$, where $d$ is the internal dimension of each particle. The ground states are uniquely characterized by vanishing populations in $\mathfrak{B}$ for every neighboring pair of particles. In the $H_{\textrm{AKLT}}$ example, $\mathfrak{B}$ corresponds to the $J$ = $2$ manifold of two neighboring spins. Similar to our protocol for AKLT states, we can depopulate $\mathfrak{B}$ by employing jump operators $c_{\mu}$ where $\operatorname{range}(c_{\mu}^{\dagger}c_{\mu})\subseteq\mathfrak{B}$. The number of jump operators can be reduced by utilizing and averaging over all symmetry rotations through $c_{\mu}\mapsto V_{g}^{\dagger}c_{\mu}V_{g}$, where $V_{g}$ is the global unitary rotation by a group element $g\in\mathcal{G}$. In order to fully depopulate the bright manifold $\mathfrak{B}$ and nothing else, the set of jump operators $\{c_{\mu}\}_{\mu=1}^{k_{\min}}$ must satisfy the necessary condition

[TABLE]

As we show below, the minimum number $k_{\min}$ of distinct jump-operators will depends on the structure of the group representation of $\mathcal{G}$.

For simplicity of discussion, let us restrict the decay channels to the rank-1 form of $c_{\mu}=|{\phi_{\mu}}\rangle\langle{\psi_{\mu}}|$. The minimum number of jump operators required can be calculated from the number of different irreducible representations (irrep) of the symmetry group $\mathcal{G}$ within the bright manifold $\mathfrak{B}$. In the case of AKLT states, $\mathfrak{B}$ consists of a single 5-dimensional irrep of the group $SO(3)$. In more general cases, the representation of $\mathcal{G}$ on $\mathfrak{B}$ may contain multiple copies of isomorphic (i.e., equivalent up to a basis change) irreps. The capability of global symmetry operations allows one decay channel to depopulate subspaces corresponding to one copy of each irrep in parallel. Hence, as we have shown earlier, one decay channel is sufficient for preparing the AKLT state. For the more general cases, however, it may be necessary to employ multiple decay channels when more than one copy of an irrep is present. Using Schur’s lemma Artin (2010), we prove that the minimum number of rank-1 decay channels is

[TABLE]

We can also construct $k_{\min}$ decay channels that satisfy the necessary condition of Eq. (19) by choosing a set of $\{\bra{\psi_{\mu}}\}$ that is supported in all irreps, with destructive interference between isomorphic irreps. The details of the proof and the construction are described in Appendix H.3.

Once we have such a minimal set of decay channels $\{c_{\mu}\}$, it remains to ascertain the uniqueness of steady states. This can be efficiently verified using our inductive proof techniques, which show that the steady states are unique as long as there are only trivial solutions to a linear equation like (11). More specifically, to prove uniqueness under open boundary condition, one simply needs to compute a $D^{2}\ell\times D^{2}\operatorname{rank}(\mathfrak{B})$ matrix $\boldsymbol{M}$ whose matrix elements are

[TABLE]

where $\boldsymbol{T}$ is the transfer matrix for the MPS, $a,b,p,q,r,c$ are indices for $D$-dimensional virtual bond space, and $\ket{z}$ enumerates the possible states in the bright manifold $\mathfrak{B}$. Here, $\ell\leq d^{2}$ is related to the maximum dimension of the irreps of the symmetry group $\mathcal{G}$ on tthe bright manifold $\mathfrak{B}$. By verifying (through exact calculations or numerics) that the steady state is unique for some small system size $n_{0}$, one can then prove uniqueness for all $n\geq n_{0}$ by showing that $\det(\boldsymbol{M}^{\dagger}\boldsymbol{M})\neq 0$. We emphasize that these proofs can be done efficiently for any given $\{A^{(s)}\}$ and decay channels $\{c_{\mu}\}$, since the calculations only involve matrices of constant dimensions, independent of system size $n$.

Our Rydberg-EIT implementation proposal can be naturally adapted for these general cases. The Rydberg-EIT scheme allows us to engineer two-body jump operators of the form $c_{\textrm{eff}}=|{s^{L1}s^{L2}}\rangle\langle{s^{R}s^{R}}|$, where $\ket{s^{Li}},\ket{s^{R}}\in\mathds{C}^{d}$ are single-spin states. Unlike the case of AKLT states, the preparation of a generic symmetric MPS may require more than one ($k_{\textrm{min}}\geq 2$) rank-1 jump operator. The implementation of multiple decay channels can be achieved, for example, by introducing extra lasers that couple (additional) hyperfine ground states to the short-lived excited state(s). By adjusting the relative strength of laser driving to each hyperfine ground state, one can engineer different EIT-dark states that acquire a dissipative interaction. This allows us to generate a set of jump operators $\{c^{\textrm{eff}}_{\mu}=|{s^{L1}_{\mu}s^{L2}_{\mu}}\rangle\langle{s^{R}_{\mu}s^{R}_{\mu}}|:\mu=1,...,k_{\text{min}}\}$ with independent $\ket{s^{R}_{\mu}s^{R}_{\mu}}$. When particles are individually addressable, one can engineer a jump operator with a more complicated right-singular vector $\bra{\psi_{\mu}}$, i.e. $\bra{\psi_{\mu}}\neq\bra{s^{R}_{\mu}}^{\otimes 2}$ for any $\bra{s^{R}_{\mu}}\in\mathds{C}^{d}$. For example, if we can engineer a unitary $U$ where $\bra{s^{R}_{\mu}}^{\otimes 2}U=\bra{\psi_{\mu}}$ and $U^{\dagger}\ket{s_{\mu}^{L1}s_{\mu}^{L2}}=\ket{\phi_{\mu}}$, then applying $U$ stroboscopically each time before turning on dissipative interaction [see Fig. 4(b)] would allow engineering of $c_{\mu}=|{\phi_{\mu}}\rangle\langle{\psi_{\mu}}|$.

V.2 Generalizing the parallelized protocol

We can also extend our strategy of parallelized connection and quantum feedback to the class of translation-invariant MPSs with internal symmetry. Recall that the idea is to prepare many segments of the desired MPS with open boundary conditions, and then connecting adjacent pairs of segments in parallel. The analysis is much simpler when the MPS is injective, which means that largest-magnitude eigenvalue of the transfer matrix $\boldsymbol{T}=\sum_{s}A^{(s)*}\otimes A^{(s)}$ is non-degenerate. When the desired MPS is injective, which is true for generic cases Perez-Garcia et al. (2007) (and also for AKLT states), we can show that the success probability of connection is typically at least $1/D^{2}$. This system-size-independent success probability implies that the scaling of preparation time of the parallelized protocol is $\mathcal{O}(\log^{2}n)$, which exponentially outperforms existing dissipative protocols that do not involve parallelization and feedback Verstraete et al. (2009). The case of non-injective MPSs is more subtle, since the associated parent Hamiltonian of the MPS has degenerate ground states. These states can be shown to be the only steady states of our protocol, but typically a mixture of them will be prepared. Nevertheless, if we place some (often reasonable) restrictions on the initial state, a pure non-injective MPS can be prepared, as we illustrate with the example of a Greenberger-Horne-Zeilinger (GHZ) state.

Injective case— We first analyze the protocol for the case of injective MPSs. Consider an arbitrary initial state of two adjacent length-$m$ chains of MPSs, which can be written as $\ket{\psi_{0}}=\mathcal{N}_{0}\sum_{b,c}C_{bc}\ket{A_{ab}^{m}}\otimes\ket{A_{cd}^{m}}$, where $C_{bc}\in\mathds{C}^{D\times D}$ is some coefficient matrix that characterizes the edge states at the interface of the chains, and $\mathcal{N}_{0}$ is a normalization constant. If the two chains are unentangled, then $C_{bc}=\alpha_{b}\beta_{c}$ for some $\vec{\alpha},\vec{\beta}\in\mathds{C}^{D}$. By turning on the jump operators acting at the interface, we can cool this state into the desired final state $\ket{\psi_{f}^{2m}}=\mathcal{N}_{f}\ket{A_{ad}^{2m}}$, for some normalization constant $\mathcal{N}_{f}$. Since the MPS is assumed to be injective, the success probability of connection can be shown to be

[TABLE]

where $\epsilon_{2}$ is the second largest eigenvalue of $\boldsymbol{T}$ (see Appendix H for more details). For random states $\vec{\alpha},\vec{\beta}\in\mathds{C}^{D}$, we have on average $p\simeq 1/D^{2}$. The maximum success probability of $p_{\text{max}}=1/D$ is obtained when $\vec{\alpha}\parallel\vec{\beta}$, i.e., when the two edge states are identical.

When our desired MPS exhibits bond-inversion symmetry $\mathcal{P}$, as the AKLT states do, we have the same issue of vanishing success probability after a quantum jump that also respects $\mathcal{P}$. To see this, consider what happens in the event of a quantum jump due to a jump operator of the form $c=|{\phi}\rangle\langle{\psi}|$. The state after discarding the two particles at the interface is $\ket{\psi_{1}}=\mathcal{N}_{1}\sum_{b^{\prime},c^{\prime}}\tilde{C}_{b^{\prime}c^{\prime}}\ket{A_{ab^{\prime}}^{m-1}}\ket{A_{c^{\prime}d}^{m-1}}$, where

[TABLE]

If $\ket{\psi}$ respects bond-inversion symmetry, i.e. $\mathcal{P}\ket{\psi}=\pm\ket{\psi}$, where $\mathcal{P}=\sum_{i,j}|{ij}\rangle\langle{ji}|$ is the swap operator, then

[TABLE]

We find that this quantity is zero regardless of the initial $C_{bc}$, entangled or unentangled, due to our requirement that $\ket{\psi}$ be orthogonal to the desired MPS $\ket{A_{cb}^{2}}$. In fact, $\operatorname{tr}(\tilde{C})$ is related to the success probability of the next connection attempt, p^{\prime}=\big{|}\braket{\psi_{f}^{2m-2}}{\psi_{1}}\big{|}^{2}\propto|\operatorname{tr}(\tilde{C})|^{2}+\mathcal{O}(\epsilon_{2}^{m})=\mathcal{O}(\epsilon_{2}^{m}), which is exponentially small for a large system size $m$.

Similar to the AKLT case, we can also try to restore the success probability by applying a global symmetry operation $U_{g}^{\otimes m-1}$ for some $g\in\mathcal{G}$ to one of the chains. In our AKLT protocol, there is a symmetry operation $U_{g}=e^{i\pi S_{y}}$ whose action on the virtual bond level $u_{g}=e^{-i\pi\sigma_{y}/2}$ yields \big{|}\operatorname{tr}(u_{g}^{\dagger}\tilde{C})\big{|}^{2}=\operatorname{tr}(\tilde{C}^{\dagger}\tilde{C}), allowing us to recover the maximum success probability of $p_{\text{max}}=1/D$ regardless of the initial state $C_{bc}$ or which quantum jump occurred. While the existence of such an operation is not known for the general case, we can at least restore the success probability to $1/D^{2}$ for many injective MPSs by applying a randomly chosen symmetry operation. This is because injectivity is typically associated with the irreducibility of the representation $u_{g}$ on the virtual bond level Sanz et al. (2009), which allows us to show that $p=1/D^{2}$. If we write $\tilde{C}=\sum_{i}\sqrt{\lambda_{i}}|{b_{i}}\rangle\langle{a_{i}}|$ in its singular value decomposition form, then

[TABLE]

If $u_{g}$ is irreducible, then by Schur’s lemma Artin (2010) we have

[TABLE]

Hence, this yields a subsequent success probability of $p\simeq 1/D^{2}$. A sufficient condition for the irreducibility of $u_{g}$ is that $U_{g}$ be irreducible and $\{A^{(s){\dagger}}A^{(s^{\prime})}:\forall s,s^{\prime}\}$ spans the whole space of $D\times D$ matrices (see Proposition 17 in Ref. Sanz et al. (2009)).

Non-injective case— For a non-injective MPS, the analysis is complicated by the presence of degenerate ground states of its parent Hamiltonian. Here we consider an illustrative example to prepare GHZ states: $\ket{\mathrm{GHZ}_{\pm}}=(\ket{0^{n}}\pm\ket{1^{n}})/\sqrt{2}$. These states have an MPS representation with $(d,D)=(2,2)$ given by the following matrices:

[TABLE]

In this representation, $\ket{\mathrm{GHZ}_{+}}\propto\ket{A^{n}_{\rightarrow\rightarrow}}=\ket{A^{n}_{\leftarrow\leftarrow}}$, and $\ket{\mathrm{GHZ}_{-}}\propto\ket{A^{n}_{\rightarrow\leftarrow}}=\ket{A^{n}_{\leftarrow\rightarrow}}$, where $\ket{\rightarrow}=(\ket{\uparrow}+\ket{\downarrow})/\sqrt{2}$ and $\ket{\leftarrow}=(\ket{\uparrow}-\ket{\downarrow})/\sqrt{2}$ are possible edge configurations. This MPS has an internal symmetry group of $\mathcal{G}=\mathbb{Z}_{2}$, which is represented by $\{\mathds{1},\sigma_{x}^{\otimes n}\}$ acting on the system. Its parent Hamiltonian is $H_{\mathrm{GHZ}}=\sum_{i}(\mathds{1}-\sigma_{z}^{(i)}\sigma_{z}^{(i+1)})$, whose ground states are doubly degenerate due to non-injectivity. The corresponding two-particle bright manifold is $\mathfrak{B}=\text{span}\{\ket{\Phi_{+}},\ket{\Phi_{-}}\}$, where $\ket{\Phi_{\pm}}=(\ket{01}\pm\ket{10})/\sqrt{2}$. The two states $\ket{\Phi_{\pm}}$ support two distinct irreducible representations of $\mathds{Z}_{2}$, which are the trivial and the sign representation, respectively. Hence, we can use just one jump operator of the form e.g. $c=\ket{00}(\kappa_{+}\bra{\Phi_{+}}+\kappa_{-}\bra{\Phi_{-}})$, with $\kappa_{\pm}\neq 0$ so that both irreps are supported (a necessary condition as shown in Appendix H.3). Then along with the global symmetry operation $\sigma_{x}^{\otimes n}$, we can depopulate the bright manifold and obtain $\text{span}\{\ket{\mathrm{GHZ}_{\pm}}\}$ as the subspace of steady states.

Now let us consider preparing $\ket{\mathrm{GHZ}_{\pm}}$ in a parallelized protocol with connections and feedback. We note that unlike in the injective case, different choices of jump operator here can lead to qualitatively different outcomes. Specifically, we consider two choices of jump operators that may result in different degrees of entanglement of the final state. First, consider an example choice of jump operator $c=|{00}\rangle\langle{01}|$ (i.e. $\kappa_{\pm}=1/\sqrt{2}$). While this along with the symmetry operation produces a dissipative dynamics that has $\ket{\mathrm{GHZ}_{\pm}}$ as the steady states, the parallelized protocol can only produce an unentangled final state of either $\ket{0^{n}}\propto\ket{\mathrm{GHZ}_{+}}+\ket{\mathrm{GHZ}_{-}}$ or $\ket{1^{n}}\propto\ket{\mathrm{GHZ}_{+}}-\ket{\mathrm{GHZ}_{-}}$ once any quantum jump occurs, regardless of states of the initial chains. Alternatively, we may choose the jump operator $c=\ket{00}(\bra{01}+i\bra{10})/\sqrt{2}$ (i.e., $\kappa_{+}=\kappa_{-}^{*}=(1+i)/2$). In this case, suppose we start with $\ket{0^{n_{0}}}+\ket{1^{n_{0}}}$ on the initial chains of length $n_{0}$, then we can produce a maximally entangled final state of $\ket{0^{n}}+\zeta\ket{1^{n}}$ even after quantum jumps, for some $\zeta\in\{\pm 1,\pm i\}$ that we can determine from recording quantum jump history. In both cases, for an arbitrary (unentangled) initial state $\ket{\psi_{0}}=\sum_{a,b=0}^{1}(\alpha_{a}\ket{a^{m}})\otimes(\beta_{b}\ket{b^{m}})$ of two chains, the success probability of connecting them is on average $1/2$ for random $\vec{\alpha},\vec{\beta}\in\mathds{C}^{2}$. This system-size-independent success probability means that the parallelized protocol for this non-injective MPS also has an efficient scaling of $\mathcal{O}(\log^{2}n)$ for the preparation time.

VI Summary and Outlook

In this work, we have shown that using only one type of nearest-neighbor decay channel and global control, an AKLT state in a large system can be efficiently prepared as the provably unique steady state of the driven-dissipative dynamics. Since the symmetry group of AKLT states is continuous, the decay channels and pulses need not be fine-tuned, and the implementation is robust against imperfections in experimental parameters. The scaling of preparation time was numerically shown to be polynomial, and can be improved exponentially using parallelization and quantum feedback. This proposal is feasible for a wide range of controlled quantum systems, and in particular optimal for a rearrangeable array of trapped neutral atoms using a Rydberg-EIT scheme. In addition, we analyzed the generalization to other translation-invariant MPSs with symmetry, and derived a bound on the minimum number of necessary decay channels that may be saturated by a construction. We also showed that the parallelized protocol can work in the general case, and provided some sufficient conditions for success, such as the injectivity of the MPS and the irreducibility of the symmetry group representation on the virtual bond level.

Finally, we note that it may be possible to generalize our symmetry-based dissipative preparation scheme to higher dimensional many-body entangled states. Many such states are described by projected entangled pair states (PEPS), a natural generalization of MPSs for arbitrary lattices, which also allow construction of frustration-free parent Hamiltonians Perez-Garcia et al. (2008) to be converted into jump operators Verstraete et al. (2009); Kraus et al. (2008); Ticozzi and Viola (2012). However, our inductive proof of uniqueness of steady states does not extend straightforwardly, since exact computation of expectation values of a generic PEPS is intractable Schuch et al. (2007). Further investigations are thus necessary to extend our strategy to higher dimensions, which can support even more interesting, long-range entangled states with symmetry-enriched topological order Mesaros and Ran (2013).

Acknowledgements.

We thank H. Pichler and F. Reiter for useful discussion. We thank Haoyang Gao for bringing to our attention a technical problem in a previous version of the manuscript. This work was supported through NSF, CUA, Vannevar Bush Faculty Fellowship, AFOSR Muri, and Moore Foundation. L.Z. was supported by NSF Graduate Research Fellowship under Grant No. DGE1144152. S.C. was supported by Kwanjeong Educational Foundation.

Appendix A Effective Liouvillian

Our proposal for preparing AKLT states uses only one type of jump operator, e.g. $c=|{00}\rangle\langle{++}|$. The dissipative dynamics due to this jump operator is $\mathcal{L}_{0}=\Gamma\sum_{i}\mathcal{D}[c^{(i)}]$, where $i$ enumerates the pair of neighboring sites $(i,i+1)$. Our key idea is to use coherent global manipulations, corresponding to operations in the symmetry group, so as to effectively realize additional jump operators. This can be achieved either by periodically applying pulsed global spin rotations (symmetry operations), or by continuously rotating the spins.

In the first, multi-pulse sequence approach, we apply pulses $V_{\theta}=(e^{i\theta S_{y}})^{\otimes n}$, each separated by an interval of $\tau$. Then in the rotating frame, we have the time-dependent Liouvillian:

[TABLE]

Suppose we choose $\theta=2\pi/\ell$ for some integer $\ell$, then this dynamics is periodic with period $\ell\tau$ , since $V_{2\pi/\ell}^{\ell}=\mathds{1}$. In the limit of fast pulses $\tau\ll 1/\Gamma$, we can use the first-order Magnus expansion to derive a simpler, effective time-independent Liouvillian that approximate the dynamics:

[TABLE]

where $\bar{c}_{\nu}=(V_{2\pi/\ell}^{\dagger})^{\nu}\,c\,V_{2\pi/\ell}^{\nu}$.

Alternatively, we may employ a continuous wave approach by introducing a time-independent Hamiltonian $H_{\textrm{CW}}=\omega\sum_{i}S_{i}^{y}$. Then in the rotating frame, we have $c(t)=e^{i\omega tS_{y}}ce^{-i\omega tS_{y}}$, and

[TABLE]

In this frame, the dynamics is periodic with period $2\pi/\omega$. Again, we compute the effective time-independent Liouvillian

[TABLE]

obtained by time-averaging (first-order Magnus expansion) as an approximation. The effective jump operator $\bar{c}_{\beta}$ in the standard form of a Liouvillian can be obtained by diagonalizing the superoperator acting on the space of density operators. More explicitly, we diagonalize a Hermitian matrix $L=U^{\dagger}\Lambda U$, whose entries $L_{ii^{\prime},jj^{\prime}}$ are given by

[TABLE]

We can read off $\Gamma_{\beta}=\Lambda_{\beta}$, and $\bar{c}_{\beta}=\sum_{k,k^{\prime}}U_{\beta,kk^{\prime}}|{k}\rangle\langle{k^{\prime}}|$. For our example choice of $c(t=0)=\ket{00}\bra{++}$, we obtain $\ell=9$ independent jump operator after evaluating the integral and diagonalizing, each with rate $\Gamma_{\beta}/\Gamma=7/32$, 3/16, 3/16, 1/8, 1/8, 1/16, 1/16, 1/64, and 1/64, respectively.

Appendix B Quantum evolution conditioned on no quantum jump

To understand how our quantum state evolves conditioned on detecting no quantum jumps, we use the stochastic wavefunction formalism for open system dynamics. Namely, we define an effective non-Hermitian Hamiltonian $H_{\textrm{eff}}=H-i\sum_{\mu}\Gamma_{\mu}c_{\mu}^{\dagger}c_{\mu}/2$, where $H$ is the Hamiltonian of the system and $c^{(i)}$ are quantum jump operators. A system evolves under $H_{\textrm{eff}}$ until it undergoes a quantum jump $\ket{\psi}\mapsto c_{\mu}\ket{\psi}$ at a rate $\bra{\psi}\Gamma_{\mu}c_{\mu}^{\dagger}c_{\mu}\ket{\psi}$. In our protocol to prepare AKLT states, $H=0$ since we work in the rotating frame. Hence, $H_{\textrm{eff}}$ is anti-Hermitian and thus diagonalizable with eigenvalues $\lambda_{\alpha}=-i\gamma_{\alpha}/2$. It is assumed that our desired states are “dark states”, which are eigenvectors of $H_{\textrm{eff}}$ with a zero imaginary part of the eigenvalue. For simplicity and illustrative purposes, let us also assume there is just one dark state $\ket{0}$ with eigenvalue $\lambda_{0}=0$, and the rest of the eigenvalues are sorted by $0<\gamma_{1}\leq\gamma_{2}\leq\cdots$. Let us decompose the initial state in the eigenbasis $\ket{\psi_{0}}=\sum_{\alpha}c_{\alpha}\ket{\alpha}=c_{0}\ket{0}+c_{1}\ket{1}+\cdots$. The evolution under $H_{\textrm{eff}}$ yields the unnormalized state

[TABLE]

The probability of undergoing no quantum jump over a time duration $T$ is

[TABLE]

Conditioned on such an event, the fidelity of the quantum state preparation is

[TABLE]

where we find that the fidelity exponentially approaches unity. This allows for effective “cooling” of the system into the desired state when there is no quantum jump for $T\gg 1/\gamma_{1}$, which occurs with probability $p_{0}\simeq\left|\braket{0}{\psi_{0}}\right|^{2}$. It is also easy to see that when there are multiple dark states, given by some projector $P_{D}$, the system is cooled into $\ket{\psi_{0}}\to P_{D}\ket{\psi_{0}}$, conditioned on no quantum jumps which occurs with probability $p_{0}=\braket{\psi_{0}}{P_{D}}{\psi_{0}}$.

Appendix C Proof of uniqueness of steady states for the AKLT example

In this appendix, we provide a detailed proof of our lemma that $\bar{\mathcal{L}}_{\textrm{MP}}$ and $\bar{\mathcal{L}}_{\textrm{CW}}$ have AKLT states as unique steady states, for $n\geq 2$ under the open boundary condition. Although the proof has already been sketched in the main text, here we provide a more detailed analysis. Let us define $\mathcal{L}_{i}=\sum_{\nu=0}^{\ell-1}\Gamma_{\nu}\mathcal{D}[\bar{c}_{\nu}^{(i)}]$ as the Liouvillian acting only on sites $i$ and $i+1$, and $\mathcal{L}_{[n]}=\sum_{i=1}^{n-1}\mathcal{L}_{i}$. Let $Q_{i,i+1}$ be the projector onto the $J=2$ manifold on spin $i$ and $i+1$. We also define

[TABLE]

as the AKLT Hamiltonian acting on $n$ spins under the open boundary condition. We also let $D_{n}=\text{span}\{\ket{G_{ab}^{n}}\}$ be the space of AKLT states on $n$ spins, and let $P_{n}$ be its corresponding projector. We prove by induction on $n$.

We start with our inductive hypothesis, which is that $\mathcal{L}_{[n]}$ only has the AKLT states $\ket{G_{ab}^{n}}$ as steady states. This means that any $n$-particle wavefunction can reach $D_{n}$ via some sequence of jump operators $\bar{c}_{\nu}$. More formally, this means that for any $n$-particle wavefunction $\ket{\psi^{n}}$, there exists a polynomial function of jump operators $f_{[n]}(\{\bar{c}_{\nu}\})$ such that

[TABLE]

We only need to show that any $(n+1)$-particle wavefunction $\ket{\psi^{n+1}}$ can reach $D_{n+1}$ via jump operators, which would yield the proof of our lemma.

From this inductive hypothesis, we first argue that there exists a polynomial $f_{[n]}(\{\bar{c}_{\nu}\})$ of jump operators acting on the first $n$ spins such that

[TABLE]

This can be done by simply looking at the the reduced density matrix $\rho_{n}$ of $\ket{\psi^{n+1}}$ on the first $n$ spins, and picking $f_{[n]}$ with respect to one of the non-trivial eigenvectors of $\rho_{n}$. We note that the AKLT Hamiltonian is frustration free, so any (ground) state minimizing the energy of $H^{[n]}_{\textrm{AKLT}}$ must also minimize the energy of $H^{[n-1]}_{\textrm{AKLT}}$, or $P_{n}P_{n-1}=P_{n}$. Hence (38) implies $f_{[n]}\ket{\psi^{n+1}}$ also contains nonzero population in the AKLT ground states among the first $n-1$ spins, i.e., $\|P_{n-1}f_{[n]}\ket{\psi^{n+1}}\|\neq 0$.

Let us then perform a general decomposition

[TABLE]

where

[TABLE]

Here $\ket{z}$ runs through the 5 states in the $J=2$ manifold on the last two spins. The remaining parts of the wavefunction $f_{[n]}\ket{\psi^{n+1}}$ are

[TABLE]

where $\ket{E_{\mu}^{n-1}}$ runs through all the excited eigenstate of the $H^{[n-1]}_{\textrm{AKLT}}$, and $s$ runs through all 9 possible 2-spin states.

We assume for the sake of contradiction that $\ket{\psi_{n+1}}$ is a state that cannot reach $D_{n+1}$ via jump operators. We now consider two cases:

Case 1: $\phi_{abz}$ = 0 for all $a,b,z$. — Then $f_{[n]}\ket{\psi^{n+1}}=\ket{\phi^{\perp}_{1}}+\ket{\phi^{\perp}_{2}}$. From our inductive hypothesis, $f_{[n]}\ket{\psi^{n+1}}$ must have nonzero population in $D_{n}$ and $D_{n-1}$, and this can only come from the $\ket{\phi^{\perp}_{1}}$ part, and thus $\phi_{abst}^{\perp}$ cannot be all vanishing. On the other hand, since we have assumed $\ket{\psi^{n+1}}$ cannot reach $D_{n+1}$, we must have $0=\braket{G_{pq}^{n+1}}{\phi^{\perp}_{1}}$. That means $Q_{i,i+1}\ket{\phi^{\perp}_{1}}\neq 0$ for some $i$. Since $Q_{i,i+1}\ket{\phi^{\perp}_{1}}=0$ for all $i\neq n-1$, this means that we must have $Q_{n-1,n}\ket{\phi^{\perp}_{1}}=\ket{\phi^{\perp}_{1}}$. Then we have $P_{n}f_{[n]}\ket{\psi^{n+1}}=P_{n}\ket{\phi^{\perp}_{1}}+P_{n}\ket{\phi^{\perp}_{2}}=0$, since $P_{n}Q_{n-1,n}=0$ and $P_{n}\ket{E_{\mu}^{n-1}}=0$. This gives a contradiction with (38).

Case 2: Some $\phi_{abz}$ $\neq$ [math]. — Consider any jump operator $\bar{c}_{\nu}^{(n)}$ acting on the spins $n$ and $n+1$. Note that

[TABLE]

because $\ket{G_{pr}^{n-1}}$ are ground states of $H^{[n]}_{\textrm{AKLT}}$ and thus orthogonal to $\ket{E_{\mu}^{n-1}}$. We also have $\braket{G_{pq}^{n+1}}{\bar{c}_{\nu}^{(n)}}{\phi^{\perp}_{1}}=0$ since $\bar{c}_{\nu}\ket{G_{st}^{2}}=0$. Then our assumption that $\ket{\psi^{n+1}}$ cannot reach $D_{n+1}$ gives the condition

[TABLE]

where $\vec{\phi}$ is a $2\times 2\times 5=20$-dimensional vector, and $\boldsymbol{M}$ is a $4\ell\times$20 dimensional matrix. The matrix elements of $\boldsymbol{M}$ are given by

[TABLE]

and can be calculated analytically, since the first factor comes from diagonalizing the transfer matrix, and the second factor is computed in a small, 9-dimensional Hilbert space of two spins. Now if $\det(\boldsymbol{M}^{\dagger}\boldsymbol{M})\neq 0$, then the matrix $\boldsymbol{M}$ has full rank, indicating that we only have the trivial solution $\phi_{abz}=0$. This would contradict this case’s assumption, and hence any $\ket{\psi_{n+1}}$ must be able to reach $D_{n+1}$ via jump operators, and the AKLT states $\ket{G_{ab}^{n+1}}$ are the unique steady states of $\mathcal{L}_{[n+1]}$.

Therefore, it only remains to compute $\det(\boldsymbol{M}^{\dagger}\boldsymbol{M})$ and show that it is nonzero for $\mathcal{L}$ in our proposals. Let us first consider $\bar{\mathcal{L}}_{\textrm{MP}}$, with $\ell=5$ corresponding to rotation pulses with $\theta=2\pi/5$. For this, we explicitly find

[TABLE]

where $x=(-3)^{n}$. It’s easy to see that this is only zero for $n=1,2$, so our inductive proof holds for $n\geq 3$. The base case of $n=2,3$ can easily be checked numerically or by exact calculations.

For $\bar{\mathcal{L}}_{\textrm{CW}}$, which has $\ell=9$, we find

[TABLE]

where again we have defined $x=(-3)^{n}$. This is also only nonzero when $n=1,2$. Since the base cases of $n=2,3$ can be checked exactly, this proves that the steady states of $\bar{\mathcal{L}}_{\textrm{MP}}$ and $\bar{\mathcal{L}}_{\textrm{CW}}$ are uniquely given by $D_{n}$ for $n\geq 2$.

Our proof method here naturally suggests a method to prepare an AKLT state with specified edge states instead of a mixture of the four $\ket{G_{ab}^{n}}$ under open boundary condition. For instance, $\ket{G_{\uparrow\uparrow}^{n}}$ can be deterministically prepared by adding two jump operators: $c_{L}$ = $|{0}\rangle\langle{-}|_{1}$ on the left edge and $c_{R}$ = $|{0}\rangle\langle{+}|_{n}$ on the right. In this case, it is easy to see that any linear combination of four ground states $\ket{G_{ab}^{n}}$ can further decay into $\ket{G_{\uparrow\uparrow}^{n}}$, which becomes the unique steady state.

Appendix D Details of our numerical simulation

We simulate the dissipative dynamics $\mathcal{L}=\sum_{\mu}\Gamma_{\mu}\mathcal{D}[c_{\mu}]$ using the stochastic wavefunction method Gardiner and Zoller (2004). In this approach, the wavefunction $\ket{\psi(t)}$ continuously evolves under the effective non-Hermitian Hamiltonian $H_{\textrm{eff}}=H-i\sum_{\mu}\Gamma_{\mu}c_{\mu}^{\dagger}c_{\mu}/2$ and stochastically undergoes quantum jumps $c_{\mu}$ at a rate $\braket{\psi}{\Gamma_{\mu}c_{\mu}^{\dagger}c_{\mu}}{\psi}$. Physical observables are extracted from an ensemble of wavefunction trajectories obtained from independent simulations. Compared to direct numerical integrations of quantum master equations, this method allows simulation of systems with a larger number of particles. In order to simulate a maximally mixed initial state (or equivalently an infinite temperature ensemble), we sample a random product state in the $S_{z}$-basis as the initial state $\ket{\psi_{\alpha}(0)}$ for the $\alpha$-th simulation. For each small time step $\delta t$, the state $\ket{\psi_{\alpha}(t)}$ evolves stochastically according to either

• •

$\ket{\psi_{\alpha}(t+\delta t)}\propto c_{\mu}\ket{\psi_{\alpha}(t)}$ with probability $\delta p_{\mu}=\braket{\psi_{\alpha}(t)}{\Gamma_{\mu}c_{\mu}^{\dagger}c_{\mu}}{\psi_{\alpha}(t)}\delta t$, or

• •

$\ket{\psi_{\alpha}(t+\delta t)}\propto e^{-iH_{\textrm{eff}}\delta t}\ket{\psi_{\alpha}(t)}$ with probability $1-\sum_{\mu}\delta p_{\mu}$.

We choose $\delta t$ so that $\delta p=\sum_{\mu}\delta p_{\mu}\ll 1$. Since this process is stochastic, we average over a sufficiently large number $N_{\text{traj}}$ of trajectories to estimate the values of observables:

[TABLE]

For the numerical data presented in the main text, we average over up to $N_{\text{traj}}=1000$ trajectories, and statistical uncertainties are estimated using the bootstrapping technique Efron (1981). For relatively large system sizes $n>8$, numerical computations of exact many-body wavefunctions are impractical. Instead, we store the wavefunction in an MPS representation, and simulate the evolution using the time-evolving block decimation algorithm Vidal (2003); *Vidal2. Since we are dissipatively preparing AKLT states that have bond dimension $D=2$, we find that restricting the maximum bond dimension of our MPS wavefunction to $D\leq 15$ is sufficient, as truncation errors are found to be $<5\times 10^{-5}$ in all simulations. This algorithm allows us to simulate systems with up to $n=25$ spins.

During our simulated evolution, we monitor two observables: (1) energy density $\braket{H_{\textrm{AKLT}}}/(n-1)$, and (2) fidelity of state preparation $\mathcal{F}=\braket{P_{G}}$ where $P_{G}$ is the projector onto AKLT ground states. Note that since our simulation in the main text is for the open boundary condition, there are four degenerate ground states that are all accepted as output; we define our fidelity $\mathcal{F}$ to be the sum of overlap with each accepted state. Another widely used measure on quantum states, trace distance Gilchrist et al. (2005), cannot be applied in this context, because it measures how close a state is to another (target) state, not to a subspace of such states. Even in a situation where the target state is a single pure state, e.g. under periodic boundary condition, the trace distance $T$ is bounded by our fidelity $\mathcal{F}$ through $1-\sqrt{\mathcal{F}}\leq T\leq\sqrt{1-\mathcal{F}}$. In our numerics where the system size goes up to $n=25$, the computational cost of calculating the density matrix and trace distance would also be prohibitively expensive. Finally, this fidelity $\mathcal{F}$ coincides with the success probability of state preparation, which is a physically meaningful metric.

Appendix E Effect of imperfect quantum jump detection

In realistic experiments, the detection of quantum jumps often entails imperfections. The presence of such imperfections affects our protocol by (i) not heralding the failure of connection of two chains (false-positive), and (ii) incorrectly heralding failure when the connection has been successful (false-negative). The former may arise due to imperfect detection efficiency, and the latter due to the dark counts in the detector. As mentioned in the main text, our parallelized protocol can still have an efficient scaling even when such imperfect quantum jump detection is accounted for. In the false-negative scenario, we can still discard affected particle pairs and continue the procedure, but we have to adopt an ideal case success probability lower than $p_{\text{max}}$ (discussed below). To minimize the occurrence of false-positives, we propose two methods that address detector inefficiency. In the following analysis, let us denote the detector efficiency by $1-\eta$, the dark count rate by $r$, and the ideal case success probability by $p$.

Success probability after false-negatives— We showed previously in Sec. III.1 that the maximum success probability of connection of $p_{\text{max}}=1/2$ can be recovered for the subsequent attempt after a failed connection (via detection of quantum jumps), if we discard the affected particles and apply a global $\pi$-rotation $U=e^{i\pi S_{y}}$ to one of the remaining chains. However, this is not the case if the quantum jump detector has only received a dark count, i.e. the failure is a false-negative, and the connection has in fact succeeded. After a dark count is registered, the experimenter will unwittingly discard the particles at the original interface, apply $U$, and retry the connection anyways. If we consider the density matrix of the remaining pair of chains, we can see that their edge states at the interface are essentially randomized in a maximally mixed state, which would intuitively yield a success probability $p\approx 1/4$ in the attempt to connect them. More precisely, if $k$ particles were discarded in each original chain of length $n$ starting from interface (i.e. $2k$ particles in the middle of the connected chain of length $2n$), the success probability of connecting the two chains of length $n-k$ can be computed to be

[TABLE]

where $\epsilon=-1/3$. In particular, when $k=1$, $p\simeq 2/9$.

Method 1 to mitigate false-positives— The first method to address false-positives from detector inefficiency is to only use jump operators of the form $c=|{++}\rangle\langle{++}|$ for the connection. Once a quantum jump occurs, the state will continue to undergo quantum jumps indefinitely, creating a much larger signal and effectively larger detection efficiency. Consequently, the detector inefficiency can be exponentially suppressed by the time $\tau_{c}$ of having the jump operators turned on. Let $\tau_{0}$ be the time-scale in which a single quantum jump would occur. In this case, the probability of diagnosing a successful connection and keeping the result is given by the probability of not detecting any quantum jumps over time $\tau_{c}$,

[TABLE]

where we defined $a=(1-p)/p$ and $\tilde{\eta}=\eta/(1-r\tau_{0})$. The fidelity is the conditional probability that the diagnosed success was a truly successful connection

[TABLE]

Note that we can only achieve a fidelity arbitrarily close to 1 if $\tilde{\eta}<1$, i.e. when the detector efficiency is larger than the dark count probability $1-\eta>r\tau_{0}$. In order to achieve a final error of $\mathcal{E}$ for a system size of $n$ from initial chains of length $n_{0}$, where $n/n_{0}$ connections are necessary, we need $1-\mathcal{F}\leq n_{0}\mathcal{E}/n$, and consequently $\tau_{c}=\mathcal{O}(\ln(n/n_{0}\mathcal{E})/\ln\tilde{\eta}^{-1})$. This is consistent with the $\tau_{c}\sim\log n$ scaling necessary in the ideal protocol. Nonetheless, our new success probability now decreases with system size $n$ as $p_{\text{succ}}\sim\mathcal{O}(n^{-\delta})$ if dark counts are non-negligible, where $\delta=\ln(1-r\tau_{0})/\ln\tilde{\eta}\approx r\tau_{0}/\ln\eta^{-1}$. Consider now the average time to prepare a chain of length $n$:

[TABLE]

where $n_{c}=2(1-p_{\text{succ}})/p_{\text{succ}}$, $T_{0}$ is the time to prepare initial length-$n_{0}$ chains, and $\tau_{r}$ is some constant time necessary to reset the edge states in the event of failure. At first sight, this indicates that our preparation time would ultimately scale polynomially instead of polylogarithmically in the infinite $n$ limit. However, in the regime of $r\tau_{0}\ll 1$, this polynomial dependence has a very small power, and its effect can be neglected if $r\tau_{c}\ll 1$. Hence, in practice, our protocol has an efficient, polylogarithmic scaling up to $n\ll n_{\text{max}}=\mathcal{O}((1/\eta)^{1/r\tau_{0}})$, beyond which it switches to a polynomial scaling. For instance, even if the single-photon detection efficiency is $1-\eta=0.2$, then assuming a dark count rate of $r=25$ Hz Hadfield (2009) and a quantum jump scattering rate of $\tau_{0}^{-1}\approx 1$ MHz, it takes an astronomically long chain of $n_{\text{max}}\sim 10^{4000}$ to reach the polynomial scaling. An example scaling under these conditions is shown in Fig. 5.

Method 2 to mitigate false-positives— The second method for addressing detector inefficiency is to slowly turn on jump operators in the vicinity of the interface. In this way, the absence of quantum jumps further confirms that the two chains have indeed been successfully connected; since only a successful connection does not lead to any subsequent quantum jumps, any false-positive diagnosis of successful connection can be corrected. More concretely, consider a $k$-step scheme where we turn on jump operators to include $k$ neighbors on each side of the original interface, one pair of neighbors at a time. At step $\ell=1,\ldots,k$, we have jump operators on for $2\ell$ particles centered at the interface, turned on for time $\tau^{c}_{\ell}$. If a quantum jump occurred and evaded detection at any step $\ell$, we assume $\tau^{c}_{\ell}$ is long enough so that the $2\ell$ particles would have formed a connected chain of length $2\ell$. At the subsequent step $\ell+1$, the length-$2\ell$ chain in the middle can be connected to the two length-$(n-\ell)$ chains on both sides if we succeed by having no quantum jump, producing a fully connected chain of length $2n$. Note that the success probability for steps $\ell>1$ is roughly $p_{2}=1/2^{4}$. From the system size scaling found in our numerical simulations presented in the main text, we expect to need $\tau^{c}_{\ell}\sim(2\ell)^{2.97}\approx(2\ell)^{3}$. Additionally, we expect the number of quantum jumps during $\tau^{c}_{\ell}$ of step $\ell$ to roughly scale as $N^{\text{jump}}_{\ell}\approx C\tau^{c}_{\ell}$ for some constant $C$. Observe that this scheme allows us to obtain a fully connected chain even in the event of initial failure(s), as long as we do not have any quantum jumps at the last step. Thus, the probability of succeeding and deciding to keep the result (due to not detecting any quantum jump) is

[TABLE]

where $T_{s}=\sum_{\ell=1}^{s}\tau^{c}_{\ell}\approx Bs^{4}$ for some constant $B$. The probability of failing at the last step but still keeping the result is

[TABLE]

The fidelity is the conditional probability of true success given that we have kept the result:

[TABLE]

where we defined $b=(1-p)/p$ and $\xi=(\eta^{C}e^{r})^{B}$. To achieve arbitrarily good fidelity, we require $\xi<1$, i.e. the dark count rate $r<C\ln\eta^{-1}$ needs to be sufficiently small. At the same time, the apparent “success” probability of keeping the result is $p_{\text{succ}}=\Pr\{\text{keep}\}\approx pe^{-rBk^{4}}$. We can carry out the same analysis as in the previous method, and a similar behavior emerges: when the dark count rate is nonzero, the efficient polylogarithmic scaling applies until a maximum chain length of $n\ll n_{\text{max}}=\mathcal{O}(\eta^{-C/r})$, beyond which a polynomial scaling of $\mathcal{O}(n^{\delta^{\prime}})$ with $\delta^{\prime}\approx r/\ln\eta^{-C}$ applies.

Appendix F Analysis of Rydberg-EIT implementation proposal

In this appendix, we derive the effective dissipative interaction between two nearby particles for our Rydberg-EIT implementation scheme introduced in the main text. Consider two particles interacting via the Rydberg shift $H_{\text{int}}=U|{rr}\rangle\langle{rr}|$. Their effective (non-Hermitian) Hamiltonian under the Rydberg-EIT scheme proposed in the main text is

[TABLE]

where $\Delta=\sqrt{\Omega^{2}+g^{2}}$, $\ket{D}=(\Omega\ket{+}-g\ket{r})/\Delta$ is the EIT-dark state, and $\ket{B}=(g\ket{+}+\Omega\ket{r})/\Delta$ is a state orthogonal to $\ket{D}$ that we call EIT-bright state. Now consider a general two-particle wavefunction $\ket{\psi}=\sum_{a}c_{aa}\ket{aa}+\sum_{a<b}c_{ab}(\ket{ab}+\ket{ba})/\sqrt{2}$, where we have restricted ourselves to working in the symmetric subspace. Then the equations of motion for the coefficients are

[TABLE]

In the limit of $U\ll g,\Omega,\gamma$, or $g\ll\Omega,\gamma,U$, and assuming we start initially with $\ket{\psi}=\ket{DD}$, we can adiabatically eliminate the fast dynamics involving coefficients $\{c_{ab}\}$ other than $c_{DD}$. This procedure can be effectively achieved by setting $\dot{c}_{ab}=0$ for ${ab}\neq DD$, allowing us to obtain

[TABLE]

Here, $\textrm{Re}[U_{DD}]$ is the interaction-induced energy shift, and $\Gamma_{DD}=-2\textrm{Im}[U_{DD}]$ is the two-body effective decay rate.

A more general version of adiabatic elimination for open system can be found in Ref. Reiter and Sørensen (2012), which allows us to obtain effective jump operators. Consider original jump operators of the form $L_{s,j}=|{s}\rangle\langle{e}|_{j}$, corresponding to the spontaneous decay from excited state $\ket{e}$ to one of the three hyperfine ground state $\ket{s}$ for $s\in\{+,0,-\}$ in atom $j$. Note that in practice, the excited state can also decay into other hyperfine ground states, which can then be repumped to the excited state using additional lasers. We denote the decay rate corresponding to $L_{s,j}$ by $\gamma_{s}$, where $\gamma_{+}+\gamma_{0}+\gamma_{-}=\gamma$. Then we can compute the effective jump operators:

[TABLE]

with rate $\Gamma_{\textrm{eff}}=\frac{\gamma_{s}}{2\gamma}\Gamma_{DD}$, and $\ket{\tilde{e}}\propto\frac{\Omega\gamma}{2\Delta}\ket{B}-i\Omega\ket{e}-\frac{g\gamma}{\Delta}\ket{D}$. Note that $\ket{\tilde{e}}$ will further decay through the original jump operator $L_{s,j}$. Assuming we are in the regime $\gamma\gg\gamma_{s}\Gamma_{DD}/2\gamma$ so that $\ket{\tilde{e}}$ is a short-lived intermediate state, we can approximate the overall effective dynamics with jump operators of the form $L^{\textrm{eff}}_{ss^{\prime}}=|{ss^{\prime}}\rangle\langle{DD}|$ for $s,s^{\prime}\in\{+,0,-\}$. While we only need a jump operator such as $L^{\textrm{eff}}_{00}=|{00}\rangle\langle{DD}|$ to ensure AKLT states are unique steady states of the engineered dissipative dynamics, the additional effective jump operators do not affect the steady states and can in fact help to more quickly depopulate the undesired states.

Appendix G Scaling of imperfection in Rydberg-EIT implementation

Our protocol prepares AKLT states with finite fidelity when experimental imperfections are taken into account. Here, we analyze how the fidelity scales as multiple chains are connected. In particular, the long-range nature of interaction in the proposed Rydberg-EIT implementation limits fidelity even in the absence of dephasing. Nevertheless, we show that by adopting the parallelized protocol, the long-range interactions only affect the initial preparation of length-$n_{0}$ chains, and such imperfection does not substantially grow in the later connection procedures. This is because the connections involve turning on the dissipative interaction on particle pairs that are spatially separated by at least $n_{0}$ particles. Since the effective decay rate scales as $\Gamma_{DD}\sim 1/R^{12}$ for the proposed implementation, the perturbative effect of long-range interaction is characterized by the very small number of $1/(n_{0}-1)^{12}$, which becomes even smaller in later rounds of connections. Hence, we neglect the effect of long-range interaction on the connections, and only consider how the induced errors on the states of initial chains propagate through the protocol. Let us assume that we initially start with individual chains of length $n_{0}$, each with bounded error $\epsilon_{0}$. At the $\ell$-th level of connections, we on average double the length $n_{\ell}\simeq 2n_{\ell-1}-n_{c}$, where $n_{c}$ is the expected number of particles discarded in each connection. The number of initial chains necessary to reach a final chain length of $n$ is $L=(n-n_{c})/(n_{0}-n_{c})$, and $L-1$ connection procedures need to be performed. Hence, the final error is bounded by

[TABLE]

where we neglect the second term which can be made small compared to the first if we choose $\tau_{c}=\mathcal{O}(\ln n)$. As we can see, the predominant source of error is due to the imperfect initial chains, whose errors add linearly.

This linear scaling of error is indeed very favorable for a many-body state preparation protocol. To put this in perspective, let us estimate how the effective temperature $T_{\textrm{eff}}$ scales in our connection procedure. We define the effective temperature through the relation

[TABLE]

where $\rho(E)$ is the density of states at energy $E$, and we assume that $T_{\textrm{eff}}$ is sufficiently small. Now we consider connecting two length-$n_{1}$ chains at effective temperature $T_{1}$ (with errors $1-\mathcal{F}_{1}$). After connection, we have a chain of length $n_{2}\approx 2n_{1}$, with bounded error $1-\mathcal{F}_{2}\lesssim 2(1-\mathcal{F}_{1})$. The corresponding effective temperature $T_{2}$ of the connected chain can be estimated from

[TABLE]

which implies

[TABLE]

where $\rho_{1}(E)$ and $\rho_{2}(E)$ denote the density of states for chains of length $n_{1}$ and $n_{2}\approx 2n_{1}$, respectively. In a generic many-body interacting system, the density of states grows exponentially in system sizes. Here, we are most interested in the density of states of low-lying excitations, e.g. the first excited band, where the scaling of $\rho(E)$ can be much weaker. Ref. Arovas et al. (1988) used a Bijl-Feynman single-mode approximation to deduce that there is a band of low-lying excited states with dispersion relation $E_{1}(k)=\frac{5}{27}(5+3\cos k)$, corresponding to magnon excitations. Therefore, we expect the number of states in the low-lying excited bands to scale at least linearly with system size, and thus $\rho_{2}(E)\geq 2\rho_{1}(E)$. Applying this to the earlier inequality, we have

[TABLE]

Hence, the effective temperature should not increase (and can potentially decrease) after each connection procedure.

Appendix H Generalization to symmetric MPSs

In this appendix, we provide additional details on how to generalize our protocol for a broader class of translation-invariant MPSs with internal symmetry. We first introduce notation and elaborate on a few useful properties (including injectivity) of translation-invariant MPSs in Sec. H.1. More detailed descriptions and proofs of these properties can be found in Ref. Perez-Garcia et al. (2007). We then discuss the meaning of internal symmetry of MPSs in Sec. H.2. Finally, we describe and analyze the generalization of our protocol in Sec. H.3, and prove the lower bound for the minimum complexity of decay channels.

H.1 Notations and relevant properties of matrix product states

Any (unnormalized) translation-invariant MPS with physical dimension $d$ (i.e. spin-$\frac{d-1}{2}$ particles) and bond dimension $D$ can be written as:

[TABLE]

where $s_{i}\in\{1,2,\ldots,d\}$ runs over the physical spin basis for the $i$-th particle, and $\ket{a},\ket{b}\in\mathds{C}^{D}$ indicate the “boundary conditions”. We denote by $\ket{A^{n}_{ab}}$ the translation-invariant MPS of $n$ particles with open boundary condition specified by $a$ and $b$, and by $\ket{A_{\circ}^{n}}$ the MPS of an $n$-particle system with periodic boundary condition, i.e. $\ket{A_{\circ}^{n}}=\sum_{a}\ket{A_{aa}^{n}}$. Under open boundary condition, there could be at most $D^{2}$ distinct states with different possible boundary conditions, e.g. four-fold degeneracy of AKLT states. In general, however, these $D^{2}$ states may not be linearly independent unless the MPS is injective (defined below).

Canonical Form— An MPS is in a canonical form if the matrices have a common block diagonal structure: $A^{(s)}=\textrm{diag}(\lambda_{1}A_{1}^{(s)},\ldots,\lambda_{B}A_{B}^{(s)})=\bigoplus_{\beta=1}^{B}\lambda_{\beta}A_{\beta}^{(s)}$, where $0<\lambda_{\beta}\leq 1$ for each block $\beta\in\{1,\ldots,B\}$. The matrices in each block must satisfy the conditions that (i) $\sum_{s}A^{(s)}_{\beta}A^{(s){\dagger}}_{\beta}=\mathds{1}$, (ii) a map defined as $\mathcal{E}_{\beta}(X)=\sum_{s}A_{\beta}^{(s)}XA_{\beta}^{(s){\dagger}}$ has $\mathds{1}$ as its only fixed point (unique eigenvector with unity eigenvalue), and finally (iii) $\sum_{s}A_{\beta}^{(s){\dagger}}\Lambda_{\beta}A_{\beta}^{(s)}=\Lambda_{\beta}$ for some diagonal positive and full-rank matrices $\Lambda_{\beta}$. From now on, we assume any MPS under consideration is written in a canonical form.

Transfer Matrix— Consider the completely positive map $\mathcal{E}(X)=\sum_{s}A^{(s)}XA^{(s){\dagger}}$, or equivalently the transfer matrix $\boldsymbol{T}=\sum_{s}A^{(s)*}\otimes A^{(s)}$. Understanding the spectrum of this transfer matrix is useful for computing the expectation value of an observable or the overlap between two quantum states, e.g. $\braket{A_{ab}^{n}}{A_{a^{\prime}b^{\prime}}^{n}}=\braket{aa^{\prime}}{\boldsymbol{T}^{n}}{bb^{\prime}}$ Perez-Garcia et al. (2007). Some known eigenvectors of $\boldsymbol{T}$ are $\frac{1}{\sqrt{D_{\beta}}}\sum_{i\in\beta}\ket{ii}$, where $D_{\beta}$ is the dimension of the $\beta$-th block, and correspond to eigenvalues $|\lambda_{\beta}|^{2}$. Denoting the other eigenvectors of $\boldsymbol{T}$ with eigenvalues $\epsilon_{\nu}$ by $\ket{\nu}$, we have

[TABLE]

Since $|\lambda_{\beta}|^{2}$ is the largest eigenvalue of each block $\beta$, typically only the first term is relevant in the limit of large $n$.

Parent Hamiltonian— For a sufficiently large $L$, the set of matrix products $\{A^{(s_{1})}\cdots A^{(s_{L})}:1\leq s_{i}\leq d\}$ spans the vector space of all matrices with the same block diagonal structure as the canonical form Perez-Garcia et al. (2007). We call $L$ the interaction length of the MPS. Without loss of generality, we can assume that $L=2$. This is because otherwise we can group $L$ sites together to get an equivalent MPS with larger physical dimension $d^{\prime}\leq d^{L}$, and a new interaction length $L^{\prime}=2$. The parent Hamiltonian of an MPS is then defined to be $H_{p}=\sum_{i}h^{(i)}$, where $h$ is any positive semi-definite operator acting on nearest neighboring sites, whose kernel is

[TABLE]

In other words, $H_{p}$ imposes a condition $h^{(i)}$ for every pair of neighboring sites $(i,i+1)$, which our MPS trivially satisfies (i.e. $h^{(i)}\ket{A_{ab}^{n}}=0$ for all $1\leq i\leq n-1$). Hence, $H_{p}$ is a frustration-free Hamiltonian of which the MPS is a zero-energy ground state. The ground-state degeneracy depends on both the boundary condition and the number of blocks in the MPS canonical form. Ref. Perez-Garcia et al. (2007) has shown that under periodic boundary condition, this degeneracy equals the number of independent blocks in the canonical form, since the ground space consists of MPSs constructed from the sub-matrices from every block.

Injectivity— Often it is useful to assume a condition that the MPS is injective, which is satisfied in the generic case except for specific, fine-tuned MPSs Perez-Garcia et al. (2007). This injectivity condition is that the transfer matrix $\boldsymbol{T}$ has only one eigenvector corresponding to its largest eigenvalue (which we normalize to 1 in the canonical form). This also implies that there is just one block in the canonical form of the MPS. In this case, Eq. (65) simplifies to

[TABLE]

where $\epsilon_{2}$ is the second largest eigenvalue of $\boldsymbol{T}$. Additionally, this implies that the parent Hamiltonian under periodic boundary condition has the MPS as its unique ground state, and that the ground state energy is gapped in the thermodynamic limit. Under open boundary condition, $D^{2}$ distinct boundary conditions give rise to $D^{2}$ linearly independent and degenerate ground states $\ket{A^{n}_{ab}}$. By appropriately modifying the parent Hamiltonian terms at the boundaries, we can break the degeneracy and make one of the $D^{2}$ states the unique ground state.

H.2 Internal Symmetries of MPSs

We say a translation-invariant MPS defined on $d$-dimensional physical spins respects an internal symmetry $\mathcal{G}$, if for some unitary representation $U:\mathcal{G}\to U(d)$, we have:

[TABLE]

That is, a global action of the symmetry operation keeps a ground state of the MPS parent Hamiltonian in the ground space under open boundary condition, or only imprints a complex phase factor under periodic boundary condition.

Assuming the symmetry group is reasonable (either a discrete or a compact connected Lie group), but without assuming injectivity, Ref. Sanz et al. (2009) showed that we can replace the action of the symmetry in the physical basis with a unitary in the virtual bond basis. More explicitly, we have

[TABLE]

Here, $u_{g}=P_{g}v_{g}$, with $v_{g}=\bigoplus_{\beta=1}^{B}v_{g}^{\beta}$ taking on the same block diagonal structure as $A^{(s)}$, and each $v_{g}^{\beta}$ a unitary in block $\beta$. $P_{g}$ is a permutation amongst the $B$ blocks. Lastly, $w_{g}=\bigoplus_{\beta}e^{i\varphi_{g}^{\beta}}\mathds{1}_{\beta}$ is a phase factor for each block. If $\mathcal{G}$ is a compact connected Lie group, Ref. Sanz et al. (2009) showed that $P_{g}=\mathds{1}$, while $g\mapsto e^{i\varphi_{g}^{\beta}}$ and $g\mapsto v_{g}^{\beta}$ are representations of $\mathcal{G}$.

For AKLT states, where $\mathcal{G}=\mathrm{SO}(3)$, the relevant representation is given by the rotations $U_{g}=\exp(i\vec{\alpha}_{g}\cdot\vec{S})$ on the spin-1 vector $\vec{S}$ for some real parameters $\vec{\alpha}_{g}=(\alpha_{g}^{x},\alpha_{g}^{y},\alpha_{g}^{z})$. In particular, we have

[TABLE]

where $u_{g}=\exp(i\vec{\alpha}^{\prime}_{g}\cdot\vec{\sigma}/2)$, with $\vec{\alpha}^{\prime}_{g}=(\alpha_{g}^{x},-\alpha_{g}^{y},\alpha_{g}^{z})$.

H.3 Finding a minimal set of decay channels using symmetry

Our goal is to find a minimal set $k_{\textrm{min}}$ of decay channels $\{c_{1},c_{2},\dots,c_{k_{\textrm{min}}}\}$ acting on neighboring pairs of particles that deterministically prepare the desired MPS, assuming global symmetry operations are available. For concreteness we will first focus on decay channels of the form $c_{\mu}=|{\phi_{\mu}}\rangle\langle{\psi_{\mu}}|$. We show there is a lower bound on $k_{\min}$ from the structure of representation of $\mathcal{G}$ on the physical particles, and provide a construction of the jump operators saturating the bound. The uniqueness of steady states under the constructed jump operators can be analytically confirmed using the same inductive proof strategy demonstrated for the case of AKLT states.

Without loss of generality, we assume that the desired states are ground states of a gapped, frustration-free parent Hamiltonian $H_{p}=\sum_{i}h^{(i)}$, where $h^{(i)}$ is a translation-invariant, nearest-neighbor projector that respects the internal symmetry $\mathcal{G}$ Perez-Garcia et al. (2007); Sanz et al. (2009). We note the projector $h$ has a block diagonal form, corresponding to different irreducible representations of $\mathcal{G}$. We call the two-particle subspace that $h$ projects onto a “bright manifold” $\mathfrak{B}\equiv\operatorname{range}(h)\subset\mathds{C}^{d^{2}}$. The ground states are uniquely characterized by vanishing populations in $\mathfrak{B}$ for every neighboring pair of particles.

The foremost necessary condition for the jump operators $\{c_{\mu}\}$ is

[TABLE]

In other words, the jump operators must be capable of depopulating the entire bright manifold after averaging over all symmetry operations. While here we have assumed that the symmetry group $\mathcal{G}$ is finite for simplicity, the following results apply to any compact group by replacing the sum over $g\in\mathcal{G}$ by an integral over the Haar measure of $\mathcal{G}$.

To find the minimum number $k_{\text{min}}$ of $\ket{\psi_{\mu}}$ (and consequently $c_{\mu}$) required, it is useful to decompose $V_{g}$ into direct sums of irreducible representations (irreps) $V_{g}=\bigoplus_{r}V_{g}^{r}$, where $r$ enumerates the irreps, each with dimension $d_{r}$. This decomposition is possible because finite-dimensional unitary representations of any group are completely reducible Artin (2010). Let $\ket{\psi_{\mu}}=\bigoplus_{r}\ket{\psi_{\mu}^{r}}$, where each $\ket{\psi_{\mu}^{r}}$ is a $d_{r}$-dimensional vector. Observe that for any $d_{r}\times d_{r^{\prime}}$ matrix $X$, we can derive the following identity using Schur’s lemma Artin (2010):

[TABLE]

where $r\cong r^{\prime}$ means $V_{g}^{r}=U_{rr^{\prime}}V_{g}^{r^{\prime}}U_{rr^{\prime}}^{\dagger}$, or $r$ is equivalent (isomorphic) to $r^{\prime}$ up to a unitary basis change. Note we can always choose a basis for the representation of $V_{g}$ that absorbs $U_{rr^{\prime}}$, so we assume $U_{rr^{\prime}}=\mathds{1}$ without loss of generality. Using the notation $\bigoplus_{r,r^{\prime}}M_{r,r^{\prime}}$ to denote the matrix whose $r$-th row, $r^{\prime}$-th column block is $M_{r,r^{\prime}}$, we can write $Q_{\mu}$ through the above identity as

[TABLE]

where the second term characterizes the possible nonzero off-diagonal blocks, which can only exist between pairs of equivalent irreps.

Since inequivalent irreps are decoupled, we for now only consider the subspace $\mathfrak{B}_{r}\subseteq\mathfrak{B}$ corresponding to $K_{r}$ copies of irreps equivalent to irrep $r$ ($\dim\mathfrak{B}_{r}=K_{r}d_{r}$). Let $Q_{\mu}^{r}=Q_{\mu}|_{\mathfrak{B}_{r}}$ as the operator $Q_{\mu}$ restricted to the subspace $\mathfrak{B}_{r}$. Observe that $\ket{\psi_{\mu}}$ restricted to this subspace is specified by the set of $K_{r}$ vectors $\{\ket{\psi_{\mu}^{s}}\in\mathds{C}^{d_{r}}\}_{s=1}^{K_{r}}$. When $K_{r}>d_{r}$, regardless of the choice of $\ket{\psi_{\mu}}$, there are $K_{r}-d_{r}$ linearly independent vectors $\vec{\beta}^{j}\in\mathds{C}^{K_{r}}$, $j=1,\ldots,K_{r}-d_{r}$, such that $\sum_{s=1}^{K_{r}}\beta^{j}_{s}\bra{\psi_{\mu}^{s}}=0$. Then any vectors of the form $\ket{\chi}=\bigoplus_{s=1}^{K_{r}}\beta_{s}^{j}\ket{v}$ are in the kernel of $Q_{\mu}^{r}$ for any $\ket{v}\in\mathds{C}^{d_{r}}$, since one can verify $Q_{\mu}^{r}\ket{\chi}=0$. Since there are $(K_{r}-d_{r})d_{r}$ linearly independent such vectors $\ket{\chi}$, we have $\operatorname{rank}(Q_{\mu}^{r})\leq d_{r}^{2}$. Hence, in order to fully depopulate $\mathfrak{B}_{r}$, we need $\text{dim}(\mathfrak{B}_{r})=\operatorname{rank}(\sum_{\mu}Q_{\mu}^{r})\leq\sum_{\mu}\operatorname{rank}(Q_{\mu}^{r})\leq k_{\text{min}}d_{r}^{2}$. Because $k_{\text{min}}$ must be an integer, we must have $k_{\text{min}}\geq\lceil K_{r}/d_{r}\rceil$, for every irrep $r$.

Note that this lower bound for $k_{\text{min}}$ can be saturated by construction as follows. First, we partition the $K_{r}$ equivalent irreps into $\lceil K_{r}/d_{r}\rceil$ groups of no more than $d_{r}$ irreps. For each group, we can assign a $\ket{\psi_{\mu}}$ that is nonzero only in the subspace corresponding to the irreps in the group. Lastly, we make all off-diagonal blocks vanish for each group $\mu\in\{1,\ldots,\lceil K_{r}/d_{r}\rceil\}$, by finding $d_{r}$ or fewer mutually orthogonal vectors $\ket{\psi_{\mu}^{r}}\in\mathds{C}^{d_{r}}$ such that $\braket{\psi_{\mu}^{r^{\prime}}}{\psi_{\mu}^{r}}=0$ for $r\neq r^{\prime}$. For a single jump operator of the form $c_{\mu}=|{\phi_{\mu}}\rangle\langle{\psi_{\mu}}|$, the state $\ket{\psi_{\mu}}$ may have supports on more than one of subspaces $\mathfrak{B}_{r}$. Therefore, the construction of a set of jump operators $\{c_{\mu}\}$ to satisfy Eq. (71) can be done in parallel for all the different $\mathfrak{B}_{r}$ corresponding to the inequivalent set of irreps, leading to the minimum number

[TABLE]

Here, $r$ enumerates inequivalent irreps of $\mathcal{G}$ in $\mathfrak{B}$, $K_{r}$ is the number of copies of $r$, and $d_{r}$ is the dimension of $r$.

We can also consider an arbitrary jump operator $c_{\mu}$ beyond the rank-1 form of $|{\phi_{\mu}}\rangle\langle{\psi_{\mu}}|$. For any operator $c_{\mu}$, we can perform singular value decomposition to write $c_{\mu}=\sum_{i_{\mu}}\sqrt{\gamma_{i_{\mu}}}|{\phi_{i_{\mu}}}\rangle\langle{\psi_{i_{\mu}}}|$, where $\braket{\phi_{i_{\mu}}}{\phi_{j_{\mu}}}=\braket{\psi_{i_{\mu}}}{\psi_{j_{\mu}}}=\delta_{i_{\mu}j_{\mu}}$. Then $c_{\mu}^{\dagger}c_{\mu}=\sum_{i_{\mu}}\gamma_{i_{\mu}}|{\psi_{i_{\mu}}}\rangle\langle{\psi_{i_{\mu}}}|$ with $\gamma_{i_{\mu}}>0$. Hence, the condition of Eq. (71) becomes a condition imposed on the set of right-singular vectors $\{\ket{\psi_{i_{\mu}}}:\forall\mu,i_{\mu}\}$, where we must have $\mathfrak{B}=\operatorname{range}(\frac{1}{|\mathcal{G}|}\sum_{g,\mu,i_{\mu}}\gamma_{i_{\mu}}V_{g}^{\dagger}|{\psi_{i_{\mu}}}\rangle\langle{\psi_{i_{\mu}}}|V_{g})$. We can thus interpret the $k_{\textrm{min}}$ found for rank-1 jump operators as the minimum number of independent $\ket{\psi_{i_{\mu}}}$’s.

While we can easily construct a minimal set of $\{c_{\mu}\}$ to satisfy the necessary condition of Eq. (71), we still need to prove the uniqueness of steady states. This can be done using our inductive proof strategy, where one simply needs to confirm that there are only trivial solutions to Eq. (44) under open boundary conditions. As discussed in Sec. V.1, this simply involves showing that a certain matrix $\boldsymbol{M}$ has full rank. Nevertheless, for non-injective MPSs, this scheme cannot break the ground-state degeneracy intrinsic to the MPS parent Hamiltonian, but it can guarantee that the ground states are the only steady states.