Exact solution of a non-Hermitian $\mathscr{PT}$-symmetric spin chain

Quantum dynamics of physical systems is governed by Hamiltonians whose eigenvalues give the energies of the states in the system. In a conventional isolated quantum system, the Hamiltonian is required to be Hermitian so that its eigenvalues are real and the time evolution is unitary, thereby guaranteeing the conservation of probabilities. This picture of quantum mechanics is quite idealized as no real quantum system exists in isolation. One way to incorporate the effect of the environment in a quantum system is to consider effective non-Hermitian Hamiltonians which naturally arises when the full Hermitian Hamiltonian describing the system and its environment is projected on the space of the system of interest [1]. These Hamiltonians generally have complex eigenspectra and they describe out-of-equilibrium systems which may gain or lose energy to the environment. The imaginary part of the energy is related to the lifetime of a state of the system.

An important class of non-Hermitian Hamiltonians, known as space-time reflection symmetric, or $\mathscr{PT}$ -symmetric Hamiltonians, was introduced in a seminal paper by Bender and Boettcher [2]. These Hamiltonians describe the borderline systems between the open and closed systems in the sense that they are not isolated like the closed systems but the effect of the environment on them is severely constrained as the loss and gain in the system are balanced out such that there is no net loss or gain in the system. Due to this constraint, $\mathscr{PT}$-symmetric Hamiltonians can have real energy levels. $\mathscr{PT}$-symmetric quantum mechanics and quantum field theories have found applications in wide range of areas including quantum optics [3–6], and condensed matter systems [7–11].

The eigenvalues of such Hamiltonians are either real, corresponding to $\mathscr{PT}$-symmetric eigenvectors, or appear as complex conjugate pairs when the corresponding eigenvectors break the $\mathscr{PT}$-symmetry. This follows from the anti-linearity of the $\mathscr{PT}$-reflection symmetry, which implies that although $[\mathscr{PT},H] = 0$, the Hamiltonian and the $\mathscr{PT}$ operator may not be simultaneously diagonalizable. Consider first the eigenvalue equation $H\vert \psi\rangle = \lambda\vert \psi\rangle$ for a state $\vert \psi\rangle$ that is symmetric, namely $\mathscr{PT}\vert \psi\rangle = \vert \psi\rangle$, then $0 = [\mathscr{PT},H]\vert \psi\rangle = (\lambda^*-\lambda)\vert \psi\rangle$ implies that the eigenvalue is real. However when the symmetry is spontaneously broken, namely $\mathscr{PT}\vert \psi\rangle = \vert \psi^{\prime}\rangle\neq \vert \psi\rangle$, then $0 = [\mathscr{PT},H]\vert \psi\rangle$ implies $H\vert \psi^{\prime}\rangle = \lambda^*\vert \psi^{\prime}\rangle$ indicating that if there exist a state $\vert \psi\rangle$ with complex energy λ, then there is another state $\vert \psi^{\prime}\rangle$ with complex conjugate energy λ. The converse proof for both of the cases also holds ⁴ . The phase transitions between the unbroken phase and the spontaneously broken $\mathscr{PT}$-phase have been experimentally realized in different systems [14–16].

Computing analytically the energy levels of the $\mathscr{PT}$-symmetric Hamiltonians is a difficult task. Even simple single-particle quantum mechanics with $\mathscr{PT}$-symmetric potentials like ix³ potential [17] or $\mathscr{PT}$- symmetric square well potential [18] can only be solved perturbatively. The proof of the reality of energy levels of certain $\mathscr{PT}$-symmetric quantum mechanics uses the conjecture of correspondence between ordinary differential equation and integrable models [19]. However, general proof of the conjecture is not known. Even though $\mathscr{PT}$-symmetric spin chains have been studied focusing on various aspects [20–22], and $\mathscr{PT}$-symmetric deformation of integrable models [23] and some other $\mathscr{PT}$-symmetric many-body systems [24–26] have been studied both numerically and analytically, there has been no complete solution showing that all the excited states can have either real energies or they form pairs having complex conjugate energies. Some non-interacting exactly solvable many-body models have been analytically studied before [27, 28].

In this paper, we introduce a many-body $\mathscr{PT}$-symmetric model which is constructed by applying complex magnetic fields to the edges of a Heisenberg spin chain. The Hamiltonian is given by

$$\begin{align} \mathcal{H} = \sum_{j = 1}^{N-1}\sum_{\alpha = \{x,y,z \}}\sigma_j^\alpha\cdot\sigma_{j+1}^\alpha+h_1\sigma_1^z+h_N\sigma_N^z, \end{align} \tag{ 1 }$$

where h₁ and h_N are the magnetic fields acting at the first and last site in the z-direction, $h_1 = h-\textrm i\gamma, h_N = h+\textrm i\gamma$. Equivalently, upon the Jordan–Wigner transformation, we obtain a fermionic version of the Hamiltonian

$$\begin{align} \mathcal{H}_f = \frac{J}{2}\sum_{i = 1}^{N-1} f_{j+1}^\dagger f_j+f_j^\dagger f_{j+1}+J\left(n_j-\frac12\right)\left(n_{j+1}-\frac12\right)+h_1n_1+h_Nn_N, \end{align} \tag{ 2 }$$

which describes a fermionic chain with nearest-neighbor interactions coupled to baths at its edges, as shown in the figure 1. This model may be experimentally realized in cold atom systems or optical systems with controlled loss and gain [29–31].

Zoom In Zoom Out Reset image size

The exact results presented in this work might be useful for providing the benchmark for various numerical techniques developed for non-Hermitian quantum mechanics [32, 33].

This system with real values of the boundary fields was solved using Bethe Ansatz [34–36] and recently the complete phase diagram and the structure of the Hilbert space were obtained [37].

In this work we consider the case where γ ≠ 0 and solve the system exactly using Bethe Ansatz. We find that the system exhibits several phases which can be broadly classified into A type and B type phases. In the B type phases the $\mathscr{PT}$-symmetry remains unbroken, whereas A type phases exhibit both $\mathscr{PT}$-symmetry broken and unbroken sectors. The existence of two sectors in the A phases is related to the existence of bound state solutions at the edges with complex energies, where the energy of the bound state at the left edge is a complex conjugate of the energy of the bound state at the right edge. Furthermore, we find that the ground state exhibited by the system depends on the orientation of the boundary fields, hence, both the A and B type phases can be further divided into sub-phases depending on the ground state exhibited by the system. Notice that γ → 0 corresponds to both the edge fields taking equal and real values, which corresponds to a special case of [37]. In this limit, we recover the results in [37] corresponding to this special case of equal fields: we find that when the values of these fields are greater than the critical value $h = 2J$, there are stable edge modes, and if the fields are below the critical value, there are exist no edge modes. In addition to this, the energies of the bound states of [37] are reproduced in this limit where γ → 0.

The paper is organized as follows: we summarize the Bethe Ansatz results in the section 2 and present the properties of the model in the section 3. We present the details of the solution in sections 4–6, and finally discuss the results in section 7.

In this section, we summarize the results obtained via the Bethe Ansatz, relegating the details of the solution to the later sections. There exist four phases that correspond to different ranges of the parameter $\xi = \Re(1/h_1) = \Re(1/h_2)$. The phases B₁ and B₂ correspond to the ranges $\xi\gt\frac{1}{2}$ and $\xi\lt-\frac{1}{2}$ respectively, whereas the phases A₁ and A₂ correspond to the ranges $0\lt\xi\lt\frac{1}{2}$ and $-\frac{1}{2}\lt\xi\lt0$ respectively. Below we summarize the results in each of these phases separately for odd and even numbers of sites chain.

2.1. Ground state

2.2. Excitations

2.2.1.  B phases.

In the B phases there exist no bound states at both edges. One can build up excitations on top of the ground state by adding an even number of spinons, bulk strings, and quartets, where the energy of all these excitations is real.

2.2.2.  A phases.

In the A phases, similar to the B phases, one can build up excitations on top of the ground state and obtain excitations which have real energies.

There exist two bound state solutions corresponding to the left and right edges. The bound state energy of the left and right boundaries are complex conjugates of each other and are given by the expressions (56), (57) and (61), (62). In order to add a bound state to either the left or the right edge, one needs to add a spinon, whose energy is given by (55). The lowest energy of the spinon corresponds to rapidity $\theta\rightarrow\infty$. Hence the state which contains a bound state at the left or the right edge has complex energy and forms a continuous branch parameterized by the rapidity θ of the spinon.

One can build excitations on top of each of these states containing one bound state, by adding an even number of spinons, wide boundary strings, bulk strings, and quartets [37]. All these excitations have complex energies due to the presence of a single bound state either at the left or right edge.

One can add both the bound states to the ground state which does not require one to add a spinon. Since the bound states have energies that are complex conjugates of each other, this state has real energy. One can construct excitations on top of this state by adding an even number of spinons, bulk strings, wide boundary strings and quartets. All of the excited states built this way have real energy.

Hence, in the A phases, there exist excited states built on top of the state with no bound states or on top of the state which contains two bound states, and obtain excitations whose energies are real. These states correspond to $\mathscr{PT}$-symmetry unbroken sector. Likewise, there exist excited states built on top of the states which contain a bound state either at the left or the right edge, and one obtains excitations which have complex energies, and hence they correspond to $\mathscr{PT}$-symmetry broken sector. The results of the odd and even number of sites chain are summarized in the tables 1 and 2 respectively.

Table 1. Energies and $\mathscr{PT}$-symmetry of the ground state and the lowest energy states containing bound states for an odd number of sites for various ranges of boundary parameter $\xi = \frac{h}{\gamma ^2+h^2}$. The B phases only consist of symmetry unbroken states with real energy eigenvalues whereas the A phases consist of both broken and unbroken $\mathscr{PT}$-symmetric states.

Phase	Parametric range	State	Spin	Energy	$\mathscr{PT}$-symmetry
B2	$\xi\lt-\frac{1}{2}$	$\left|\frac{1}{2}\right\rangle$	$\frac{1}{2}$	E0	Unbroken
A2	$-\frac{1}{2}\lt\xi\lt0$	$\left|\frac{1}{2}\right\rangle$	$\frac{1}{2}$	E0	Unbroken
		$\left|-\frac{1}{2}\right\rangle_{L/R}$	$-\frac{1}{2}$	$E_0+m+m^*$
		$\left|\pm\frac{1}{2}\right\rangle_{\theta,L}$	$\pm\frac{1}{2}$	$E_0+E_\theta+m$	Broken
		$\left|\pm\frac{1}{2}\right\rangle_{\theta,R}$	$\pm\frac{1}{2}$	$E_0+E_\theta+m^*$
A1	$0\lt\xi\lt\frac{1}{2}$	$\left|-\frac{1}{2}\right\rangle$	$-\frac{1}{2}$	E0	Unbroken
		$\left|\frac{1}{2}\right\rangle_{L/R}$	$\frac{1}{2}$	$E_0+m+m^*$
		$\left|\pm\frac{1}{2}\right\rangle_{\theta,L}$	$\pm\frac{1}{2}$	$E_0+E_\theta+m$	Broken
		$\left|\pm\frac{1}{2}\right\rangle_{\theta,R}$	$\pm\frac{1}{2}$	$E_0+E_\theta+m^*$
B1	$\xi\gt\frac{1}{2}$	$\left|-\frac{1}{2}\right\rangle$	$-\frac{1}{2}$	E0	Unbroken

Table 2. Energies and $\mathscr{PT}$-symmetry of the ground states and the lowest energy states with bound states for even number of sites for various ranges of boundary parameter $\xi = \frac{h}{\gamma ^2+h^2}$. Similar to the odd number of sites chain, the B phases only consist of symmetry unbroken states with real energy eigenvalues whereas the A phases consist of both broken and unbroken $\mathscr{PT}$-symmetric states.

Phase	Parametric range	State	Spin	Energy	$\mathscr{PT}$-symmetry
B2	$\xi\lt-\frac{1}{2}$	$\vert 0\rangle_\theta$	0	$E_0+E_\theta$	Unbroken
		$\vert 1\rangle_\theta$	1	$E_0+E_\theta$
A2	$-\frac{1}{2}\lt\xi\lt0$	$\vert 1\rangle_\theta$	1	$E_0+E_\theta$	Unbroken
		$\vert 0\rangle_\theta$	0	$E_0+E_\theta$
		$\vert 1\rangle_{\theta,L,R}$	1	$E_0+E_\theta+m+m^*$
		$\vert 0\rangle_{\theta,L,R}$	0	$E_0+E_\theta+m^*+m$
		$\vert 0\rangle_L$	0	$E_0+m$	Broken
		$\vert 0\rangle_R$	0	$E_0+m^*$
A1	$0\lt\xi\lt\frac{1}{2}$	$\vert -1\rangle_\theta$	−1	$E_0+E_\theta$	Unbroken
		$\vert 0\rangle_\theta$	0	$E_0+E_\theta$
		$\vert -1\rangle_{\theta,L,R}$	−1	$E_0+E_\theta+m+m^*$
		$\vert 0\rangle_{\theta,L,R}$	0	$E_0+E_\theta+m^*+m$
		$\vert 0\rangle_L$	0	$E_0+m$	Broken
		$\vert 0\rangle_R$	0	$E_0+m^*$
B1	$\xi\gt\frac{1}{2}$	$\vert -1\rangle_\theta$	−1	$E_0+E_\theta$	Unbroken
		$\vert 0\rangle_\theta$	0	$E_0+E_\theta$

The isotropic Heisenberg $XXX_{\frac{1}{2}}$ Hamiltonian is a linear endomorphism in the product space $\displaystyle\otimes_{j = 1}^N h_j$ where the local Hilbert space of each quantum spin variable is the two-dimensional complex vector space $h_n = \mathbb{C}^2$

$$\begin{align} \mathcal{H} = \sum_{j = 1}^{N}\sum_{\alpha = \{x,y,z \}}\sigma_j^\alpha\cdot\sigma_{j+1}^\alpha, \end{align} \tag{ 3 }$$

where $\sigma_j^\alpha$ are the Pauli matrices acting on vector space h_j . Imposing the periodic boundary condition, $\sigma_1^\alpha = \sigma_{N+1}^\alpha$, Bethe solved the eigenvalue problem $H\psi = E\psi$ exactly using a method now called Bethe Ansatz method [38]. Alcaraz et al solved the model with open boundary conditions applying real boundary fields [34]. Here we consider the complex deformation as mentioned in Hamiltonian equation (1) by allowing the boundary fields to be complex i.e. $h_{1/N}\in\mathbb{C}$. Our strategy is to use the Algebraic Bethe Ansatz along with Cherednik–Sklyanin reflection algebra to diagonalize the Hamiltonian [35]. The complex boundary terms lower the SU(2) symmetry to U(1) symmetry $[\mathcal{H},\sum_i S_i^z] = 0$ . Thus, the total z-component of the spin M is still a good quantum number for the model. The boundary terms, however, break the $\mathbb{Z}_2$ flip symmetry. The bulk is symmetric under the $\mathbb{Z}_2$ flip, but the entire system with boundary remains invariant only under the following transformation

$$\begin{align} \mathcal{H}\to\prod_{i = 1}^N \sigma_i^x\mathcal{H}\sigma_i^x, \quad h_1\to-h_1 \quad \text{and}\quad h_N\to-h_N. \end{align} \tag{ 4 }$$

We wish to solve the eigenvalue problem for the Hamiltonian given by equation (1) $\mathcal{H}\in \text{End}(\stackrel{N}{\otimes}h_n)$ with complex boundary fields using Algebraic Bethe Ansatz.

In $V_a\otimes h_n ~\exists$ a L-matrix

$$\begin{align} L_{a,n}(\lambda) = \left(\lambda-\theta_n\right)I_{a,n}+P_{a,n} , \end{align} \tag{ 5 }$$

where $V_a = \mathbb{C}^2$ is an auxiliary space and θ_n is the site-dependent inhomogeneity parameter. The L-matrix satisfies the fundamental commutation relation

$$\begin{align} R_{a,b}(\lambda-\mu)L_{n,a}(\lambda)L_{n,b}(\mu) = L_{n,b}(\mu)L_{n,a}(\lambda)R_{a,b}(\lambda-\mu) , \end{align} \tag{ 6 }$$

where the R-matrix acting as an intertwiner $\in V_a\otimes V_b$ given by

$$\begin{align} R_{a,b}(\lambda) = \lambda I_{a,b}+ P_{a,b} \end{align} \tag{ 7 }$$

is a matrix that solves the Yang–Baxter equation

$$\begin{align} R_{12}(\lambda-\lambda^{\prime})R_{13}(\lambda)R_{23}(\lambda^{\prime}) = R_{23}(\lambda^{\prime})R_{13}(\lambda)R_{12}(\lambda-\lambda^{\prime}), \end{align} \tag{ 8 }$$

and satisfies the unitary condition

$$\begin{align} R_{12}(\lambda)R_{21}(-\lambda)\propto \mathbb{I}. \end{align} \tag{ 9 }$$

To describe the open-boundary condition, apart from the RLL relation equation (6), the R-matrix also has to satisfy Sklyanin's reflection algebra given by

$$\begin{align} R_{1,2} (\lambda-\mu) K_{1}^{-}(\lambda) R_{2,1}(\lambda+\mu) K_{2}^{-}(\mu) = K_{2}^{-}(\mu) R_{1,2}(\lambda+\mu) K_{1}^{-}(\lambda) R_{2,1}(\lambda-\mu). \end{align} \tag{ 10 }$$

The diagonal $K^-$ matrix that satisfies the reflection equation (10) is [35, 39]

$$\begin{align} K^-(\lambda) = \begin{pmatrix} \xi_-+\lambda & 0 \\ 0 & \xi_-\lambda \end{pmatrix}, \end{align} \tag{ 11 }$$

where $\xi_-$ is a free complex parameter.

There exists a dual reflection algebra given by

$$\begin{align} & R_{1,2}(-\lambda+\mu) K_{1}^{+}(\lambda) R_{2,1}(-\lambda-\mu-2 ) K_{2}^{+}(\mu) \nonumber\\ & \quad = K_{2}^{+}(\mu) R_{1,2}(-\lambda-\mu-2 ) K_{1}^{+}(\lambda) R_{2,1}(-\lambda+\mu). \end{align} \tag{ 12 }$$

The $K^+$ matrix that satisfies the dual reflection equation (12) is [35, 39]

$$\begin{align} K^+(\lambda) = \begin{pmatrix} \xi_++\lambda+1 & 0 \\ 0 & \xi_+-\lambda-1 \end{pmatrix}, \end{align} \tag{ 13 }$$

where $\xi_+$ is a free parameter. One can define a transfer matrix shown below

$$\begin{align} t(\lambda) = \text{tr}_a\left(K_a^+(\lambda)T_{a}(\lambda) K_{a}^{-}(\lambda) \hat{T}_{a}(\lambda) \right), \end{align} \tag{ 14 }$$

where the trace is over the auxiliary space V_a and the operators $T_a(\lambda)$ and $\hat T_a(\lambda)$ are given by

$$\begin{align} T_a(\lambda)& = L_{a,N}(\lambda)L_{a,N-1}(\lambda)\cdots L_{a,2}(\lambda)L_{a,1}(\lambda), \end{align} \tag{ 15 }$$

$$\begin{align} \hat T_a(\lambda)& = L_{a,1}(\lambda)L_{a,2}(\lambda)\cdots L_{a,N-1}(\lambda)L_{a,N}(\lambda). \end{align} \tag{ 16 }$$

It is straightforward to prove the commutative property of the transfer matrix for any two distinct spectral parameters i.e. $[t(\lambda_i), t(\lambda_j)] = 0~\forall~ \lambda_i,\lambda_j\in\mathbb{C}$ and its eigenvalues are holomorphic functions of $\lambda\in\mathbb{C}$. Moreover, the commutative property guarantees that the transfer matrix behaves as a generating function for a family of infinite towers of commuting operators, among which is the Hamiltonian. Simultaneous diagonalization of the operators $t(\lambda_i)$ allows us to obtain complete spectral characteristics of the quantum integrable model. In this model, Hamiltonian is related to the first derivative of the logarithm of $t(\lambda)$ at λ = 0 and $\theta_i = 0~ \forall ~i$ i.e.

$$\begin{align} \mathcal{H} & = \left.\frac{\partial \ln t(\lambda)}{\partial \lambda}\right|_{\lambda = 0,\left\{\theta_{j} = 0\right\}}-N \nonumber\\ & = \sum_{j = 1}^{N-1}\sum_{\alpha = 1}^3 \sigma_{j}^\alpha \sigma_{j+1}^\alpha+\frac{1}{\xi_-} \sigma_{1}^{z}+\frac{1}{\xi_+}\sigma_{N}^{z}, \end{align} \tag{ 17 }$$

if one identifies the free complex parameters $\xi_\pm$ as $\xi_- = \frac{1}{h_1}$ and $\xi_+ = \frac{1}{h_N}$.

The eigenvalues $\Lambda(\lambda)$ of the transfer matrix equation (14) satisfy the Baxter's T − Q relation [40]

$$\begin{align} \begin{aligned} \Lambda(\lambda)& = \frac{2(\lambda+1)^{2 N+1}}{2 \lambda+1}(\lambda+\xi_-)\left( \lambda+\xi_+\right)\frac{Q(\lambda-1)}{Q(\lambda)} \\ &\quad+\frac{2 \lambda^{2 N+1}}{2\lambda+1}(\lambda+1-\xi_-)\left(\lambda+1-\xi_+\right) \frac{Q(\lambda+1)}{Q(\lambda)} , \end{aligned} \end{align} \tag{ 18 }$$

where the Q-function is given by

$$\begin{align} Q(\lambda) = \prod_{\ell = 1}^M (\lambda-\lambda_\ell)(\lambda+\lambda_\ell+1). \end{align} \tag{ 19 }$$

Regularity of the T–Q equation gives the BAEs corresponding to the reference state where all spins are up

$$\begin{align} \left(\frac{\lambda_j+1}{\lambda_j}\right)^{2N} \frac{(\lambda_j+\xi_-)(\lambda_j+\xi_+)}{(\lambda_j+1-\xi_-)(\lambda_j+1-\xi_+)}\prod_{\ell\neq j}^M\frac{(\lambda_j+\lambda_\ell)(\lambda_j-\lambda_\ell-1)}{(\lambda_j-\lambda_\ell+1)(\lambda_j+\lambda_\ell+2)} = 1. \end{align} \tag{ 20 }$$

Taking the derivative of the logarithm of the eigenvalue equation (18), we obtain

$$\begin{align} E = \frac{\partial\ln\Lambda(\lambda)}{\partial\lambda}\big\vert_{\lambda = 0}-N = \sum_{j = 1}^M \frac{2}{\lambda_j(\lambda_j+1)} +N-1+\frac{1}{\xi_-}+\frac{1}{\xi_+}. \end{align} \tag{ 21 }$$

Changing the variables $\lambda_j = \textrm i\mu_j-\frac{1}{2}$, the energy equation (21) becomes

$$\begin{align} E = -\sum_{j = 1}^M \frac{2}{\mu_j^2+\frac{1}{4}}+N-1+\frac{1}{\xi_-}+\frac{1}{\xi_+}. \end{align} \tag{ 22 }$$

Similarly, the Bethe Ansatz equation (20) can be written as

$$\begin{align} \left(\frac{\mu_{j}-\frac{\textrm i}{2}}{\mu_{j}+\frac{\textrm i}{2}}\right)^{2 N} \left(\frac{\mu_{j}+\textrm i\left(\frac{1}{2}-\xi_-\right)}{\mu_{j}-\textrm i\left(\frac{1}{2}-\xi_-\right)}\right)\left(\frac{\mu_{j}+\textrm i\left(\frac{1}{2}-\xi_+\right)}{\mu_{j}-\textrm i\left(\frac{1}{2}-\xi_+\right)}\right) = \prod_{j \neq \ell = 1}^{M}\left(\frac{\mu_{j}-\mu_{\ell}-\textrm i}{\mu_{j}-\mu_{\ell}+\textrm i}\right)\left(\frac{\mu_{j}+\mu_{\ell}-\textrm i}{\mu_{j}+\mu_{\ell}+\textrm i}\right).\nonumber\\ \end{align} \tag{ 23 }$$

From the Bethe equation equation (23), we can see that if $\mu_j\in\mathbb{C}$ is a solution, then −µ_j is also the solution. In other words, the distribution of roots of Bethe equations is symmetric about the origin.

Taking logarithm on both sides of equation (23) and using $\ln\left(\frac{i-z}{i+z} \right) = 2i\tan^{-1}(z)$, we get

$$\begin{align} & (2N+1)\tan^{-1}(2\mu_j)-\tan^{-1}\left(\frac{\mu_j}{\frac{1}{2}-\xi_-}\right)-\tan^{-1}\left(\frac{\mu_j}{\frac{1}{2}-\xi_+}\right)\nonumber\\ & \qquad = \sum_{\ell = 1}^M\left[\tan^{-1}(\mu_j-\mu_\ell)+\tan^{-1}(\mu_j+\mu_\ell) \right]+\pi I_j. \end{align} \tag{ 24 }$$

The last term comes from the fact that exp function over $\mathbb{C}$ is not injective and hence its inverse function ln is multivalued. Here $I_j \in \mathbb{Z}$ acts as the quantum number of the model.

To analyze equation (24) in the thermodynamic limit, we define the density of Bethe roots as

$$\begin{align} \rho(\mu_j) = \frac{1}{\mu_{j+1}-\mu_j}. \end{align} \tag{ 25 }$$

Such that we convert the sums over j in equations (22) and (24) into integral over µ as

$$\begin{align} E = -\int_\mathcal{C}\mathrm{d}\mu~ \rho(\mu)\frac{2}{\mu^2+\frac{1}{4}}+N-1+\frac{1}{\xi_-}+\frac{1}{\xi_+}, \end{align} \tag{ 26 }$$

and

$$\begin{align} & (2N+1)\tan^{-1}(2\mu_j)-\tan^{-1}\left(\frac{\mu_j}{\frac{1}{2}-\xi_-}\right)-\tan^{-1}\left(\frac{\mu_j}{\frac{1}{2}-\xi_+}\right)\nonumber\\ & \qquad = \int_\mathcal{C} \mathrm{d}\mu^{\prime} \rho(\mu^{\prime})\left[\tan^{-1}(\mu_j-\mu^{\prime})+\tan^{-1}(\mu_j+\mu^{\prime}) \right]+\pi I_j, \end{align} \tag{ 27 }$$

where the integral is over the locus of M-coupled algebraic equations equation (23). If the parameters $h_{1/N}\in \mathbb{R}$ are real, then the locus would simply span the real line, and hence the limit of the integral would be $(-\infty,\infty)$. However when the parameters are allowed to be complex $h_{1/N}\in \mathbb{C}$, then the integration contour may lie in different regions of the complex plane depending on the boundary parameters as shown in figure 2.

Zoom In Zoom Out Reset image size

We extract the density of roots in the ground state density $\rho_0(\mu)$ by subtracting equation (27) written for µ_j from the same equation written for $\mu_{j+1}$ and expanding in the difference $\Delta\mu = \mu_{j+1}-\mu_j$. This gives

$$\begin{align} 2\rho_0(\mu) = f(\mu)-\int_\mathcal{C}\mathrm{d}\mu~ K(\mu-\mu^{\prime})\rho_0(\mu^{\prime})-\int_\mathcal{C}\mathrm{d}\mu~ K(\mu+\mu^{\prime})\rho_0(\mu^{\prime})+\mathcal{O}\left(\frac{1}{N}\right). \end{align} \tag{ 28 }$$

Since we are interested in thermodynamics limit $N\to \infty$, the higher order terms are negligible. Here,

$$\begin{align} f(\mu)& = (2N+1)a_{\frac12}(\mu)-a_{\frac{1}{2}-\xi_-}(\mu)-a_{\frac{1}{2}-\xi_+}(\mu), \end{align} \tag{ 29 }$$

$$\begin{align} \kern-7.5pt a_\gamma(\mu)& = \frac{1}{\pi}\frac{\gamma}{\mu^2+\gamma^2}, \quad\text{and} \end{align} \tag{ 30 }$$

$$\begin{align} K(\mu)& = \frac{1}{\pi \left(\mu ^2+1\right)} = a_1(\mu). \end{align} \tag{ 31 }$$

We need to solve the integral equation equation (28) in order to find the eigenvalues given by equation (26). However, it is very difficult to solve the coupled algebraic Bethe Ansatz equation (23) exactly to find its locus, over which we need to perform the integration to find the root density. In order to get around this problem, we choose $\mathscr{PT}$-symmetric boundary fields $h_1 = \frac{1}{\xi_-} = \frac{1}{\xi+i\chi}$ and $h_N = \frac{1}{\xi_+} = \frac{1}{\xi-i\chi}$ where $\{\xi,\chi\}\in\mathbb{R}$. Let us recall the action of time reversal and parity in a discrete system; the time reversal operation $\mathscr{T}$ is such that $\mathscr{T} \textrm i \mathscr{T} = -\textrm i$ and the effect of the parity on a system of N spins is such that $\mathscr{P} \sigma_{l}^{\alpha} \mathscr{P} = \sigma_{N+1-l}^{\alpha}$. Thus, a $\mathscr{PT}$-symmetric XXX model has the Hamiltonian of the form

$$\begin{align} \mathcal{H} = \sum_{j = 1}^{N-1}\sum_{\alpha = 1}^3\sigma_j^\alpha\sigma_{j+1}^\alpha+\frac{1}{\xi+\textrm i\chi} \sigma_1^z+\frac{1}{\xi-\textrm i\chi}\sigma_N^z. \end{align} \tag{ 32 }$$

One can immediately see that the two $\tan^{-1}$ functions in the first line of equation (24) combine to have a real argument and hence the solution of the logarithmic Bethe Ansatz is real just as in the real case. One can numerically look at the solution of equation (24) with the $\mathscr{PT}$-symmetric boundary fields.

Since the roots lie in the real line as shown in figure 3, we can take the integration over the entire real line in equation (28) and solve for the density distribution of Bethe roots.

Zoom In Zoom Out Reset image size

5.1. Boundary strings solution

It is well known that the solutions of the Bethe Ansatz equations for arbitrary M, in the thermodynamics limit $N\to\infty$, contain spin complexes of various lengths $\ell$. An $\ell$-spin complex consists of an aggregate of $\ell$ complex rapidities µ_j , which all have a common real part and whose imaginary parts differ by unity. The spin complex is thus described by $\ell$ complex numbers,

$$\begin{align} \mu^{\ell, j} = \mu^\ell+\frac{\textrm i}{2}(\ell+1-2 j), \quad j = 1,2, \ldots, \ell \end{align} \tag{ 33 }$$

where $\mu^\ell \in \mathbb{R}$. These solutions describe the complex bound states of magnetic excitations.

From the Bethe Ansatz equation equation (23), it is easy to see that when $\xi\lt\frac{1}{2}$, there are two additional solutions in the thermodynamic limit. We call these boundary string solutions. The left boundary string is of the form

$$\begin{align} \mu_{L} = \pm \textrm i\left(\frac{1}{2}+(\xi+\textrm i\chi)\right), \end{align} \tag{ 34 }$$

and the right boundary string is of the form

$$\begin{align} \mu_{R} = \pm \textrm i\left(\frac{1}{2}+(\xi-i\chi)\right). \end{align} \tag{ 35 }$$

The parametric region $|\xi|\lt\frac{1}{2}$ where boundary string solutions exist are henceforth called A phases and the region $|\xi|\gt\frac{1}{2}$ where boundary solution do not exist are called B phases. We solve the Bethe Ansatz equation in different boundary parametric regimes for odd and even numbers of sites separately and list all the elementary excitations.

5.2. A phases

5.2.1.  A₁ phase.

Consider the phase A₁ where the real part of the boundary parameter takes the value $0\lt\xi\lt\frac{1}{2}$. The Bethe equation corresponding to the reference state with all spins down is obtained from equation (23) after applying the transformation equation (4) i.e.

$$\begin{align} & \left(\frac{\mu_{j}-\frac{\textrm i}{2}}{\mu_{j}+\frac{\textrm i}{2}}\right)^{2 N} \left(\frac{\mu_{j}-\chi+\textrm i \left(\frac 1 2+\xi\right)}{\mu_{j}+\chi-\textrm i\left(\frac{1}{2 }+\xi\right)}\right)\left(\frac{\mu_{j}+\chi+\textrm i\left(\frac{1}{2 }+\xi\right)}{\mu_{j}-\chi-\textrm i\left(\frac{1}{2 }+\xi\right)}\right)\nonumber\\ & \quad = \prod_{j \neq \ell = 1}^{M}\left(\frac{\mu_{j}-\mu_{\ell}-\textrm i}{\mu_{j}-\mu_{\ell}+\textrm i}\right)\left(\frac{\mu_{j}+\mu_{\ell}-\textrm i}{\mu_{j}+\mu_{\ell}+\textrm i}\right). \end{align} \tag{ 36 }$$

The eigenvalue of the Hamiltonian is given by

$$\begin{align} E = -\sum_{j = 1}^M \frac{2}{\mu_j^2+\frac{1}{4}}+N-1-\frac{1}{\xi+\textrm i\chi}-\frac{1}{\xi-\textrm i\chi} = -\sum_{j = 1}^M \frac{2}{\mu_j^2+\frac{1}{4}}+N-1-\frac{2 \xi }{\xi ^2+\chi ^2}. \end{align} \tag{ 37 }$$

Due to the $\mathscr{PT}$-symmetry that the energy is real even though the Hamiltonian is not Hermitian.

5.2.1.1.  A₁ phase: odd number of sites.

The ground state root density is obtained from equation (28), which is given by

$$\begin{align} 2\rho_{\left|-\frac12\right\rangle}(\mu) & = (2N+1)a_{\frac12}(\mu)-a_{\frac{1}{2}+\xi+\textrm i\chi}(\mu)-a_{\frac{1}{2}+\xi-\textrm i\chi}(\mu)\nonumber\\ & \quad -\sum_{\upsilon = \pm}\int_{-\infty}^\infty\mathrm{d}\mu~ a_1(\mu+\upsilon\mu^{\prime})\rho_{\left|-\frac12\right\rangle}(\mu^{\prime})-\delta(\mu). \end{align} \tag{ 38 }$$

Where the delta function is added to remove the solution $\mu_j = 0$ solution which results in a vanishing wavefunction [36].

The solution of the above equation in the Fourier space takes the following form

$$\begin{align} \tilde\rho_{\left|-\frac12\right\rangle}(\omega) = \frac{(2 N+1) \textrm e^{-\frac{| \omega | }{2}}-\textrm e^{- \left(\frac12+ \xi +\textrm i\chi\right) | \omega | }-\textrm e^{- \left( \frac12+\xi-\textrm i\chi\right) | \omega | }-1}{2 \left(\textrm e^{-| \omega | }+1\right)} .\end{align} \tag{ 39 }$$

The total number of Bethe roots is given by

$$\begin{align} M_{\left|-\frac12\right\rangle} = \tilde\rho_{\left|-\frac12\right\rangle}(\omega = 0) = \frac{N-1}{2}. \end{align} \tag{ 40 }$$

We can compute the z-component of the total spin, where

$$\begin{align} S^z_{\left|-\frac12\right\rangle} = -\left(\frac{N}{2}-M_{\left|-\frac12\right\rangle}\right) = -\frac{1}{2}. \end{align} \tag{ 41 }$$

To compute the energy using equation (26), we need to compute the integral

$$\begin{align} -\int_{-\infty}^\infty \mathrm{d}\mu \frac{2\rho_{\left|-\frac12\right\rangle}(\mu)}{\mu^2+\frac{1}{4}} & = \pi-(2 N+1) \ln (4)+\psi ^{(0)}\left(\frac{\xi+\textrm i\chi}{2}+1\right)-\psi ^{(0)}\left(\frac{\xi+\textrm i\chi}{2}+\frac{1}{2}\right)\nonumber\\ & \quad +\psi ^{(0)}\left(\frac{\xi-\textrm i\chi}{2}+1\right)-\psi ^{(0)}\left(\frac{\xi-\textrm i\chi}{2}+\frac{1}{2}\right), \end{align} \tag{ 42 }$$

where $\psi^0(z) = \frac{\textrm d}{\textrm dz}\left(\ln\Gamma(z)\right)$ is the diGamma function.

Thus, the energy of this state is

$$\begin{align} E_{\left|-\frac12\right\rangle} = E_0 & = N-1+\pi-(2 N+1) \ln (4)-\frac{1}{\xi+\textrm i\chi}-\frac{1}{\xi-\textrm i\chi}+\psi ^{(0)}\left(\frac{\xi+\textrm i\chi}{2}+1\right)\nonumber\\ & \quad -\psi ^{(0)}\left(\frac{\xi+\textrm i\chi}{2}+\frac{1}{2}\right)+\psi ^{(0)}\left(\frac{\xi-\textrm i\chi}{2}+1\right)-\psi ^{(0)}\left(\frac{\xi -\textrm i\chi}{2}+\frac{1}{2}\right)\in\mathbb{R}. \end{align} \tag{ 43 }$$

Adding the boundary string solution $\mu_L = \pm \textrm i\left(\frac{1}{2}+(\xi+\textrm i\chi)\right)$, we get

$$\begin{align} 2\rho_{\vert 0\rangle_L}(\mu)& = (2N+1)a_{\frac12}(\mu)-a_{\frac{1}{2}+\xi+\textrm i\chi}(\mu)-a_{\frac{1}{2}+\xi-\textrm i\chi}(\mu)-a_{\frac{1}{2}-\xi-\textrm i\chi}(\mu)-a_{\frac{3}{2}+\xi +\textrm i \chi }(\mu)\nonumber\\ & \quad -2\int_{-\infty}^\infty\mathrm{d}\mu~ a_1(\mu-\mu^{\prime})\rho_{\vert 0\rangle_L}(\mu^{\prime})-\delta(\mu). \end{align} \tag{ 44 }$$

In the Fourier space, we obtain the solution

$$\begin{align} \tilde\rho_{{\left|0\right\rangle}_L}(\omega) = \tilde{\rho}_{\left|-\frac12\right\rangle}(\omega)+\Delta\tilde\rho_{L}(\omega) \quad \text{where}~ \Delta\tilde\rho_{L}(\omega) = -\frac{\textrm e^{-\left(\frac{1}{2}-\xi -\textrm i \chi \right)| \omega | }+\textrm e^{-\left(\frac{3}{2}+\xi +\textrm i \chi \right)| \omega | }}{2 \left(\textrm e^{-| \omega | }+1\right)}. \end{align} \tag{ 45 }$$

The total number of roots is given by

$$\begin{align} M_L = 1+\tilde\rho_{0_L}(0) = 1+\frac{N-1}{2}-\frac{1}{2} = \frac{N}{2}. \end{align} \tag{ 46 }$$

Hence, the z-component of the total spin is

$$\begin{align} S^z_{L} = 0. \end{align} \tag{ 47 }$$

Notice that the number of roots is an integer only if the number of sites is even. Since we have a chain with an odd number of sites, in order for one to add a boundary string to the above state, a propagating hole (spinon) needs to be added. Adding a spinon with rapidity θ, we obtain

$$\begin{align} 2\rho_{{\left|-\frac12\right\rangle}_{\theta,L}}(\mu) & = (2N+1) a_{\frac12}(\mu)- a_{\frac{1}{2}+\xi+\textrm i\chi}(\mu)-a_{\frac{1}{2}+\xi-\textrm i\chi}(\mu)-a_{\frac{1}{2}-\xi-\textrm i\chi}(\mu)-a_{\frac{3}{2}+\xi +\textrm i \chi }(\mu)\nonumber\\ & \quad -2\int_{-\infty}^\infty\mathrm{d}\mu~ a_1(\mu-\mu^{\prime})\rho_{{\left|-\frac12\right\rangle}_{\theta,L}}(\mu^{\prime})-\delta(\mu)-\delta(\mu-\theta)-\delta(\mu+\theta). \end{align} \tag{ 48 }$$

The solution of the above equation in the Fourier space is immediate

$$\begin{align} \tilde\rho_{{\left|-\frac12\right\rangle}_{\theta,L}}(\omega) & = \tilde{\rho}_{\left|-\frac12\right\rangle}(\omega)+\Delta\tilde\rho_{L}(\omega)+\Delta\tilde\rho_\theta(\omega)\nonumber\\ \text{where}~ \Delta\tilde\rho_{\theta}(\omega) & = \frac{-\textrm e^{-\textrm i \theta \omega }-\textrm e^{\textrm i \theta \omega }}{2 \left(\textrm e^{-| \omega | }+1\right)} = -\frac{\cos(\theta\omega)}{\left(\textrm e^{-| \omega | }+1\right)}. \end{align} \tag{ 49 }$$

The total number of real roots is given by

$$\begin{align} M_{{\left|-\frac12\right\rangle}_{\theta,L}} = 1+\tilde{\rho}_{{\left|-\frac12\right\rangle}_{\theta,L}}(0) = \frac{N-1}{2}. \end{align} \tag{ 50 }$$

Thus, the z-component of the total spin is

$$\begin{align} S_{\theta,L}^z = -\left(\frac{N}{2}-\left(\frac{N}{2}-\frac{1}{2}\right)\right) = -\frac12. \end{align} \tag{ 51 }$$

We can compute the energy of the system using equation (26), and we obtain

$$\begin{align} E_{{\left|-\frac12\right\rangle}_{\theta,L}} & = -\int_{-\infty}^\infty\mathrm{d}\mu~ \rho_{{\left|-\frac12\right\rangle}_{\theta,L}}(\mu)\frac{2}{\mu^2+\frac{1}{4}}+N-1+\frac{1}{\xi+\textrm i\chi}\nonumber\\ & \quad +\frac{1}{\xi-\textrm i\chi}-\frac{2}{\frac{1}{4}+\left(\textrm i\left(\frac{1}{2}+(\xi+\textrm i\chi)\right)\right)^2}. \end{align} \tag{ 52 }$$

Using equation (43), the energy given by equation (52) can be simplified as

$$\begin{align} E_{{\left|-\frac12\right\rangle}_{\theta,L}} = E_0+\frac{2\pi}{\sin(\pi(\xi+\textrm i\chi))}+\frac{2\pi}{\cosh(\pi\theta)}. \end{align} \tag{ 53 }$$

The first term is the energy of the ground state $\left|-\frac12\right\rangle$, the second term

$$\begin{align} m = \frac{2\pi}{\sin(\pi(\xi+\textrm i\chi))} \end{align} \tag{ 54 }$$

is the energy of the bound state at the left boundary and the third term

$$\begin{align} E_\theta = \frac{2\pi}{\cosh(\pi\theta)} \end{align} \tag{ 55 }$$

is the energy of the spinon propagating with rapidity θ. Notice that the first and the third term are real but the second term is complex.

The real part of the energy is

$$\begin{align} \Re \left(E_{{\left|-\frac12\right\rangle}_{\theta,L}}\right) = E_0+\frac{2\pi}{\cos(\pi\theta)}+\frac{4 \pi \sin (\pi \xi ) \cosh (\pi \chi )}{\cosh (2 \pi \chi )-\cos (2 \pi \xi )}. \end{align} \tag{ 56 }$$

The energy of the spinon is strictly positive with energy range $0\lt E_\theta\lt2\pi$. Likewise, the real part of the energy of the bound state is also strictly positive with a range $E_L\gt0$. When χ = 0, the lowest value of the real part of the energy of the bound state is $E_{L}^- = 2\pi$, but as $\chi\to\infty$, the lowest value of the real part of the energy vanishes.

The imaginary part of the energy is

$$\begin{align} \Im \left(E_{{\left|-\frac12\right\rangle}_{\theta,L}}\right) = -\frac{4 \pi \cos (\pi \xi ) \sinh (\pi \chi )}{\cosh (2 \pi \chi )-\cos (2 \pi \xi )}. \end{align} \tag{ 57 }$$

Since $0\lt\xi\lt\frac{1}{2}$ and χ > 0, the imaginary part of the equation is always negative. Exactly at $\xi = \frac12$, the imaginary part of the energy vanishes for all values of χ, and at χ = 0 and $\chi = \infty$, the imaginary part of the energy vanishes for all values of ξ. The negativity of the imaginary part of the energy suggests that there is loss at the left boundary of the system.

We can obtain the state $\left|-\frac{1}{2}\right\rangle_{\theta,R}$ by adding the right boundary string and a spinon with rapidity θ. This state is described by following the density distribution

$$\begin{align} \tilde\rho_{{\left|-\frac12\right\rangle}_{\theta,R}}(\omega) & = \tilde{\rho}_{\left|-\frac12\right\rangle}(\omega)+\Delta\tilde\rho_{R}(\omega)+\Delta\tilde\rho_\theta(\omega)\nonumber\\ \text{where}~ \; \Delta\tilde\rho_{R}(\omega) & = -\frac{\textrm e^{-\left(\frac{1}{2}-\xi +\textrm i \chi \right)| \omega | }+\textrm e^{-\left(\frac{3}{2}+\xi -\textrm i \chi \right)| \omega | }}{2 \left(\textrm e^{-| \omega | }+1\right)}, \end{align} \tag{ 58 }$$

and its energy is given by

$$\begin{align} E_{{\left|-\frac12\right\rangle}_{\theta,R}} = E_0+\frac{2\pi}{\sin(\pi(\xi-\textrm i\chi))}+\frac{2\pi}{\cosh(\pi\theta)}. \end{align} \tag{ 59 }$$

Notice that this is exactly the complex conjugate of the energy equation (53) of the state $\vert -\frac{1}{2}\rangle_{\theta,L}$ constructed above. Once again, the first term is the energy of the ground state configuration, the third term is the energy of the spinon, and the second term

$$\begin{align} m^* = \frac{2\pi}{\sin(\pi(\xi-\textrm i\chi))} \end{align} \tag{ 60 }$$

is the energy of the boundary mode situated at the right end of the chain.

The real part of the energy is

$$\begin{align} \Re \left(E_{{\left|-\frac12\right\rangle}_{\theta,R}}\right) = E_0+\frac{2\pi}{\cos(\pi\theta)}+\frac{4 \pi \sin (\pi \xi ) \cosh (\pi \chi )}{\cosh (2 \pi \chi )-\cos (2 \pi \xi )}, \end{align} \tag{ 61 }$$

and the imaginary part of the energy is

$$\begin{align} \Im \left(E_{{\left|-\frac12\right\rangle}_{\theta,R}}\right) = \frac{4 \pi \cos (\pi \xi ) \sinh (\pi \chi )}{\cosh (2 \pi \chi )-\cos (2 \pi \xi )}. \end{align} \tag{ 62 }$$

Notice that the imaginary part of the energy of the bound state at the right boundary is exactly equal to the negative of the imaginary part of the energy of the bound state at the left boundary. The imaginary part of the energy is strictly positive and hence there is gain at the right boundary of the system.

Now, we consider the Bethe Ansatz equation corresponding to the reference state with all spins up given by equation (23).

The eigenvalues of the Hamiltonian are given by

$$\begin{align} E = -\sum_{j = 1}^M \frac{2}{\mu_j^2+\frac{1}{4}}+N-1+\frac{1}{\xi+\textrm i\chi}+\frac{1}{\xi-\textrm i\chi} = -\sum_{j = 1}^M \frac{2}{\mu_j^2+\frac{1}{4}}+N-1+\frac{2 \xi }{\xi ^2+\chi ^2}. \end{align} \tag{ 63 }$$

Taking the logarithm of the Bethe Ansatz equation, we extract the density distribution of the Bethe roots using the same process as above. We obtain the following density distribution in the Fourier space

$$\begin{align} \tilde\rho_{\left|\frac12\right\rangle}(\omega) = \frac{(2 N+1) \textrm e^{-\frac{| \omega | }{2}}-\textrm e^{- \left(\frac12- \xi -\textrm i\chi\right) | \omega | }-\textrm e^{- \left( \frac12-\xi+\textrm i\chi\right) | \omega | }-1}{2 \left(\textrm e^{-| \omega | }+1\right)}. \end{align} \tag{ 64 }$$

The total number of Bethe roots is given by

$$\begin{align} M_{\left|\frac12\right\rangle} = \tilde\rho_{\left|\frac12\right\rangle}(\omega = 0) = \frac{N-1}{2}. \end{align} \tag{ 65 }$$

We can compute the z-component of the total spin, where

$$\begin{align} S^z = \left(\frac{N}{2}-M_{\left|-\frac12\right\rangle}\right) = \frac{1}{2}. \end{align} \tag{ 66 }$$

The energy of the state is given by

$$\begin{align} E_{\left|\frac{1}{2}\right\rangle} = -\int_{-\infty}^\infty\mathrm{d}\mu~ \rho_{\left|\frac{1}{2}\right\rangle}(\mu)\frac{2}{\mu^2+\frac{1}{4}}+N-1+\frac{1}{\xi+\textrm i\chi}+\frac{1}{\xi-\textrm i\chi}. \end{align} \tag{ 67 }$$

Writing

$$\begin{align} \tilde\rho_{\left|\frac12\right\rangle}(\omega) = \tilde\rho_{-\left|\frac12\right\rangle}(\omega)+\tilde\rho_*(\omega), \end{align} \tag{ 68 }$$

where

$$\begin{align} \tilde\rho_*(\omega) = \frac{\textrm e^{\left(\frac12+ \xi+\textrm i\chi\right) | \omega | }-\textrm e^{- \left(\frac12- \xi -\textrm i\chi\right) | \omega | }+\textrm e^{- \left( \frac12+\xi-\textrm i\chi\right) | \omega | }-\textrm e^{- \left( \frac12-\xi+\textrm i\chi\right) | \omega | }}{2 \left(\textrm e^{-| \omega | }+1\right)}, \end{align} \tag{ 69 }$$

and using equations (68) and (69), we can write the energy of the state $\left|\frac12\right\rangle$ given by equation (67) as

$$\begin{align} E_{\left|\frac{1}{2}\right\rangle} = E_{\left|-\frac{1}{2}\right\rangle}+\frac{2}{\xi+\textrm i\chi}+\frac{2}{\xi-\textrm i\chi}-\int_{-\infty}^\infty\mathrm{d}\mu~ \rho_*(\mu)\frac{2}{\mu^2+\frac{1}{4}}. \end{align} \tag{ 70 }$$

By computing the integral

$$\begin{align} -\int_{-\infty}^\infty\mathrm{d}\mu~ \rho_*(\mu)\frac{2}{\mu^2+\frac{1}{4}} = 2 \pi \csc (\pi (\xi +\textrm i \chi ))-\frac{2}{\xi +\textrm i \chi }+2 \pi \csc (\pi (\xi -\textrm i \chi ))-\frac{2}{\xi -i \chi }, \end{align} \tag{ 71 }$$

we obtain

$$\begin{align} E_{\left|\frac{1}{2}\right\rangle} = E_{\left|-\frac{1}{2}\right\rangle}+\frac{2\pi}{\sin(\pi(\xi+\textrm i\chi))}+\frac{2\pi}{\sin(\pi(\xi-\textrm i\chi))}. \end{align} \tag{ 72 }$$

The last two terms are the energies of the bound state at the left boundary and right boundary respectively. Thus, this state contains a bound state at both the left and the right boundaries. We will henceforth represent this state as $\left|\frac12\right\rangle_{L,R}$. We can rewrite the energy as

$$\begin{align} E_{{\left|\frac{1}{2}\right\rangle}_{L,R}} = E_{\left|-\frac{1}{2}\right\rangle}+\frac{8 \pi \sin (\pi \xi ) \cosh (\pi \chi )}{\cosh (2 \pi \chi )-\cos (2 \pi \xi )}, \end{align} \tag{ 73 }$$

which shows that the energy of the state is real and strictly positive.

For the BAE with reference state with all spins up, the boundary strings are of the form

$$\begin{align} \mu_{L^{\prime}} = \pm \textrm i\left(\frac{1}{2}-(\xi+\textrm i\chi)\right) \end{align} \tag{ 74 }$$

and

$$\begin{align} \mu_{R^{\prime}} = \pm \textrm i\left(\frac{1}{2}-(\xi-\textrm i\chi)\right). \end{align} \tag{ 75 }$$

Now, we can construct a state ${\left|\frac{1}{2}\right\rangle}_{\theta,L}$ by adding a spinon and the boundary string $\mu_{R^{\prime}} = \pm \textrm i\left(\frac{1}{2}-(\xi-\textrm i\chi)\right)$ to equation (64) to obtain

$$\begin{align} \tilde{\rho}_{{\left|\frac{1}{2}\right\rangle}_{\theta,L}}(\omega) & = \tilde{\rho}_{\left|\frac{1}{2}\right\rangle}(\omega)+\Delta\tilde\rho_\theta(\omega)+\Delta\tilde\rho_{R^{\prime}}(\omega) \nonumber\\ \text{where}~ \Delta\tilde\rho_{R^{\prime}}(\omega) & = -\frac{\textrm e^{-\left(\frac{1}{2}+\xi -\textrm i \chi \right)| \omega | }+\textrm e^{-\left(\frac{3}{2}-\xi +\textrm i \chi \right)| \omega | }}{2 \left(\textrm e^{-| \omega | }+1\right)}. \end{align} \tag{ 76 }$$

The total number of roots is given by

$$\begin{align} M_{{\left|\frac12\right\rangle}_{\theta,L}} = 1+\tilde{\rho}_{{\left|\frac12\right\rangle}_{\theta,L}}(0) = \frac{N-1}{2}. \end{align} \tag{ 77 }$$

Thus, the z-component of the total spin is

$$\begin{align} S_{\theta,L}^z = \left(\frac{N}{2}-\left(\frac{N}{2}-\frac{1}{2}\right)\right) = \frac12. \end{align} \tag{ 78 }$$

The energy of this state is given by

$$\begin{align} E_{{\left|\frac{1}{2}\right\rangle}_{\theta,L}} = E_{\left|\frac12\right\rangle}+\frac{2\pi}{\cosh(\pi\theta)}-\frac{2\pi}{\sin(\pi(\xi+\textrm i\chi))}. \end{align} \tag{ 79 }$$

Using equation (72), we get

$$\begin{align} E_{{\left|\frac{1}{2}\right\rangle}_{\theta,L}} = E_{\left|-\frac12\right\rangle}+\frac{2\pi}{\cosh(\pi\theta)}+\frac{2\pi}{\sin(\pi(\xi-\textrm i\chi))}. \end{align} \tag{ 80 }$$

This state contains a bound state at the left boundary and is degenerate with the state $\left|-\frac12\right\rangle_{\theta,L}$.

If we add a spinon with rapidity θ and the boundary string $\mu_{L^{\prime}} = \pm \textrm i\left(\frac{1}{2}-(\xi+\textrm i\chi)\right)$ to the state $\left|\frac12\right\rangle$, we obtain the state ${\left|\frac{1}{2}\right\rangle}_{\theta,R}$ which contains a bound state at the right boundary and is degenerate to the state ${\left|-\frac{1}{2}\right\rangle}_{\theta,R}$. All other excited states can be constructed by adding bulk excitations with purely real energy, such as an even number of spinons, bulk strings, and quartets. Thus, adding these bulk excitations on top of a $\mathscr{PT}$ symmetric elementary excitations with real energies like $\left|-\frac{1}{2}\right\rangle$ and $\left|\frac{1}{2}\right\rangle_{L/R}$, we construct excited states with real energies. Whereas, adding the bulk excitations on top of $\mathscr{PT}$ broken states $\left|\pm\frac{1}{2}\right\rangle_{\theta,R}$ and $\left|\pm\frac{1}{2}\right\rangle_{\theta,L}$, we create excited states in $\mathscr{PT}$ broken phase. In the thermodynamic limit, for each excited state built on top of state with bound state on the right edge, there exists a corresponding excited state built on top of state with bound state on the left edge whose energies are complex conjugates of each other.

5.2.1.2.  A₁ phase: even number of sites.

Recall that the ground state configuration of the spin chain with an odd number of sites has the root density

$$\begin{align} {\tilde\rho_{\left|-\frac12\right\rangle}}(\omega) = \frac{(2 N+1) \textrm e^{-\frac{| \omega | }{2}}-\textrm e^{- \left(\frac12+ \xi +\textrm i\chi\right) | \omega | }-\textrm e^{- \left( \frac12+\xi-\textrm i\chi\right) | \omega | }-1}{2 \left(\textrm e^{-| \omega | }+1\right)}, \end{align} \tag{ 81 }$$

thus, the number of roots is $M = \frac{N-1}{2}$. If N is even, the number of roots is not an integer. Thus, to consider the ground state configuration of the spin chain with an even number of sites, we need to add a spinon to equation (81) i.e. consider the following distribution

$$\begin{align} \tilde\rho_{{\left|-1\right\rangle}_\theta}(\omega) = \frac{(2 N+1) \textrm e^{-\frac{| \omega | }{2}}-\textrm e^{- \left(\frac12+ \xi +\textrm i\chi\right) | \omega | }-\textrm e^{- \left( \frac12+\xi-\textrm i\chi\right) | \omega | }-1}{2 \left(\textrm e^{-| \omega | }+1\right)}+\Delta\tilde\rho_\theta(\omega), \end{align} \tag{ 82 }$$

where $\Delta\tilde\rho_\theta(\omega) = -\frac{\cos(\omega\theta)}{(1+\textrm e^{-|\omega|}}$. The total number of roots is given by

$$\begin{align} M_{{\left|-1\right\rangle}_\theta} = \tilde\rho_{{\left|-1\right\rangle}_\theta}(\omega = 0) = \frac{N-2}{2}. \end{align} \tag{ 83 }$$

Clearly, the number of roots is an integer only for the spin chain with an even number of sites. The z-component of the total spin of this state is

$$\begin{align} S^z = -\left(\frac{N}{2}-M_{{\left|-1\right\rangle}_\theta}\right) = -1. \end{align} \tag{ 84 }$$

Using equation (26), we can compute the energy of this state and we obtain

$$\begin{align} E_{{\left|-1\right\rangle}_\theta} = E_0+\frac{2\pi}{\cosh(\pi\theta)}. \end{align} \tag{ 85 }$$

The lowest energy state for a spin chain with an even number of sites is parameterized by the rapidity θ of the spinon. Since cosh is a monotonically non-decreasing function, the ground state is obtained in the limit $\theta\to \infty$.

Next, considering the reference state with all spins down, we take the state with all real roots and add the left boundary string solution $\mu_L = \pm \textrm i\left(\frac{1}{2}+(\xi+\textrm i\chi)\right)$, which gives the root density of the form

$$\begin{align} \tilde\rho_{{\left|0\right\rangle}_L}(\omega)& = \frac{(2 N+1) \textrm e^{-\frac{| \omega | }{2}}-\textrm e^{- \left(\frac12+ \xi +\textrm i\chi\right) | \omega | }-\textrm e^{- \left( \frac12+\xi-\textrm i\chi\right) | \omega | }-1}{2 \left(\textrm e^{-| \omega | }+1\right)}+\Delta\tilde\rho_{L}(\omega),\\ & \quad \text{where }~ \Delta\tilde\rho_{L}(\omega) = -\frac{\textrm e^{-\left(\frac{1}{2}-\xi -\textrm i \chi \right)| \omega | }+\textrm e^{-\left(\frac{3}{2}+\xi +\textrm i \chi \right)| \omega | }}{2 \left(\textrm e^{-| \omega | }+1\right)}\nonumber . \end{align} \tag{ 86 }$$

The total number of roots is given by

$$\begin{align} M_{\vert 0\rangle_L} = 1+\tilde\rho_{0_L}(0) = 1+\frac{N-1}{2}-\frac{1}{2} = \frac{N}{2}. \end{align} \tag{ 87 }$$

Hence, the z-component of the total spin is

$$\begin{align} S^z_{\vert 0\rangle_L} = 0. \end{align} \tag{ 88 }$$

The energy of this state henceforth represented at $\vert 0\rangle_L$ is

$$\begin{align} E_{\vert 0\rangle_L} = E_0+\frac{2\pi}{\sin(\pi(\xi+\textrm i\chi))}. \end{align} \tag{ 89 }$$

Likewise, we can obtain the state $\vert 0\rangle_R$ with energy

$$\begin{align} E_{\vert 0\rangle_R} = E_0+\frac{2\pi}{\sin(\pi(\xi-\textrm i\chi))}, \end{align} \tag{ 90 }$$

by adding right boundary string solution $\mu_R = \pm \textrm i\left(\frac{1}{2}+(\xi-\textrm i\chi)\right)$ to the state containing all real roots with respect to the reference state with all spins down.

Let us add both string solutions and a spinon to create a state with root distribution

$$\begin{align} \tilde{\rho}_{\vert 0\rangle_{\theta,L,R}}(\omega) & = \frac{(2 N+1) \textrm e^{-\frac{| \omega | }{2}}-\textrm e^{- \left(\frac12+ \xi +\textrm i\chi\right) | \omega | }-\textrm e^{- \left( \frac12+\xi-\textrm i\chi\right) | \omega | }-1}{2 \left(\textrm e^{-| \omega | }+1\right)}\nonumber\\ & \quad +\Delta\tilde\rho_{L}(\omega)+\Delta\tilde\rho_{R}(\omega)+\Delta\tilde\rho_{\theta}(\omega). \end{align} \tag{ 91 }$$

The total number of roots is given by

$$\begin{align} M_{\vert 0\rangle_{\theta,L,R}} = 2+\rho_{\vert 0\rangle_{\theta,L,R}} = \frac{N}{2}, \end{align} \tag{ 92 }$$

And hence the z-component of the spin is

$$\begin{align} S^z_{\vert 0\rangle_{\theta,L,R}} = 0. \end{align} \tag{ 93 }$$

The energy of the state is given by

$$\begin{align} E_{\vert 0\rangle_{\theta,L,R}} = E_0+\frac{2\pi}{\sin(\pi(x+\textrm i\chi))}+\frac{2\pi}{\sin(\pi(\xi-\textrm i\chi))}+\frac{2\pi}{\cosh(\pi\theta)}. \end{align} \tag{ 94 }$$

Now, consider the state with all real roots corresponding to the reference state with all spins up and a spinon described by the density

$$\begin{align} \tilde\rho_{{\left|1\right\rangle}_\theta}(\omega) = \frac{(2 N+1) \textrm e^{-\frac{| \omega | }{2}}-\textrm e^{- \left(\frac12- \xi -\textrm i\chi\right) | \omega | }-\textrm e^{- \left( \frac12-\xi+\textrm i\chi\right) | \omega | }-1}{2 \left(\textrm e^{-| \omega | }+1\right)}+\Delta\tilde\rho_\theta(\omega). \end{align} \tag{ 95 }$$

The total number of roots is given by

$$\begin{align} M_{{\left|1\right\rangle}_\theta} = \tilde\rho_{{\left|1\right\rangle}_\theta}(\omega = 0) = \frac{N-2}{2}, \end{align} \tag{ 96 }$$

which is an integer only for the even number of sites. And the z-component of the spin of this state is

$$\begin{align} S^z = \left(\frac{N}{2}-M_{{\left|1\right\rangle}_\theta}\right) = 1. \end{align} \tag{ 97 }$$

We can compute the energy of this state, and we obtain

$$\begin{align} E_{\vert 1\rangle_\theta} = E_0+\frac{2\pi}{\sin(\pi(x+\textrm i\chi))}+\frac{2\pi}{\sin(\pi(\xi-\textrm i\chi))}+\frac{2\pi}{\cosh(\pi\theta)}. \end{align} \tag{ 98 }$$

This state contains a spinon and bound states at both boundaries. This state represented as $\vert 1\rangle_{\theta,L,R}$ henceforth is degenerate with the state $\vert 0\rangle_{\theta,L,R}$.

We can add both boundary strings $\mu_{L^{\prime}} = \pm \textrm i\left(\frac{1}{2}-(\xi+\textrm i\chi)\right)$ and $\mu_{R^{\prime}} = \mu_{L^{\prime}} = \pm \textrm i\left(\frac{1}{2}-(\xi-\textrm i\chi)\right)$ and a spinon to the state with all real roots with respect to the reference state with all spins up. We obtain a state with the following root density

$$\begin{align} \tilde{\rho}_{\vert 0\rangle_\theta}(\omega) & = \frac{(2 N+1) \textrm e^{-\frac{| \omega | }{2}}-\textrm e^{- \left(\frac12- \xi -\textrm i\chi\right) | \omega | }-\textrm e^{- \left( \frac12-\xi+\textrm i\chi\right) | \omega | }-1}{2 \left(\textrm e^{-| \omega | }+1\right)}\nonumber\\ & \quad+\Delta\tilde\rho_\theta(\omega)+\Delta\tilde\rho_\theta(\omega)+\Delta\tilde{\rho}_{R^{\prime}}(\omega)+\Delta\tilde\rho_{L^{\prime}}(\omega). \end{align} \tag{ 99 }$$

The total number of roots in this state is

$$\begin{align} M_{\vert 0\rangle_\theta} = 2+\tilde{\rho}_{\vert 0\rangle_\theta}(0) = \frac{N}{2}, \end{align} \tag{ 100 }$$

and the z-component of the spin of this state is equal to

$$\begin{align} S_{\vert 0\rangle_\theta} = 0. \end{align} \tag{ 101 }$$

The energy of this state is given by

$$\begin{align} E_{\vert 0\rangle_\theta} = E_0+\frac{2\pi}{\cosh(\pi\theta)}, \end{align} \tag{ 102 }$$

and is degenerate with the state $\vert -1\rangle_\theta$.

Note that the state $\vert 0\rangle_L$ with energy given by equation (89) can also be created by starting with all the real roots with respect to the reference state with all spins up and adding the right boundary string $\mu_R^{\prime} = \pm \textrm i\left(\frac{1}{2}-(\xi-\textrm i\chi)\right)$ and the state $\vert 0\rangle_R$ with energy given by equation (90) can likewise be constructed by adding left boundary string $\mu_R^{\prime}=\pm \textrm i\left(\frac{1}{2}-(\xi-\textrm i\chi)\right)$ on top of all the real roots with respect to the reference state with all up spins. All the other excited states in this phase can be constructed by adding the bulk excitations (spinons, bulk strings, and quartets) on top of these elementary excitations. Since the bulk excitations have purely real energies, adding them on top of the state in $\mathscr{PT}$-symmetric phases results in states in excited states in $\mathscr{PT}$-symmetric phases and similarly adding these excitations on top of the state in spontaneously broken $\mathscr{PT}$ phase results in excited states in spontaneously broken $\mathscr{PT}$ phase as explained earlier in detail.

5.2.2.  A₂ phase.

Consider the phase A₂ where the real part of the boundary parameter takes the value $0\gt\xi\gt-\frac{1}{2}$. In this phase, all the states can be constructed from the phase A₁ by using the property equation (4).

We can take all the states in phase A₁ and apply the transformation that all the spins up and spins downs are interchanged and both the real and imaginary parts of the boundary parameters change the signs.

$$\begin{align} \vert \uparrow\rangle\leftrightarrow \vert \downarrow\rangle;\quad \xi\to -\xi \quad \text{and} \quad \chi\to -\chi \end{align} \tag{ 103 }$$

to construct all the states in phase A₂.

5.3. B phases

5.3.1.  B₁ phase.

Let us consider the parameter regime $\xi\gt\frac{1}{2}$. The BAE with respect to the reference state with all spins up is

$$\begin{align} & \left(\frac{\mu_{j}-\frac{\textrm i}{2}}{\mu_{j}+\frac{i}{2}}\right)^{2 N} \left(\frac{\mu_{j}-\textrm i\left((\xi+\textrm i\chi)-\frac{1}{2}\right)}{\mu_{j}+\textrm i\left((\xi+\textrm i\chi)-\frac{1}{2}\right)}\right)\left(\frac{\mu_{j}-\textrm i\left((\xi-\textrm i\chi)-\frac{1}{2}\right)}{\mu_{j}+\textrm i\left((\xi-\textrm i\chi)-\frac{1}{2}\right)}\right)\nonumber\\ & \qquad = \prod_{j \neq \ell = 1}^{M}\left(\frac{\mu_{j}-\mu_{\ell}-\textrm i}{\mu_{j}-\mu_{\ell}+\textrm i}\right)\left(\frac{\mu_{j}+\mu_{\ell}-\textrm i}{\mu_{j}+\mu_{\ell}+\textrm i}\right), \end{align} \tag{ 104 }$$

which can be written as

$$\begin{align} \left(\frac{\mu_{j}-\frac{\textrm i}{2}}{\mu_{j}+\frac{\textrm i}{2}}\right)^{2 N} & = \left(\frac{\mu_j-\chi+\textrm i\left(\xi-\frac{1}{2}\right)}{\mu_j-\chi-\textrm i\left(\xi-\frac{1}{2}\right)} \right) \left(\frac{\mu_j+\chi+\textrm i\left(\xi-\frac{1}{2}\right)}{\mu_j+\chi-\textrm i\left(\xi-\frac{1}{2}\right)} \right)\nonumber\\ & \quad \times \prod_{j \neq \ell = 1}^{M}\left(\frac{\mu_{j}-\mu_{\ell}-\textrm i}{\mu_{j}-\mu_{\ell}+\textrm i}\right)\left(\frac{\mu_{j}+\mu_{\ell}-\textrm i}{\mu_{j}+\mu_{\ell}+\textrm i}\right). \end{align} \tag{ 105 }$$

The complex solution of the form $\mu_j = a+\textrm ib$ with b > 0 makes the numerator of the LHS vanish, but the complex µ_j for which RHS vanishes are of the from $\mu_j = \pm\chi-\textrm i\left(\xi-\frac{1}{2}\right)$ i.e. the imaginary part is negative. Thus, there are no boundary string solutions in this region.

5.3.1.1.  B₁ phase: odd number of sites.

The ground state is constructed from the reference state with all spins down. The root density is of the form

$$\begin{align} \tilde\rho_{\left|-\frac12\right\rangle}(\omega) = \frac{(2 N+1) \textrm e^{-\frac{| \omega | }{2}}-\textrm e^{- \left(\frac12+ \xi +\textrm i\chi\right) | \omega | }-\textrm e^{- \left( \frac12+\xi-\textrm i\chi\right) | \omega | }-1}{2 \left(\textrm e^{-| \omega | }+1\right)}. \end{align} \tag{ 106 }$$

The total number of Bethe roots is given by

$$\begin{align} M_{\left|-\frac12\right\rangle} = \tilde\rho_{\left|-\frac12\right\rangle}(\omega = 0) = \frac{N-1}{2}. \end{align} \tag{ 107 }$$

We can compute the z-component of the total spin, where

$$\begin{align} S^z = -\left(\frac{N}{2}-M_{\left|-\frac12\right\rangle}\right) = -\frac{1}{2}. \end{align} \tag{ 108 }$$

The energy of the state is $E = E_0$. All the other excited states in this phase can be constructed by adding the bulk excitations (spinons, bulk strings, and quartets) on top of the ground state. Since the ground state energy is real in this phase, all the excited states constructed by adding bulk excitations on top of this state have real energies. Hence all the states in this phase are in $\mathscr{PT}$ unbroken phase.

5.3.1.2.  B₁ phase: even number of sites.

Starting from the reference state with all spins down, the total number of roots is integer only when the number of sites is odd. Thus, we need to add a hole and consider the state with root density

$$\begin{align} \tilde\rho_{{\left|-1\right\rangle}_\theta}(\omega) = \frac{(2 N+1) \textrm e^{-\frac{| \omega | }{2}}-\textrm e^{- \left(\frac12+ \xi +\textrm i\chi\right) | \omega | }-\textrm e^{- \left( \frac12+\xi-\textrm i\chi\right) | \omega | }-1}{2 \left(\textrm e^{-| \omega | }+1\right)}+\Delta\tilde\rho_\theta(\omega). \end{align} \tag{ 109 }$$

Now, the total number of roots is given by

$$\begin{align} M_{{\left|-1\right\rangle}_\theta} = \tilde\rho_{{\left|-1\right\rangle}_\theta} = \frac{N-2}{2}. \end{align} \tag{ 110 }$$

Thus, the z-component of the spin is

$$\begin{align} S^z_{{\left|-1\right\rangle}_\theta} = -\left(\frac{N}{2}-\frac{N-2}{2}\right) = -1. \end{align} \tag{ 111 }$$

The energy of this state is

$$\begin{align} E_{{\left|-1\right\rangle}_\theta} = E_0+\frac{2\pi}{\cosh(\pi\theta)}. \end{align} \tag{ 112 }$$

We can create another state $\vert 0\rangle_\theta$ by adding a spinon on top of the state with all real roots with respect to the reference state with all spins up. The root density becomes

$$\begin{align} \tilde\rho_{{\left|0\right\rangle}_\theta}(\omega) = \frac{(2 N+1) \textrm e^{-\frac{| \omega | }{2}}+\textrm e^{- \left( \xi +\textrm i\chi-\frac12\right) | \omega | }+\textrm e^{- \left( \xi-\textrm i\chi-\frac12\right) | \omega | }-1}{2 \left(\textrm e^{-| \omega | }+1\right)}+\Delta\tilde\rho_\theta(\omega). \end{align} \tag{ 113 }$$

The total number of roots is given by

$$\begin{align} M_{{\left|0\right\rangle}_\theta} = \rho_{{\left|0\right\rangle}_\theta}(0) = \frac{N}{2}, \end{align} \tag{ 114 }$$

and hence the z-component of the spin is

$$\begin{align} S_{{\left|0\right\rangle}_\theta}^z = 0. \end{align} \tag{ 115 }$$

Once again using equation (26), the energy of the state is obtained as

$$\begin{align} E_{{\left|0\right\rangle}_\theta} = E_0+\frac{2\pi}{\cosh(\pi\theta)}. \end{align} \tag{ 116 }$$

The minimal energy state is obtained in the limit $\theta\to\infty$. Thus, the two-fold degenerate ground state is represented by $\vert 0\rangle_{\theta\to\infty}$ and $\vert -1\rangle_{\theta\to\infty}$.

5.3.2.  B₂ phase.

Let us consider the parameter regime $\xi\lt-\frac{1}{2}$.

5.3.2.1.  B₂ phase: odd number of sites.

The ground state is constructed from the reference state with all spins up. The root density is of the form

$$\begin{align} \tilde\rho_{\left|\frac12\right\rangle}(\omega) = \frac{(2 N+1) \textrm e^{-\frac{| \omega | }{2}}-\textrm e^{- \left(\frac12- \xi -\textrm i\chi\right) | \omega | }-\textrm e^{- \left( \frac12-\xi+\textrm i\chi\right) | \omega | }-1}{2 \left(\textrm e^{-| \omega | }+1\right)}. \end{align} \tag{ 117 }$$

The total number of Bethe roots is given by

$$\begin{align} M_{\left|\frac12\right\rangle} = \tilde\rho_{\left|\frac12\right\rangle}(\omega = 0) = \frac{N-1}{2}, \end{align} \tag{ 118 }$$

from which we can compute the z-component of the total spin, where

$$\begin{align} S^z = \left(\frac{N}{2}-M_{\left|\frac12\right\rangle}\right) = \frac{1}{2}. \end{align} \tag{ 119 }$$

The energy of the state is $E = E_0$. All the other excited states in this phase can be constructed by adding the bulk excitations (spinons, bulk strings, and quartets) on top of these doubly degenerate ground states. Since both the ground states have completely real energy in this phase, all the excited states constructed by adding bulk excitations on top of these states have real energies. Hence all the states in this phase are in $\mathscr{PT}$ unbroken phase.

5.3.2.2.  B₂ phase: even number of sites.

Using the transformation equation (103), we can construct the state in this phase from the B₁ phase with an even number of sites. Here, the ground state is two-fold degenerate $\vert 0\rangle_{\theta\to\infty}$ and $\vert 1\rangle_{\theta\to\infty}$. The former state is created by adding a spinon on the top of the real roots with respect to the reference state with all spins down, and the latter is created by adding a spinon on the top of the real roots with respect to the reference state with all spins up. The former has $\frac{N}{2}$ number of roots and $S^z = 0$, and the latter has $\frac{N-2}{2}$ number of the real roots and spin $S^z = 1$.

Following [34, 37, 41], we study the pure boundary effect by computing the wavefunction in the 1-magnon sector.

The Bethe Ansatz equation in the 1-magnon sector becomes

$$\begin{align} \left(\frac{\mu-\frac{\textrm i}{2}}{\mu+\frac{\textrm i}{2}}\right)^{2N}\frac{\mu+\textrm i\left(\frac{1}{2}-\xi_-\right)}{\mu-\textrm i\left(\frac{1}{2}-\xi_-\right)}\frac{\mu+\textrm i\left(\frac{1}{2}-\xi_+\right)}{\mu-\textrm i\left(\frac{1}{2}-\xi_+\right)} = 1. \end{align} \tag{ 120 }$$

And the wavefunction reads

$$\begin{align} F(x) = \left(1+\frac{\left(\mu -\frac{\textrm i}{2}\right) }{\mu +\frac{\textrm i}{2}}\left(\frac{1}{\xi_-}-1\right)\right) \left(\frac{\mu +\frac{\textrm i}{2}}{\mu -\frac{\textrm i}{2}}\right)^x-\left(1+\frac{\left(\mu +\frac{\textrm i}{2}\right) }{\mu -\frac{\textrm i}{2}}\left(\frac{1}{\xi _-}-1\right)\right) \left(\frac{\mu +\frac{\textrm i}{2}}{\mu -\frac{\textrm i}{2}}\right)^{-x}.\nonumber\\ \end{align} \tag{ 121 }$$

For $\mu = \pm \textrm i(\frac{1}{2}-\xi_-)$, the above wavefunction reduces to the non-normalized wavefunction of the left bound mode

$$\begin{align} F_L(x) = \pm\frac{\left(2 \xi _-1\right) }{\xi _-^2}\left(\frac{\xi _-1}{\xi _-}\right)^{-x}. \end{align} \tag{ 122 }$$

For $\mu = \pm \textrm i(\frac{1}{2}+\xi_-)$, we get the non-normalized wavefunction of the right bound mode

$$\begin{align} F_R(x) = \pm\frac{1-2 \xi _+}{\left(\xi _+-1\right) \xi _+}\left(\frac{\xi _+-1}{\xi _+}\right)^{-(N-x)}. \end{align} \tag{ 123 }$$

Properly normalizing the modulus of the squares of wavefunctions given by equations (122) and (123), we can plot the modulus squared wavefunction for both real and complex boundary fields. As shown in figure 4, we see that the wavefunction is exponentially localized sharply at the boundary in the case of the real boundary fields and as the imaginary part of the field is increased, the localized wavefunction starts to broaden thereby making the demarcation of the bulk and boundary flimsy.

Zoom In Zoom Out Reset image size

We considered the integrable Heisenberg spin chain with complex boundary fields and diagonalized by using the Bethe Ansatz method. For generic complex boundary fields, the roots of the Bethe equations lie in a complex plane and it is not, in general, possible to solve the resulting integral equation. However, for the choice of $\mathscr{PT}$-symmetric boundary fields, the roots of Bethe equations lie on the real line which enables us to solve the resulting integral equation for the density of the roots. This choice corresponds to the case where there is a balanced loss and gain in the model. For this choice of the boundary fields, we have shown that the system exhibits two types of phases named A and B. In the B type phase all the eigenstates have real energies in the thermodynamic limit, hence the $\mathscr{PT}$ symmetry is unbroken. The A type phase exhibits two sectors, where one sector comprises of eigenstates with real energies and the other sector comprises of eigenstates with complex energies, corresponding to $\mathscr{PT}$ symmetry unbroken and broken sectors respectively. Furthermore, we find that the ground state exhibited by the system changes depending on the orientation of the boundary fields. Each of the A and B type phases can be divided into two sub-phases named A₁, A₂ and B₁, B₂ respectively, depending on the ground state exhibited by the system.

The existence of two sectors in the A type phase is related to the existence of localized bound states at the left and right edges, whose energies are complex conjugates of each other. We computed the wavefunction in one magnon sector and we find that the exponentially localized boundary wavefunction starts to broaden as we increase the complex part of the boundary field. In addition to this, as the complex part of the boundary fields is increased, the energy of the bound state becomes lesser than the maximum energy of a single spinon which is 2π. Taking into account the broadening of the bound state wavefunction, and the merging of the energy of the bound state into the band of a single spinon branch, it is natural to interpret that for large values of χ the boundary string no longer corresponds to a bound state exponentially localized at the edge.

Recently, it was shown [37] that in the Heisenberg spin chain with real boundary fields, the Hilbert space comprises of a certain number of towers which depends on the number of bound states exhibited by the system. It was also shown that as the boundary parameters are changed, the bound states leak into the bulk, which results in an eigenstate phase transition where the number of towers of the Hilbert space changes. It is natural that the Heisenberg chain with complex boundary fields that we considered here may also exhibit such towers of excited states in the A type phases, where there exists exponentially localized bound states at the edges. If so, it would be interesting to understand the fate of these towers when one moves from A type to B type phases and also for large values of χ where the bound states may no longer exist.

No new data were created or analysed in this study.