Preparing Greenberger-Horne-Zeilinger and $W$ states on a long-range Ising spin model by global controls

TL;DR

This paper develops methods to prepare GHZ and W entangled states on a long-range Ising spin model using global controls, providing analytical solutions for GHZ and a scalable approach for W states.

Contribution

It introduces analytical solutions for GHZ state preparation and a scalable numerical method for W states in long-range Ising models using global controls.

Findings

Analytical solutions for GHZ state preparation for any system size.

A standard, scalable method for W state preparation demonstrated numerically.

Number of parameters increases linearly with system size.

Abstract

Entanglement, a unique quantum resource with no classical counterpart, remains at the heart of quantum information. The Greenberger-Horne-Zeilinger (GHZ) and $W$ states are two inequivalent classes of multipartite entangled states which can not be transformed into each other by means of local operations and classic communication. In this paper, we present the methods to prepare the GHZ and $W$ states via global controls on a long-range Ising spin model. For the GHZ state, general solutions are analytically obtained for an arbitrary-size spin system, while for the $W$ state, we find a standard way to prepare the $W$ state that is analytically illustrated in three- and four-spin systems and numerically demonstrated for larger-size systems. The number of parameters required in the numerical search increases only linearly with the size of the system.

Figures

Equations

$$|\mathrm{GHZ}⟩=\frac{1}{\sqrt{2}}\left(|000\cdots 0⟩+|111\cdots 1⟩\right)$$

$$|W⟩=\frac{1}{\sqrt{n}}(|100\cdots 0⟩+|010\cdots 0⟩$$

$$+|001\cdots 0⟩+\cdots +|000\cdots 1⟩).$$

$$H_{\text{Ising}}=H_{zz}=\sum _{k<m}^n{\sigma }_k^z{\sigma }_m^z,$$

$$H_x=\sum _{k=1}^n{\sigma }_k^x,H_y=\sum _{k=1}^n{\sigma }_k^y,$$

$$H\left(t\right)=H_{zz}+f\left(t\right)H_x+g\left(t\right)H_y,$$

$$U(M,{\tau }_m,{\beta }_m^x,{\beta }_m^y)=\prod _{m=1}^M{e}^{-i{H}_{zz}{\tau }_m}{e}^{-i{\beta }_m^x{H}_x}{e}^{-i{\beta }_m^y{H}_y}.$$

$$F\left[U\right]=\left|⟨T|U(M,{\tau }_m,{\beta }_m^x,{\beta }_m^y)|00\cdots 0⟩\right|.$$

$$\begin{array}{rl}& ZZ(\tau )=e^{-i\tau H_{zz}},\\ & X(\beta )=e^{-i\beta H_x},\\ & Y(\beta )=e^{-i\beta H_y}.\end{array}$$

$$Y(\pi /4)-ZZ(\pi /4)-X(\pi /4).$$

$$Y\left(\frac{\pi }{4}\right)|0\cdots 00⟩=\frac{1}{\sqrt{2^n}}\sum _{k=0}^{2^n-1}|k⟩.$$

$$H_{\text{Ising}}|k⟩=\left\{\frac{n(n-1)}{2}-2s(k)\left[n-s(k)\right]\right\}|k⟩,$$

$$\begin{array}{rl}& ZZ\left(\frac{\pi }{4}\right)Y\left(\frac{\pi }{4}\right)|0\cdots 00⟩\\ & =\frac{1}{\sqrt{2^n}}\sum _{k=0}^{2^n-1}{e}^{-i\left\{\frac{n(n-1)}{2}-2s(k)\left[n-s(k)\right]\right\}\frac{\pi }{4}}|k⟩\\ & =e^{-i\frac{n(n-1)\pi }{8}}\frac{1}{\sqrt{2^n}}\sum _{k=0}^{2^n-1}{i}^{s(k)\left[n-s(k)\right]}|k⟩.\end{array}$$

$$\begin{array}{rl}& X\left(-\frac{\pi }{4}\right)|\text{GHZ}⟩\\ & =\frac{1}{{\sqrt{2}}^{n+1}}\left[(|0⟩+i|1⟩{)}^{\otimes n}+(i|0⟩+|1⟩{)}^{\otimes n}\right]\\ & =\frac{1}{{\sqrt{2}}^{n+1}}\sum _{k=0}^{2^n-1}\left(i^{s(k)}+i^{n-s(k)}\right)|k⟩.\end{array}$$

$$\begin{array}{rl}f(s(k))& =\left(i^{s(k)}+i^{n-s(k)}\right)/i^{s(k)[n-s(k)]}\\ & =i^{s(k)-s(k)[n-s(k)]}+i^{n-s(k)-s(k)[n-s(k)]}\\ & =i^{s(k)[1-n+s(k)]}+i^{[1-s(k)][n-s(k)]}.\end{array}$$

$$f(s(k)+1)=i^{[s(k)+1][2-n+s(k)]}+i^{s(k)[1-n+s(k)]}.$$

$$f(s(k))=f(s(k)+1),$$

$$|00\cdots 0⟩=\frac{1}{\sqrt{2^n}}{\left(|+⟩+|-⟩\right)}^{\otimes n}=\frac{1}{\sqrt{2^n}}\sum _{k=0}^{2^n-1}|k{⟩}_x.$$

$$XX\left(\frac{\pi }{4}\right)|k{⟩}_x=e^{-i\frac{n(n-1)\pi }{8}}{i}^{s_x(k)[n-s_x(k)]}|k{⟩}_x,$$

$$\begin{array}{rl}& XX\left(\frac{\pi }{4}\right)|00\cdots 0⟩\\ & =\frac{1}{{\sqrt{2}}^{n+1}}\sum _{k=0}^{2^n-1}\left[1+(-1{)}^{n-s_x(k)}{i}^{n+1}\right]|k{⟩}_x\\ & =\frac{1}{\sqrt{2}}\left(|00\cdots 0⟩+i^{n+1}|11\cdots 1⟩\right).\end{array}$$

$$Z\left(-\frac{(n+1)\pi }{4n}\right)=Y\left(\pi /4\right)-X\left(-\frac{(n+1)\pi }{4n}\right)-Y\left(-\pi /4\right).$$

$$Y\left(\pi /4\right)-ZZ\left(\pi /4\right)-X\left(-\frac{(n+1)\pi }{4n}\right)-Y\left(-\pi /4\right).$$

$$|{\phi }_m⟩=\frac{1}{\sqrt{C_n^{m-1}}}\sum _{s(k)=m-1}|k⟩.$$

$${\mathcal{X}}_{\pm }:|x_k^{\pm }⟩=c\left(|{\phi }_k⟩\pm |{\phi }_{n+2-k}⟩\right),$$

$${\mathcal{Y}}_{\pm }:|y_k^{\pm }⟩=c\left[|{\phi }_k⟩\pm i^n{\left(-1\right)}^{k-1}|{\phi }_{n+2-k}⟩\right],$$

$$Y\left(\frac{\pi }{4}\right)|00\cdots 0⟩=\sum _{k=0}^{\frac{n-1}{2}}\sqrt{\frac{C_n^k}{2^{n-1}}}|x_{k+1}^{+}⟩;$$

$$Y\left(-\frac{\pi }{4}\right)|W⟩=\sum _{k=0}^{\frac{n-1}{2}}\sqrt{\frac{C_n^k}{2^{n-1}n}}{\left(-1\right)}^k\left(n-2k\right)|x_{k+1}^{+}⟩.$$

$$\begin{gathered} Y\left(\frac{\pi }{4}\right)|00\cdots 0⟩ \\ =\sum _{k=0}^{\frac{n}{2}-1}\sqrt{\frac{C_n^k}{2^{n-1}}}|x_{k+1}^{+}⟩+\sqrt{\frac{C_n^{n/2}}{2^n}}|x_{n/2+1}^{+}⟩, \end{gathered}$$

$$X\left(-\frac{\pi }{4}\right)|W⟩=\sum _{k=0}^{\frac{n}{2}-1}\sqrt{\frac{C_n^k}{2^{n-1}n}}\left[i^{k-1}\left(2k-n\right)\right]|y_{k+1}^{-}⟩.$$

$$Y\left(\frac{\pi }{4}\right)-U_x-\{\begin{array}{ll}{U}_y-X\left(\frac{\pi }{4}\right)& n\text{ is even},\\ Y\left(\frac{\pi }{4}\right)& n\text{ is odd},\end{array}$$

$$\begin{array}{rl}& Y(\pi /4)-ZZ\left(\left[\pi -\arccos \left(\frac{1}{3}\right)\right]/4\right)\\ & -X\left(\arccos \left(\frac{1}{3}\right)/4\right)-ZZ\left(\left[\pi -\arccos \left(\frac{1}{3}\right)\right]/4\right)\\ & -Y(\pi /4).\end{array}$$

Full text

I Introduction

Entanglement is one of the most intriguing features of quantum physics, which remains at the heart of applications such as quantum computation qu-com , quantum teleportation qu-tele , and quantum cryptography CrApp ; CrPhoton . Meanwhile, it plays a crucial role in a variety of phenomena, e.g., the fractional quantum Hall effect qu-ha and quantum phase transitions qu-tr . Therefore, preparation of the entangled states is of significance in many-particle physics. It is well-known that there are two important and different types of entangled states, which cannot be transformed into each other under local operations and classical communication LOCC , i.e., the Greenberger-Horne-Zeilinger (GHZ) state GHZstate

[TABLE]

and the $W$ state Wstate

[TABLE]

The study of these two states has attracted much interest Property , and different methods to prepare these two states have been proposed recently, e.g. via dissipative preparation of entangled states 11 ; 12 .

One of the most commonly used methods to prepare entangled states is via quantum circuits, which has been implemented in experiments on various kinds of quantum systems, such as nuclear magnetic resonance (NMR) QC-NMR , ion traps QC-ion and cavity QED QC-QED . However, this method usually requires individual addressability of qubits, which makes the experiment difficult in large systems with many qubits. To overcome this problem, one of the practical ways is to employ global controls that act on all the spins.

The Ising spin model is one of the most ubiquitous in many physical systems such as optical lattices jack , NMR systems Vandersypen , ion traps porras and polar molecules Micheli . It also plays an important role in both condensed-matter physics ising-cmp and quantum information theory ising-qit . Recently, the study of global control methods to generate entangled states has attracted a lot of attention in the Ising-type spin-spin interaction systems global control GHZ ; cluster ; ben ; gao . The global control method relaxes the demanding experimental requirement to address and operate a single spin.

In this paper, we study how to prepare GHZ and $W$ states on a long-range Ising spin model via global controls. The rest of this paper is organized as follows. In Sec. II we introduce the studied spin model under global controls and our quantum control problem. In Sec. III we present the general solutions to obtain the GHZ state. In Sec. IV, we establish a standard way to generate the $W$ state. Finally, a brief summary with a discussion is presented in Sec. V.

II long-range Ising spin model and control problem

An $n$-spin long-range Ising spin model has the following Hamiltonian:

[TABLE]

where $\sigma_{k}^{\alpha}\left(\alpha=x,y,z\right)$ are the spin-1/2 Pauli matrices acting on the $k$th qubit. Here, we adopt the following two available global controls

[TABLE]

i.e., applying transverse magnetic fields on all of the spins along the $x$ or $y$ direction.

The total time-dependent Hamiltonian is

[TABLE]

where $f(t)$ and $g(t)$ are time-dependent functions. We assume here that the global fields are strong enough that the pulses can be regarded as the form of a $\delta$ function. This condition is easily satisfied in some physical systems such as NMR. This implies that the propagator can be regarded as a product of the time evolution operators under Hamiltonians in Eqs. (3) and (4):

[TABLE]

Our task is to find a pulse sequence with a set of suitable parameters $\{M,\tau_{m},\beta_{m}^{x},\beta_{m}^{y}\}$ to maximize the fidelity

[TABLE]

Here, $|T\rangle=|\text{GHZ}\rangle$ or $|W\rangle$. In the following sections, let

[TABLE]

III GHZ state

In this section, we discuss the cases in odd- and even-spin systems, respectively.

III.1 Odd-spin systems

Inspired by the solution in three-spin system gao , we found the general sequence for preparing the GHZ state in odd-spin systems as

[TABLE]

We will prove this in the following. A $Y(\pi/4)$ pulse applied to the initial state $|0\cdots 00\rangle$ yields a uniform superposition:

[TABLE]

Here, $k$ is a binary number indicating the index of the basis vector. Any vector from the computational basis is an eigenvector of the interaction Hamiltonian:

[TABLE]

where $s(k)$ is number of spin-up spins in the state $|k\rangle$. Therefore,

[TABLE]

Applying $X(-\pi/4)$ on the GHZ state, we obtain

[TABLE]

Except for the global phase, the ratio of the coefficient of the state $|k\rangle$ in Eq. (13) to the coefficient of Eq. (12) is a function of $s(k)$:

[TABLE]

Hence,

[TABLE]

Note that $i^{[1-s(k)][n-s(k)]}/i^{[s(k)+1][2-n+s(k)]}=i^{2[n-2s(k)-1]}=1$ when $n$ is odd. Therefore, we have

[TABLE]

which implies that $f(s(k))$ is independent of $s(k)$ and thus the state in Eq. (13) is equal to that of Eq. (12) up to an overall phase. Therefore the sequence $Y(\pi/4)-ZZ(\pi/4)-X(\pi/4)$ creates the GHZ state in odd systems.

III.2 Even-spin systems

For even systems, preparing GHZ states can be regarded as creating the total spin coherence between the state $|00\cdots 0\rangle$ and $|11\cdots 1\rangle$ with a specific phase difference slichter . The $ZZ$ coupling alone is zero coherence, which cannot change the coherence order of a state. However, the $ZZ$ coupling can be transformed into $H_{xx}=\sum_{k<m}^{n}\sigma_{k}^{x}\sigma_{m}^{x}=e^{-iH_{y}\pi/4}H_{zz}e^{iH_{y}\pi/4}$, which is a sum of zero-coherence and double-coherence terms. With this Hamiltonian, it is more convenient to expand the initial state as a sum of tensor products of $\sigma^{x}$’s eigenstates $\left\{|+\rangle,|-\rangle\right\}$. Let $|+\rangle=|0\rangle_{x}$ and $|-\rangle=|1\rangle_{x}$, and then $|k\rangle_{x}$ is the binary representation in the $\sigma^{x}$’s basis. The initial state can be expanded as

[TABLE]

Like for Eq. (12), we have

[TABLE]

where $s_{x}(k)$ is the number of plus signs in state $|k\rangle_{x}$ and $XX\left(t\right)=e^{-iH_{xx}t}$ is the evolution operator under the Hamiltonian $H_{xx}$. With the same method used in Eq. (16), it can be proved that $i^{s_{x}(k)[n-s_{x}(k)]}/\left[1+(-1)^{n-s_{x}(k)}i^{n+1}\right]$ is independent of $s_{x}(k)$ when $n$ is even. Combining this result with Eqs. (17) and (18), regardless of the global phase, we have

[TABLE]

This is the GHZ-type state. In order to obtain the exact state $|\text{GHZ}\rangle$ in Eq. (1), one needs to further apply a global $Z$ rotation:

[TABLE]

In addition, as we discussed, the operator $XX(t)$ can be implemented as $Y(\pi/4)-ZZ(t)-Y(-\pi/4)$. Accordingly the sequence to prepare the GHZ state in even-spin systems is

[TABLE]

IV $W$ state

In this section, we study how to generate the $W$ state. Since the Hamiltonian (3) and the global control Hamiltonians (4) are symmetric with respect to permutation of the qubits, they can be block diagonalized under the following symmetry-adapted basis:

[TABLE]

Here, $C^{m}_{n}$ is the number of $m$ combinations over $n$ elements; the sum is taken over all the computational basis states $|k\rangle$ with $s(k)$ spin-up spins, $m=1,\cdots,n+1$. The initial state $|0\cdots 00\rangle$ and the target state $|W\rangle$ are within subspace spanned by basis (22), implying that the calculation can be analyzed in this subspace. Moreover, $H_{zz}$ and $H_{x}$ commute with the “$X$-parity” operator $X=\prod_{i=1}^{n}\sigma_{i}^{x}$; therefore, the representations of $H_{zz}$ and $H_{x}$ in the symmetry-adapted basis can be further simultaneously block diagonalized in the two eigenspaces of $X$ with eigenvalues $\pm 1$:

[TABLE]

where $c$ is the normalizing constant and $k=1,2,\cdots,\left\lfloor n/2\right\rfloor+1$ ($\left\lfloor x\right\rfloor$ is the largest integer less than or equal to $x$). Likewise, the representations of $H_{zz}$ and $H_{y}$ in the symmetry-adapted basis can be block-diagonalized under the eigenstates of the “$Y$-parity” operator $Y=\prod_{i=1}^{n}\sigma_{i}^{y}$:

[TABLE]

where $c$ is the normalizing constant and $k=1,2,\cdots,\left\lfloor n/2\right\rfloor+1$.

Using different operations ($X,Y$, and $ZZ$ operations), we can shuttle between these subspaces. According to Eqs. (22), (23), and (24), the following properties can be obtained when $n\geqslant 2$ (the corresponding transformations are represented in Fig. 1): (1) When $n$ is odd [Fig. 1 (a)]

(i) $Y(\pi/4)$ rotates $\left|00\cdots 0\right\rangle$ (point $I$) into the subspace $\mathcal{X}_{+}$ (point $A$),

[TABLE]

and (ii) $Y(-\pi/4)$ rotate the $W$ state (point $J$) into the subspace $\mathcal{X}_{+}$ (point $B$),

[TABLE]

(2) When $n$ is even. [Fig. 1 (b)], (i) $Y(\pi/4)$ rotates $\left|00\cdots 0\right\rangle$ (point $I$) into the subspace $\mathcal{X}_{+}$ (point $A$),

[TABLE]

and (ii) $X(-\pi/4)$ rotates the $W$ state (point $J$) into the subspace $\mathcal{Y}_{-}$ (point $D$),

[TABLE]

We also note here that the subspaces $\mathcal{X}_{\pm}$ overlap the subspaces $\mathcal{Y}_{\pm}$ when $n$ is even. Utilizing these properties, a routine way to generate the $W$ state can be summarized as:

[TABLE]

where $U_{x}=\prod_{i}{ZZ(\theta_{zzi})X(\theta_{xi})}$ and $U_{y}=\prod_{i}{ZZ(\theta_{zzi})Y(\theta_{yi})}$. The number of operations needed here is uniformly bounded shandeluo . Since the dimension of the subspace increases linearly with the size $n$ of the system, we expect that the number of operations required also has a linear dependence on $n$. In addition, the system is pure state controllable in the subspaces $\mathcal{X}_{\pm}$ and $\mathcal{Y}_{\pm}$, respectively; a solution through the decomposition (29) thus always exists. The proof is included in the Appendix. Therefore, a solution with high fidelity can always be found by setting the error threshold sufficiently small in principle.

As examples, we show how to prepare the $W$ state on three- and four-spin Ising models. Here, we denote the matrix representations of different Hamiltonians in different subspaces by $H_{\alpha}^{\Xi}$, where $\alpha=x$, $y$ or $zz$ and $\Xi=\mathcal{X}_{\pm}$ or $\mathcal{Y}_{\pm}$ represents the corresponding subspaces. We use the right-hand side of Eqs. (25) and (27) as the initial states, and the right-hand side of Eqs. (26) and (28) as the target states for odd and even $n$, respectively.

For a three-spin case, following Eq. (29), we can obtain the $W$ state with the sequence:

[TABLE]

which can be represented on the Bloch sphere in Fig. 2 (a). The detailed analysis is identical to that in gao .

As for the four-spin case, according to Eq. (27),

[TABLE]

We should find the sequence $U_{x}$ to rotate this state into the subspace $\mathcal{Y}_{-}$. Since $\left|x^{+}_{2}\right\rangle$ is the only common basis vector for the subspaces $\mathcal{X}_{+}$ and $\mathcal{Y}_{-}$, suppose $U_{x}=X\left(\theta_{1}\right)-ZZ\left(\theta_{2}\right)-X\left(\theta_{3}\right)$, and let

[TABLE]

where we have $\theta_{2}=\pi/4$, $\theta_{1}=0$, and $\theta_{3}=\pi/16$ regarding simplicity.

According to Eq. (28), $X\left(-\pi/4\right)\left|W\right\rangle=\left(i\left|y^{-}_{1}\right\rangle-\left|y^{-}_{2}\right\rangle\right)/\sqrt{2}$; we next need to find $U_{y}$, which rotates $\left|y_{1}^{-}\right\rangle$ to $(i\left|y^{-}_{1}\right\rangle-\left|y^{-}_{2}\right\rangle)/\sqrt{2}$. Since

[TABLE]

and

[TABLE]

we have

[TABLE]

when $\theta_{1}=\pi/12$ and $\theta_{2}=3\pi/8$, we obtain $\left(i\left|y^{-}_{1}\right\rangle-\left|y^{-}_{2}\right\rangle\right)/\sqrt{2}$. The process , as illustrated in Fig. 2 (b), can be expressed as

[TABLE]

Therefore, we achieve the target $W$ state with the following control sequence:

[TABLE]

For large systems, it is difficult to find the analytic solutions. However, since the dimension of the subspaces increases linearly with the system size, with the numerical optimization algorithm it is possible to generate the $W$ states via the global-control sequence in principle. Applying the approach to a numerical search, one can find the solution in much larger systems (Fig. 3). We tried 100 numerical searches with different initial parameters for each $n$ and different numbers of parameters. Blue squares in Fig. 3 represent the minimum number of parameters in 100 sets of data to keep the fidelity above 0.999 by the numerical optimization algorithm. The red line represents the linear fitting. It is clear that the number of parameters increases almost linearly with the system size; that is, using only $O(n)$ operations we can prepare the $W$ states with a high fidelity in a long-range Ising model using only global controls.

V Conclusion

In summary, we presented a universal method to prepare the GHZ state and a standard procedure to prepare the $W$ state on a long-range Ising model using only global controls. Actually the solutions to preparing the GHZ state also obey the rules of the partition of subspaces described for the $W$ state: the GHZ state is within the subspace $\mathcal{X}_{+}$. One can follow a procedure similar to that used in finding solutions for the $W$ state to obtain different solutions for preparing the GHZ state if needed. Remarkably, the dimension of the irreducible subspace we employed increases only linearly with the size of the system due to the symmetry, so the scheme for preparing the $W$ state is also feasible in larger spin systems. In addition, it is also applicable to a wide range of physical implementations, and will contribute to quantum control for implementing quantum information processing in the future.

ACKNOWLEDGMENTS

We thank Prof. D. D’Alessandro, I. Hinks, and H. Le for helpful discussions. This work is supported by the National Key Basic Research Program under Grant No. 2013CB921800 and No. 2014CB848700, the National Science Fund for Distinguished Young Scholars under Grant No. 11425523, the National Natural Science Foundation of China under Grant No. 11375167 and No. 11661161018, the Strategic Priority Research Program (B) of the CAS under Grant No. XDB01030400.

Appendix A Controllability and irreducibility

In this appendix, we give the matrix representations of the relevant Hamiltonians in the corresponding subspaces and proofs for controllability of the system and irreducibility of the subspaces $\mathcal{X}_{\pm}$ and $\mathcal{Y}_{\pm}$.

Let $H_{x}^{s}$, $H_{y}^{s}$ and $H_{zz}^{s}$ denote the matrix representations of Hamiltonians $H_{y}$, $H_{y}$ and $H_{zz}$ in the basis (22), respectively, and $\left\{a_{ij}\right\}$ and $\left\{b_{ij}\right\}$ denote the matrix elements of $H_{x}^{s}$ and $H_{y}^{s}$, respectively. All the elements of $H_{x}^{s}$ and $H_{y}^{s}$ are zero except those on the two lines next to the diagonal line:

[TABLE]

for $k=1,\cdots,n$. The matrix representation of $H_{zz}$ is diagonal in this basis:

[TABLE]

where $\lambda_{k}=2\left(k-1-\frac{n}{2}\right)^{2}-\frac{n}{2}$ for $k=1,\cdots,n+1$.

From Fig. 1 we know that it is only necessary to know the matrix representations of $H_{x}$ in the subspace $\mathcal{X}_{+}$ when $n$ is either even or odd, $H_{y}$ in the subspace $\mathcal{Y}_{-}$ when $n$ is even and $H_{zz}$. According to Eqs. (38) and (39), we can write down the specific matrix representations of the Hamiltonians in the corresponding subspaces for $n\geqslant 2$.

When $n$ is even, the matrix elements of $H_{x}$, $\left\{a_{k}^{+}\right\}$, in the subspace $\mathcal{X}_{+}$ are zero except those on the lines next to the diagonal line:

[TABLE]

When $n$ is odd, the matrix elements of $H_{x}$, $\left\{\tilde{a}_{k}^{+}\right\}$ (we use a tilde to denote the case when $n$ is odd), in the subspace $\mathcal{X}_{+}$ are zero except those on the lines next to the diagonal line and the last element at the bottom right corner:

[TABLE]

When $n$ is even, the matrix elements of $H_{y}$, $\left\{b_{k}^{-}\right\}$, in the subspace $\mathcal{Y}_{-}$ are zero except those on the lines next to the diagonal line:

[TABLE]

where $1\leqslant k\leqslant\frac{n}{2}-1$. Likewise, $H_{zz}$ is diagonal in $\mathcal{X}_{\pm}$ and $\mathcal{Y}_{\pm}$, and the matrix representations in these subspaces are:

[TABLE]

where ${\Xi}=\mathcal{X}_{+},\mathcal{X}_{-},\mathcal{Y}_{+},\mathcal{Y}_{-}$. All the diagonal elements are different from each other:

[TABLE]

When $n$ is even, $H_{zz}^{\mathcal{X}_{+}}$ and $H_{zz}^{\mathcal{Y}_{+}}$ include an additional diagonal element, $\lambda_{n/2+1}=-n/2$. Both the analytic solutions and the numerical search rely on the matrix representations of these Hamiltonians.

We are concerned with pure-state controllability of the system. One can prove this by looking at the connectivity graph where the vertices represent the eigenstates of $H_{zz}^{\mathcal{X}_{+}}$ and the edges connect vertices $j$ and $k$ if the element $H_{x}^{\mathcal{X}_{+}}(j,k)$ is nonzero shandeluo ; turinici . For example, in the subspace $\mathcal{X}_{+}$ (Fig. 4), since the element $H_{x}^{\mathcal{X}_{+}}(k,k+1)$ is nonzero, it transfers magnitude between the $k$th and $(k+1)$th eigenstates of $H_{zz}^{\mathcal{X}_{+}}$. Therefore, in the graph there is an edge connecting two adjacent vertices. Moreover, the difference between the adjacent eigenvalues $\lambda_{k}$ and $\lambda_{k+1}$ of $H_{zz}^{\mathcal{X}_{+}}$ is $2(n+1)-4k$, which means there is no degenerate transition. The same result holds for subspaces $\mathcal{X}_{-}$ and $\mathcal{Y}_{\pm}$. Therefore, the system is pure state controllable in the subspaces $\mathcal{X}_{\pm}$ and $\mathcal{Y}_{\pm}$ according to Theorem 3.7.1 in shandeluo . This proof also guarantees the existence of the solutions in Eq. (29).

Next, we prove the subspaces $\mathcal{X}_{\pm}$ and $\mathcal{Y}_{\pm}$ are irreducible. Take $\mathcal{X}_{+}$ with odd $n$. In order to partition $\mathcal{X}_{+}$ further, one must find another operator $A$ such that $\left[H_{x},A\right]=\left[H_{zz},A\right]=0$. Since the eigenstates of $H_{zz}$ are nondegenerate according to (44), the matrix representation of $A$ must be diagonal:

[TABLE]

Since $\left[H_{x},A\right]\left|s_{k}\right\rangle=\left(H_{x}A-AH_{x}\right)\left|s_{k}\right\rangle=0$, we can always obtain $\lambda^{A}_{1}=\lambda^{A}_{2}=\cdots=\lambda^{A}_{\left\lfloor\frac{n+1}{2}\right\rfloor}$ according to (40) and (41). Therefore, $A$ is a trivial matrix (a product of a constant and a unit matrix). Thus the subspace cannot be reduced any more. Similar results can be proved for the subspaces $\mathcal{X}_{-}$ and $\mathcal{Y}_{\pm}$ and for the case with even $n$.