Tensor product approach to quantum control

TL;DR

This paper demonstrates that tensor product algorithms can efficiently control large quantum systems, enabling the simulation and optimization of quantum states with reduced computational resources, which is crucial for advancing quantum computing.

Contribution

It introduces a tensor product approach combined with the tAMEn algorithm for optimal quantum control, overcoming dimensionality challenges in large quantum systems.

Findings

Successfully controlled a 41-spin quantum system on a single workstation.

Achieved significant savings in computational time and memory.

Enabled potential development of quantum computers with 50-100 qubits.

Abstract

In this proof-of-concept paper we show that tensor product approach is efficient for control of large quantum systems, such as Heisenberg spin wires, which are essential for emerging quantum computing technologies. We compute optimal control sequences using GRAPE method, applying the recently developed tAMEn algorithm to calculate evolution of quantum states represented in the tensor train format to reduce storage. Using tensor product algorithms we can overcome the curse of dimensionality and compute the optimal control pulse for a 41 spin system on a single workstation with fully controlled accuracy and huge savings of computational time and memory. The use of tensor product algorithms opens new approaches for development of quantum computers with 50 to 100 qubits.

Equations

$$\dot{\rho }=-i\text{mathaccent}28766\text{mathaccent}28766H(t)\rho +\text{mathaccent}28766\text{mathaccent}28766R\rho ,$$

$$\text{mathaccent}28766\text{mathaccent}28766H(t)=I\otimes \text{mathaccent}28766H(t)-\text{mathaccent}28766H(t)\otimes I.$$

$$\text{mathaccent}28766H(t)=\text{mathaccent}28766H_0+\sum _{k=1}^K{c}_k(t)\text{mathaccent}28766H_k.$$

$$\dot{\rho }=-i\left(\text{mathaccent}28766{\text{mathaccent}28766H}_0+\sum _{k=1}^K\sum _{n=1}^N{c}_{n,k}\text{mathaccent}28766{\text{mathaccent}28766H}_k{\phi }_n(t)\right)\rho ,\rho (0)={\rho }_0.$$

$$F(c)=\Re ⟨{\rho }_T⟩\rho (T).$$

$$\frac{\partial F}{\partial c_{n,k}}=\Re ⟨{\rho }_T⟩\frac{\partial \rho (T)}{\partial c_{n,k}}.$$

$$\frac{\partial }{\partial c_{n,k}}\frac{\partial \rho }{\partial t}=-i\left({\phi }_n(t)\text{mathaccent}28766{\text{mathaccent}28766H}_k\rho +\text{mathaccent}28766\text{mathaccent}28766H(t)\frac{\partial \rho }{\partial c_{n,k}}\right).$$

$$\frac{\partial }{\partial t}\left[\begin{array}{c}\frac{\partial \rho }{\partial c_{n,k}}\\ \rho \end{array}\right]=-i\left[\begin{array}{cc}\text{mathaccent}28766\text{mathaccent}28766H(t)& {\phi }_n(t)\text{mathaccent}28766{\text{mathaccent}28766H}_k\\ 0& \text{mathaccent}28766\text{mathaccent}28766H(t)\end{array}\right]\left[\begin{array}{c}\frac{\partial \rho }{\partial c_{n,k}}\\ \rho \end{array}\right],{\left[\begin{array}{c}\frac{\partial \rho }{\partial c_{n,k}}\\ \rho \end{array}\right]}_{t=0}=\left[\begin{array}{c}0\\ {\rho }_0\end{array}\right].$$

$$\rho (\delta t)={\rho }_0-i\delta t\left(\text{mathaccent}28766{\text{mathaccent}28766H}_0+\sum _{k=1}^K\sum _{n=1}^N{c}_{n,k}\text{mathaccent}28766{\text{mathaccent}28766H}_k{\phi }_n(t)\right){\rho }_0+\mathcal{O}(\delta t^2)$$

$${\theta }_n(t)=\{\begin{array}{ll}1,& \text{ifleavevmodenobreak tn−1<t⩽tn},\\ 0,& \text{otherwise},\end{array}$$

$$\frac{\partial }{\partial t}\frac{\partial \rho }{\partial c_{n,k}}=-i\text{mathaccent}28766\text{mathaccent}28766H\frac{\partial \rho }{\partial c_{n,k}}\text{for}t\in (0,t_{n-1}],\text{s.t.}\frac{\partial \rho }{\partial c_{n,k}}(0)=0,$$

$$\frac{\partial \rho }{\partial c_{n,k}}(T)=\exp \left[-i{\int }_{t_n}^T\text{mathaccent}28766\text{mathaccent}28766H(s)\mathrm{d}s\right]\frac{\partial \rho }{\partial c_{n,k}}(t_n)=\prod _{m=n+1}^N\underset{Q_m}{\underbrace{\exp \left[-i\tau \text{mathaccent}28766\text{mathaccent}28766H(t_m)\right]}}\frac{\partial \rho }{\partial c_{n,k}}(t_n),$$

$$\frac{\partial F}{\partial c_{n,k}}=\Re ⟨{\rho }_T⟩\prod _{m=n+1}^N{Q}_m\frac{\partial \rho }{\partial c_{n,k}}(t_n)=\Re ⟨\prod _{m=n+1}^N{Q}_m^{†}{\rho }_T⟩\frac{\partial \rho }{\partial c_{n,k}}(t_n).$$

$$\lambda (t_n)=\prod _{m=n+1}^N{Q}_m^{†}{\rho }_T$$

$$\frac{\partial \lambda }{\partial t}=-i\text{mathaccent}28766\text{mathaccent}28766H(t)\lambda \text{for}t\in [t_n,T),\text{s.t.}\lambda (T)={\rho }_T.$$

$$\rho =[\rho (i_1,\dots ,i_d;j_1,\dots ,j_d)].$$

$$\rho (i_1,\dots ,i_d;\text{leavevmode}{j}_1,\dots ,j_d)\approx \sum _{{\alpha }_0,\dots ,{\alpha }_d=1}^{r_0,\dots ,r_d}{\rho }_{{\alpha }_0,{\alpha }_1}^{(1)}(i_1,j_1)\cdots {\rho }_{{\alpha }_{d-1},{\alpha }_d}^{(d)}(i_d,j_d).$$

$$\rho (t)\approx \sum _{m=1}^M\rho ({\tau }_m)L_m(t),t\in (0,\tau ],0<{\tau }_1<{\tau }_2<\cdots <{\tau }_M=\tau .$$

$$(\underset{A}{\underbrace{S\otimes I+I\otimes i\text{mathaccent}28766\text{mathaccent}28766H}})\overline{\rho }=\underset{f}{\underbrace{(S{1}_M)\otimes \rho (0)}}.$$

$$\begin{array}{rl}\overline{\rho }(i_1,\dots ,i_d;j_1,\dots ,j_d,m)& =\sum _{\alpha }{\rho }_{{\alpha }_0,{\alpha }_1}^{(1)}(i_1,j_1)\cdots {\rho }_{{\alpha }_{d-1},{\alpha }_d}^{(d)}(i_d,j_d){\rho }_{{\alpha }_d,{\alpha }_{d+1}}^{(d+1)}(m),\\ \overline{\rho }& =\sum _{\alpha }{\rho }_{{\alpha }_0,{\alpha }_1}^{(1)}\otimes \cdots \otimes {\rho }_{{\alpha }_{d-1},{\alpha }_d}^{(d)}\otimes {\rho }_{{\alpha }_d,{\alpha }_{d+1}}^{(d+1)}.\end{array}$$

$$P_{\ne k}=\left(\sum _{{\alpha }_0\dots {\alpha }_{k-2}}{\rho }_{{\alpha }_0,{\alpha }_1}^{(1)}\otimes \cdots \otimes {\rho }_{{\alpha }_{k-2}}^{(k-1)}\right)\otimes I\otimes \left(\sum _{{\alpha }_{k+1}\dots {\alpha }_{d+1}}{\rho }_{{\alpha }_{k+1}}^{(k+1)}\otimes \cdots \otimes {\rho }_{{\alpha }_d,{\alpha }_{d+1}}^{(d+1)}\right),$$

$$\text{mathaccent}28766H=\sum _{k=1}^{d-1}{J}_x\text{mathaccent}28766{\sigma }_x^{(k)}\text{mathaccent}28766{\sigma }_x^{(k+1)}+J_y\text{mathaccent}28766{\sigma }_y^{(k)}\text{mathaccent}28766{\sigma }_y^{(k+1)}+J_z\text{mathaccent}28766{\sigma }_z^{(k)}\text{mathaccent}28766{\sigma }_z^{(k+1)},$$

$$\text{mathaccent}28766{\sigma }_{\{x,y,z\}}^{(k)}=I\otimes \cdots \otimes I\otimes \underset{\text{k–th place}}{\underbrace{\text{mathaccent}28766{\sigma }_{\{x,y,z\}}}}\otimes I\otimes \dots \otimes I.$$

$${\psi }_0=|\uparrow \downarrow \downarrow \cdots \downarrow \downarrow ⟩,{\psi }_T=|\downarrow \downarrow \cdots \downarrow \downarrow \uparrow ⟩,$$

Full text

1 Introduction

Quantum control plays a central part in both established technologies such as nuclear magnetic resonance (NMR) and emerging technologies such as quantum computing [16]. The ‘hardware’ of a quantum computer is a quantum system (e.g. a spin system), which we manipulate with time–dependent fields (e.g. laser pulses). Our goal is to design a pulse sequence that steers the system from an initial state $\ket{\psi_{0}}$ to a desired state $\ket{\psi_{T}}$ during time $T$ — a ‘firmware’, so to say. The main challenge for numerical methods is the notoriously known curse of dimensionality — the dimension of the state space growing exponentially with the number of particles $d.$ Even for simplest spin–$\frac{1}{2}$ systems, the wavefunction $\ket{\psi},$ the density matrix $\rho,$ and the Hamiltonian $\mathaccent 28766{H}$ require at least $\mathcal{O}(2^{d})$ storage. To find the optimal trajectory we need to compute $\rho(t)$ in $100$s to $1000$s points of time. This leads to gigabyte–scale storage for $d\simeq 20$ spins and becomes unfeasible for $d\gtrsim 30.$

The main approaches used currently are the sparse format for $\mathaccent 28766{H}$, $\ket{\psi}$’s, $\rho$’s (only large entries are stored, zeroes and small elements are discarded), and dimensionality reduction — the equations of system dynamics are projected into a subspace of low–lying eigenstates [44] or its heuristically chosen equivalent [27, 4].

Tensor product algorithms are a novel approach to high–dimensional problems. Based on the idea of separation of variables, they approximate high–dimensional data with a tensor product of low–dimensional factors. For example, when interaction of particles is a superposition of a few near–neighbour interactions, the corresponding quantum state is a superposition of a few unentangled states, i.e it has a low Schmidt rank. The corresponding tensor format also has a low rank, providing good compression and low storage costs. By maintaining the compressed format for all operations we can compute the system’s trajectory and optimise it without any uncontrollable truncations or heuristic assumptions. Model examples of such systems include Ising and Heisenberg spin chains [35], more realistic applications include quantum wires [3] and backbones of simple protein molecules [34].

In quantum physics a tensor product format was first proposed in 1970s as the renormalization group formalism [43], followed in 1990s by MPS/MPO [14, 26] and DMRG [42] algorithms for finding the ground state of a quantum spin chain. Similar formats were applied in statistics and can be traced back to works of Hitchcock in 1930s [22, 21] and Tucker in 1960s [40]. Tensor formats in three dimensions were rigorously studied in the numerical linear algebra community as a generalisation of the low–rank decomposition of matrices [30, 15, 18, 33]. Generalisations to higher dimensions eventually led to re-discovery of the MPS/DMRG framework under the name of the tensor train (TT) format [31, 29]. In the last decade the TT format and the more general Hierarchical Tucker format [20] have been successfully applied to a variety of problems, such as finding several low–lying states of a quantum system [8], performing superfast Fourier transform [9], solving high–dimensional linear systems [10] and first–order differential equations [7]. First applications of tensor product methods to optimal control problems are considered in [6, 1].

The recently proposed tAMEn algorithm [7] allows to compute evolution of a quantum system under a given (time–dependent) Hamiltonian with high accuracy. Although the trajectory lies in extremely high–dimensional state space, tAMEn keeps the computation costs feasible by performing all calculations in the TT format and adapting the TT ranks to meet the required accuracy. In this paper we employ the tAMEn algorithm as a building block to develop a version of a classical GRAPE algorithm [25] and use it to control a spin chain with $d=41$ spins. This shows that tensor product algorithms can be used to design control sequences for quantum computers with $50$ to $100$ qubits, which can be expected within a decade [19].

2 Optimal quantum control

2.1 Dynamic optimisation problem

The dynamics of a quantum system with Hamiltonian $\mathaccent 28766{H}$ is described by Liouville–von Neumann equation in superoperator form

[TABLE]

where $\rho$ is the quantum density matrix stretched into a vector form (referred to later as quantum state), and

[TABLE]

Traditionally $\mathaccent 28766{H}(t)$ is split into a constant drift and a time–dependent control term:

[TABLE]

Control operators $\mathaccent 28766{H}_{k}$ are usually defined by model or instrument used; we manipulate the system by changing control pulse functions $c_{k}(t).$ Now choosing a suitable basis set $\{\varphi_{n}(t)\}_{n=1}^{N}$ and dropping the relaxation part $\mathaccent 28766{\mathaccent 28766{R}}$ for simplicity, we obtain an ODE for system dynamics:

[TABLE]

Our goal is to choose control parameters $c=[c_{n,k}]$ to maximise fidelity (or overlap)

[TABLE]

2.2 First–order optimisation framework

Classical optimisation methods require a gradient of $F$ w.r.t. $c_{n,k}$ which reads

[TABLE]

To find the gradient, we differentiate (1) w.r.t. control parameters and obtain

[TABLE]

Changing the order of derivatives in the left hand side, we obtain a system of coupled ODEs for the density matrix and its gradient:

[TABLE]

To justify the initial condition for the gradient, write a first–order approximation for $\rho(\delta t)$ with an infinitesimal time $\delta t$:

[TABLE]

and differentiate it w.r.t. $c_{n,k}$ to see $\frac{\partial\rho}{\partial c_{n,k}}(\delta t)=\mathcal{O}(\delta t)\to 0.$

The system (4) can be integrated on $[0,T]$ using a suitable method (e.g. Runge–Kutta scheme [32]) for all $k,n$ to produce both the fidelity (2) and its gradient (3).

2.3 GRAPE algorithm

A classical gradient ascend pulse engineering (GRAPE) method [25] is a simple version of the approach described above. Control pulses $c_{k}(t)$ are assumed piecewise–constant on intervals $(t_{n-1},t_{n}]$ between the nodes of the grid $\{t_{n}\}_{n=0}^{N}$ with $t_{0}=0$ and $t_{N}=T.$ We will use a uniform grid $t_{k}=\tau k$ with step–size $\tau=T/N.$ The basis functions $\varphi_{n}(t)$ can be taken as Heaviside functions

[TABLE]

thus making $c_{n,k}$ the amplitude of control pulse $c_{k}(t)$ on $(t_{n-1},t_{n}].$ The time–locality of basis functions $\theta_{n}(t)$ allows us to compute (3) more efficiently. Note that $\theta_{n}(t)=0$ for $t\leqslant t_{n-1}$ hence ODE (4) for the gradient component simplifies to the homogeneous ODE

[TABLE]

which of course gives $\frac{\partial\rho}{\partial c_{n,k}}(t_{n-1})=0.$ Therefore, to compute (3) we can propagate the state $\rho(t)$ on $(0,t_{n-1}]$ and start solving the coupled system (4) only from $t_{n-1}.$ When we reach $t_{n}$ and compute $\frac{\partial\rho}{\partial c_{n,k}}(t_{n}),$ we can note that $\theta_{n}(t)=0$ for $t>t_{n}$ hence (5) holds on $t\in(t_{n},T]$ as well. Since $\mathaccent 28766{\mathaccent 28766{H}}(t)$ is piecewise–constant,

[TABLE]

and plugging this in (3) gives

[TABLE]

The state

[TABLE]

is the solution of ODE with the boundary condition at the right end, colloquially referred to as ‘back–propagation’:

[TABLE]

Hence the coupled system (4) needs to be solved only on $[t_{n-1},t_{n}]$ to produce the $n$-th components of the gradient. By pre-computing all $\lambda(t_{n})$ for $n=1,\ldots,N,$ we can implement the GRAPE iteration with complexity scaling linearly in $N.$

For $n=N,\ldots,1,$ propagate (7) backward: $\lambda(t_{n})\>\longrightarrow\>\lambda(t_{n-1})$ 2. 2.

For $n=1,\ldots,N,$ propagate (4) forward: $\begin{bmatrix}0\\ \rho(t_{n-1})\end{bmatrix}\>\longrightarrow\>\begin{bmatrix}\frac{\partial\rho}{\partial c_{n,k}}(t_{n})\\ \rho(t_{n})\end{bmatrix}$ 3. 3.

Compute (6) and update the control amplitudes: $c_{n,k}\>\longrightarrow\>c_{n,k}-\epsilon\frac{\partial F}{\partial c_{n,k}}.$

Here $\epsilon$ is the step size for the gradient ascend. It should be carefully selected, since a too small value will make the algorithm converging very slowly while a very big one will make the algorithm divergent. In our implementation the step size is chosen adaptively: it is increased each time the step is successful (i.e. results in better fidelity) and decreased each time we have an overstep (i.e. fidelity decreases and the step has to be rejected). Note that line–search strategies for this optimisation problem are rather pointless, because the gradient (6) and the fidelity function (2) are computed simultaneously from (4), and have therefore the same complexity.

In the classic GRAPE algorithm, the propagators $Q_{m}$ in (6) are calculated as matrix exponentials [28] and applied to the states; in more advanced versions of the algorithm the action of matrix exponentials on the states is computed on–fly maintaining the sparse format for $\lambda$’s and $\rho$’s [36].

Instead of computing propagators, we use the tAMEn algorithm to solve the first–order ODEs (4) and (7) with Hamiltonians and states represented in the TT format. The next section focuses on mathematical details of this approach.

3 Tensor train format and the tAMEn algorithm

Introducing individual state indices for each site, we can treat $\rho$ as a multi–index (high–dimensional) array which is referred to as tensor in numerical linear algebra:

[TABLE]

We should admit that our use of the term tensor does not imply proper differentiation of upper and lower indices or Einstein summation convention [11]. We use the Tensor Train (TT) decomposition [29] to separate pairs of $i_{k},j_{k}$, belonging to different sites, by a product approximation of the form

[TABLE]

Here, $\rho^{(k)}$, $k=1,\ldots,d,$ are called TT blocks, and the summation ranges $r_{1},\ldots,r_{d}$ are called TT ranks. The grouping of $i_{k}$ and $j_{k}$ is such that unentangled (e.g. pure) states are represented by (8) exactly with elementary TT ranks $r_{0}=\cdots=r_{d}=1$. The TT representation is equivalent to matrix product states (MPS) [14, 26] with open boundary conditions $r_{0}=r_{d}=1$. Entangled states require $r_{k}>1$ for at least one intermediate $k=1,\ldots,d-1.$ However, we aim at representing states with weak correlations, so that $r_{1},\ldots,r_{d-1}$ can be kept bounded by a moderate constant $r_{k}\leqslant r\ll 2^{d}.$ This yields an $\mathcal{O}(dr^{2})$ storage cost of the TT format, linear in the system size. For example, ground states of non–critical one–dimensional spin chains admit such bounded TT ranks due to the area law [12].

We use a linear piecewise Chebyshëv scheme [7] to discretise ODEs (1), (4) and (7) in time. Recall that we solve ODEs on time intervals $t\in(t_{n-1},t_{n}]$, where control pulses $c_{k}(t)$ and hence $\mathaccent 28766{\mathaccent 28766{H}}(t)$ are constant and only $\rho(t)$ depends on time. Without loss of generality we can assume that $t\in(0,\tau]$ and drop the time index $n$ for the rest of this section. We choose a basis of Lagrange polynomials $\{L_{m}(t)\}_{m=1}^{M}$ centered at Chebyshëv nodes $\{\tau_{m}\}_{m=1}^{M}$ and represent the state as

[TABLE]

Since the length of the interval $\tau$ is usually quite small to provide better resolution of control pulses, the basis size $M$ can be very moderate. We collect the nodal values $\rho(\tau_{m})$ into a vector $\overline{\rho}=[\rho(\tau_{m})]_{m=1}^{M}$ of length $4^{d}M,$ which we aim to find. The time derivative operator is discretised as a differentiation matrix $S=[L^{\prime}_{\ell}(\tau_{m})]_{m,\ell=1}^{M}$ in accordance to the spectral approximation theory [39]. Applying this discretisation scheme to ODE (1) we obtain the following linear system

[TABLE]

Here $\overline{\rho}$ is the unknown vector of $\rho(\tau_{m})$’s, $\mathrm{1}_{M}$ is the vector of size $M$ full of $1$’s, and $I$ denotes identity matrices of appropriate size. Note that the Kronecker products are implicitly realised by the TT representation (10). Therefore, we never actually compute them in (9), but just collect the corresponding factors into TT blocks. The same applies to the states. In particular, we consider unentangled states $\rho_{0}$ and $\rho_{T},$ and construct their TT blocks explicitly. Simple mixed states can be made by concatenation of TT blocks instead of a direct summation of a large number of density matrices [34]. The TT ranks can be truncated down to optimal values for a desired accuracy using the singular value decomposition (SVD) [29].

Algebraic equations (9) can be solved in the TT approximation by an alternating iteration, such as the Alternating Minimal Energy (AMEn) algorithm [10], which is an enhanced version of Alternating Least Squares (ALS) [23]. We expand the TT representation (8) by an extra block which carries the time index $m$ and write:

[TABLE]

Now the last block $\rho^{(d+1)}=[\rho^{(d+1)}_{\alpha_{d},\alpha_{d+1}}(m)]$ carries degrees of freedom in time, the rightmost TT rank is $r_{d+1}=1,$ and $r_{d}$ can now be larger than $1.$

The basic ALS algorithm would solve (9) by sweeping over TT blocks, i.e. representing $\overline{\rho}$ by (10) and reducing the equation to only one block $\rho^{(k)}$ in each step. This can be written using a frame matrix

[TABLE]

which is of size $4^{d}M\times 4r_{k-1}r_{k}$ for $k\leqslant d,$ and of size $4^{d}M\times r_{k-1}Mr_{k}$ for $k=d+1.$ The frame matrix realises a linear map from the elements of a single TT block into the elements of the whole vector given (10), that is $\overline{\rho}=P_{\neq k}\rho^{(k)},$ assuming that elements of $\rho^{(k)}$ are stretched in a vector. The ALS method solves the reduced Galerkin systems $(P_{\neq k}^{\dagger}AP_{\neq k})\rho^{(k)}=P_{\neq k}^{\dagger}f$ subsequently for $k=1,\cdots,d+1.$ This system can be assembled efficiently [23] due to the TT formats of $P_{\neq k}$, $A$ and $f$ and solved using standard methods.

The AMEn method [10] improves the convergence of ALS by expanding $\rho^{(k)}$ (and hence all $P_{\neq q}$ for $q>k$) with a TT approximation of the residual $f-A\overline{\rho},$ where $\overline{\rho}$ is formed by (10) with an updated block $\rho^{(k)}$ plugged in. This also allows to adapt the TT ranks to ensure a desired accuracy of the TT approximation. The tAMEn (time–dependent AMEn) algorithm [7] utilises the special meaning of the last TT block, carrying the time variable $m,$ in order to preserve conservation laws (e.g. the Frobenius norm), and to estimate the time discretisation error.

For the interval lengths $\tau$ chosen in our numerical experiments, we have found that $M=8$ Chebyshëv nodes in each interval are sufficient to resolve the time derivative with the relative accuracy of $10^{-5}$ or better. In each step the tAMEn algorithm produces a set of TT blocks $\rho^{(1)},\ldots,\rho^{(d+1)},$ comprising the discrete–time solution $\overline{\rho}$ on the interval $(0,\tau].$ After $\overline{\rho}$ is obtained, we compute $\rho(\tau)$ in the TT format by taking the slice of the last TT block $\rho^{(d+1)}(M)$ and hence removing the time index $m$ from consideration. Since $t=\tau$ represents the end point of the interval $(t_{n-1},t_{n}]$ the obtained state represents $\rho(t_{n})$ and can be used now as an initial state in the next GRAPE step over the time interval $(t_{n},t_{n+1}].$

The tAMEn algorithm is agnostic to a particular form of $\mathaccent 28766{\mathaccent 28766{H}}$ (the only assumption is that it must admit a TT decomposition), and can be applied to any tensor–structured differential equation $\dot{x}=Ax$ such as the auxiliary matrix formalism (4) and the back–propagation (7).

4 Numerical experiments

We consider a Heisenberg chain of spin–$\frac{1}{2}$ particles, for which a Hamiltonian is a simple sum of nearest neighbour interactions

[TABLE]

where the Pauli matrices $\mathaccent 28766{\sigma}_{\{x,y,z\}}^{(k)}$ [37] act only on spin in position $k$ of the chain,

[TABLE]

When $J_{x}=J_{y}=J_{z},$ this system is called the XXX Heisenberg model, when $J_{x}=J_{y}\neq J_{z}$ it is called the XXZ model, and for $J_{x}\neq J_{y}\neq J_{z}$ it is called the XYZ model. It is commonly mentioned that linear Heisenberg chains can be diagonalised exactly using the Bethe ansatz [2], however the eigenstates can only be written via roots of a system of algebraic equations of degree $d$ [37], that are not computable in closed form neither using analytic methods nor numerically with reasonable accuracy for $d\gtrsim 20.$ For simpler XX models ($J_{x}=J_{y},$ $J_{z}=0$) eigenvectors are available in closed form and can be used for dimensionality reduction, allowing chains with $d\simeq 200$ spins to be controlled [44]. We could not find examples of optimal control pulse computed numerically for XXX, XXZ and XYZ Heisenberg models with the number of spins $d\gtrsim 20.$ We decided therefore to test our tensor product approach for XXX and XXZ Heisenberg spin-$\frac{1}{2}$ chains with $d=11,$ $d=21$ and $d=41.$

The initial and target states are taken as

[TABLE]

so the task is to move the $\ket{\uparrow}$ state from the first to the last position in the chain. The control operator $H_{c}(t)=c(t)\mathaccent 28766{\sigma}_{z}^{(1)}$ is the magnetic field acting on the first spin only.

In theory, Heisenberg chains are fully controllable [41], which means the fidelity (2) can be made infinitely close to $1$ as $T\to\infty.$ In practice, however, the available time $T$ for the pulse is limited, and the final state $\rho(T)$ will not be fully focused, leaving some infidelity $1-|\Braket{\rho_{T}}{\rho(T)}|>0.$ We optimise the pulse sequence to reduce the infidelity as much as possible.

The results are shown on Fig. 1. Computing a short high–fidelity pulse for large spin system appears to be a challenging task. Using 4 cores on a Xeon E5-2650 CPU and 10 GB of memory, we reduced infidelity to $10^{-3}$ for $d=21$ and to $10^{-2}$ for $d=41$ within a few days of computation. This seems fast compared to high–fidelity calculations reported in [38], where infidelity $10^{-4}$ was reached after weeks/month of optimisation for a system with $d=4$ qubits with $4$ levels each.

The efficiency of the method depends primarily on reasonable choice of the pulse length $T,$ which is usually of the order of the inverse of the natural frequency of the system [24]. The time $T$ typically should increase with the number of spins in the chain due to the finite rate of the information exchange among them. We see in our experiments that larger intervals $T$ are required to reduce infidelity for larger chains; also, XXZ chains demand larger $T$ than XXX chains. We keep resolution $\tau=10^{-2}$ for the control pulse constant for all experiments, which is the main reason why computational time per iteration grows with $T.$ The TT ranks grow only mildly with $d$ and remain moderate $r_{k}\lesssim 30$ for $d=41.$ This resulted in about $70\cdot 10^{3}$ unknowns in the TT decomposition of each $\rho(t_{k})$ for $d=21$ and $140\cdot 10^{3}$ unknowns for $d=41$. For comparison, a single diagonal of $\rho(t_{k})$ has $2^{21}\approx 2\cdot 10^{6}$ entries for $d=21$ and $2^{41}=2\cdot 10^{12}$ entries for $d=41.$

5 Conclusions and future work

Tensor product algorithms open new possibilities to control long spin wires and multi–qubit gates in quantum computers. We explored a proof-of-concept example, in which we considered a simple XXX or XXZ Heisenberg chain, optimised control pulse using a classical GRAPE method [25], while representing all states n the TT format [29] and propagating them using the tAMEn algorithm [7]. This combination of relatively simple algorithms allowed us to reach fidelity of $99\%$ for a chain of $d=41$ spins using a single workstation for calculations. We used only controllable approximation techniques and avoided any heuristic or random truncations of the state space. Our algorithm is deterministic and flexible, i.e. can be applied to any linear or quasi–linear quantum system, for which we can expect the states to have moderate entanglement (as measured by Schmidt ranks).

We have not yet reached the fidelity of $\sim 99.99\%$ required for topological error correction [38], mostly because of slow first–order convergence of GRAPE, but we hope to make it possible using a second–order optimisation algorithm, such as Newton–Raphson [17] or a quasi–Newton e.g. BFGS [5, 13]. We have not imposed any constraints on pulse shape, which resulted in a few outbreaks seen on Fig. 1, but implementation of the box constraints $|c(t)|<C$ is a straightforward step which we postpone for future work.

Acknowledgements

Authors appreciate financial support of EPSRC which makes this work possible — grants EP/P033954/1 (D.Q. and D.S.) and EP/M019004/1 (S.D.) Numerical calculations for this paper were performed on a dedicated server, provided to D.S. by the University of Brighton (Rising Stars grant).