Time optimal quantum state transfer in a fully-connected quantum computer

As investment into quantum computing platforms continues to grow and the scale of such platforms rapidly expands [1], transferring a quantum state between two distant qubits becomes an important problem. If the physical interactions among the qubits are short-ranged, then the speed of information propagation is always finite according to the Lieb–Robinson bound [2]. As a result, the time for performing quantum state transfer (QST) within such systems has to be at least proportional to the distance between the two qubits. Recently however, long-range interacting quantum systems have become increasingly popular for both quantum computing and quantum simulation, with these systems including trapped ions [3], Rydberg atoms [4], polar molecules [5], defect centers in solids [6], atoms coupled to cavities and photonic crystals [7]. The physical interactions in these systems typically decay as a power law $1/r^{\alpha}$ with r being the inter-particle distance. For sufficiently small α, these long-range interactions can be utilized to provide speedups for quantum information processing [8–11] or realize novel quantum many-body dynamics and phases [12–18]. However, the exact degree of possible speedup for performing QST or other important quantum information processing tasks is unknown [10].

More concretely, with the recent development of optimal Lieb–Robinson-type bounds for a wide-range of long-range interacting systems [19–25], it has been shown that the speed of QST remains finite as long as $\alpha\gt2D+1$, where D is the dimension of the system. However, for systems where the interaction decays with $\alpha\lt D$, no such optimal Lieb–Robinson-type bound has been found. Nonetheless, these systems are important as they could provide a dramatic speedup in quantum information processing. For example, a fully-connected quantum computer [8] (where α = 0) could allow exponentially faster creation of topologically ordered states when compared to quantum computers with short-range interacting qubits [9]. With respect to the QST, the best known Lieb–Robinson-type bound [26] predicts a time of at least $O(\log(N)/N)$ for α = 0, but the fastest known QST protocol in this case requires a minimum time of $O(1/\sqrt{N})$ [26].

While it remains challenging to find an optimal Lieb–Robinson-type bound for generic strongly long-range interacting systems, here we take a different but novel approach to bound the speed limit of QST, using the framework of the quantum brachistochrone (QB) [27–32]. This approach requires us to study only a subclass of Hamiltonians relevant for a fully-connected quantum computer but allows us to obtain an exactly tight bound for the speed of QST. To the best of our knowledge, no such bound has been previously obtained for long-range interacting many-body systems. Our work, therefore, bridges an important gap between the pursuit of optimal Lieb-Robinson-type bounds and the study of quantum speed limits, as the former often only predicts the correct scaling in the large system size limit [2, 22–25] while the latter is precise but is usually only obtained for small quantum systems [27–35]. In addition, our work extends the existing study of QST in spin systems [36–42] that is largely limited to time-independent Hamiltonians and thus does not predict a speed limit for general, time-dependent Hamiltonians studied here. Our work is also complementary to other optimal control methods for speeding up quantum information processing tasks such as the shortcut to adiabaticity [43–47], which seeks to speedup adiabatic processes while preserving robustness, but does not guarantee a time optimal path in general.

The QB method is an important tool in quantum optimal control theory that aims to find the fastest Hamiltonian to achieve a certain unitary or a particular state evolution. While the method can provide an exact speed limit for quantum information processing tasks, it is typically only used for small and simple quantum systems (such as one or two qubits) [29–32]. More importantly, the application of QB typically requires the control parameters of the Hamiltonian to satisfy equality constraints. This is not suitable for studying the speed limit of QST or Lieb–Robinson-type bounds, which usually involves only an upper bound on the strengths of two-body interactions in the Hamiltonian. Therefore, a major contribution of this work is to generalize the standard QB method to incorporate inequality constraints. While there exist other methods in optimal control theory, such as those using Pontryagin's maximum principle [48], that can also deal with inequality constraints [49], our QB based method is technically simpler and more intuitive.

The remaining content of this paper is organized as follows: In section 1, we introduce the particular problem we are trying to solve, i.e. finding the speed limit of QST for a subclass of Hamiltonians realizable by a fully-connected quantum computer. In section 2, we introduce the standard QB method and our extension of it to incorporate inequality constraints. In section 3, we analytically solve the equations from our generalized QB method and obtain an exact speed limit for the QST defined in section 1. The details of the solution are presented in the appendix A. Finally, in section 4, we discuss the implications of our results along with future directions our work can point toward.

In this paper, we consider the problem of (perfect) QST between two qubits in an N-qubit system [36]. Specifically, the goal of QST is to transfer an arbitrary unknown state $c_0 |0\rangle + c_1|1\rangle$ from one qubit to another qubit. Without loss of generality, we may label the source qubit as our qubit 1 and the target qubit as our qubit N. We may further assume that the target qubit is initially in state $|0\rangle$ before the QST and the source qubit is in state $|0\rangle$ after the QST. This assumption can be easily made true via single-qubit gates. Finally, we will always assume all remaining qubits $\{2,3,\ldots, N-1\}\equiv \mathcal{A}$ to be in some state $|\psi_{\mathcal{A}}\rangle$ before and after the QST. With these assumptions in place, the process of QST can be expressed as:

$$\begin{align} \left(c_0 |0\rangle + c_1|1\rangle\right)\otimes |\psi_{\mathcal{A}}\rangle \otimes |0\rangle \rightarrow |0\rangle \otimes |\psi_{\mathcal{A}}\rangle \otimes \left(c_0 |0\rangle + c_1|1\rangle\right). \end{align} \tag{ 1 }$$

Next, we define the class of Hamiltonians considered to achieve the QST described above. For most controllable quantum systems, the interactions among the qubits are two-body, and we can write a general form of the Hamiltonian as:

$$\begin{align} H = \sum_{i\ne j}h_{ij}\left(t\right) + \sum_{i} h_i\left(t\right) \end{align} \tag{ 2 }$$

Due to the allowed time dependence, this Hamiltonian can represent any quantum circuit made of single-qubit and two-qubit gates, and can thus achieve universal quantum information processing. For most experimental systems, the two-qubit interaction strength, quantified by $\lVert h_{ij}\rVert$ where $\lVert\cdot\rVert$ denotes the operator norm of h_ij , is limited by the hardware design, while the strength of single-qubit terms $\{\lVert h_i\rVert\}$ can usually be made much larger via strong applied drives or fields. As a result, it is reasonable to assume that $\{\lVert h_i\rVert\}$ are not constrained while $\lVert h_{ij}\rVert$ is upper bounded by $J_0/r_{ij}^{\alpha}$, where J₀ denotes the maximum interaction strength and $\alpha\unicode{x2A7E}0$ captures how fast interactions decay in the inter-particle distance r_ij between qubits i and j. This assumption is made in most studies of Lieb–Robinson-type bounds [19–25].

In the Heisenberg picture, a Lieb–Robinson-type bound is typically written

$$\begin{align} \lVert \left[A_i\left(t\right),B_j\left(0\right)\right]\rVert \unicode{x2A7D} f\left(t,r_{ij}\right) \end{align} \tag{ 3 }$$

where the operators A_i and B_j act on qubits i and j, respectively, at t = 0. The function $f(t,r_{ij})$ depends on the exact form of the bound and typically increases monotonically as t increases.

To see why equation (3) can lower bound the time it takes to perform the QST defined in equation (1), we set $B_j = \sigma_1^x$ and $A_i = \sigma_N^y$ at t = 0, and assume qubit 1 is initially in the state $|0\rangle$. Then it is simple to see $\langle [A_i(t),B_j(0)]\rangle = 2i$ for a QST assumed to occur in time t. As a result, we must have $f(t,r_{1N})\unicode{x2A7E} ||[A_i(t),B_j(0)]||\unicode{x2A7E} 2$, which can only be satisfied if $t\gt t_{\text{LR}}$ where t_LR depends on $f(t,r_{ij})$ and $r_{1N}$. If the globally minimum time to perform QST is denoted by t_min, we call the bound equation (3) (qualitatively) optimal when $t_{\text{LR}}\unicode{x2A7E} c t_{\text{min}}$ where the constant c is independent of $r_{1N}$ and N. In other words, the speed limit predicted by the Lieb–Robinson-type bound is only loose by a constant factor.

Optimal Lieb–Robinson-type bounds have recently been found for a wide-range of quantum information processing tasks [22, 23, 25] when $\alpha\gt D$, but remain largely undiscovered for the case where $\alpha\lt D$. For $\alpha\lt D$, we often call interactions 'strongly long-ranged' [26] since the total interaction energy per qubit diverges in N. In particular, the case of α = 0 corresponds to a fully-connected universal quantum computer where any two qubits can interact at equal strength. Such a quantum computer can be experimentally realized by trapped ions [8], superconducting qubits [50], atoms coupled to cavities or waveguides [7, 51], etc. In a fully-connected quantum computer, the physical distance between any two qubits is irrelevant, and the QST time should only depend on N. The best known Lieb–Robinson-type bound [26] predicts a minimum QST time $t_{\text{min}} \unicode{x2A7E} O(\log(N)/N)$. However, this bound is believed to be non-optimal as the fastest QST protocol known for α = 0 takes time $t = O(1/\sqrt{N})$.

It is worth mentioning that the fastest QST protocols known for long-range interacting systems so far are conditioned on the intermediate qubits $\{2,3,\ldots,N-1\}$ being in some specific state, usually the state $|\psi_{\mathcal{A}}\rangle = |0\rangle ^{\otimes {N-2}}$ [9, 10, 24–26]. An unconditional QST, which is similar to a remote two-qubit gate, does not have this requirement of the intermediate qubits being in some specific state, but appears to take a much longer time than a conditional QST [10]. We therefore only focus on conditional QST in this paper as our goal is to find the fastest QST protocol for a fully-connected quantum computer. In particular, we will assume from this point that all intermediate qubits are in the $|0\rangle$ state before and after the QST. This condition can be easily realized via initial state preparation on a quantum computer.

In order to find an exact speed limit of QST, we need to further constrain the general Hamiltonian expressed in equation (2). A reasonable constraint is to require the Hamiltonian to conserve the total $\sigma^z\equiv \sum_{i = 1}^N \sigma_i^z$. This constraint is compatible with the QST process defined in equation (1), and leads the Hamiltonian in equation (2) to simplify to

$$\begin{align} H\left(t\right) = \sum_{i\neq j}\left(J_{ij}\left(t\right) \sigma_i^{+}\sigma_j^{-}+\mathrm{h.c.}+U_{ij}\left(t\right)\sigma_i^z\sigma_j^z\right)+\sum_{j}B_j\left(t\right) \sigma_j^z. \end{align} \tag{ 4 }$$

We emphasize that this constraint is in fact natural for most quantum computing or simulation platforms where the physical interactions among qubits are flip–flop, Heisenberg, or XXZ interactions. Such a constraint is also assumed in most previous studies on QST in spin systems [37–42]. For platforms with Ising interactions, this constraint can also be approximately achieved by applying a large transverse field in the z direction, effectively reducing the Ising interactions to flip-flop interactions [12, 17]. We also note that the ZZ-interaction terms $\{U_{ij}(t)\sigma_i^z\sigma_j^z\}$ are not important as they play the same role as the field terms $\{B_j(t) \sigma_j^z\}$ for our QST process because we only ever have one qubit in state $|1\rangle$. Thus we will only need to constrain the flip–flop interaction strength. For a fully-connected quantum computer, the constraint we apply to equation (4) will simply be

$$\begin{align} |J_{ij}\left(t\right)| \unicode{x2A7D} J_0. \end{align} \tag{ 5 }$$

A bound on the speed of QST in this case has been derived by mapping equation (4) to non-interacting bosons/fermions, which predicts $t_{\text{min}} \unicode{x2A7E} 1/(J_0\sqrt{N-1})$ [26]. This bound is qualitatively optimal as [26] includes a QST protocol with $t = \pi/(J_0\sqrt{N-2})$ ($N\unicode{x2A7E}3$ assumed).

The main results of this paper are:

• 1.  
  For any $N\unicode{x2A7E}2$, we find a family of spin Hamiltonians of the form in equation (4) that achieves perfect QST in time

  $$\begin{align} t = \frac{\pi}{J_0\sqrt{2N}}. \end{align} \tag{ 6 }$$
• 2.  
  We prove rigorously that this is the globally minimum QST time assuming the Hamiltonian in equation (4) is invariant under the permutations of the qubits $2,3,\ldots,N$.

The additional assumption on permutation symmetry is natural due to the constraint in equation (5) and the setup of our QST. We do not believe a faster QST protocol could exist without this assumption and will provide justifications and possible ways to prove them in section 4. This assumption is also not needed for N = 3.

The permutation symmetry allows us to restrict the QST to within a subspace spanned by the following four basis states:

$$\begin{align} |\phi_0\rangle &\equiv |0\rangle \otimes |0\rangle ^{\otimes \left(N-2\right)} \otimes |0\rangle, \end{align} \tag{ 7 }$$

$$\begin{align} |\phi_1\rangle &\equiv |1\rangle \otimes |0\rangle ^{\otimes \left(N-2\right)} \otimes |0\rangle, \end{align} \tag{ 8 }$$

$$\begin{align} |\phi_2\rangle &\equiv |0\rangle \otimes |W_{N-2}\rangle \otimes |0\rangle, \end{align} \tag{ 9 }$$

$$\begin{align} |\phi_3\rangle &\equiv |0\rangle \otimes |0\rangle ^{\otimes \left(N-2\right)} \otimes |1\rangle, \end{align} \tag{ 10 }$$

where $|W_{N-2}\rangle \equiv \frac{1}{\sqrt{N-2}} (|10\cdots 0\rangle + \cdots + |0\cdots 1\rangle)$ is the W-state for qubits $2,3,\ldots,N-1$. Since $|\phi_0\rangle$ is an eigenstate of the Hamiltonian in equation (4), the non-trivial dynamics of the Hamiltonian can be captured by an effective 3-level Hamiltonian in the subspace spanned by $\{|\phi_1\rangle,|\phi_2\rangle,|\phi_3\rangle\}$:

$$\begin{align} H_{\text{eff}}\left(t\right) = \begin{pmatrix} \delta_1\left(t\right)& J_{1\mathcal{A}}\left(t\right)&J_{1 N}\left(t\right)\\ J^*_{1\mathcal{A}}\left(t\right)&\delta_{\mathcal{A}}\left(t\right)& J_{\mathcal{A} N}\left(t\right)\\ J^*_{1 N}\left(t\right)& J^*_{\mathcal{A}N}\left(t\right)&\delta_N\left(t\right) \end{pmatrix}, \end{align} \tag{ 11 }$$

where the constraint given by equation (5) now becomes

$$\begin{align} |J_{1N}\left(t\right)| \unicode{x2A7D} J_0, \quad |J_{1\mathcal{A}}\left(t\right)|,|J_{\mathcal{A} N}\left(t\right)|\unicode{x2A7D} \sqrt{N-2} J_0. \end{align} \tag{ 12 }$$

while $\delta_1(t),\delta_{\mathcal{A}}(t),\delta_N(t)$ are unconstrained. Now our goal is to find the optimal H_eff satisfying the constraints in equation (12) that takes the state $|\phi_1\rangle$ to $|\phi_3\rangle$ in minimal time. This is a much simpler problem but still highly non-trivial due to the time dependence of all matrix elements in H_eff, such that we are optimizing over infinitely many possible 3-level Hamiltonians. As it is impossible to find the analytical solution to the evolution operator of equation (11), one cannot obtain the time optimal Hamiltonian directly.

The QB method [27–32] offers us a way to deal with this type of problem, but it does not permit inequality constraints as given in equation (12). As a result, we will need to first generalize the QB method to handle our constraints, as detailed in the next section.

We introduce here an extension to the QB method first developed in [27] that allows for the consideration of problems containing both inequality and equality constraints on the Hamiltonian. This extension is not only necessary to find the time optimal Hamiltonian of the form equation (11) under the constraints in equation (12), but may also be helpful in dealing with more general Hamiltonians in equation (2) with bounded interaction strengths.

We note that [49] has also developed an extension to the standard QB method that allows for inequality constraints, but it uses Pontraygin's maximum principle [48] and the formalism is quite involved. Moreover, it does not intuitively demonstrate the physical difference between a QB problem with equality and inequality constraints. Here we utilize a distinct approach based on the method of slack variables, which intuitively speaking, divides the optimization problem into two scenarios: one where the inequality constraint is simply removed and the other where the inequality constraint becomes an equality constraint.

Formally, the method of slack variables allows us to change all inequality constraints on the Hamiltonian to equality constraints in the form of $f_j(H,s_j) = 0$ for the jth constraint, where s_j is a slack variable. For example, the constraint $|J_{1N}(t)| \unicode{x2A7D} J_0$ can be rewritten as $|J_{1N}(t)|^2+|s(t)|^2 = |J_0|^2$ with a slack variable s(t) added.

We can now introduce our generalized QB method, which also includes a brief review of the standard QB method. Suppose we want to find the minimal evolution time T for a particular Hamiltonian H(t) to generate a target unitary operator U₀. We can formulate this problem as an optimization problem for minimizing the functional $T = \int_0^T dt$ over the Hamiltonian H(t) and the corresponding evolution operator U(t) under the constraints of the Hamiltonian specified by $\{f_{j}(H,s_j)\}$ and the Schrodinger's equation $dU(t)/dt = -i H(t)U(t)$. We can then define an action integral [27, 29, 49]

$$\begin{align} S& \equiv \int_0^T dt \left[1+L_S+L_C\right] \end{align} \tag{ 13 }$$

$$\begin{align} L_S & \equiv \mathrm{Tr}\left[F\left(t\right) \left(i\frac{dU\left(t\right)}{dt}U\left(t\right)^\dagger- H\left(t\right)\right)\right] \end{align} \tag{ 14 }$$

$$\begin{align} L_C &\equiv \sum_j\lambda_j\left(t\right) f_j\left(H,s_j\right) \end{align} \tag{ 15 }$$

where both F(t) (which is a matrix/operator) and $\{\lambda_j(t)\}$ (which are scalars) are Lagrange multipliers. Since H(t) is Hermitian, we can assume F(t) is Hermitian and $\{\lambda_j(t)\}$ are real without loss of generality.

If T is minimized by a time optimal Hamiltonian H(t), according to the variational principle, the variations of S with respect to H(t), U(t), F(t), $\{\lambda_j(t)\}$, $\{s_j(t)\}$, and T should all vanish. In particular, the variation of S with respect to H(t) leads to the equation that defines F(t):

$$\begin{align} F\left(t\right) = \sum_j \lambda_j\left(t\right) \frac{\partial f_j\left(H,s_j\right)}{\partial H}. \end{align} \tag{ 16 }$$

The variation of S with respect to U(t) leads to the so-called QB equation:

$$\begin{align} i\frac{dF\left(t\right)}{dt} = \left[F\left(t\right),H\left(t\right)\right]. \end{align} \tag{ 17 }$$

The variation of S with respect to T and U(T) leads to a normalization condition:

$$\begin{align} \mathrm{Tr}\left[F\left(t\right)H\left(t\right)\right] = 1. \end{align} \tag{ 18 }$$

The variation of S with respect to $\lambda_j(t)$ gives the constraint equation:

$$\begin{align} f_j\left(H,s_j\right) = 0. \end{align} \tag{ 19 }$$

And finally, the variation of S with respect to the slack variables $s_j(t)$ leads to:

$$\begin{align} \lambda_j\left(t\right)\frac{\partial f_j\left(H,s_j\right)}{\partial s_j} = 0. \end{align} \tag{ 20 }$$

This last equation is our key addition to the standard QB method. Let us try to understand intuitively what it means. According to equation (20), we need to have either $\lambda_j(t) = 0$ or $\frac{\partial f^j(H,s_j)}{\partial s_j} = 0$. Notably, the case where $\lambda_j(t) = 0$ is equivalent to removing the constraint $f_j(H,s_j)$ from the action, while the other possible case where $\frac{\partial f^j(H,s_j)}{\partial s_j} = 0$ leads to $s_j(t) = 0$. As an example, we consider the inequality constraint $|J_{1N}(t)|\unicode{x2A7D} J_0$, which is captured by $f(H,s)\equiv|J_{1N}(t)|^2+|s(t)|^2 - |J_0|^2 = 0$. It is easy to see that $\frac{\partial f(H,s)}{\partial s} = 0$ leads to $s(t) = 0$, and we end up with an equality constraint given by $|J_{1N}| = J_0$. On the other hand, when the slack variable $s(t)\ne0$, we end up with $|J_{1N}(t)|\lt J_0$, and we treat $|J_{1N}(t)|$ as effectively unconstrained while solving the QB equation and enforce the inequality only after finding the relevant solution.

However, because both $\lambda_j(t)$ and $s_j(t)$ are time-dependent, it appears that we have to consider the above two cases for every value of time t for equation (20). Fortunately, we can avoid this complication by exploiting the fact that all variables in the action S are smooth functions of t as long as the Hamiltonian H(t) is a smooth function of t, which is indeed true for any physical (experimentally realizable) Hamiltonian. As a result, at least one of the following two conditions will always hold: (1) $\lambda_j(t) = 0$ for all t, and (2) $s_j(t) = 0$ for all t. It is easy to check if neither of the conditions hold, at least one of $\lambda_j(t)$ and $s_j(t)$ has to be discontinuous at some value of t. To avoid double counting the case where $\lambda_j(t) = s_j(t) = 0$ for all t, we will implicitly assume that $s_j(t)$ is not always zero for case (1). As we show below for our concrete QB problem defined in section I, this continuity analysis is crucial to making the inequality constraints manageable.

With the aforementioned recipe for treating inequality constraints in place, our next step is to solve equations (16)–(20) together with the boundary condition that U(T) fulfills the QST. The solution obtained will form the necessary conditions for H(t) to be time optimal. We note that this will lead us to obtain a family $\mathcal{H}$ of H(t) that are all at least locally time optimal Hamiltonians, and we need to compare all members of $\mathcal{H}$ to find the globally time optimal Hamiltonian among them in order to identify the speed limit of QST.

To do so, let us further simplify our effective Hamiltonian H_eff by going into an interaction picture of its diagonal part (which only add phases to the basis states)

$$\begin{align} H_{\text{eff}}^I\left(t\right)& = \begin{pmatrix} 0& \tilde{J}_{1\mathcal{A}}\left(t\right)&\tilde{J}_{1 N}\left(t\right)\\ \tilde{J}^*_{1\mathcal{A}}\left(t\right)& 0 & \tilde{J}_{\mathcal{A} N}\left(t\right)\\ \tilde{J}^*_{1 N}\left(t\right)& \tilde{J}^*_{\mathcal{A}N}\left(t\right)& 0 \end{pmatrix}, \end{align} \tag{ 21 }$$

with the constraints:

$$\begin{align} |\tilde{J}_{1\mathcal{A}}\left(t\right)|& = |J_{1\mathcal{A}}\left(t\right)|\unicode{x2A7D} \sqrt{N-2}J_0 \end{align} \tag{ 22 }$$

$$\begin{align} |\tilde{J}_{\mathcal{A}N}\left(t\right)|& = |J_{\mathcal{A}N}\left(t\right)|\unicode{x2A7D} \sqrt{N-2}J_0 \end{align} \tag{ 23 }$$

$$\begin{align} |\tilde{J}_{1N}\left(t\right)|& = |J_{1N}\left(t\right)|\unicode{x2A7D} J_0. \end{align} \tag{ 24 }$$

Since $H_{\text{eff}}^I(t)$ has zero diagonals, we need to add two more constraints

$$\begin{align} \mathrm{Tr}\left[H\left(t\right) \begin{pmatrix} 1&0&0\\ 0&0&0\\ 0&0&-1 \end{pmatrix}\right]& = 0 \end{align} \tag{ 25 }$$

$$\begin{align} \mathrm{Tr}\left[H\left(t\right) \begin{pmatrix} 1&0&0\\ 0&-2&0\\ 0&0&1 \end{pmatrix}\right]& = 0 \end{align} \tag{ 26 }$$

while a third constraint of $\mathrm{Tr}[H(t)] = 0$ is not needed as removing the constraint only contributes to a global phase [27].

Our goal now is to find the time optimal Hamiltonian $H_{\text{eff}}^I(t)$ that fulfills the QST, which requires evolving the state $|\phi_1\rangle$ to $|\phi_3\rangle$ (defined in equation (7)) up to a phase. The evolution operator in the interaction picture thus needs to take the following form at time T (end of the state transfer):

$$\begin{align} U\left(T\right) = \begin{pmatrix} 0&0&e^{i\phi}\\ \cos \theta e^{-i \alpha}&-\sin \theta e^{-i\beta}&0\\ \sin \theta e^{i\beta}&\cos \theta e^{i \alpha}&0 \end{pmatrix} \end{align} \tag{ 27 }$$

where $\theta,\alpha,\beta$ are arbitrary angles. The phase φ is not important as one can apply a single-qubit gate (which is assumed to take negligible time) on qubit N at the end of the evolution to make sure the phase acquired by $|\phi_0\rangle$ matches the phase acquired by $|\phi_1\rangle$ (which evolves to $|\phi_3\rangle$).

Introducing five Lagrange multipliers $\lambda_{1\mathcal{A}}$, $\lambda_{\mathcal{A}N}$, $\lambda_{1N}$, λ₁, and λ₂ for the above five constraint equations (equations (22)–(26)), respectively, the F(t) operator in equation (16) takes the form:

$$\begin{align} F\left(t\right)& = \begin{pmatrix} \lambda_1+\lambda_2&&\lambda_{1\mathcal{A}} J_{1\mathcal{A}}&&\lambda_{1N} J_{1N}\\ \lambda_{1\mathcal{A}} J_{1\mathcal{A}}^*&&-2\lambda_2&&\lambda_{\mathcal{A}N} J_{\mathcal{A}N}\\ \lambda_{1N} J_{1N}^*&&\lambda_{\mathcal{A}N} J_{\mathcal{A}N}^*&&-\lambda_1+\lambda_2 \end{pmatrix} \end{align} \tag{ 28 }$$

where to simplify the notation we suppress the time-dependence of all terms in equation (28) and exchange $\tilde{J}_{1N}$, $\tilde{J}_{1\mathcal{A}}$, and $\tilde{J}_{\mathcal{A}N}$ for $J_{1N}$, $J_{1\mathcal{A}}$, and $J_{\mathcal{A}N}$ as they obey the same constraints.

The normalization condition equation (18) now reduces to:

$$\begin{align} \lambda_{1\mathcal{A}}|J_{1\mathcal{A}}|^2+\lambda_{\mathcal{A}N}|J_{\mathcal{A}N}|^2+\lambda_{1N}|J_{1N}|^2 = 1. \end{align} \tag{ 29 }$$

And to deal with the inequality constraints, we require from equation (20) that

$$\begin{align} \lambda_{1\mathcal{A}}\left(t\right)s_{1\mathcal{A}}\left(t\right)& = 0 \end{align} \tag{ 30 }$$

$$\begin{align} \lambda_{\mathcal{A}N}\left(t\right)s_{\mathcal{A}N}\left(t\right)& = 0 \end{align} \tag{ 31 }$$

$$\begin{align} \lambda_{1N}\left(t\right)s_{1N}\left(t\right)& = 0. \end{align} \tag{ 32 }$$

As we discussed towards the end of section 2, each equation above requires us to investigate two different cases, resulting in the eight different cases listed in table 1. For each case, we shall solve the QB equation (equation (17)) together with the above equations (27)–(29). We defer the detailed solutions of each case to appendix A, and provide a summary of them in table 1.

Table 1. Summary of the solutions to the QB equation (equation (17)) and accompanying equations (27)–(29) for eight different cases resulting from the three inequality constraints equation (30)–(32) The cases where no minimum time exists is due to an empty set of solutions to these equations. If the solution set is non-empty, we analytically solve the Hamiltonians obeying such solutions and find the minimum QST time. In case 7, $\overline{J}_{1N}\equiv \frac{1}{T}\int_0^T J_{1N}(t)dt$.

#	Case	Minimum QST Time T
1	$\lambda_{1\mathcal{A}} = \lambda_{\mathcal{A}N} = \lambda_{1N} = 0$	No minimum
2	$s_{1\mathcal{A}} = \lambda_{\mathcal{A}N} = \lambda_{1N} = 0$	No minimum
3	$s_{\mathcal{A}N} = \lambda_{1\mathcal{A}} = \lambda_{1N} = 0$	No minimum
4	$s_{1\mathcal{A}} = s_{1N} = \lambda_{\mathcal{A}N} = 0$	No minimum
5	$s_{\mathcal{A}N} = s_{1N} = \lambda_{1\mathcal{A}} = 0$	No minimum
6	$s_{1N} = \lambda_{1\mathcal{A}} = \lambda_{\mathcal{A}N} = 0$	$\pi/(2J_0)$
7	$s_{1\mathcal{A}} = s_{\mathcal{A}N} = \lambda_{1N} = 0$	$\pi/\sqrt{2(N-2)J_0^2+4|\overline{J}_{1N}|^2}$
8	$s_{1\mathcal{A}} = s_{\mathcal{A}N} = s_{1N} = 0$	$\pi/(J_0\sqrt{2N})$

In each case of table 1, whenever a member of the set of slack variables $\{s_{1\mathcal{A}}(t),s_{\mathcal{A}N}(t),s_{1N}(t)\}$ is set to zero for all time, the corresponding interaction strength $\{J_{1\mathcal{A}}(t),J_{\mathcal{A}N}(t),J_{1N}(t)\}$ is maximized for all time. On the other hand, if a member of the set of Lagrange multipliers $\{\lambda_{1\mathcal{A}}(t),\lambda_{\mathcal{A}N}(t),\lambda_{1N}(t)\}$ is set to zero for all time, then the corresponding interaction strength cannot be maximized at all time (otherwise it reduces to the previous case). And since each interaction strength must vary smoothly in time, the interaction strength can only be maximized at a finite number (including zero) of specific time points in this case.

We provide some intuitive explanations of the results in table 1: case 1 will not have a minimum in the QST time because all interaction strengths $\{J_{1\mathcal{A}}(t),J_{\mathcal{A}N}(t),J_{1N}(t)\}$ are not maximized for at least some finite period of time. It is therefore always possible to speed up the QST by maximizing at least one of $\{J_{1\mathcal{A}}(t),J_{\mathcal{A}N}(t),J_{1N}(t)\}$ while increasing the other two proportionally for a finite period of time. Cases 2–5 share a common structure that the interaction strengths are not mirror symmetric, i.e. $J_{1\mathcal{A}}(t)\ne J_{A\mathcal{N}}(t)$ for some finite period of time. It makes intuitive sense that a non-mirror-symmetric Hamiltonian is likely not time optimal, although the actual proof of this is much more involved (see appendix A). Case 6 in fact requires $J_{1\mathcal{A}} = J_{\mathcal{A}N} = 0$ after solving the QB equation and therefore the effective Hamiltonian reduces to a two-level Hamiltonian with a direct coupling strength of $|J_{1N}(t)| = J_0$ between basis states $|\phi_1\rangle$ and $|\phi_3\rangle$, where it is easy to show that the minimum QST time is simply given by $\pi/(2J_0)$. The optimal Hamiltonians in case 7 and case 8 are similar, with the only difference being in case 7 that one cannot maximize $|J_{1N}(t)|$ for all time. Case 7 reduces to case 8 upon setting $|J_{1N}(t)| = J_0$.

Case 8 offers the fastest QST protocol of all cases, and we can in fact construct a simple time-independent Hamiltonian $H_{\text{eff}}^{{\,}\text{opt}}$ in the Schrodinger picture that satisfies case 8:

$$\begin{align} H_{\text{eff}}^{{\,}\text{opt}} = \begin{pmatrix} 0 & \sqrt{N-2} J_0 & J_0 \\ \sqrt{N-2} J_0 & -3 J_0 & \sqrt{N-2} J_0\\ J_0 & \sqrt{N-2} J_0 & 0 \end{pmatrix}. \end{align} \tag{ 33 }$$

Note that there are other equivalently time optimal Hamiltonians satisfying case 8, such as $-H_{\text{eff}}^{{\,}\text{opt}}$. A family of spin Hamiltonians in the form of equation (4) can be mapped to $H_{\text{eff}}^{{\,}\text{opt}}$. One simple example is:

$$\begin{align} H_{\text{opt}} = J_0 \left[\sum_{i\neq j}\left(\sigma_i^{+}\sigma_j^{-}+\mathrm{h.c.}\right)-\frac{N}{2}\left(\sigma_1^z + \sigma_N^z\right)\right]. \end{align} \tag{ 34 }$$

This Hamiltonian can be readily implemented in many experimental platforms with uniform, all-to-all interactions, such as trapped ions [3], neutral atoms in an optical lattice coupled to a cavity [51], circuit QED systems [50], or atoms trapped near a photonic waveguide [7]. Individual addressing of the qubits is required as one needs to apply a different field along z on the source and target qubits. This Hamiltonian was also previously found in the study of perfect QST in fully-connected graphs [52].

For experimental realization of the time-optimal QST Hamiltonian in equation (34), one potential concern is that the experimental Hamiltonian is usually subject to noise or imperfect control. If the fidelity of the QST deteriorates with an increasing number of qubits, then the speedup of the QST brought by additional qubits may not be useful in practice. Fortunately, as shown in appendix B, we find that noise on the interaction and/or the field strengths in equation (34) actually have vanishing effects on the fidelity of QST as N increases, assuming that the noise has a finite amplitude. This means high-fidelity, time-optimal QST can be practically achieved for a large number of qubits.

A drawback of the above optimal spin Hamiltonian H_opt is that the field strength diverges in N, which makes it harder to control the field strength accurately for a larger N. Here we also construct an alternative spin Hamiltonian that maps to the same optimal effective Hamiltonian $H_{\text{eff}}^{{\,}\text{opt}}$ but has a constant field strength applied to the source and target qubits:

$$\begin{align} H_{\text{opt}}^{^{\prime}} = & J_0 \left[\sigma_1^+\sigma_N^- + \left(\sum_{i = 2}^{N-1}\sigma_i^{+}\right)\left(\sigma_1^{-}+\sigma_N^{-}\right)+\mathrm{h.c.} \right] -\frac{3}{2}J_0 \left(\sigma_1^z + \sigma_N^z\right).\nonumber \end{align} \tag{ 35 }$$

The main drawback of this equally optimal Hamiltonian is that one needs to remove the uniform couplings among the qubits $2,3,\ldots,N-1$, which is more challenging to achieve experimentally. Nonetheless, we anticipate that noise on the interaction and/or the field strengths in equation (35) will have vanishing effects on the fidelity of QST as N increases similarly to in the case of noise in equation (B2) due to the similarity between the two protocols. To our best knowledge, this Hamiltonian has not been previously found.

At $T = \pi/(J_0\sqrt{2N})$, one can verify that $e^{-i H_{\text{opt}} T}|\phi_0\rangle \propto|\phi_0\rangle$ and $e^{-i H_{\text{opt}} T}|\phi_1\rangle \propto |\phi_3\rangle$, with some phases that only depend on N. The same applies to $H_{\text{opt}}^{^{\prime}}$. The desired QST described by equation (1) with $|\psi_{\mathcal{A}}\rangle = |0\rangle ^{\otimes {N-2}}$ can then be achieved after performing a single-qubit gate (assumed to take negligible time) on the target qubit to cancel the unwanted phases.

In this paper, we aimed to find the globally time optimal Hamiltonian for realizing QST in a subclass of Hamiltonians relevant for a fully-connected computer. By developing a new QB method that can deal with inequality constraints, we have been able to find the exact minimum time needed for QST to be implemented, with an additional assumption on the permutation symmetry of the Hamiltonian over the intermediate qubits. The QB method also allows us to construct a family of explicit, time-independent Hamiltonians that can be experimentally realized to achieve the exact speed limit of QST. This is in contrast to existing QST protocols for long-range interacting systems [9, 10, 24, 26] that are only qualitatively optimal.

Concretely, we recall that for spin-systems where the physical interactions decay as $1/r^\alpha$ with r being the inter-particle distance, qualitatively optimal bounds on the time to perform QST have only been discovered for the case of $\alpha\gt D$ with protocols existing that perform QST in a time that scales as $O(\log{r})$ for $D\lt\alpha\lt2D$ [10]. However, in the cases of $0\lt\alpha\lt D$, the best known lower bound on the time to perform QST scales as $O(\log{(N)}/N^{1-\frac{\alpha}{D}})$ and is conjectured to be loose. Our work demonstrates that this best known Lieb–Robinson bound is certainly loose with respect to the important subclass of Hamiltonians considered here, and we conjecture that a qualitatively optimal bound for $0\lt\alpha\lt D$ should scale as $1/N^{\frac{1}{2}-\frac{\alpha}{D}}$. Our results also improve the speed of the fastest known QST protocols [26] for $0\lt\alpha\lt D$ by a factor of at least $\sqrt{2}$. For noisy intermediate scale quantum computers, even a constant factor of speedup can noticeably improve the overall fidelity of the quantum circuit [1]. Important to this point, we have also shown that our time-optimal QST Hamiltonian is very robust to control noises. In fact, noises of finite amplitude have vanishing effects on the fidelity of QST as the number of qubit increases.

In the future, we expect the QB method we develop here, with its ability to intuitively deal with inequality constraints, to be helpful in many other quantum optimal control problems where the interaction/drive strengths are only upper bounded. For a Hamiltonian with M inequality constraints, we need to deal with 2^M cases. But as in the problem we studied, if certain symmetries are enforced, the number of independent constraints can be drastically reduced.

An immediate open question following this work is whether one can relax the permutation symmetry assumption we enforced on the Hamiltonian. We expect a positive answer based on how the spatial inversion (mirror) symmetry naturally emerges in the time optimal Hamiltonian from our solutions to the relevant QB equations (i.e. the cases in table 1 violating the mirror symmetry all lead to no time optimal Hamiltonians). The permutation symmetry would likely require a much more cumbersome calculation. One possible way we are currently investigating to carry out such a calculation is to develop a software algorithm that can solve the QB equations automatically. Such an algorithm will also allow solutions to a variety of quantum optimal control problems for much larger systems.

Another future direction related to the current work is to find the speed limit of unconditional QST protocols or two-qubit gates for a fully-connected quantum computer or more generally, strongly long-range interacting systems. We expect our generalized QB method to be useful in this endeavor at least for small system sizes, and finding time optimal two-qubit gates is important for NISQ quantum computers even if only a small number of qubits are interacting with each other [53]. The alternative approach to finding such a speed limit is through the derivations of optimal Lieb-Robinson-type bounds using the Frobenius norm instead of the operator norm, and such bounds can also be used to determine the speed limit of quantum information scrambling [54, 55]. This approach is likely capable of providing the correct scaling of the speed limits and can be complemented by the QB method-based approach we develop here to find more accurate speed limits.

Finally, it is interesting to see whether QB or other methods in quantum optimal theory can be applied to bound the speed limits for multi-qubit gates or the generation of many-body entangled states. Alternatively, it is also worth exploring whether the lack of scalability in most quantum optimal control methods could be compensated by the recently developed Lieb–Robinson-type bounds which often apply to thermodynamic systems but lack the precision needed for finite quantum systems.

We acknowledge funding support from the NSF RAISE-TAQS program under Grant No. CCF-1839232 and the W M Keck Foundation.

All data that support the findings of this study are included within the article (and any supplementary files).

In this appendix, we provide the details of how we solve the QB equation (equation (17)) together with the equations (27)–(29) for all cases listed in table 1.

Defining $N_{\mathcal{A}}\equiv N-2$ and taking $|J_{1\mathcal{A}}(t)|\unicode{x2A7D} J_0$ and $|J_{\mathcal{A}N}(t)|\unicode{x2A7D} J_0$, we note that equation (17) reduces to:

$$\begin{align} \frac{d\lambda_1}{dt}& = \frac{d\lambda_2}{dt} = 0\nonumber\\ i\frac{d\left(\lambda_{1\mathcal{A}}J_{1\mathcal{A}}\right)}{dt}& = \left(-\lambda_1+3\lambda_2\right)J_{1\mathcal{A}}+\left(-\lambda_{1N}+\lambda_{\mathcal{A}N}\right)J_{1N}J_{\mathcal{A}N}^*\nonumber\\ i\frac{d\left(\lambda_{\mathcal{A}N}J_{\mathcal{A}N}\right)}{dt}& = \left(-\lambda_1+3\lambda_2\right)J_{\mathcal{A}N}+\left(-\lambda_{1\mathcal{A}}+\lambda_{1N}\right)J_{1N}J_{1\mathcal{A}}^*\nonumber\\ i\frac{d\left(\lambda_{1N}J_{1N}\right)}{dt}& = -2\lambda_1J_{1N}-N_{\mathcal{A}}\left(\lambda_{1\mathcal{A}}-\lambda_{\mathcal{A}N}\right)J_{1\mathcal{A}}J_{\mathcal{A}N}.\nonumber \end{align} \tag{ A1 }$$

The boundary condition equation (27) becomes:

$$\begin{align} \lambda_1& = 0\nonumber\\ 6\lambda_2\cos^2\theta& = \sqrt{N_\mathcal{A}}|\lambda_{\mathcal{A}N}J_{\mathcal{A}N}|\sin2\theta\cos\left(\alpha+\beta-\text{arg}\left[J_{\mathcal{A}N}\right]\right)\nonumber\\ \lambda_{1\mathcal{A}}J_{1\mathcal{A}}& = \lambda_{\mathcal{A}N}\left(J_{\mathcal{A}N} \cos^2\theta e^{i 2\alpha}-J_{\mathcal{A}N}^* \sin^2\theta e^{-i 2\beta}\right) +3 \frac{\lambda_2}{\sqrt{N_\mathcal{A}}} \sin \theta\cos\theta e^{i\left(\alpha-\beta\right)}\nonumber\\ \lambda_{\mathcal{A}N}J_{\mathcal{A}N}& = e^{-i\phi}\left(\frac{\lambda_{1N}}{\sqrt{N_\mathcal{A}}}J_{1N}^*\cos\theta e^{-i\alpha}-\lambda_{1\mathcal{A}}J_{1\mathcal{A}}^*\sin\theta e^{i \beta}\right)\nonumber\\ \frac{\lambda_{1N}}{\sqrt{N_\mathcal{A}}}J_{1N}& = e^{-i \phi} \left(\frac{\lambda_{1N}}{\sqrt{N_\mathcal{A}}}J_{1N}^*\sin\theta e^{-i\beta}+\lambda_{1\mathcal{A}}J_{1\mathcal{A}}^*\cos\theta e^{i \alpha}\right).\nonumber \end{align} \tag{ A2 }$$

To proceed, we will look at each of the eight cases listed in table 1 one after the other.

A.1.  $\lambda_{1\mathcal{A}} = \lambda_{\mathcal{A}N} = \lambda_{1N} = 0$

A.2.  $s_{1\mathcal{A}}(t) = \lambda_{\mathcal{A}N}(t) = \lambda_{1N}(t) = 0$

In this case, we replace $J_{1\mathcal{A}}(t)$ with $J_0e^{-i \phi_{1\mathcal{A}}(t)}$. The equations (A1) and (A2) now become:

$$\begin{align} 1& = 2N_{\mathcal{A}}\lambda_{1\mathcal{A}}J_0\nonumber\\ \lambda_2 \cos^2\theta\nonumber& = 0\nonumber\\ \lambda_{1\mathcal{A}}J_0 e^{-i \phi_{1\mathcal{A}}}& = \frac{3 \lambda_2 \sin 2\theta e^{i\left(\alpha-\beta\right)}}{2\sqrt{N_{\mathcal{A}}}}\nonumber\\ 0& = \sqrt{N_{\mathcal{A}}}\lambda_{1\mathcal{A}}J_0\sin\theta e^{-i\left(\phi-\beta-\phi_{1\mathcal{A}}\right)}\nonumber\\ 0& = \sqrt{N_{\mathcal{A}}}\lambda_{1\mathcal{A}}J_0\cos\theta e^{-i\left(\phi-\alpha-\phi_{1\mathcal{A}}\right)}. \end{align} \tag{ A3 }$$

These equations are self-contradictory and again indicate that there is no locally time optimal Hamiltonian in this case.

A.3.  $s_{\mathcal{A}N}(t) = \lambda_{1\mathcal{A}}(t) = \lambda_{1N}(t) = 0$

This case is very similar to case 2. Replacing $J_{\mathcal{A}N}(t)$ with $J_0e^{-i \phi_{\mathcal{A}N}(t)}$, equations (A1) and (A2) become:

$$\begin{align} 1& = 2N_{\mathcal{A}}\lambda_{\mathcal{A}N}J_0\nonumber\\ \lambda_2 \cos^2\theta\nonumber& = \sqrt{N_{\mathcal{A}}}\lambda_{\mathcal{A}N}J_0\sin2\theta\cos\left(\alpha+\beta-\phi_{\mathcal{A}N}\right)\nonumber\\ 0& = \lambda_{\mathcal{A}N}J_0\left( \cos^2\theta e^{i \left(2\alpha-\phi_{\mathcal{A}N}\right)}-\sin^2\theta e^{-i \left(2\beta-\phi_{\mathcal{A}N}\right)}\right) +\frac{3 \lambda_2 \sin \theta\cos\theta e^{i\left(\alpha-\beta\right)}}{\sqrt{N_{\mathcal{A}}}}\nonumber\\ 0& = \sqrt{N_{\mathcal{A}}}\lambda_{\mathcal{A}N}J_0e^{-i\phi_{\mathcal{A}N}}. \end{align} \tag{ A4 }$$

These equations are again self-contradictory and indicate that there is no locally time optimal Hamiltonian in this case.

A.4.  $s_{1\mathcal{A}}(t) = s_{1N}(t) = \lambda_{\mathcal{A}N}(t) = 0$

Replacing $J_{1\mathcal{A}}(t)$ with $J_0 e^{-i\phi_{1\mathcal{A}}(t)}$ and $J_{1N}(t)$ with $J_0 e^{-i\phi_{1N}(t)}$ in equation (A1), the equation reduces to:

$$\begin{align} \frac{d\lambda_{1\mathcal{A}}}{dt}& = \lambda_{1N}|J_{\mathcal{A}N}|\sin\left(\phi_{1N}-\phi_{1\mathcal{A}}-\text{arg}\left(J_{\mathcal{A}N}\right)\right)\nonumber\\ \lambda_{1\mathcal{A}}\frac{d\phi_{1\mathcal{A}}}{dt}& = -3\lambda_2-\lambda_{1N}|J_{\mathcal{A}N}|\cos\left(\phi_{1N}-\phi_{1\mathcal{A}}-\text{arg}\left(J_{\mathcal{A}N}\right)\right)\nonumber\\ \frac{d\lambda_{1N}}{dt}& = -N_{\mathcal{A}}\lambda_{1\mathcal{A}}|J_{\mathcal{A}N}|\sin\left(\phi_{1N}-\phi_{1\mathcal{A}}-\text{arg}\left(J_{\mathcal{A}N}\right)\right)\nonumber\\ \lambda_{1N}\frac{d\phi_{1N}}{dt}& = -N_{\mathcal{A}}\lambda_{1\mathcal{A}}|J_{\mathcal{A}N}|\cos\left(\phi_{1N}-\phi_{1\mathcal{A}}-\text{arg}\left(J_{\mathcal{A}N}\right)\right)\nonumber\\ 0& = 3J_{\mathcal{A}N}\lambda_2-\left(\lambda_{1\mathcal{A}}-\lambda_{1N}\right)J_0^2e^{-i\left(\phi_{1N}-\phi_{1\mathcal{A}}\right)}.\nonumber \end{align} \tag{ A5 }$$

The same replacements in equation (A2) lead to:

$$\begin{align} 1& = 2\left(N_{\mathcal{A}}\lambda_{1\mathcal{A}}+\lambda_{1N}\right)J_0^2\nonumber\\ \lambda_2\cos^2\theta& = 0\nonumber\\ \lambda_{1\mathcal{A}}J_0 e^{-i\phi_{1\mathcal{A}}}& = \frac{3\lambda_2\sin2\theta e^{i\left(\alpha-\beta\right)}}{2\sqrt{N_{\mathcal{A}}}}\nonumber\\ 0& = \lambda_{1N}\cos\theta e^{i\left(\phi_{1N}-\alpha\right)} -\sqrt{N_{\mathcal{A}}}\lambda_{1\mathcal{A}}\sin\theta e^{i\left(\phi_{1\mathcal{A}}+\beta\right)}\nonumber\\ \lambda_{1N}e^{-i\left(\phi_{1N}-\phi\right)}& = \sqrt{N_{\mathcal{A}}}\lambda_{1\mathcal{A}}\cos\theta e^{i\left(\alpha+\phi_{1\mathcal{A}}\right)}+\lambda_{1N}\sin\theta e^{i\left(\phi_{1N}-\beta\right)}\nonumber \end{align} \tag{ A6 }$$

implying that either $\lambda_2 = 0$ or $\cos \theta = 0$ such that we will always have $\lambda_{1\mathcal{A}}(T) = 0$.

The first case implies (from equation (A5)) that $\lambda_{1\mathcal{A}} = \lambda_{1N}$ which leads to a contradiction in equation (A6). Therefore, any generated unitary of the form in equation (27) for this case must have $\cos\theta = 0$ such that equation (A6) further reduces to:

$$\begin{align} 1& = 2\left(N_{\mathcal{A}}\lambda_{1\mathcal{A}}+\lambda_{1N}\right)J_0^2\nonumber\\ 0& = \lambda_{1\mathcal{A}}\left(T\right) = \lambda_{1\mathcal{A}}\left(0\right)\nonumber\\ \lambda_{1N}e^{-i\phi_{1N}}& = \pm\lambda_{1N}e^{i\left(\phi_{1N}-\phi-\beta\right)}.\nonumber \end{align} \tag{ A7 }$$

By taking the derivative with respect to t of the first equality of equation (A7), we find:

$$\begin{align} \frac{d\lambda_{1N}}{dt} = -N_{\mathcal{A}}\frac{d\lambda_{1\mathcal{A}}}{dt}. \end{align} \tag{ A8 }$$

Combining equations (A5) and (A8):

$$\begin{align} \frac{d\lambda_{1\mathcal{A}}}{dt}& = \lambda_{1\mathcal{A}}|J_{\mathcal{A}N}|\sin\left(\phi_{1N}-\phi_{1\mathcal{A}}-\text{arg}\left(J_{\mathcal{A}N}\right)\right)\nonumber\\ & = \lambda_{1N}|J_{\mathcal{A}N}|\sin\left(\phi_{1N}-\phi_{1\mathcal{A}}-\text{arg}\left(J_{\mathcal{A}N}\right)\right)\nonumber \end{align} \tag{ A9 }$$

such that equation (A9) can be satisfied only if $\lambda_{1\mathcal{A}} = \lambda_{1N}$, $|J_{\mathcal{A}N}(t)| = 0$, or $\sin(\phi_{1N}-\phi_{1\mathcal{A}}-\text{arg}(J_{\mathcal{A}N})) = 0$, which implies $\frac{d\lambda_{1\mathcal{A}}}{dt} = \frac{d\lambda_{1N}}{dt} = 0$. Combined with equation (A7), this requires $\lambda_{1\mathcal{A}} = 0$, $\lambda_{1N}\neq0$ such that equation (A5) further reduces to:

$$\begin{align} 0& = \frac{d\phi_{1N}}{dt} = \sin\left(\phi_{1N}-\phi_{1\mathcal{A}}-\text{arg}\left(J_{\mathcal{A}N}\right)\right)\nonumber\\ \frac{d \phi_{1\mathcal{A}}}{dt}& = -\frac{d}{dt}\text{arg}\left(J_{\mathcal{A}N}\right)\nonumber\\ \cos\left(n\pi\right)& = \cos\left(\phi_{1N}-\phi_{1\mathcal{A}}-\text{arg}\left(J_{\mathcal{A}N}\right)\right),\,n\in\mathbb{Z}\nonumber\\ 3\lambda_2& = -\lambda_{1N}|J_{\mathcal{A}N}|\cos\left(n\pi\right)\nonumber\\ |J_{\mathcal{A}N}|^2& = J_0^2.\nonumber \end{align} \tag{ A10 }$$

Thus, equation (A10) requires us to replace $J_{\mathcal{A}N}$ with $\pm J_0e^{-i\phi_{\mathcal{A}N}}$ such that we can transform the Hamiltonian of the form in equation (21) to a basis where it is real-valued

$$\begin{align} H^{^{\prime}}\left(t\right)& = e^{-i\int_0^t d\tau H_{\phi}\left(\tau\right)} \left(H+H_{\phi}\right)e^{i\int_0^t d\tau H_{\phi}\left(\tau\right)}\nonumber\\ & = \begin{pmatrix} c_{1\mathcal{A}}&\sqrt{N_{\mathcal{A}}}J_0&J_0\\ \sqrt{N_{\mathcal{A}}}J_0&0&\pm\sqrt{N_{\mathcal{A}}}J_0\\ J_0&\pm\sqrt{N_{\mathcal{A}}}J_0&c_{1\mathcal{A}} \end{pmatrix}\nonumber \end{align} \tag{ A11 }$$

where $c_{1\mathcal{A}} = -\frac{d\phi_{1\mathcal{A}}}{dt}$ is potentially time-dependent and H_φ is given by:

$$\begin{align} H_{\phi}& = \begin{pmatrix} c_{1\mathcal{A}}&0&0\\ 0&0&0\\ 0&0&c_{1\mathcal{A}} \end{pmatrix}. \end{align} \tag{ A12 }$$

Note that the above basis transformation does not affect the universality class of equation (27). Then, a second iteration of the QB method can be performed on equation (A11) with respect to

$$\begin{align} F^{\prime}\left(t\right) = \begin{pmatrix} \lambda_1^{\prime}&\sqrt{N_{\mathcal{A}}}\lambda_{1\mathcal{A}}^{\prime}&\lambda_{1N}^{\prime}\\ \sqrt{N_{\mathcal{A}}}\lambda_{1\mathcal{A}}^{\prime}&0&\pm\sqrt{N_{\mathcal{A}}}\lambda_{1\mathcal{A}}^{\prime}\\ \lambda_{1N}^{\prime}&\pm\sqrt{N_{\mathcal{A}}}\lambda_{1\mathcal{A}}^{\prime}&-\lambda_1^{\prime} \end{pmatrix} \end{align} \tag{ A13 }$$

such that equations (17) and (18) lead to:

$$\begin{align} 1& = 2\left(2N_{\mathcal{A}}\lambda_{1\mathcal{A}}^{\prime}+\lambda_{1N}^{\prime}\right)J_0\nonumber\\ 0& = \lambda^{\prime}_1 = \frac{d\lambda_{1\mathcal{A}}^{\prime}}{dt} = \frac{d\lambda_{1N}^{\prime}}{dt}\nonumber\\ 0& = \lambda_{1\mathcal{A}}^{\prime}\left(c_{1\mathcal{A}}\mp J_0\right)\pm\lambda_{1N}^{\prime} J_0.\nonumber \end{align} \tag{ A14 }$$

Equation (A14) can only be satisfied when $c_{1\mathcal{A}}$ is a constant such that the unitary operator generated by equation (A11) can be expressed as:

$$\begin{align} U_{\pm}\left(t\right)& = \begin{pmatrix} u^{\pm}_{11}&u^{\pm}_{12}&u^{\pm}_{13}\nonumber\\ u^{\pm}_{12}&u^{\pm}_{22}&\pm u^{\pm}_{12}\nonumber\\ u^{\pm}_{13}&\pm u^{\pm}_{12}&u^{\pm}_{11} \end{pmatrix}\nonumber\\ u^{\pm}_{11}& = \frac{e^{it\left(c_{1\mathcal{A}}\pm J_0\right)}}{2} +\frac{e^{i\frac{t}{2}\left(c_{1\mathcal{A}}\mp J_0\right)}}{2}\left(\cos\frac{t\omega_{\pm}}{2}+i\frac{c_{1\mathcal{A}}\mp J_0}{\omega_{\pm}}\sin\frac{t\omega_{\pm}}{2}\right)\nonumber\\ u^{\pm}_{12}& = -ie^{i\frac{t}{2}\left(c_{1\mathcal{A}}\mp J_0\right)}\frac{2\Omega_{1}}{\omega_{\pm}}\sin\frac{t\omega_{\pm}}{2}\nonumber\\ u^{\pm}_{13}& = \pm\left(u_{11}- e^{it\left(c_{1\mathcal{A}}\pm J_0\right)}\right)\nonumber\\ u^{\pm}_{22}& = e^{i\frac{t}{2}\left(c_{1\mathcal{A}}\mp J_0\right)}\left(\cos\frac{t\omega_{\pm}}{2}-i\frac{c_{1\mathcal{A}}\mp J_0}{\omega_{\pm}}\sin\frac{t\omega_{\pm}}{2}\right)\nonumber\\ \omega_{\pm}& = \sqrt{8N_{\mathcal{A}}J^2_{0}+\left(c_{1\mathcal{A}}\mp J_0\right)^2}.\nonumber \end{align} \tag{ A15 }$$

In order for equation (A15) to be equivalent to a state transfer unitary operator of the form in equation (27), we require that $u^{\pm}_{12} = 0$ such that $t = \frac{2m\pi}{\omega_{\pm}}$ where $m\in\mathbb{N}$ in order to respect the ambiguity of the phase associated with evolving with respect to t. Further, we must require that $u^{\pm}_{11} = 0$, such that $\frac{m(c_{1\mathcal{A}}\pm3J_0)}{\omega} = \eta$ where $\eta = 2n$ when $m = \text{odd}$ and $\eta = (2n+1)$ when $m = \text{even}$ with $n\in\mathbb{Z}$. Now, to determine what value of $c_{1\mathcal{A}}\pm 3J_0$ and m will minimize t, we will solve the related problem

$$\begin{align} \text{minimize } &g\left(x,y\right) = \frac{x}{\sqrt{q^2+y^2}}\nonumber\\ \text{given } &\frac{y-r}{\sqrt{q^2+y^2}} = \gamma = \frac{z}{x},\;|z\pm x|\unicode{x2A7E} 1\;\text{and}\; x\unicode{x2A7E} 1\nonumber \end{align} \tag{ A16 }$$

where q and r are some constants. Then, we will find that

$$\begin{align} y = \frac{r\pm |\gamma|\sqrt{r^2+\left(1-\gamma\right)^2q^2}}{1-\gamma^2} \end{align} \tag{ A17 }$$

such that equation (A17) is clearly maximal when $\gamma\rightarrow 1$ which is nearest met when $\gamma = \frac{x\pm 1}{x}$. Further, without loss of generality, we can take $\gamma\unicode{x2A7E}0$ and $r\unicode{x2A7E}0$ such that

$$\begin{align} y = \begin{cases} -\frac{rx^2+\left(x+1\right)\sqrt{r^2x^2-\left(2x+1\right)q^2}}{2x+1},&\gamma = \frac{x+1}{x}\nonumber\\ -\frac{rx^2+\left(x-1\right)\sqrt{r^2x^2+\left(2x-1\right)q^2}}{2x-1},&\gamma = \frac{x-1}{x}. \end{cases} \end{align} \tag{ A18 }$$

We now take the derivative of $g(x,y)$ from equation (A16) with respect to x after substituting in the expressions of equation (A18)

$$\begin{align} \frac{d g\left(x\right)}{dx}& = \begin{cases} \frac{-\left(r^3 x+rs^2 \left(x+2\right)\left(2x+1\right)-\left(r^2-s^2\left(2x+1\right)\right)\tilde{r}_+\right)}{\left(2x+1\right)\left(r^2\left(x+1\right)^2+\tilde{r}_+^2+2rx\left(1+x\right)\tilde{r}_+\right)\tilde{r}_+\tilde{s}_+},& \gamma = \frac{x+1}{x}\nonumber\\ \frac{-\left(r^3 x-rs^2 \left(x-2\right)\left(2x-1\right)-\left(r^2+s^2\left(2x-1\right)\right)\tilde{r}_-\right)}{\left(2x-1\right)\left(r^2\left(x-1\right)^2-\tilde{r}_-^2+2rx\left(1+x\right)\tilde{r}_-\right)\tilde{r}_-\tilde{s}_-},& \gamma = \frac{x-1}{x} \end{cases} \nonumber\\ \tilde{r}_+& = \sqrt{r^2x^2-s^2\left(2x+1\right)}\nonumber\\ \tilde{s}_+& = \sqrt{s^2+\frac{x^2}{\left(2x+1\right)^2}\left(rx+\left(x+1\right)\tilde{r}_+\right)^2}\nonumber\\ \tilde{r}_-& = \sqrt{r^2x^2+s^2\left(2x-1\right)}\nonumber\\ \tilde{s}_-& = \sqrt{s^2+\frac{x^2}{(2x-1)^2}(rx+(x-1)\tilde{r}_-)^2}. \end{align} \tag{ A19 }$$

By inspection, equation (A19) is a decreasing function of x > 0 when $\gamma = \frac{x+1}{x}$ and an increasing function of x > 0 when $\gamma = \frac{x-1}{x}$. Given that both of these cases have g(x) asymptotically approach the same value, the smallest possible value of g(x) occurs for x = 1.

Applying the result of equation (A16), in order for equation (A15) to be equivalent to a state transfer unitary operator of the form in equation (27), we require that $c_{1\mathcal{A}}\pm3J_0 = 0$ and m = 1 such that the target unitary is generated in a minimum time $T_{\text{min}} = \frac{\pi}{\sqrt{2N_{\mathcal{A}}}J_0}$. However, in this case we require $|J_{1\mathcal{A}}(t)| = J_0$ for all time, which requires $s_{1\mathcal{A}} = 0$, and hence this case reduces to case 8.

A.5.  $s_{\mathcal{A}N}(t) = s_{1N}(t) = \lambda_{1\mathcal{A}}(t) = 0$

A.6.  $s_{1N}(t) = \lambda_{1\mathcal{A}}(t) = \lambda_{\mathcal{A}N}(t) = 0$

Replacing $J_{1N}(t)$ with $J_0e^{-i \phi_{1N}(t)}$, equations (A1) and (A2) become:

$$\begin{align} 1& = 2\lambda_{1N}J^2_{0}\nonumber\\ 0& = \frac{d\phi_{1N}}{dt} = \frac{d\lambda_2}{dt}\nonumber\\ 0& = 3J_{1\mathcal{A}}\lambda_2+J_{\mathcal{A}N}^*\lambda_{1N}J_0e^{-i \phi_{1N}}\nonumber\\ 0& = 3J_{\mathcal{A}N}\lambda_2+J_{1\mathcal{A}}^*\lambda_{1N}J_0e^{-i \phi_{1N}}.\nonumber \end{align} \tag{ A20 }$$

Given $\lambda_{1N}J_0\neq0$, equation (A20) can be satisfied if $J_{1\mathcal{A}} = J_{\mathcal{A}N} = 0$ such that $\phi_{1N}$ is an arbitrary constant, or if $\phi_{1N} = m\pi,\,m\in\mathbb{Z}$, $|J_{1\mathcal{A}}(t)| = |J_{\mathcal{A}N}(t)|$, and $\text{arg}(J_{1\mathcal{A}}) = \text{arg}(J_{\mathcal{A}N})+n\pi,\,n\in\mathbb{Z}$.

In the case where $J_{1\mathcal{A}}(t) = J_{\mathcal{A}N}(t) = 0$, the unitary generated by equation (21) can be expressed as:

$$\begin{align} U\left(t\right) = \begin{pmatrix} \cos J_0 t&0& -ie^{-i\phi_{1N}}\sin J_0 t\\ 0&1&0\\ -ie^{-i\phi_{1N}}\sin J_0 t&0&\cos J_0t \end{pmatrix} \end{align} \tag{ A21 }$$

which is equivalent to a state transfer unitary operator of the form in equation (27) in a minimum time of $T_{\text{min}} = \frac{\pi}{2J_0}$.

In the case where $\phi_{1N} = m\pi,\,m\in\mathbb{Z}$, $|J_{1\mathcal{A}}(t)| = |J_{\mathcal{A}N}(t)|$, and $\text{arg}(J_{1\mathcal{A}}) = \text{arg}(J_{\mathcal{A}N})+n\pi,\,n\in\mathbb{Z}$, the Hamiltonian operator of the form in equation (21) can be transformed to a basis where it is real-valued

$$\begin{align} H^{^{\prime}}\left(t\right)& = e^{-i\int_0^t d\tau H_{\phi}\left(\tau\right)} \left(H+H_{\phi}\right)e^{i\int_0^t d\tau H_{\phi}\left(\tau\right)}\nonumber\\ & = \begin{pmatrix} -\frac{d}{dt}\text{arg}\left(J_{1\mathcal{A}}\right)&\sqrt{N_{\mathcal{A}}}|J_{1\mathcal{A}}\left(t\right)|&\pm J_0\\ \sqrt{N_{\mathcal{A}}}|J_{1\mathcal{A}}\left(t\right)|&0&\sqrt{N_{\mathcal{A}}}|J_{1\mathcal{A}}\left(t\right)|\\ \pm J_0&\sqrt{N_{\mathcal{A}}}|J_{1\mathcal{A}}\left(t\right)|&-\frac{d}{dt}\text{arg}\left(J_{1\mathcal{A}}\right) \end{pmatrix}\nonumber\\ H_{\phi}& = \begin{pmatrix} -\frac{d}{dt}\text{arg}\left(J_{1\mathcal{A}}\right)&0&0\\ 0&0&0\\ 0&0&-\frac{d}{dt}\text{arg}\left(J_{1\mathcal{A}}\right) \end{pmatrix}\nonumber \end{align} \tag{ A22 }$$

where this basis transformation does not affect the constraints on the Hamiltonian. Therefore, a second iteration of the QB method can be performed on equation (A22) with respect to an $F^{\prime}(t)$ operator

$$\begin{align} F^{\prime}\left(t\right) = \begin{pmatrix} \lambda_1^{\prime}&\lambda_3^{\prime}&\lambda_{1N}^{\prime}\\ \lambda_3^{\prime}&0&-\lambda_3^{\prime}\\ \lambda_{1N}^{\prime}&-\lambda_3^{\prime}&-\lambda_1^{\prime} \end{pmatrix} \end{align} \tag{ A23 }$$

where equations (17) and (18) become:

$$\begin{align} 1& = \pm2\lambda_{1N}^{\prime} J_0\nonumber\\ 0& = \frac{d\lambda_1^{\prime}}{dt} = \frac{d\lambda_3^{\prime}}{dt} = \frac{d\lambda_{1N}^{\prime}}{dt}\nonumber\\ 0& = \sqrt{N_{\mathcal{A}}}|J_{1\mathcal{A}}|\left(\lambda_1^{\prime}+\lambda_{1N}^{\prime}\right)-\left(\frac{d}{dt}\text{arg}\left(J_{1\mathcal{A}}\right)\pm J_0\right)\lambda_3^{\prime}\nonumber\\ 0& = \sqrt{N_{\mathcal{A}}}|J_{1\mathcal{A}}|\left(\lambda_1^{\prime}-\lambda_{1N}^{\prime}\right)-\left(\frac{d}{dt}\text{arg}\left(J_{1\mathcal{A}}\right)\pm J_0\right)\lambda_3^{\prime}\nonumber\\ 0& = \sqrt{N_{\mathcal{A}}}|J_{1\mathcal{A}}|\lambda_3^{\prime}\pm J_0\lambda_1^{\prime}.\nonumber \end{align} \tag{ A24 }$$

In order to satisfy equation (A24), we must require that $|J_{1\mathcal{A}}| = \lambda_1^{\prime} = 0$. Given the phase term $\text{arg}(J_{1\mathcal{A}})$ will have no physical source, we must further require that $\frac{d}{dt}\text{arg}(J_{1\mathcal{A}}) = 0$. Then, the unitary generated by equation (A22) can be expressed by equation (A21), which is equivalent to a state transfer unitary operator of the form in equation (27) in a minimum time of $T_{\text{min}} = \frac{\pi}{2J_0}$.

A.7.  $s_{1\mathcal{A}}(t) = s_{\mathcal{A}N}(t) = \lambda_{1N}(t) = 0$

Replacing $J_{1\mathcal{A}}(t)$ with $J_0 e^{-i\phi_{1\mathcal{A}}(t)}$ and $J_{\mathcal{A}N}(t)$ with $J_0 e^{-i\phi_{\mathcal{A}N}(t)}$, equations (A1) and (A2) reduce to:

$$\begin{align} 0& = N_{\mathcal{A}}\left(\lambda_{1\mathcal{A}}-\lambda_{\mathcal{A}N}\right)J^2_{0}\nonumber\\ 1& = 2N_{\mathcal{A}}\left(\lambda_{1\mathcal{A}}+\lambda_{\mathcal{A}N}\right)J^2_{0}\nonumber\\ \frac{d\lambda_{1\mathcal{A}}}{dt}& = -\lambda_{\mathcal{A}N}|J_{1N}|\sin\left(\text{arg}\left(J_{1N}\right)-\phi_{1\mathcal{A}}-\phi_{\mathcal{A}N}\right)\nonumber\\ \lambda_{1\mathcal{A}}\frac{d\phi_{1\mathcal{A}}}{dt}& = -3\lambda_2+\lambda_{\mathcal{A}N}|J_{1N}|\cos\left(\text{arg}\left(J_{1N}\right)-\phi_{1\mathcal{A}}-\phi_{\mathcal{A}N}\right)\nonumber\\ \frac{d\lambda_{\mathcal{A}N}}{dt}& = \lambda_{1\mathcal{A}}|J_{1N}|\sin\left(\text{arg}\left(J_{1N}\right)-\phi_{1\mathcal{A}}-\phi_{\mathcal{A}N}\right)\nonumber\\ \lambda_{\mathcal{A}N}\frac{d\phi_{\mathcal{A}N}}{dt}& = 3\lambda_2-\lambda_{1\mathcal{A}}|J_{1N}|\cos\left(\text{arg}\left(J_{1N}\right)-\phi_{1\mathcal{A}}-\phi_{\mathcal{A}N}\right).\nonumber \end{align} \tag{ A25 }$$

In order to satisfy equation (A25), we must require that $\lambda_{1\mathcal{A}}(t) = \lambda_{\mathcal{A}N}(t)$ and $\frac{d\phi_{1\mathcal{A}}}{dt} = -\frac{d\phi_{\mathcal{A}N}}{dt}$, such that it can only be satisfied if $|J_{1N}| = 0$ or if $\sin(\text{arg}(J_{1N})-\phi_{1\mathcal{A}}-\phi_{\mathcal{A}N}) = 0$

In the case where $|J_{1N}(t)| = 0$, we can transform the Hamiltonian of the form in equation (21) to a basis where it is constant, real-valued

$$\begin{align} H^{^{\prime}}\left(t\right)& = e^{-i t H_{\phi}} \left(H+H_{\phi}\right)e^{i t H_{\phi}}\nonumber\\ & = \begin{pmatrix} c_{1\mathcal{A}}&\sqrt{N_{\mathcal{A}}}J_0&0\\ \sqrt{N_{\mathcal{A}}}J_0&0&\sqrt{N_{\mathcal{A}}}J_0\\ 0&\sqrt{N_{\mathcal{A}}}J_0&c_{1\mathcal{A}} \end{pmatrix}\nonumber \end{align} \tag{ A26 }$$

where $c_{1\mathcal{A}} = -\frac{d\phi_{1\mathcal{A}}}{dt} = 12N_{\mathcal{A}}\lambda_2J_0^2$ and H_φ is given in equation (A12). Given that $c_{1\mathcal{A}}$ is proportional to λ₂ it can be taken as a free parameter. Now, the unitary operator generated by equation (A26) can be expressed

$$\begin{align} U\left(t\right) = & \begin{pmatrix} u_{11}&u_{12}&u_{13}\\ u_{12}&u_{22}&u_{12}\\ u_{13}&u_{12}&u_{11} \end{pmatrix}\nonumber\\ u_{11}& = \frac{e^{-it c_{1\mathcal{A}}}}{2}+\frac{e^{-i\frac{t}{2}c_{1\mathcal{A}}}}{2}\left(\cos\frac{t\omega}{2}-i\frac{c_{1\mathcal{A}}}{\omega}\sin\frac{t\omega}{2}\right)\nonumber\\ u_{12}& = -ie^{-i\frac{t}{2}c_{1\mathcal{A}}}\frac{\sqrt{N_{\mathcal{A}}}J_0}{\omega}\sin\frac{t\omega}{2}\nonumber\\ u_{13}& = u_{11}- e^{-itc_{1\mathcal{A}}}\nonumber\\ u_{22}& = e^{-i\frac{t}{2}c_{1\mathcal{A}}}\left(\cos\frac{t\omega}{2}+i\frac{c_{1\mathcal{A}}}{\omega}\sin\frac{t\omega}{2}\right)\nonumber\\ \omega& = \sqrt{8N_{\mathcal{A}}J^2_{0}+c^2_{1\mathcal{A}}}.\nonumber \end{align} \tag{ A27 }$$

In order for equation (A27) to be equivalent to a state transfer unitary operator of the form in equation (27), we require that $u_{12} = 0$ such that $t = \frac{2m\pi}{\omega}$ where $m\in\mathbb{N}$ in order to respect the ambiguity of the phase associated with evolving with respect to t. Further, we must require that $u_{11} = 0$, such that $\frac{mc_{1\mathcal{A}}}{\omega} = 0$ where $\eta = 2 n$ when $m = \text{odd}$ and $\eta = (2n+1)$ when $m = \text{even}$ with $n\in\mathbb{Z}$.

Applying the result of equation (A16), we find that in order for equation (A27) to be equivalent to a state transfer unitary operator of the form in equation (27), we require that $c_{1\mathcal{A}} = 0$ and m = 1 such that the target unitary is generated in a minimum time $T_{\text{min}} = \frac{\pi}{\sqrt{2(N-2)J_0^2}}$.

In the case where $\sin(\text{arg}(J_{1N})-\phi_{1\mathcal{A}}-\phi_{\mathcal{A}N}) = 0$, we can transform the Hamiltonian of the form in equation (21) to a basis where it is real-valued

$$\begin{align} H^{^{\prime}}\left(t\right)& = e^{-i t H_{\phi}} \left(H+H_{\phi}\right)e^{i t H_{\phi}}\nonumber\\ & = \begin{pmatrix} c_{1\mathcal{A}}-|J_{1N}|&\sqrt{N_{\mathcal{A}}}J_0&|J_{1N}|\\ \sqrt{N_{\mathcal{A}}}J_0&0&\sqrt{N_{\mathcal{A}}}J_0\\ |J_{1N}|&\sqrt{N_{\mathcal{A}}}J_0&c_{1\mathcal{A}}-|J_{1N}| \end{pmatrix}\nonumber\\ H_{\phi}& = \begin{pmatrix} c_{1\mathcal{A}}-|J_{1N}|&0&0\\ 0&0&0\\ 0&0&c_{1\mathcal{A}}-|J_{1N}| \end{pmatrix}\nonumber \end{align} \tag{ A28 }$$

where $c_{1\mathcal{A}} = -\frac{d\phi_{1\mathcal{A}}}{dt} = 12N_{\mathcal{A}}\lambda_2J_0^2$. Given that $c_{1\mathcal{A}}$ is proportional to λ₂ it can be taken as a free-parameter, and because the terms dependent on $|J_{1N}(t)|$ commute with the rest of the Hamiltonian, the unitary operator generated by equation (A28) can be expressed

$$\begin{align} U\left(t\right)& = \begin{pmatrix} u_{11}&u_{12}&u_{13}\\ u_{12}&u_{22}&u_{12}\\ u_{13}&u_{12}&u_{11} \end{pmatrix}\nonumber\\ u_{11}& = \frac{e^{-it\left(c_{1\mathcal{A}}-2|\overline{J}_{1N}|\right)}}{2}+\frac{e^{-i\frac{t}{2}c_{1\mathcal{A}}}}{2}\left(\cos\frac{t\omega}{2}+i\frac{c_{1\mathcal{A}}}{\omega}\sin\frac{t\omega}{2}\right)\nonumber\\ u_{12}& = -ie^{-i\frac{t}{2}c_{1\mathcal{A}}}\frac{\sqrt{N_{\mathcal{A}}}J_0}{\omega}\sin\frac{t\omega}{2}\nonumber\\ u_{13}& = u_{11}- e^{it\left(c_{1\mathcal{A}}-2|\overline{J}_{1N}|\right)}\nonumber\\ u_{22}& = e^{-i\frac{t}{2}c_{1\mathcal{A}}}\left(\cos\frac{t\omega}{2}+i\frac{c_{1\mathcal{A}}}{\omega}\sin\frac{t\omega}{2}\right)\nonumber\\ \omega& = \sqrt{8N_{\mathcal{A}}J^2_{0}+c^2_{1\mathcal{A}}}.\nonumber \end{align} \tag{ A29 }$$

In order for equation (A29) to be equivalent to a state transfer unitary operator of the form in equation (27), we require that $u_{12} = 0$ such that $t = \frac{2m\pi}{\omega}$ where $m\in\mathbb{N}$ in order to respect the ambiguity of the phase associated with evolving with respect to t. Further, we must require that $u_{11} = 0$, such that $\frac{m(c_{1\mathcal{A}}-4 |\overline{J}_{1N}|)}{\omega} = 0$ where $\eta = 2 n$ when $m = \text{odd}$ and $\eta = (2n+1)$ when $m = \text{even}$ with $n\in\mathbb{Z}$.

Applying the result of equation (A16), we find that in order for equation (A29) to be equivalent to a state transfer unitary operator of the form in equation (27), we require that $c_{1\mathcal{A}} = 4 |\overline{J}_{1N}|\unicode{x2A7D}4J_0$ and m = 1 such that the target unitary is generated in a minimum time $T_{\text{min}} = \frac{\pi}{\sqrt{2}\sqrt{N_{\mathcal{A}}J_0^2+2\overline{J}_{1N}^2}}\unicode{x2A7E}\frac{\pi}{\sqrt{2N}J_0}$ which saturates only when $|J_{1N}(t)| = J_0$. Hence, at its saturation limit, this case reduces to case 8.

A.8.  $s_{1\mathcal{A}}(t) = s_{\mathcal{A}N}(t) = s_{1N}(t) = 0$

In equations (34) and (35), two experimentally realizable Hamiltonians that can implement QST in the minimal time are presented. In practice, these Hamiltonians will experience some noise in their implementation. Here we provide an analysis of the performance of QST using equation (34) when the control on the interaction strengths and single-qubit fields is assumed to be noisy. We expect similar results for the QST Hamiltonian in equation (35).

We add noise terms to the Hamiltonian in equation (34) as:

$$\begin{align} H_{\text{noisy}} = H_{\text{opt}} + \sum_{i\neq j}\epsilon_{ij} \left(\sigma_i^{+}\sigma_j^{-}+\mathrm{h.c.}\right) + \epsilon_1\sigma_1^z + \epsilon_N \sigma_N^z \end{align} \tag{ B1 }$$

where $\{\epsilon_{ij}\}$ are noisy coefficients sampled from a normal distribution of zero mean and standard deviation σ_c , while $\{\epsilon_1,\epsilon_N\}$ are sampled from a normal distribution of zero mean and standard deviation σ_f .

Since this noisy Hamiltonian does not change the state of the system if all qubits are in $|0\rangle$, the fidelity of the QST for an arbitrary state of the source qubit is lower bounded by the fidelity of performing QST when the source qubit is in $|1\rangle$. We can express this latter fidelity as

$$\begin{align} F = \langle \phi_1 | e^{i H_{\text{opt}} t} e^{-i H_{\text{noisy}} t} | \phi_1 \rangle \end{align} \tag{ B2 }$$

where $|\phi_1\rangle$ is defined in equation (7) and $t = \frac{\pi}{\sqrt{2N}}$.

In the following, we will show that the infidelity $1-F$ should scale linearly in σ_c and σ_f , and more importantly, as $1/\sqrt{N}$. As a result, as long as the noise amplitudes σ_c and σ_f are finite, the infidelity should vanish in the large N limit. Alternatively, we can argue that the infidelity can remain finite and small as long as $\sigma_c, \sigma_f \ll \sqrt{N}$. Therefore, we believe our time optimal QST protocol is robust against control noises.

We first provide some intuitive analysis to explain the above mentioned scaling. If only the noise terms $\epsilon_1\sigma_1^z + \epsilon_N \sigma_N^z$ are present in equation (B1), a simple first-order perturbation theory shows that

$$\begin{align} 1-F \approx |\langle \phi_1 | \left(\epsilon_1 - \epsilon_N\right) t |\phi_1 \rangle|. \end{align} \tag{ B3 }$$

Since $t = O(1/\sqrt{N})$, we have $1-F = O(\sigma_f/\sqrt{N})$.

Next, we consider the case where only the noise on the interactions $\sum_{i\neq j}\epsilon_{ij} (\sigma_i^{+}\sigma_j^{-}+\mathrm{h.c.}$) are present in equation (B1). Intuitively, these additional interactions are disordered and should have little effect on the QST. The time optimal QST achieved by equation (34) is through the collective couplings between the qubits $1,N$ and the qubits $2,3,\ldots,N$. We therefore expect the noise on the interactions to mainly perturb the collective coupling strength $J_{1\mathcal{A}}$ and $J_{\mathcal{A}N}$ in equation (11), with an amplitude of $\sigma_c$. Since $J_{1\mathcal{A}} = J_{\mathcal{A}N} = O(\sqrt{N})$ for equation (34), as long as $\sigma_c \ll \sqrt{N}$, the effective Hamiltonian in equation (11) is only weakly perturbed and a small infidelity is therefore expected.

We have also performed exact numerical calculations to support the above intuitive analysis. By averaging $1-F$ calculated using equation (B2) over 1000 different random realizations of the noisy Hamiltonian in equation  (B1), we show the dependence of this statistically averaged infidelity as a function of N in figure  1 and as a function of σ_c or σ_f in figure  2. We clearly see that the dependence of the infidelity on N is approximately given by $O(1/\sqrt{N})$, while the infidelity scales linearly with either σ_c or σ_f . For a system of N = 500 qubits, the QST fidelity is above $99.5\%$ even for a fairly large noise amplitude of $\sigma_{c,f} = 0.1$, i.e. $10\%$ of the interaction strength.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size