Eigenvalue structure of a Bose–Einstein condensate in a $\mathcal {PT}$-symmetric double well

Non-Hermitian Hamiltonians obeying parity-time ($\mathcal {PT}$) symmetry can have entirely real eigenvalues, which is referred to as the case of unbroken $\mathcal {PT}$ symmetry [1]. Much effort has been put into the development of a more general formulation of quantum mechanics, where the requirement for Hermiticity is replaced by the weaker requirement of $\mathcal {PT}$ symmetry [2–6]. This includes the construction of a unitary time development and a positive definite scalar product in the case of unbroken $\mathcal {PT}$ symmetry [2, 7]. Also it is possible to embed $\mathcal {PT}$ symmetry in the more general approach of pseudo-Hermiticity [8–10].

$\mathcal {PT}$ symmetry is not only the subject of theoretical but also of experimental studies [11–14]. In particular the observation of $\mathcal {PT}$ symmetry in optical wave guide systems was a great success [11, 12, 15]. The experimental realization in a true quantum system, however, is still missing although suggestions already exist [16] for systems which have been extensively studied [17, 18]. In contrast to the above-named generalization of quantum mechanics, in experimental realizations one acts on the assumption that the complete system is Hermitian. Here the non-Hermiticity enters the picture as an effective description of an open quantum system by using it to model gain and loss contributions. The system we will investigate falls into this category. Following a suggestion by Klaiman et al [15] we study a Bose–Einstein condensate (BEC) where particles are coherently injected and removed from the system. Both a coherent influx and outflux of atoms have been experimentally realized. Electron beams have been successfully used to ionize atoms of the BEC and extract them with an electrostatic field [19]. Coherent incoupling of atoms into a BEC from an additional source condensate has been realized by exploiting various electronic excitations of the incoupled atoms [20].

While most of the work done in the field of $\mathcal {PT}$ symmetry treats the linear Schrödinger equation, a BEC in mean-field approximation is described by the nonlinear Gross–Pitaevskii equation (GPE). It is known that this equation provides accurate results for temperatures considerably smaller than the critical temperature of the condensate [21, 22], however, for dynamically unstable states the mean-field description breaks down [23–25]. The nonlinearity leads to remarkable effects, such as the coexistence of $\mathcal {PT}$-symmetric and $\mathcal {PT}$-broken states in certain parameter regimes [26–29].

In this paper we analyse the spectral properties of a BEC with contact interaction described by the dimensionless GPE,

$$\begin{equation} \left[-\frac{\partial ^2}{\partial x^2} +V(x) -g\left|\psi (x,t)\right|^2\right] \psi (x,t) = \mathrm{i}\frac{\partial }{\partial t}\psi (x,t), \end{equation} \tag{ 1 }$$

in the $\mathcal {PT}$-symmetric double-well potential

$$\begin{equation} V(x)=\case{1}{4} x^2 + v_0 \,\mathrm{e}^{-\sigma x^2} + \mathrm{i}\gamma x\,\mathrm{e}^{-\rho x^2}, \end{equation} \tag{ 2 }$$

shown in figure 1. Imaginary potentials are sinks or sources of probability density depending on the sign of the imaginary part. The physical interpretation is a coherent injection or removal of particles into or from the condensate. The non-Hermitian GPE can also be obtained in the mean-field limit of an open Bose–Hubbard model described by a master equation [30, 31]. The symmetric real part of the potential forms a double well consisting of a harmonic trap and a Gaussian barrier with height v₀ and width parameter σ. The parameter ρ of the antisymmetric imaginary part is chosen in such a way that its extrema coincide with the minima of the double well. Positive values of g describe an attractive contact interaction whereas negative values relate to a repulsive interaction. The separation ψ(x, t) = e^−iμtψ(x) leads to the time-independent GPE with the (nonlinear) eigenvalue μ which can be identified with the chemical potential and the (nonlinear) eigenstate ψ(x). In all following calculations v₀ = 4, σ = 0.5 and ρ ≈ 0.12 are fixed, while the gain/loss parameter γ and the strength of the nonlinearity g are varied to study the behaviour of the system. For the sake of clarity only a one-dimensional system is considered. An extension to three dimensions with additional harmonic potentials in y and z direction does not show further effects [28].

Zoom In Zoom Out Reset image size

Exact numerical results are gained by integrating the one-dimensional wave function outwards and by fulfilling the boundary conditions via a root search.

Before studying a BEC in the $\mathcal {PT}$-symmetric double well it is necessary to expand the concept of linear $\mathcal {PT}$ symmetry to nonlinear systems, which will be done in section 2. We apply an analytic continuation to the GPE from complex to bicomplex numbers in section 3 and present a full analysis of the four-dimensional bicomplex eigenvalue spectrum, the eigenfunctions and the symmetries of the system. The continuation allows the detailed investigation of the exceptional points of the system and their mathematical structure in section 4. A comparison of the eigenvalues and the exceptional points with those of a linear 4 × 4 matrix model is drawn in section 5. Finally it is shown in section 6 that the matrix model has the interesting feature of four interacting eigenmodes in the limit of vanishing interaction, and that this behaviour can also be found in the GPE.

A linear system is called $\mathcal {PT}$-symmetric if its Hamiltonian commutes with the combined action of the parity and time reversal operator $[\hat{H}, \mathcal {PT}]=0$. With $\hat{H}=\hat{p}^2/2m+V(\hat{x})$ and the operators $\mathcal {P}:\ \hat{x}\rightarrow -\hat{x},\ \hat{p}\rightarrow -\hat{p}$ and $\mathcal {T}:\ \hat{p}\rightarrow -\hat{p},\ \mathrm{i}\rightarrow -\mathrm{i}$ one gains the necessary condition $V(\hat{x})=V^*(-\hat{x})$, i.e. the real part of the potential must be symmetric and the imaginary part antisymmetric. Linear $\mathcal {PT}$-symmetric systems have the following properties [2]:

If an eigenstate of $\hat{H}$ is also an eigenstate of the operator $\mathcal {PT}$, i.e. $\mathcal {PT}\psi =\exp (\mathrm{i}\varphi )\psi , \varphi \in \mathbb {R}$, then the corresponding energy eigenvalue E is real. Also the inverse statement holds. The eigenspace to a real eigenvalue can always be constructed out of $\mathcal {PT}$-symmetric states.

With these two properties the following important conclusion can be derived. If and only if all eigenstates of the Hamiltonian can be written as eigenstates of the $\mathcal {PT}$ operator, the spectrum is entirely real. One refers to this case as unbroken $\mathcal {PT}$ symmetry, otherwise the $\mathcal {PT}$ symmetry is broken. Complex eigenvalues appear as complex conjugate pairs and their wave functions can be mapped on each other by application of the $\mathcal {PT}$ operator.

In nonlinear systems like BECs, we cannot use the simple commutator relation $[\hat{H},\mathcal {PT}]=0$ to define $\mathcal {PT}$ symmetry. In the following we expand this condition to nonlinear systems described by a Gross–Pitaevskii-like equation,

$$\begin{equation} \hat{H}_\mathrm{lin} \psi + f\left(\psi \right)\psi =\mathrm{i}\frac{\partial }{\partial t}\psi , \end{equation} \tag{ 3 }$$

where $\hat{H}_\mathrm{lin}=-\hat{p}^2/2m+V(\hat{x})$ and f(ψ) is a general nonlinear part. We restrict the discussion to nonlinearities which are invariant under the change of a global phase,

$$\begin{equation} f(\mathrm{e}^{\mathrm{i}\chi }\psi )=f(\psi), \quad \chi \in \mathbb {R}. \end{equation} \tag{ 4 }$$

This is necessary for the time-independent nonlinear Schrödinger equation to have the usual form. For nonlinear systems we replace the commutator relation of linear $\mathcal {PT}$-symmetric systems with the requirement

$$\begin{equation} \mathcal {PT}[(\hat{H}_\mathrm{lin} + f(\psi ))\psi ] \stackrel{!}{=} [\hat{H}_\mathrm{lin} + f(\mathcal {PT}\psi )]\mathcal {PT}\psi . \end{equation} \tag{ 5 }$$

We will see that this condition suffices to regain the properties of linear $\mathcal {PT}$ symmetry. If the linear part $\hat{H}_\mathrm{lin}$ is $\mathcal {PT}$-symmetric, (5) is reduced to a simple condition for the nonlinear part f(ψ),

$$\begin{equation} \mathcal {PT}f(\psi ) = f(\mathcal {PT}\psi ). \end{equation} \tag{ 6 }$$

In this work we are mainly interested in the properties of stationary states. Therefore we look at the time-independent nonlinear Schrödinger equation,

$$\begin{equation} \hat{H}_\mathrm{lin} \psi + f\left(\psi \right)\psi = \mu \psi . \end{equation} \tag{ 7 }$$

Application of the $\mathcal {PT}$ operator leads to

$$\begin{equation} \hat{H}_\mathrm{lin} \mathcal {PT}\psi + f\left(\mathcal {PT}\psi \right)\mathcal {PT}\psi = \mu ^*\mathcal {PT}\psi , \end{equation} \tag{ 8 }$$

where (6) and $[\hat{H}_\mathrm{lin},\mathcal {PT}]=0$ were used. Equations (7) and (8) show immediately that the energy eigenvalues μ occur in complex conjugate pairs with the eigenstates ψ and $\mathcal {PT}\psi$, respectively. Also the most striking property of $\mathcal {PT}$ symmetry, namely the concurrence of $\mathcal {PT}$-symmetric states and real eigenvalues, is true for such nonlinear systems. This can be seen by evaluating (8) for $\mathcal {PT}$-symmetric states $\mathcal {PT}\psi =\exp (\mathrm{i}\varphi )\psi$,

$$\begin{equation} \hat{H}_\mathrm{lin}\psi + f(\mathrm{e}^{\mathrm{i}\varphi }\psi )\psi = \mu ^*\psi . \end{equation} \tag{ 9 }$$

As stated in (4) only phase independent nonlinearities are considered, and thus for $\mathcal {PT}$-symmetric states the energy eigenvalue must be real μ = μ*.

Again this proof is also valid in the inverse direction. For non-degenerate real eigenvalues (7) and (8) demand that the eigenfunction is $\mathcal {PT}$-symmetric because ψ and $\mathcal {PT}\psi$ fulfil the same eigenvalue equation. However if a real eigenvalue is degenerate, it is in general not possible to choose $\mathcal {PT}$-symmetric eigenstates because the superposition principle is only valid in linear systems.

These results can be summarized as follows. In nonlinear non-degenerate systems of type (3) with a $\mathcal {PT}$-symmetric linear part and a nonlinear part which fulfils the conditions (4) and (6)

• the eigenvalues are either real or occur in complex conjugate pairs,
• the eigenvalues are real if and only if the eigenstate itself is $\mathcal {PT}$-symmetric,
• if ψ is an eigenstate to μ then $\mathcal {PT}\psi$ is eigenstate to μ*.

It is worth mentioning that for $\mathcal {PT}$-symmetric wave functions ψ the conditions (4) and (6) also guarantee that f(ψ) is $\mathcal {PT}$-symmetric, i.e. the real part of f is symmetric and the imaginary part is antisymmetric. Since for a given state ψ the nonlinearity has to be seen as a contribution to the potential V, this establishes the link to linear $\mathcal {PT}$ symmetry.

These general considerations are now applied to the Gross–Pitaevskii nonlinearity,

$$\begin{equation} f(\psi ,{\boldsymbol r}) = \int \mathrm{d}^3{\boldsymbol r}^{\prime }\, V({\boldsymbol r},{\boldsymbol r}^{\prime }) |\psi ({\boldsymbol r}^{\prime })|^2. \end{equation} \tag{ 10 }$$

Because the wave function only appears as square modulus, the nonlinearity is not changed by an arbitrary phase and thus condition (4) is always fulfilled independent of the interaction type. The second condition (6) carries over to a condition for the interaction potential $V({\boldsymbol r},{\boldsymbol r}^{\prime })$,

$$\begin{eqnarray} &&\int \mathrm{d}^3{\boldsymbol r}^{\prime }\, V^*(-{\boldsymbol r},{\boldsymbol r}^{\prime }) |\psi ({\boldsymbol r}^{\prime })|^2 = \int \mathrm{d}^3 {\boldsymbol r}^{\prime }\, V({\boldsymbol r},{\boldsymbol r}^{\prime }) |\psi ^*(-{\boldsymbol r}^{\prime })|^2,\nonumber \\ &&V({\boldsymbol r},{\boldsymbol r}^{\prime }) = V^*(-{\boldsymbol r},-{\boldsymbol r}^{\prime }). \end{eqnarray} \tag{ 11 }$$

The most common interaction potentials, namely the contact,

$$\begin{equation} V_\mathrm{c}({\boldsymbol r},{\boldsymbol r}^{\prime }) \propto \delta ^3({\boldsymbol r}-{\boldsymbol r}^{\prime }), \end{equation} \tag{ 12 }$$

monopolar [32],

$$\begin{equation} V_\mathrm{m}({\boldsymbol r},{\boldsymbol r}^{\prime }) \propto \frac{1}{|{\boldsymbol r}-{\boldsymbol r}^{\prime }|}, \end{equation} \tag{ 13 }$$

and dipolar interaction [33],

$$\begin{equation} V_\mathrm{d}({\boldsymbol r},{\boldsymbol r}^{\prime }) \propto \frac{1-3\cos ^2\vartheta }{|{\boldsymbol r}-{\boldsymbol r}^{\prime }|^3}, \end{equation} \tag{ 14 }$$

fulfil this requirement and therefore are possible candidates for $\mathcal {PT}$-symmetric BECs.

The main goal of this paper is to gain deeper insight into the mathematical structure of the eigenvalue spectrum of a BEC in the $\mathcal {PT}$-symmetric double well (2). We limit the following investigations to the behaviour of the two $\mathcal {PT}$-symmetric and the two $\mathcal {PT}$-broken bound states with the lowest real eigenvalues. This is justified because calculations of higher excited states show a similar behaviour but at significantly higher gain/loss contributions. In the parameter range considered the higher excited states are only merely affected by the imaginary part of the potential.

3.1. Eigenvalues and eigenfunctions

The eigenvalues and eigenfunctions have already been presented in [28]. The basic results necessary to understand the following discussions are briefly recapitulated. The eigenvalue spectrum of the ground state and first excited state is shown in figure 2 for different values of the nonlinearity parameter g. The linear case g = 0 shows the usual properties known from $\mathcal {PT}$-symmetric systems. The two $\mathcal {PT}$-symmetric solutions with real eigenvalues coalesce in a tangent bifurcation and for larger values of the gain/loss parameter γ two $\mathcal {PT}$-broken solutions with complex conjugate eigenvalues arise.

Zoom In Zoom Out Reset image size

For nonvanishing contact interactions we again expect a spectrum that contains solutions with real or complex conjugate eigenvalues. Figure 2 confirms this prediction. There is, however, a crucial difference to the linear case. The $\mathcal {PT}$-symmetric solutions coalesce and vanish at a critical value γ = γ_c2 in a tangent bifurcation while the $\mathcal {PT}$-broken states exist already at smaller values of γ and emerge from the ground state at γ_c1 in a pitchfork bifurcation. The physical reason for this behaviour is the stronger localization of the $\mathcal {PT}$-broken states which due to the attractive contact interaction leads to an energy reduction relative to the $\mathcal {PT}$-symmetric states. It is worth noting that for a repulsive interaction, i.e. g < 0, the $\mathcal {PT}$-broken states are shifted upwards and emerge from the first excited state.

The eigenvalue structure can be divided into the following regions. In the interval γ < γ_c1 only $\mathcal {PT}$-symmetric solutions exist, whereas in the regime γ > γ_c2 only the two $\mathcal {PT}$-broken solutions are present. In the intermediate region the $\mathcal {PT}$-symmetric and $\mathcal {PT}$-broken states coexist. This region grows for stronger nonlinearities g and for g ≳ 0.24 the $\mathcal {PT}$-broken states exist even at γ = 0. The existence of symmetry breaking states at γ = 0 is known as macroscopic quantum self-trapping [34].

As shown in section 2 the eigenvalues determine the symmetries of the corresponding wave functions. We immediately see in figure 3 that real eigenvalues coincide with $\mathcal {PT}$-symmetric wave functions. The states with complex eigenvalues have asymmetric, and thus $\mathcal {PT}$-broken, wave functions. They can be mapped on each other by application of the $\mathcal {PT}$ operator. The square modulus of the $\mathcal {PT}$-symmetric solutions is symmetric whereas it is asymmetric for the $\mathcal {PT}$-broken solutions thus breaking the symmetry of the system.

Zoom In Zoom Out Reset image size

3.2. Analytic continuation

The eigenvalue spectrum has the unusual property that the number of solutions is not conserved when the states undergo bifurcations even though complex solutions are also taken into account. Previous works showed that the nonanalytic nonlinearity of the GPE, which is proportional to the square modulus of the wave function |ψ|², is responsible for the change in the number of solutions at the bifurcation points [28, 35, 36].

It is possible to apply an analytic continuation in such a way that the number of solutions is conserved. This allows for a better understanding of the mathematical structure of the solutions, and in particular of the bifurcation scenarios and exceptional points. Also the analytic continuation is used to compare our solutions to a linear matrix model.

To apply the analytic continuation [35], the complex GPE is split into its real and imaginary parts,

$$\begin{equation} \left[{-}\Delta + V_\mathrm{R}-g \left(\psi _\mathrm{R}^2 +\psi _\mathrm{I}^2\right)\right] \psi _\mathrm{R} - V_\mathrm{I} \psi _\mathrm{I} = \mu _\mathrm{R}\psi _\mathrm{R}-\mu _\mathrm{I}\psi _\mathrm{I}, \end{equation} \tag{ 15a }$$

$$\begin{equation} \left[{-}\Delta +V_\mathrm{R}-g \left(\psi _\mathrm{R}^2 +\psi _\mathrm{I}^2\right)\right] \psi _\mathrm{I} + V_\mathrm{I} \psi _\mathrm{R} = \mu _\mathrm{R}\psi _\mathrm{I}+\mu _\mathrm{I}\psi _\mathrm{R}, \end{equation} \tag{ 15b }$$

with the real and imaginary parts denoted by R and I, respectively.

If we allow for complex values for the real and imaginary part of the wave function and the chemical potential the two coupled differential equations are analytic and the number of solutions is conserved. Since we have to clearly distinguish between the usual imaginary unit i and the complex extension in the calculations, a new imaginary unit j² = −1 is introduced to avoid any confusion. Then the complexified real and imaginary parts of the wave function read

$$\begin{equation} \psi _\mathrm{R} \rightarrow \psi _1 + \mathrm{j}\psi _\mathrm{j}, \end{equation} \tag{ 16a }$$

$$\begin{equation} \psi _\mathrm{I} \rightarrow \psi _\mathrm{i}+ \mathrm{j}\psi _\mathrm{k}, \end{equation} \tag{ 16b }$$

and analogously for μ_R and μ_I. Actually this is identical to extending the wave function in the original GPE to a bicomplex number with four components,

$$\begin{equation} \psi = \psi _1 + \mathrm{j}\psi _\mathrm{j}+ \mathrm{i}\psi _\mathrm{i}+ \mathrm{k}\psi _\mathrm{k}, \end{equation} \tag{ 17 }$$

with k = ij [36]. Bicomplex numbers are four-dimensional hypercomplex numbers with a commutative multiplication which allow an elegant description of the analytic continuation that complexifies real and imaginary part.

We now have to distinguish between the complex conjugation with respect to i and j, represented by the operators $\mathcal {T}_\mathrm{i}$ and $\mathcal {T}_\mathrm{j}$. The operator $\mathcal {T}_\mathrm{i}$ again has the physical interpretation of time reversal and its action in coordinate space is reduced to i → −i. Analogously we define $\mathcal {T}_\mathrm{j}$ as the complex conjugation j → −j which, however, has no inherent physical interpretation. The system now has two symmetries, the usual $\mathcal {PT}_\mathrm{i}$ symmetry and a new symmetry $\mathcal {T}_\mathrm{j}$. This leads to the following properties. For $\mathcal {PT}_\mathrm{i}$-symmetric wave functions the μ_i and μ_k components must vanish. If one of these components has a finite value then this value must occur with positive and negative sign and the corresponding $\mathcal {PT}_\mathrm{i}$-broken wave functions can be mapped on each other by application of the $\mathcal {PT}_\mathrm{i}$ operator. The same holds for the $\mathcal {T}_\mathrm{j}$ operator for the components μ_k and μ_j and additionally for the combined operator $\mathcal {PT}_\mathrm{i}\mathcal {T}_\mathrm{j}$ for the components μ_j and μ_i. This is summarized in table 1.

Table 1. A wave function with a specific symmetry can only have two nonvanishing components given by the second column. The components in the third column must vanish. If the wave function breaks a symmetry, the components in the third column occur in pairs with the same absolute value but different sign, and the corresponding wave functions can be mapped on each other by application of the respective symmetry operator.

Symmetry	Nonvanishing	Vanishing
$\mathcal {PT}_\mathrm{i}$	μ1, μj	μi, μk
$\mathcal {T}_\mathrm{j}$	μ1, μi	μk, μj
$\mathcal {PT}_\mathrm{i}\mathcal {T}_\mathrm{j}$	μ1, μk	μj, μi

3.3. Analytically continued spectrum

Using bicomplex wave functions the chemical potential μ of the eigenstates also becomes a bicomplex number with four components

$$\begin{equation} \mu =(\mu _1+\mathrm{j}\mu _\mathrm{j})+\mathrm{i}(\mu _\mathrm{i}+\mathrm{j}\mu _\mathrm{k}). \end{equation} \tag{ 18 }$$

Figure 4 shows that the GPE now has four solutions independent of the parameter γ and thus the number of solutions is constant, as expected. The already known solutions occur unchanged in the spectrum.

Zoom In Zoom Out Reset image size

For the $\mathcal {PT}_\mathrm{i}$-symmetric states with real eigenvalues existing for γ < γ_c2 the value of μ₁ is identical to the real part of μ without the analytic continuation. The other three components vanish μ_j = μ_i = μ_k = 0. These solutions coalesce as before in an tangent bifurcation at γ_c2. However, they do not vanish at this point. They are continued to the regime γ > γ_c2 by two states that only exist within the complex extension. These two additional states have nonzero μ₁ and μ_j components, i.e. their chemical potentials have complexified real parts μ = μ₁ + jμ_j.

At γ_c1 the branches with the already known $\mathcal {PT}_\mathrm{i}$-broken solutions emerge from the ground state in a pitchfork bifurcation. The chemical potential of these states has the form μ = μ₁ + iμ_i. On these branches the chemical potential is a usual complex number because the components μ_j and μ_k are zero. Since the number of solutions is also not changed at γ_c1, two additional branches arise in the interval 0 < γ < γ_c1 which continue the $\mathcal {PT}_\mathrm{i}$-broken solutions to values of γ smaller than γ_c1. The eigenvalues on these branches have the form μ = μ₁ + kμ_k.

The wave functions on the new branches ψ_{k ±} and ψ_{j ±} are shown in figure 5. To discuss the symmetries of the eigenstates one should point out the actions of the symmetry operators on the wave functions' components. The $\mathcal {PT}_\mathrm{i}$ operator reflects all components at x = 0 and additionally the ψ_i and ψ_k components at the x axis. The $\mathcal {T}_\mathrm{j}$ operator only reflects the ψ_j and ψ_k components at the x axis. The combined action is represented by the $\mathcal {PT}_\mathrm{i}\mathcal {T}_\mathrm{j}$ operator which reflects all components at x = 0 and the ψ_i and ψ_j components at the x axis.

Zoom In Zoom Out Reset image size

The branches ψ_{k ±} with nonzero μ₁ and μ_k components have, according to table 1, wave functions with $\mathcal {PT}_\mathrm{i}\mathcal {T}_\mathrm{j}$ symmetry and broken $\mathcal {PT}_\mathrm{i}$ and $\mathcal {T}_\mathrm{j}$ symmetry. This can be verified by looking at the corresponding wave functions in figures 5(a) and (b). The two wave functions are mapped on themselves by application of the $\mathcal {PT}_\mathrm{i}\mathcal {T}_\mathrm{j}$ and are mapped on each other by application of the $\mathcal {PT}_\mathrm{i}$ or $\mathcal {T}_\mathrm{j}$ operator. Therefore we can conclude that the branches with complex eigenvalues are continued to the interval γ < γ_c1 using the analytic continuation, and the property is conserved that they are $\mathcal {PT}_\mathrm{i}$-broken. However, there is a significant difference. The $\mathcal {PT}_\mathrm{i}$-broken branches ψ_{i ±} in the interval γ > γ_c1 have asymmetric eigenstates, whereas the wave functions of the continued branches have symmetric and antisymmetric components although the complete wave functions are $\mathcal {PT}_\mathrm{i}$-broken. The square modulus appearing in the GPE

$$\begin{equation} \left|\psi \right|^2=\psi _\mathrm{R}^2+\psi _\mathrm{I}^2 =(\psi _1+\mathrm{j}\psi _\mathrm{j})^2 + (\psi _\mathrm{i}+\mathrm{j}\psi _\mathrm{k})^2 \end{equation} \tag{ 19 }$$

has a symmetric one-component and an antisymmetric j-component for the branches ψ_{k ±}. Thus the square modulus is asymmetric, which is consistent with the wave functions being $\mathcal {PT}_\mathrm{i}$-broken.

The figures 5(c) and (d) show the wave functions on the branches which arise at the pitchfork bifurcation at γ_c2. Because only the μ₁ and μ_j components are nonzero, comparison with table 1 shows that the eigenstates must obey $\mathcal {PT}_\mathrm{i}$ symmetry, but are $\mathcal {T}_\mathrm{j}$-broken and hence also $\mathcal {PT}_\mathrm{i}\mathcal {T}_\mathrm{j}$-broken. Therefore the wave functions are invariant under the action of the $\mathcal {PT}_\mathrm{i}$ operator and are mapped on each other by the $\mathcal {T}_\mathrm{j}$ and $\mathcal {PT}_\mathrm{i}\mathcal {T}_\mathrm{j}$ operator. This means $\mathcal {PT}_\mathrm{i}$-symmetric solutions exist for all values of γ if the analytic continuation is applied. For all $\mathcal {PT}_\mathrm{i}$-symmetric solutions the square modulus is symmetric, which can be easily checked in (19).

The properties of the additional symmetries can also be applied to the already known $\mathcal {PT}_\mathrm{i}$-symmetric and broken states but this only leads to trivial results because the j and k components of these states are always zero.

It must be emphasized that although the analytic continuation allows us to gain a deeper insight into the mathematical structure, the additional states have no physical meaning. In the following the $\mathcal {T}_\mathrm{j}$ symmetry is not further discussed and we again write $\mathcal {PT}$ instead of $\mathcal {PT}_\mathrm{i}$.

At an exceptional point two or more eigenvalues and the corresponding eigenvectors of a non-Hermitian quantum system are degenerate. An exceptional point at which n eigenvalues and eigenvectors coalesce is called an nth-order exceptional point (EPn). The characteristic signature of an EPn is the cyclic permutation of the n eigenvalues and eigenvectors while encircling the exceptional point in the complex parameter space.

Two exceptional points are found in the eigenvalue structure of a BEC in a $\mathcal {PT}$-symmetric double-delta-trap [26, 37] which shows a similar behaviour as the system studied in this paper. Using the delta-potential the point at γ_c2 at which the ground state and the excited state coalesce was found to be a second-order exceptional point (EP2). A second exceptional point arises due to the contact interaction at γ_c1 where the $\mathcal {PT}$-broken states emerge from the ground state. Here three eigenvalues and eigenstates coalesce and the exceptional point has the characteristic behaviour of an EP3. In this section we show that the same exceptional points can also be found in the spatially extended double-well potential. However, in contrast to the above mentioned results here we distinguish between the four components of the continuation and do not use a simplified approach with a projection to two components.

Before treating the exceptional points of the full problem with complex potential and nonlinearity the linear limit (g = 0) is investigated.

4.1. Exceptional point in the linear case

In the linear case the spectrum has only one exceptional point γ_c at which the ground and excited state coalesce and complex solutions arise. The characteristic properties of exceptional points are studied by encircling the exceptional point with a complex parameter.

In the case of the critical point at γ_c we expand the real parameter γ to a complex number. One cycle with radius r is parameterized by

$$\begin{equation} \gamma =\gamma _\mathrm{c} + r\,\mathrm{e}^{\mathrm{i}\vartheta },\qquad \vartheta =0,\ldots ,2\pi . \end{equation} \tag{ 20 }$$

The imaginary part of the perturbation Im γ = rsin ϑ must be interpreted as a real contribution to the potential of the system. Inserting (20) into the potential (2) shows that an additional real term

$$\begin{equation} V_\mathrm{p}(x)=-\mathrm{Im}\, \gamma \, x\,\mathrm{e}^{-\rho x^2} \end{equation} \tag{ 21 }$$

occurs. Thus the imaginary part of γ breaks the $\mathcal {PT}$ symmetry of the potential, because the symmetric real part of the potential obtains an antisymmetric contribution proportional to Im γ.

The behaviour of the energy eigenvalues in the vicinity of the exceptional point is shown in figure 6(a),(c). The two eigenvalues involved in the exceptional point are represented by two surfaces (Riemann sheets). Following a path around the exceptional point leads to a transition from the upper to the lower surface and vice versa. The cycle with fixed radius r = 7 × 10⁻⁴ shown in figure 6(b),(d) starts on the branches with $\mathcal {PT}$-broken solutions with complex eigenvalues at ϑ = 0. Here the eigenvalues have the same real part but a positive or negative imaginary part. At ϑ = π the parameter γ is real again, but on the branches with $\mathcal {PT}$-symmetric solutions on the other side of the exceptional point. After one cycle each state has turned into the respective other state, and after two cycles the initial configuration is restored. This is consistent with the square root behaviour of an EP2.

Zoom In Zoom Out Reset image size

4.2. Exceptional points in the nonlinear case

Encircling the exceptional points is apparently only possible if the number of states is conserved and therefore in the nonlinear case the analytic continuation has to be used. The eigenvalue spectrum shows two exceptional points at γ_c1 and γ_c2. We first discuss the exceptional point at γ_c2 where the ground and excited state coalesce. One cycle around the exceptional point is parameterized with a complexified real part of the gain/loss parameter,

$$\begin{equation} \gamma =\gamma _\mathrm{c2}-r\,\mathrm{e}^{\,\mathrm{j}\vartheta }, \qquad \vartheta =0,\ldots ,2\pi . \end{equation} \tag{ 22 }$$

The path starts at γ = γ_c2 − r in the regime of entirely real eigenvalues. After half a cycle the parameter γ = γ_c2 + r is real again, but on the other side of the exceptional point. At this value the two additional states with nonzero μ₁ and μ_j components exist. On the whole cycle the two components μ_i and μ_k are zero. As can be seen in figure 7 the two eigenvalues and respective eigenvectors interchange after one cycle, and after two cycles the initial configuration is restored. This is the characteristic square root behaviour of an EP2.

Zoom In Zoom Out Reset image size

At the exceptional point at γ_c1 three eigenvalues and eigenvectors become equal. This exceptional point is encircled on two different paths. On the one hand a complexified real part of the nonlinearity parameter is used,

$$\begin{equation} g = g_\mathrm{c1} + r\,\mathrm{e}^{\,\mathrm{j}\vartheta }, \qquad \vartheta =0,\ldots ,2\pi , \end{equation} \tag{ 23 }$$

where the exceptional point resides at the real value g_c1. On the other hand a complexified imaginary part of the gain/loss parameter is used,

$$\begin{equation} \gamma =\gamma _\mathrm{c1}+\mathrm{i}r\,\mathrm{e}^{\,\mathrm{j}\vartheta } \equiv \gamma _\mathrm{c1}+\mathrm{i}\varepsilon , \qquad \vartheta =0,\ldots ,2\pi , \end{equation} \tag{ 24 }$$

which again breaks the symmetry of the real part of the potential. The symmetry-breaking contribution is controlled by the asymmetry parameter ε.

The evolution of the eigenvalues while encircling the exceptional point with the perturbed nonlinearity (23) is shown in the figures 8(a) and (b). Figure 8(a) shows the component μ_i of the eigenvalues in the vicinity of the exceptional point. Each surface corresponds to one of the three eigenvalues. For real values of the parameter g the surfaces can be identified with states of the known eigenvalue spectrum. The middle surface represents the $\mathcal {PT}$-symmetric ground state ψ_g with real eigenvalues. The top (bottom) surface corresponds to the $\mathcal {PT}$-broken state ψ_{i +} (ψ_{i −}) with positive (negative) imaginary part. One recognizes that a path starting on the ground state does not leave its surface, i.e. the ground state does not interchange with another state while encircling the exceptional point. The behaviour of the $\mathcal {PT}$-broken states corresponds to that of an EP2. After one cycle the states are interchanged, and after two cycles the initial configuration is restored. It is known that for certain perturbations an EP3 can show such an EP2 signature [38].

Zoom In Zoom Out Reset image size

To confirm that this point is indeed an EP3 it is necessary to find a perturbation for which the three eigenvalues and eigenvectors interchange cyclically. This behaviour is observed on the path parameterized by (24). Figures 8(c) and (d) show the eigenvalue component μ_k in the proximity of the exceptional point for the whole parameterized region (c) and for a cycle with specific radius (d). Starting on the bottom Riemann sheet, after one counterclockwise cycle the path lies on the middle sheet, after two circles on the top sheet, and finally after three circles on the bottom sheet again. This corresponds to the expected cubic root behaviour of an EP3.

The results presented so far have been obtained with the nonlinear GPE (1). We now want to show that quantitatively very similar results can be obtained from a simple linear matrix model. Our starting point is a nonlinear model system introduced by Graefe et al [29, 39], which shows a very good qualitative agreement with the properties of our system, viz. a two mode approximation of a BEC in a double well. The chemical potential of this approximation of the BEC is given as

$$\begin{equation} \tilde{\mu } = 2(\tilde{\varepsilon } -\mathrm{i}\tilde{\gamma })s_z + 4\tilde{g} s_z^2 + 2v s_x, \end{equation} \tag{ 25 }$$

with the vector ${\boldsymbol s}=(s_x,s_y,s_z)$ the stationary solution (fixed point) of the nonlinear and non-Hermitian Bloch equations

$$\begin{equation} \dot{s}_x = -2\tilde{\varepsilon } s_y - 4\tilde{g} s_y s_z + 4\tilde{\gamma } s_x s_z, \end{equation} \tag{ 26a }$$

$$\begin{equation} \dot{s}_y = 2\tilde{\varepsilon } s_x + 4\tilde{g} s_x s_z -2v s_z + 4\tilde{\gamma } s_y s_z, \end{equation} \tag{ 26b }$$

$$\begin{equation} \dot{s}_z = 2v s_y - \tilde{\gamma } (1 - 4s_z^2). \end{equation} \tag{ 26c }$$

The parameter v defines the coupling between the two wells, and $\tilde{\gamma }$ corresponds to the gain/loss parameter γ of the spatially extended potential. The strength of the nonlinearity g and the symmetry breaking potential ε are represented by $\tilde{g}$ and $\tilde{\varepsilon }$ in the model. We have derived a linear matrix model, whose eigenvalues exactly agree with the solution of the nonlinear model. The 4 × 4 matrix reads

$$\begin{equation} H = \left(\begin{array}{@{}cccc@{}}0 & 1 & 0 & 0 \\\ms v^2-\tilde{\gamma }^2+\tilde{\varepsilon }^2 & 0 & 1 & 0 \\\ms -\dsty\frac{\mathrm{i}\tilde{g}^3 \tilde{\gamma }\tilde{\varepsilon }}{\tilde{g}^2+\tilde{\gamma }^2} & -4\mathrm{i}\tilde{\gamma }\tilde{\varepsilon } + \dsty\dsty\frac{4\tilde{\gamma }^2-2\tilde{g}^2}{\tilde{g}^2 +\tilde{\gamma }^2}\tilde{\varepsilon }^2 & \tilde{g}+\dsty\frac{2\mathrm{i}\tilde{g} \tilde{\gamma }\tilde{\varepsilon }}{\tilde{g}^2+\tilde{\gamma }^2} & 1 \\\ms h_{41} & -\dsty\frac{\mathrm{i}\tilde{g}^3 \tilde{\gamma }\tilde{\varepsilon }}{\tilde{g}^2+\tilde{\gamma }^2} & h_{43} & \tilde{g} \end{array} \right) \end{equation} \tag{ 27 }$$

with the matrix elements

$$\begin{eqnarray*} &&\fl h_{41} \,{=}\, \frac{-\mathrm{i}\tilde{\gamma }\tilde{\varepsilon }}{\tilde{g}^2\! +\! \tilde{\gamma }^2} [\tilde{g}^4+(2\tilde{g}^2+4\tilde{\gamma }^2) (\tilde{\gamma }^2-v^2)]+\frac{2\tilde{\gamma }^2-\tilde{g}^2}{\tilde{g}^2+\tilde{\gamma }^2} [(4\tilde{\gamma }^2+\tilde{g}^2-2v^2)\tilde{\varepsilon }^2 + 2\mathrm{i}\tilde{\gamma }\tilde{\varepsilon }^3 - 2\tilde{\varepsilon }^4],\\ &&\fl h_{43} = \frac{v^2\tilde{\gamma }^2}{\tilde{g}^2+\tilde{\gamma }^2} -\tilde{\gamma }^2 + \frac{2\tilde{g}^2-3\tilde{\gamma }^2}{\tilde{g}^2+\tilde{\gamma }^2} \tilde{\varepsilon }^2. \end{eqnarray*}$$

Details of the derivation are given in the appendix.

The eigenvalues of the linear matrix model (27) are related not only to the $\mathcal {PT}$-symmetric real and the symmetry breaking complex solutions of the GPE (1) but also to states obtained with the analytically continued equations (15). However, in this section we do not compare individual bicomplex components of the chemical potential but use the projections

$$\begin{equation} \mu _+ = (\mu _1 + \mu _\mathrm{k}) + \mathrm{j}\,(\mu _\mathrm{j}- \mu _\mathrm{i}), \end{equation} \tag{ 28a }$$

$$\begin{equation} \mu _- = (\mu _1 - \mu _\mathrm{k}) + \mathrm{j}\,(\mu _\mathrm{j}+ \mu _\mathrm{i}) \end{equation} \tag{ 28b }$$

of the bicomplex solutions on complex numbers, which are formally obtained by the replacement i → ∓j in (18). For vanishing asymmetry parameter $\tilde{\varepsilon }=0$ the eigenvalues $\tilde{\mu }_+$ and $\tilde{\mu }_-$ coincide, and the matrix (27) simplifies to

$$\begin{equation} H = \left(\begin{array}{@{}cccc@{}}0 & 1 & 0 & 0 \\\ms v^2-\tilde{\gamma }^2 & 0 & 1 & 0 \\\ms 0 & 0 & \tilde{g} & 1 \\\ms 0 & 0 & \dsty\frac{v^2\tilde{\gamma }^2}{\tilde{g}^2+\tilde{\gamma }^2} -\tilde{\gamma }^2 & \tilde{g}\end{array} \right), \end{equation} \tag{ 29 }$$

with the eigenvalues

$$\begin{equation} \tilde{\mu } = \left\lbrace \begin{array}{@{}ll@{}}\pm \sqrt{v^2-\tilde{\gamma }^2},\\\ms \tilde{g}\pm \tilde{\gamma }\sqrt{\dsty\frac{v^2}{\tilde{g}^2 +\tilde{\gamma }^2}-1}.\end{array}\right. \end{equation} \tag{ 30 }$$

To compare the matrix model (27) with the numerically exact solutions of the GPE (1) the parameters $\tilde{g}$, $\tilde{\gamma }$, and $\tilde{\varepsilon }$ of the model are replaced with appropriately scaled quantities, $\tilde{g}=g/g_0$, $\tilde{\gamma }=\gamma /\gamma _0$, $\tilde{\varepsilon }=\varepsilon /\gamma _0$, and a shift $\tilde{\mu }=\mu -\mu _0$ of the chemical potential is introduced. Figure 9 compares the two models for the nonlinearities g = 0.1 and 0.2. The values of g₀ = −5.64, γ₀ = 0.953 and v = 0.0426, obtained by adjusting the eigenvalues (30) of the matrix model to the exact results, are independent of the nonlinearity, however, the energy shift μ₀ differs. The eigenvalues are in excellent agreement and thus the matrix model predicts the correct eigenvalues for a wide range of gain/loss and nonlinearity parameters.

Zoom In Zoom Out Reset image size

The matrix model (27) also provides the correct properties of the exceptional points. The eigenvalues μ_± while encircling the exceptional point at γ_c1 along the path given in (24) are shown in figures 10(a) and (b) for the nonlinearity g = 0.2. The eigenvalues are shifted by the mean value of the chemical potential obtained with the four states (including the isolated state not shown in figure 10). The solid and dashed lines present the numerically exact results and the eigenvalues of the matrix model, respectively. Both figures clearly exhibit the permutation of the three states after one cycle in the parameter space indicating the EP3. Note that in the matrix model (27) the exceptional point at γ_c1 is an EP3, whereas in the matrix model in [29] an EP2 is found, i.e. only two instead of three eigenvectors coincide. Heiss et al [37] constructed a three-dimensional matrix model with interacting levels that shows an EP3 at γ_c1 and has three of the eigenvalues given by (30). The four-dimensional matrix model (27) has the advantage that both the EP3 at the pitchfork bifurcation and the EP2 at the tangent bifurcation can be described, and the model even allows us to study interactions between all four eigenmodes in the next section.

Zoom In Zoom Out Reset image size

The degeneracy of two states in an EP2 can be achieved in systems which are controlled by at least two external parameters, and has been observed in a quite large variety of systems (see [40] and references therein). The degeneracy of more than two eigenvalues and eigenvectors is possible, in principle, however, an EPn in general requires the adjustment of (n² + n − 2)/2 parameters [41]. Therefore, the occurrence of higher order exceptional points appears to be more rare. EP3s have been observed in dipolar condensates [36] and in the $\mathcal {PT}$-symmetric condensate discussed here (see sections 4 and 5). Signatures of three interacting eigenmodes related to an EP3 have also been found, e.g. in resonance spectra of the hydrogen atom in crossed magnetic and electric fields [42] and two dielectric microdiscs [43]. Optical wave guides can be arranged such that they contain EP3s [38]. Higher order exceptional points occur in a non-Hermitian $\mathcal {PT}$-symmetric Bose–Hubbard model [44]. However, to the best of our knowledge more than three interacting eigenmodes have not yet been observed in any non-Hermitian quantum system describing a single particle or a BEC with a mean-field approach.

The pitchfork bifurcation at γ_c1 and the tangent bifurcation at γ_c2 merge in the limit g → 0. At $\tilde{g}=0$ and $\tilde{\gamma }=v$ the four eigenvalues in (30) are degenerate. The matrix model (29) at those parameters has the form of a Jordan block of an EP4, viz.

$$\begin{equation} H = J[H] = \left(\begin{array}{@{}cccc@{}}0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0\end{array} \right), \end{equation} \tag{ 31 }$$

which means there is only one eigenvector belonging to the four degenerate eigenvalues $\tilde{\mu }=0$. The structure of the Jordan block in (31) strongly motivates the search for signatures of an EP4, such as the permutation of four states after one cycle on an appropriately chosen parameter path, in both the matrix model and in the system described by the GPE (1).

To reveal the EP4 we have to take into account that the limit g → 0 in the GPE is nontrivial [37]. For $\tilde{g}=0$ the matrix H in (27) has two doubly degenerate branches of eigenvalues $\tilde{\mu }=\pm \sqrt{v^2-(\tilde{\gamma } +\mathrm{i}\tilde{\varepsilon })^2}$. We choose a small but nonzero nonlinearity parameter $\tilde{g}$, which lifts that degeneracy but also causes the branching singularity to split into the tangent and the pitchfork bifurcation. To search for signatures of the EP4 we encircle both branching points simultaneously on paths with sufficiently large radii.

For a path in the complex $\tilde{\gamma }$ plane (with centre $\tilde{\gamma }=v$ and constant $\tilde{\varepsilon }=0$) and for a path in the complex $\tilde{\varepsilon }$ plane (with centre $\tilde{\varepsilon }=0$ and constant $\tilde{\gamma }=v$) the eigenvalues of the matrix H show qualitatively the same behaviour: the eigenvalues move in pairs on two paths where two states permute after one cycle in the parameter space, i.e. they show the typical behaviour of two EP2s related to two tangent bifurcations.

However, the behaviour changes completely for simultaneous perturbations in $\tilde{\gamma }$ and $\tilde{\varepsilon }$. The eigenvalues μ_± (28) of the matrix H with constant parameters g = 0.02 and ε = 0.0002j, and the parameter γ along the path γ = γ₀v + 0.0004 e^jϑ are shown in figures 11(a) and (b) and compared to the numerical results of the GPE. The chemical potential is shifted by the mean value of all four states. The solid lines for the numerically exact results of the GPE agree at least qualitatively with the dashed lines indicating the eigenvalues of the matrix model. In both the matrix model and the GPE the four eigenvalues move on a single path in such a way that they permute after one cycle (ϑ = 0, ..., 2π) in the parameter space, which provides clear evidence for the existence of an EP4 in this system.

Zoom In Zoom Out Reset image size

We have generalized the concept of linear $\mathcal {PT}$ symmetry to nonlinear Gross–Pitaevskii-like systems of the type $\hat{H}_\mathrm{lin}\psi +f(\psi )\psi =\mathrm{i}\dot{\psi }$. If the spectrum is non-degenerate, the linear part $\hat{H}_\mathrm{lin}$ is $\mathcal {PT}$-symmetric and the function f(ψ) fulfils the two conditions,

$$\begin{eqnarray*} f(\mathrm{e}^{\mathrm{i}\chi }\psi ) &=& f(\psi ), \quad \chi \in \mathbb {R},\\ \mathcal {PT}f(\psi ) &=& f(\mathcal {PT}\psi ), \end{eqnarray*}$$

the usual properties of linear $\mathcal {PT}$ symmetry hold, in particular $\mathcal {PT}$-symmetric states and real eigenvalues coincide. These conditions imply that $\mathcal {PT}$-symmetric wave functions render the nonlinear part and thus the Hamiltonian $\mathcal {PT}$-symmetric.

The eigenvalue structure of one special system fulfilling the above-named conditions was studied, namely a BEC with contact interaction in a $\mathcal {PT}$-symmetric double well. The eigenvalue structure shows a tangent and a pitchfork bifurcation at which the number of solutions changes. For further investigation of the bifurcations an analytic continuation to bicomplex numbers was introduced by complexification of the real and imaginary parts of the wave functions. It was shown that this continuation leads to a conserved number of solutions. Additionally it dictates new symmetries to the system and completes the eigenvalue picture. The tangent bifurcation at which the two $\mathcal {PT}$-symmetric solutions coalesce was found to be a second-order exceptional point (EP2). At the pitchfork bifurcation the two $\mathcal {PT}$-broken solutions and one $\mathcal {PT}$-symmetric solution coalesce. Due to studies in a $\mathcal {PT}$-symmetric double-delta-trap potential [37] we expect a third-order exceptional point (EP3) at this bifurcation. To see the true nature of this exceptional point, perturbations breaking the symmetry of the real part of the potential were necessary. Using such perturbations revealed that this exceptional point is indeed an EP3.

We derived a linear 4 × 4 matrix model which eigenvalues coincide with the solutions of the non-Hermitian Bose–Hubbard model in mean-field limit by Graefe et al [29, 39]. It was shown that the eigenvalues of the matrix model are in excellent agreement with the numerical exact solutions of the GPE equation in the spatially extended double-well potential. Also both the EP2 at the tangent bifurcation and the EP3 at the pitchfork bifurcation are described by the matrix model. In the limit of vanishing interaction the matrix model has the Jordan block of an EP4. On an appropriate path the permutation of four states was observed in the matrix model and in the numerical results, thus providing strong evidence for the existence also of an EP4.

Since we showed that $\mathcal {PT}$-symmetric effects should be observable for a certain class of interaction types including the long-range dipolar interaction, such investigations are a starting point for future work as well as for calculations in more complex geometries.

We want to construct a linear matrix model for $\mathcal {PT}$-symmetric condensates with asymmetric perturbation where the eigenvalues of a 4 × 4 matrix H agree with the chemical potential μ in (25) with ${\boldsymbol s}=(s_x,s_y,s_z)$ the fixed points of the nonlinear and non-Hermitian Bloch equations (26) introduced in [39]. (We omit the tildes on the parameters μ, γ, g, and ε in the appendix.) The Bloch equations yield

$$\begin{equation} \mu = g - 2\mathrm{i}\gamma s_z + \frac{\varepsilon}{2s_z} \end{equation} \tag{ A.1 }$$

with s_z the roots of the polynomial

$$\begin{equation} \fl s_z^4 + \underbrace{\frac{g\varepsilon }{g^2+\gamma ^2}}_{a_3} s_z^3 + \underbrace{\frac{\varepsilon ^2+v^2-g^2-\gamma ^2}{g^2+\gamma ^2}}_{a_2} s_z^2 + \underbrace{\frac{-g\varepsilon }{4(g^2+\gamma ^2)}}_{a_1} s_z + \underbrace{\frac{-\varepsilon ^2}{16(g^2+\gamma ^2)}}_{a_0} = 0, \end{equation} \tag{ A.2 }$$

and

$$\begin{eqnarray*} a_3 &=& -(s_{z_1}+s_{z_2}+s_{z_3}+s_{z_4}),\quad a_2 = \sum _{i<j} s_{z_i} s_{z_j},\\ a_1 &=& -\sum _{i<j<k} s_{z_i} s_{z_j} s_{z_k},\quad a_0 = s_{z_1} s_{z_2} s_{z_3} s_{z_4}. \end{eqnarray*}$$

The decisive step is now the construction of the characteristic polynomial

$$\begin{equation} \chi (\mu ) = \mu ^4 + c_3 \mu ^3 + c_2 \mu ^2 + c_1 \mu + c_0 \end{equation} \tag{ A.3 }$$

of the matrix H. The eigenvalues of H coincide with the solutions μ₁ to μ₄ of (A.1) when the conditions

$$\begin{eqnarray*} c_3 &=& -{\rm Tr}\, H = {-}(\mu _1+\mu _2+\mu _3+\mu _4),\quad c_2 = \sum _{i<j} \mu _i \mu _j,\\ c_1 &=& {-}\sum _{i<j<k} \mu _i \mu _j \mu _k,\quad c_0 = \det H = \mu _1 \mu _2 \mu _3 \mu _4 \end{eqnarray*}$$

are fulfilled. For the coefficient c₃ = −(μ₁ + μ₂ + μ₃ + μ₄) we obtain

$$\begin{eqnarray*} c_3 &=& -4 g + 2\mathrm{i}\gamma (s_{z_1}+s_{z_2}+s_{z_3}+s_{z_4}) - \frac{\varepsilon }{2} \left(\frac{1}{s_{z_1}}+\frac{1}{s_{z_2}}+\frac{1}{s_{z_3}} +\frac{1}{s_{z_4}}\right) \\ &\Rightarrow & s_{z_1} s_{z_2} s_{z_3} s_{z_4} c_3 = a_0 c_3 = -2 a_0 (2 g + \mathrm{i}a_3 \gamma ) + \frac{\varepsilon }{2} a_1, \end{eqnarray*}$$

resulting in

$$\begin{equation} c_3 = -2 (2g + \mathrm{i}a_3 \gamma ) + \frac{\varepsilon }{2} \frac{a_1}{a_0} = -2g - \frac{2\mathrm{i}g \gamma \varepsilon }{g^2+\gamma ^2}. \end{equation} \tag{ A.4 }$$

The coefficients c₀ to c₂ are derived in a similar way, though more tediously, and read

$$\begin{eqnarray} &&\fl c_2 = 2(\gamma ^2-v^2) + g^2 + \frac{g^2v^2}{g^2+\gamma ^2} + \frac{2\mathrm{i}\gamma (3g^2+2\gamma ^2)}{g^2+\gamma ^2}\varepsilon - \frac{g^2+2\gamma ^2}{g^2+\gamma ^2}\varepsilon ^2, \end{eqnarray} \tag{ A.5 }$$

$$\begin{eqnarray} &&\fl c_1 = -2g(\gamma ^2-v^2) - \frac{2\mathrm{i}g \gamma (3\gamma ^2+g^2-v^2)}{g^2+\gamma ^2}\varepsilon + \frac{6g \gamma ^2}{g^2+\gamma ^2}\varepsilon ^2 + \frac{2\mathrm{i}g \gamma }{g^2+\gamma ^2}\varepsilon ^3, \end{eqnarray} \tag{ A.6 }$$

$$\begin{eqnarray} &&\fl c_0= (\gamma ^2-v^2)\frac{(g^2+\gamma ^2)^2-\gamma ^2v^2}{g^2+\gamma ^2} + 4\mathrm{i}\gamma (\gamma ^2-v^2)\varepsilon \nonumber \\ &&-\, \frac{2\gamma ^2(g^2+3\gamma ^2-v^2)}{g^2+\gamma ^2}\varepsilon ^2 - \frac{4\mathrm{i}\gamma ^3}{g^2+\gamma ^2}\varepsilon ^3 + \frac{\gamma ^2}{g^2+\gamma ^2}\varepsilon ^4. \end{eqnarray} \tag{ A.7 }$$

It can now be checked directly that the characteristic polynomial of the matrix H given in (27) exactly coincides with (A.3), and thus the eigenvalues of H agree with the solution of the nonlinear model of [39].