Robust control of long-distance entanglement in disordered spin chains

Many elementary tasks of quantum information processing can be performed on small scales with existing technology. For example, state tomography [1] on a single qubit is routinely done in many laboratories, and the number 15 has been factorized with a quantum device [2]. Realizing these tasks on a larger scale is one of the most pressing challenges today in research on engineered quantum systems: characterizing the state of many qubits is an actively pursued problem even on the theoretical side [3], and factorizing 77 or 187 is still impossible with our available technology. Similarly, entangled states of two qubits can be prepared with many systems [4, 5], but most setups are not scalable; the number of entangled photons is limited by the increasingly low probabilities of spontaneous events [6], and satisfactory scaling has so far been demonstrated for trapped ions only [7].

A central difference between trapped ions and other systems is that ions interact via long-range interactions, which facilitates the creation of strongly entangled states. Most other systems, however, are limited by rapidly decreasing interactions, and this disadvantage easily compensates for the added value of the long coherence times [8] that can be found e.g. in impurities of solid state lattices like nitrogen-vacancy (NV) centres. These unfavourable interaction properties can be improved if auxiliary quantum systems are available to mediate interactions [9–11], and the establishment of long-distance entanglement via chains of auxiliary spins seems to be a very promising route [12, 13].

Most approaches, however, rely on perfectly ordered chains [14–16], whereas any implantation of auxiliary spins is likely to result in a slightly disordered chain with non-uniform interactions. Our goal is to devise temporally shaped control fields that permit the establishment of long-distance entanglement independently of the specific realization of such disorder.

Compensation of disorder through suitably designed control fields is well established and has been demonstrated abundantly. In particular numerical pulse design [17, 18] has proven very successful. In our current study, however, numerical approaches suffer from the inherent growth of complexity of composite quantum systems with the number of constituents. We therefore resort to a description in terms of matrix product states (MPSs); that is, a state $| {\rm{\Psi }}\rangle$ of N spins is parametrized in terms of matrices ${A}_{{i}_{j}}$ [19, 20] via the prescription

$$\begin{eqnarray}&&| {\rm{\Psi }}\rangle =\displaystyle \sum _{{i}_{1}\cdots {i}_{N}}\mathrm{tr}[{A}_{{i}_{1}}\cdots {A}_{{i}_{N}}]| {i}_{1}\cdots {i}_{N}\rangle .\end{eqnarray} \tag{ 1 }$$

The dimensions of these matrices (bond dimension) limit the overall entanglement of states that can be described with this ansatz, but a bond dimension between 20 and 30 is enough for the present purposes. MPSs permit treating systems of several hundreds of spins [21, 22]. Based on the MPSs description and the underlying variational ansatz [19], advanced numerical algorithms for the simulation of large systems with weakly entangled states have been developed [23, 24], and such efficient descriptions provide an extremely promising starting point for the control of large quantum systems [25]. In general, MPSs can simulate the time evolution of many-body systems efficiently and accurately only for a short period of time due to the increase of entanglement between any two blocks of components [26, 27]. Despite their limitation to describe weakly entangled states, however, MPSs are a viable option for our purpose since we target the creation of strong entanglement among few distant spins. Whereas MPSs fail to describe strongly entangled states of N spins without losing their favourable scaling in N, they are perfectly capable of describing the state of an N spin system in which only a subset $M\ll N$ is strongly entangled. We will therefore strive for a control strategy that ensures that the N-spin system is weakly entangled during the entire time-window while M-body entanglement is being enhanced.

Despite the favourable scaling of the computational effort with N, simulating a system with $N\gg 1$ spins in terms of MPSs is a numerically expensive endeavour. Typical pulse shaping algorithms, however, rely on an iterative refinement that requires many repeated propagations [17, 18], which pushes the problem from hard to practically impossible. In order to be able to treat sufficiently large spin chains, we therefore resort to a variation of Lyapunov control [28] that permits us to identify a good pulse with a single propagation only. Normally, Lyapunov control is based on the identification of the control field that maximizes the increment of the selected goal at each instance of time. The present goal is the creation of entanglement, and since entanglement is independent of local spin orientations, we can always choose our target functional to be independent of single-spin dynamics.

Suppose initially the system is in a completely separable state and there is no direct interaction between the end spins; in this case the present time-local control scheme will not identify any control Hamiltonian that results in an increase of pairwise entanglement between site 1 and N (defined more rigorously later in equation (2)); only at a later stage when some entanglement has been built up, will a suitable control Hamiltonian be identified. We will therefore define a time-dependent target that is such that one can always find a suitable control Hamiltonian, and that will eventually coincide with the desired entanglement measure. If we consider a chain with nearest neighbour interaction only, then the only goal that is initially achievable is entanglement between neighbouring spins, say spin 1 and 2. Once this goal is achieved, one may strive for the creation of entanglement between spin 1 and 3. Such an entanglement swapping scheme can be realized with a sequence of N − 1 time-intervals; pairwise entanglement between site 1 and j is targeted in the $j-1\mathrm{st}$ interval, and the end of this interval is reached once a satisfactory value for the target has been reached.

Such target functionals should favour pairwise entanglement between two selected spins, and, additionally, penalize entanglement shared by any other spin in order to ensure validity of the description in terms of MPSs. The entanglement between spin i and the rest of the system can be characterized in terms of the purity $S({\varrho }_{i})=1-\mathrm{tr}{\varrho }_{i}^{2}$ of the reduced density matrix ${\varrho }_{i}.$ A target functional that is maximized if spins i and j form a Bell state can be chosen as $S({\varrho }_{i})+S({\varrho }_{j})-\mu S({\varrho }_{i,j}).$ The first two terms favour entanglement shared by spins i and j; the third term takes into account that mixed states tend to be weakly entangled or separable [29]. For $\mu =2,$ this is a lower bound [30] to the concurrence [31] of mixed states, which is particularly good for weakly mixed states [32]; for the use of the control target, however, we found that lower values of μ are favourable, and we will use $\mu =1/5$ later on. With an additional penalty for entanglement shared by spins different from i or j, we arrive at the target

$$\begin{eqnarray}&&{\tau }_{{ij}}=S({\varrho }_{i})+S({\varrho }_{j})-\mu S({\varrho }_{i,j})-\displaystyle \sum _{k\ne i,j}{\alpha }_{k}S({\varrho }_{k}),\end{eqnarray} \tag{ 2 }$$

where the non-negative scalars ${\alpha }_{k}$ permit choosing the emphasis on the penalty for spin k.

For the specific physical situations to be considered, we will limit ourselves to single-spin control Hamiltonians, as realistically available means of control like microwave or laser-fields induce such single-spin dynamics. In regular Lyapunov control, the control Hamiltonian is constructed based on the time-derivative of the target functional, but since the present functionals are invariant under single-spin dynamics, ${\dot{\tau }}_{{ij}}$ does not permit constructing an optimal control Hamiltonian. It is, however, possible to consider the curvature $\ddot{\tau }$ rather than the increase $\dot{\tau }$ of a target functional τ to read off an instantaneously optimal control Hamiltonian [28, 33]. The curvatures $\ddot{\tau }$ are defined in terms of the first two temporal derivatives

$$\begin{eqnarray}&&{\dot{\varrho }}_{j}=-i{\mathrm{tr}}_{\bar{j}}[H,| {\rm{\Psi }}\rangle \langle {\rm{\Psi }}| ]\ {\rm{and}}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{\ddot{\varrho }}_{j}=-i{\mathrm{tr}}_{\bar{j}}[\dot{H},| {\rm{\Psi }}\rangle \langle {\rm{\Psi }}| ]-{\mathrm{tr}}_{\bar{j}}[H,[H,| {\rm{\Psi }}\rangle \langle {\rm{\Psi }}| ]],\end{eqnarray} \tag{ 4 }$$

of the reduced density matrices ${\varrho }_{j},$ which, in turn depend on the state $| {\rm{\Psi }}\rangle$ of the full chain, and the chain's Hamiltonian $H={H}_{s}+{H}_{c}(t)$ [34] comprised of the static system Hamiltonian H_s and the to-be-designed time-dependent control Hamiltonian ${H}_{c}(t);$ '${\mathrm{tr}}_{\bar{j}}$' denotes the partial trace over all spins but spin j. Since the tunable control Hamiltonian does not contain any interaction terms, it can be written as

$$\begin{eqnarray}&&{H}_{c}(t)=\displaystyle \sum _{i=1}^{N}\displaystyle \sum _{\theta }\;{g}_{i}^{(\theta )}(t)\ {\sigma }_{i}^{(\theta )},\end{eqnarray} \tag{ 5 }$$

in terms of the usual Pauli matrices ${\sigma }_{i}^{(\theta )}$ for spin half systems where $\theta =\{x,y,z\},$ or different operators that correspond to implementable Hamiltonians.

Equations (3) and (4) suggest that ${\ddot{\tau }}_{{ij}}$ depends bi-linearly on ${H}_{c}(t),$ and that it depends linearly on ${\dot{H}}_{c}(t).$ Similarly to why ${\dot{\tau }}_{{ij}}$ does not depend of ${H}_{c}(t),$ however, one may see that ${\ddot{\tau }}_{{ij}}$ depends on H_c(t) only linearly, and that it does not depend on ${\dot{H}}_{c}(t)$ at all; bilinear terms in H_c(t) and linear terms in ${\dot{H}}_{c}(t)$ correspond to local unitary dynamics, i.e. dynamics that ${\tau }_{{ij}}$ is invariant under. Given the linear dependence in ${H}_{c}(t),$ one can thus express ${\ddot{\tau }}_{{ij}}$ as

$$\begin{eqnarray}&&{\ddot{\tau }}_{{ij}}=\displaystyle \sum _{p=1}^{N}\displaystyle \sum _{\theta }\displaystyle \frac{\partial {\ddot{\tau }}_{{ij}}({\rm{\Psi }})}{\partial {g}_{p}^{(\theta )}}{g}_{p}^{(\theta )}+{{\ddot{\tau }}_{{ij}}| }_{{g}_{p}^{\theta }=0\ \forall \theta }.\end{eqnarray} \tag{ 6 }$$

The maximium of ${\ddot{\tau }}_{{ij}}$ under the constraint that the magnitude ${\displaystyle \sum }_{k}{({g}_{i}^{(\theta )})}^{2}$ of a local control Hamiltonian is limited by some maximally admitted value is obtained for

$$\begin{eqnarray}&&{{g}_{i}^{(\theta )}| }_{\mathrm{opt}}={Z}_{i}\displaystyle \frac{\partial {\ddot{\tau }}_{1j}({\rm{\Psi }})}{\partial {g}_{i}^{(\theta )}},\end{eqnarray} \tag{ 7 }$$

with the normalization constant Z_i chosen such that the control does not exceed its admitted strength. This optimal choice can be constructed as analytic function of the system state $| {\rm{\Psi }}\rangle$ and the system Hamiltonian, so that at any instance during the propagation the optimal choice of H_c is available [28, 34]. For the actual propagation, it is practical not to work with time-dependent Hamiltonians, but rather choose the control parameter's constant during some short time interval Δ and update this choice after each multiple of this period [23].

Since τ is based on single-spin and two-spin reduces density matrices only, and the Hamiltonian contains only single-spin and two-spin terms, the optimal control Hamiltonian is characterized completely in terms of up to three-spin reduced density matrix, which can be constructed efficiently from the MPSs. In fact, for μ = 0, ${\tau }_{{ij}},$ defined in equation (2), is defined in terms of single-spin reduced density matrices only, so that the optimal control Hamiltonian can be constructed exclusively in terms of two-spin reduced density matrices.

To be specific, let us consider the example of a spin chain with an Ising interaction

$$\begin{eqnarray}&&{H}_{s}(\vec{J})=\displaystyle \sum _{i=1}^{N-1}{J}_{i}{\sigma }_{i}^{z}{\sigma }_{i+1}^{z},\end{eqnarray} \tag{ 8 }$$

that is characterized by a vector $\vec{J}$ that contains the coupling constants for the N − 1 nearest neighbour interactions. For this specific Hamiltonian, the ${\sigma }_{z}$ components of the optimal control Hamiltonian vanish (since they commute with H_s), and the ${\sigma }_{x}$ and ${\sigma }_{y}$ components are obtained from

$$\begin{eqnarray*}\displaystyle \frac{\partial \ddot{S}({\rho }_{i})}{\partial {g}_{i}^{\theta }} & = & {J}_{i-1}({\rho }_{i}^{\theta }{\rho }_{i-1,i}^{z,z}-{\rho }_{i}^{z}{\rho }_{i-1,i}^{z,\theta })+{J}_{i}({\rho }_{i}^{\theta }{\rho }_{i,i+1}^{z,z}-{\rho }_{i}^{z}{\rho }_{i,i+1}^{\theta ,z})\\ \displaystyle \frac{\partial \ddot{S}({\rho }_{i})}{\partial {g}_{i+1}^{\theta }} & = & {J}_{i}({\rho }_{i}^{\phi }{\rho }_{i,i+1}^{\theta ,\phi }-{\rho }_{i}^{\theta }{\rho }_{i,i+1}^{\phi ,\phi }),\\ \displaystyle \frac{\partial \ddot{S}({\rho }_{i})}{\partial {g}_{i-1}^{\theta }} & = & {J}_{i-1}({\rho }_{i}^{\phi }{\rho }_{i-1,i}^{\phi ,\theta }-{\rho }_{i}^{\theta }{\rho }_{i-1,i}^{\phi ,\phi }),\\ \displaystyle \frac{\partial \ddot{S}({\rho }_{12})}{\partial {g}_{1}^{\theta }} & = & 0,\\ \displaystyle \frac{\partial \ddot{S}({\rho }_{12})}{\partial {g}_{2}^{\theta }} & = & {J}_{2}\displaystyle \sum _{{\rm{\Delta }}=0}^{3}({\rho }_{1,2}^{{\rm{\Delta }},\theta }{\rho }_{1,2,3}^{{\rm{\Delta }},z,z}-{\rho }_{1,2}^{{\rm{\Delta }},z}{\rho }_{1,2,3}^{{\rm{\Delta }},\theta ,z}),\\ \displaystyle \frac{\partial \ddot{S}({\rho }_{12})}{\partial {g}_{3}^{\theta }} & = & {J}_{2}\displaystyle \sum _{{\rm{\Delta }}=0}^{3}({\rho }_{1,2}^{{\rm{\Delta }},\phi }{\rho }_{1,2,3}^{{\rm{\Delta }},\theta ,\phi }-{\rho }_{1,2}^{{\rm{\Delta }},\theta }{\rho }_{1,2,3}^{{\rm{\Delta }},\phi ,\phi }),\end{eqnarray*}$$

and (for $j\gt 2$) from

$$\begin{eqnarray*}\displaystyle \frac{\partial \ddot{S}({\rho }_{1j})}{\partial {g}_{1}^{\theta }} & = & {J}_{1}\displaystyle \sum _{{\rm{\Delta }}=0}^{3}({\rho }_{1,j}^{\theta ,{\rm{\Delta }}}{\rho }_{1,2,j}^{z,z,{\rm{\Delta }}}-{\rho }_{1,j}^{z,\delta }{\rho }_{1,2,j}^{\theta ,z,{\rm{\Delta }}}),\\ \displaystyle \frac{\partial \ddot{S}({\rho }_{1j})}{\partial {g}_{2}^{\theta }} & = & {J}_{1}\displaystyle \sum _{{\rm{\Delta }}=0}^{3}({\rho }_{1,j}^{\phi ,{\rm{\Delta }}}{\rho }_{1,2,j}^{\theta ,\phi ,{\rm{\Delta }}}-{\rho }_{1,j}^{\theta ,{\rm{\Delta }}}{\rho }_{1,2,j}^{\phi ,\phi ,{\rm{\Delta }}}),\\ \displaystyle \frac{\partial \ddot{S}({\rho }_{1j})}{\partial {g}_{j+1}^{\theta }} & = & {J}_{j}\displaystyle \sum _{{\rm{\Delta }}=0}^{3}({\rho }_{1,j}^{{\rm{\Delta }},\phi }{\rho }_{1,j,j+1}^{{\rm{\Delta }},\theta ,\phi }-{\rho }_{1,j}^{{\rm{\Delta }},\theta }{\rho }_{1,j,j+1}^{{\rm{\Delta }},\phi ,\phi }),\\ \displaystyle \frac{\partial \ddot{S}({\rho }_{1j})}{\partial {g}_{j-1}^{\theta }} & = & {J}_{j-1}\displaystyle \sum _{{\rm{\Delta }}=0}^{3}({\rho }_{1,j}^{{\rm{\Delta }},\phi }{\rho }_{1,j-1,j}^{{\rm{\Delta }},\phi ,\theta }-{\rho }_{1,j}^{{\rm{\Delta }},\theta }{\rho }_{1,j-1,j}^{{\rm{\Delta }},\phi ,\phi }),\\ \displaystyle \frac{\partial \ddot{S}({\rho }_{1j})}{\partial {g}_{j}^{\theta }} & = & \displaystyle \sum _{{\rm{\Delta }}=0}^{3}{J}_{j-1}({\rho }_{1,j}^{{\rm{\Delta }},\theta }{\rho }_{1,j-1,j}^{{\rm{\Delta }},z,z}-{\rho }_{1,j}^{{\rm{\Delta }},z}{\rho }_{1,j-1,j}^{{\rm{\Delta }},z,\theta })+{J}_{j}({\rho }_{1,j}^{{\rm{\Delta }},\theta }{\rho }_{1,j,j+1}^{{\rm{\Delta }},z,z}-{\rho }_{1,j}^{{\rm{\Delta }},z}{\rho }_{1,j,j+1}^{{\rm{\Delta }},\theta ,z}),\end{eqnarray*}$$

where $\{\theta ,\phi \}=\{x,y\},$ and $\theta \ne \phi .$ ${\rho }_{i}^{{\theta }_{i}}=\mathrm{tr}({\sigma }_{i}^{{\theta }_{i}}\rho ),$ ${\rho }_{i,j}^{{\theta }_{i},{\theta }_{j}}=\mathrm{tr}({\sigma }_{i}^{{\theta }_{i}}{\sigma }_{j}^{{\theta }_{j}}\rho )$ and ${\rho }_{i,j,k}^{{\theta }_{i},{\theta }_{j},{\theta }_{k}}=\mathrm{tr}({\sigma }_{i}^{{\theta }_{i}}{\sigma }_{j}^{{\theta }_{j}}{\sigma }_{k}^{{\theta }_{k}}\rho )$ are the expectation values of the spin operators ${\sigma }_{{\theta }_{i}},$ ${\sigma }_{{\theta }_{i}}\otimes {\sigma }_{{\theta }_{j}}$ and ${\sigma }_{{\theta }_{i}}\otimes {\sigma }_{{\theta }_{j}}\otimes {\sigma }_{{\theta }_{k}}$ for spins $i,j$ and k. In the case ofspins at the end of the chain, i.e. i = 1, i = N or j = N, it is implied that ${J}_{0}={J}_{N}=0$ since there is no corresponding interaction partner.

4.1. Ideal Ising chain

Before discussing disordered chains, let us first demonstrate the functionality in terms of an ideal, ordered chain. Figure 1 depicts the sequential increase and decrease of the different control targets ${\tau }_{1j}$ with μ = 0 and ${\alpha }_{k}=1$ for a chain of N = 10 spins with uniform interactions, i.e. ${J}_{i}=J$ for all i. The control Hamiltonians are limited by $\beta =\sqrt{{({g}_{i}^{x})}^{2}+{({g}_{i}^{y})}^{2}}=70J$ and they remain constant over periods of ${\rm{\Delta }}=1/(1000J)$ ³ . The chain is initialized in a separable state $| +{\rangle }^{\otimes N},$ where $| +\rangle$ is the eigenstate of ${\sigma }^{x}$ with eigenvalue $+1,$ and the target ${\tau }_{1j}$ is replaced by ${\tau }_{1,j+1}$ if ${\tau }_{1j}$ saturates.

Zoom In Zoom Out Reset image size

The value of ${\tau }_{1j}$ is reduced if spins other than 1 and j participate in any entanglement. This is merely due to the necessity to keep many-body entanglement sufficiently small for an efficient simulation, but the original goal of creating long-distant entanglement is not jeopardized by an unintentional creation of additional entanglement. One should therefore characterize the performance by the pairwise entanglement, i.e. the entanglement of the reduced density matrix ϱ_1j of spins 1 and j, rather than ${\tau }_{1j}.$ We characterize the pairwise entanglement via Wootters' convex roof construction of concurrence ${c}_{1j}=\sqrt{{\mu }_{1}}-{\displaystyle \sum }_{i=2}^{4}\sqrt{{\mu }_{i}}$ with the decreasingly ordered eigenvalues ${\mu }_{i}$ of $({\sigma }_{y}\otimes {\sigma }_{y}){\varrho }_{1j}^{*}({\sigma }_{y}\otimes {\sigma }_{y}){\varrho }_{1j}$ [31].

Indeed, one observes a sequential growth and decline of the different selections of pairwise entanglement c_1j similar to the behaviour of ${\tau }_{1j}$ depicted in figure 1. There is an essentially negligible decrease of the peak height as j increases, and substantial entanglement of ${c}_{1\ 10}\simeq 0.999$ is established between the two spins at the end of the chain. Tests with longer chains gave ${c}_{1N}\simeq 0.997$ for N = 20 and ${c}_{1N}\simeq 0.994$ for N = 40. There is thus only a negligible decay of the achievable entanglement with increasing systems size.

Figure 2 depicts as an example the time-dependent control (equation (7)) of spins 5–6 that achieve the swapping from ${\tau }_{15}$ to ${\tau }_{16}$ in an ideal Ising chain. The construction in terms of finite time steps tends to result in unnecessary high-frequency components; we explicitly verified that such components can be dropped without sizeable reduction of performance and figure 2 depicts such a 'smoothened' pulse⁴ . As one can see, the control is mostly limited to short time windows, with extended time windows of free dynamics in-between; that is, the control tends to align the spins such that their subsequent dynamics exploits the intrinsic interaction in an optimal fashion.

Zoom In Zoom Out Reset image size

4.2. Disordered chain

The performance of the present control methods is by no means specific to uniform chains. Repeating this procedure for disordered chains with J_i = r_iJ where r_i are random numbers drawn from a uniform distribution $[0.9,1.1]$ resulted in similar values of ${c}_{1N},$ and in the following we can address the central question of whether this method permits identification of a control pulse that works independently of the specific realization of the coupling constants J_i.

Since there is typically no unique optimal (or close to optimal) pulse for a given set of coupling constants, it is possible to find a pulse that performs well for different interaction landscapes. We can therefore consider an ensemble of chains with different realizations of coupling constants, and construct the optimal control Hamiltonians via the ensemble average. If the utilized ensemble is sufficiently large, one can expect the resulting pulse to perform irrespective of the details of the actual system. Figure 3 shows the entanglement dynamics around the peaks of c_1j resulting from a control sequence constructed with an ensemble average over N = 50 different realizations. Only a minor tribute is paid to the disorder, as the maximally reached value of c_1N is 0.956; that is, there is a loss of about $4\%.$ To address the question of whether this pulse is applicable to this specific ensemble only, or, if it will perform equally well for any other random realization of coupling constants, we can apply the pulse to a test ensemble of another 50 randomly chosen spin chains. As the dots in figure 3 show, the behaviour on this test ensemble is hardly different from that of the original ensemble, which substantiates that the control induces dynamics that is largely insensitive to variations in the interaction landscape. Despite the fact that our method does not involve iterative refinement of the pulse, it yields good results even for more strongly disordered chains as we explicitly verified for an ensemble with coupling constants drawn from the interval $[0.8,1.2];$ even though the maximal amplitude of the static noise is comparable in size to the typical interaction constant, one obtains substantial entanglement ${c}_{1N}\simeq 0.865$ between the two end spins for N = 10.

Zoom In Zoom Out Reset image size

4.3. Reduced control

So far, we have considered the control of all spins. For most of the spins, however, control was introduced only to ensure that MPSs are a good description. We may, therefore also relax the control and choose ${\alpha }_{k}=0$ for many spins.

Control of the first spin is essential; and during the $j-1\mathrm{st}$ interval, in which the entanglement shared with spin 1 is swapped between the $j-1\mathrm{st}$ and the jth spin, control on spin j − 1 and j is also essential. It seems plausible that control of spins that are far away from any of these three spins is less important than control of spins that are close to one of the essential spins. We have investigated the performance of control on a reduced number of spins, and found that control of many spins can be forfeited. Quite surprisingly, control of spin 2 is not necessary for good performance after the third interval. Control of $j-2\mathrm{nd}$ and $j+1\mathrm{st},$ however, is important during the swapping procedure from spin j − 1 to spin j. Spins that had participated in a swapping operation will no longer be relevant; it is therefore not surprising that control on spins $\lt j-2$ can readily be given up. However, spins that will participate in a swapping operation in the near future need to be controlled and we found that control on spin $j+2$ is required for good performance. One may thus reduce the control of spin 1 and spins j − 2 through to $j+2$ with j increasing by 1 after each swapping operation. With control on these six spins only, one obtains very good performance as depicted in figure 4. The loss of entanglement during the swapping operations is essentially negligible and substantial entanglement ${c}_{1\ 40}\simeq 0.985$ is established over a chain of 40 spins.

Zoom In Zoom Out Reset image size

Given the control on a reduced number of spins, a thorough check of the accuracy of the numerical propagation is in order. We have therefore simulated the dynamics with MPSs of different bond dimension. The entanglement that builds up during the dynamics suggests a necessary minimal bond dimension of 8. We worked with a bond dimension of 10, and explicitly confirmed that an increase to a bond dimension of 20 does not result in any discernible change in dynamics, which confirms the validity of the MPSs' description with low bond dimension.

It is interesting to notice that our approach is quite different from the notion of 'perfect state transfer' (PST) or 'almost perfect state transfer' (APST), where two parties employ a spin chain with perfect coupling strength as quantum wire. The protocol is initialized with the preparation of the first spin in the to-be-transferred state. After some period of evolution induced by the system Hamiltonian, the final spin will have this state with non-zero fidelity [35–37]. This protocol certainly also permits the creation of some distant entanglement, but our control scheme has the following advantages: i) It is robust against disorder in the coupling of the spins and against the lengths of the spin chain, whereas (A)PST greatly depends on the coupling and the length of the spin wire [36, 37]. ii) The amount of entanglement that can be achieved with the present control scheme is substantially higher than that of (A)PST, and iii) the time cost to achieve such high entanglement is much lower. For example, the maximally achieved entanglement between site 1 and site 10 within time cost $4000/J$ in a perfect spin chain by means of (A)PST is 0.95 [35] as compared to 0.999 for a perfect spin chain and 0.958 for a disordered spin chain that can be created within $12.76/J$ with control. In particular with increasing system size, the advantage of the present method becomes apparent: entanglement established over 80 sites within time cost $4000/J$ in a perfect spin chain with (A)PST not exceeding 0.5 [35], whereas using our control strategy permits us to create entanglement between site 1 and site 80 amounting to 0.990 with the time cost roughly equal to $130/J.$

The creation of long-distance entanglement is also by no means limited to bipartite entanglement, but may also be employed for the creation of entanglement between three or more distant spins. The demonstration of the usability of MPSs for the control of large spin system, in particular, also suggests that other goals like the implementation of multi-qubit quantum gates involving distant spins can be realized in a similar fashion. In all such situations, the applicability of MPSs for a given situation can be ensured through a suitably extended target functional that makes sure that many-body entanglement remains sufficiently low. As simple control strategies like Lyapunov control might fail to identify goals whose realization requires many elementary interactions, it helps to define a sequence of intermediate goals. In the present case we did so by sudden changes of the target functional, but also smooth, continuous modulations of targets (which itself might become the object of optimization) is conceivable. Such well-designed dynamical goals together with advanced numerical techniques like MPSs promise to help us make the step from small scale proof-of-principle demonstrations towards large scales.

The numerical work was strongly supported by Mari Carmen Bañuls. JC acknowledges the hospitality by the Institute for Interdisciplinary Information Sciences Tsinghua University and the Institute of Physics, Chinese Academy of Sciences. Computational resources in Rechenzentrum Garching and financial support the European Research Council within the project 'Optimal dynamical control of quantum entanglement' is gratefully acknowledged.