Spinor solitons and their $\mathcal{PT}$-symmetric offspring

TL;DR

This paper develops a stable PT-symmetric extension of spinor models in (1+1) dimensions, demonstrating integrability, conservation laws, and stable soliton solutions despite symmetry-breaking perturbations.

Contribution

It introduces a novel PT-symmetric extension of spinor models that avoids instability and shows that these models possess exact, stable soliton solutions with varying conservation properties.

Findings

The PT-symmetric Thirring model is integrable with infinitely many conserved quantities.

The PT-symmetric Gross-Neveu model conserves energy and momentum but not charge.

All three models exhibit stable soliton solutions, with stability near the existence boundary.

Abstract

Although the spinor field in (1+1) dimensions has the right structure to model a dispersive bimodal system with gain and loss, the plain addition of gain to one component of the field and loss to the other one results in an unstable dispersion relation. In this paper, we advocate a different recipe for the $\mathcal{PT}$-symmetric extension of spinor models --- the recipe that does not produce instability of the linear Dirac equation. Having exemplified the physical origins of the $\mathcal P$- and $\mathcal T$-breaking terms, we consider the extensions of three U(1)-invariant spinor models with cubic nonlinearity. Of these, the \PT-symmetric extension of the Thirring model is shown to be completely integrable and possess infinitely many conserved quantities. The \PT-symmetric Gross-Neveu equation conserves energy and momentum but does not conserve charge. The third model is introduced for the purpose of comparison with the previous two; its \PT-symmetric extension has no conservation laws at all. Despite this dramatic difference in the integrability and conservation properties, all three \PT-symmetric models are shown to have exact soliton solutions. Similar to the solitons of the extended Thirring and Gross-Neveu equations, the solitons of the new model are found to be stable --- except for a narrow band of frequencies adjacent to the soliton existence boundary. The persistence under the $\mathcal P$- and $\mathcal T$-breaking perturbations as well as the prevalence of stability highlight a remarkable sturdiness of spinor solitons in (1+1) dimensions.

Figures

Equations

$$i{\gamma }^{\mu }{\partial }_{\mu }\psi +\psi =0.$$

$$\psi=\left(\begin{array}[]{c}u\\ v\end{array}\right),$$

$$\gamma^{0}=\left(\begin{array}[]{cc}0&1\\ 1&0\end{array}\right),\quad\gamma^{1}=\left(\begin{array}[]{cc}0&1\\ -1&0\end{array}\right).$$

$$i(u_t-u_x)+v=0,i(v_t+v_x)+u=0,$$

$$i{u}_{\xi }+v=0,i{v}_{\eta }+u=0,$$

$$\eta =\frac{t+x}{2},\xi =\frac{t-x}{2}.$$

$$\eta \to e^{\beta }\eta ,\xi \to e^{-\beta }\xi ,$$

$$u\to e^{-\beta /2}u,v\to e^{\beta /2}v.$$

$$\mathcal{P}:x^1\to -x^1,u\to v,v\to u,$$

$$\mathcal{T}:x^0\to -x^0,u\to v^{\ast },v\to u^{\ast },$$

$$|u(x,0){|}^2+|v(x,0){|}^2<C(-\infty <x<\infty )$$

$$i{u}_t+u_{xx}+v=0,i{v}_t+v_{xx}+u=0,$$

$$i{u}_t+u_{xx}+v$$

$$i{v}_t+v_{xx}+u$$

$$\mathcal{T}:t\to -t,u\to u^{\ast },v\to v^{\ast }.$$

$$i{u}_t-i{u}_x+v$$

$$i{v}_t+i{v}_x+u$$

$${\omega }^2=1+(k-i\gamma {)}^2.$$

$$i{u}_t-i{u}_x+v$$

$$i{v}_t+i{v}_x+u$$

$$(\omega +\gamma k{)}^2=1+k^2,$$

$$i{u}_t-i{u}_x+v$$

$$i{v}_t+i{v}_x+u$$

$${\omega }^2=1+(k+\gamma {)}^2.$$

$$i(u_t-u_x)+v$$

$$i(v_t+v_x)+u$$

$${\omega }^2=1-{\gamma }^2+k^2.$$

$$u_1=\frac{u-iv}{2},u_2=\frac{v-iu}{2},$$

$$i{u}_{1t}+u_{2x}+u_2$$

$$i{u}_{2t}-u_{1x}+u_1$$

$$i(u_t-u_x)+v=\gamma v,i(v_t+v_x)+u=\gamma u,$$

$$i{u}_{\xi }+v+(Au{v}^{\ast }+Bv{u}^{\ast })v=0,$$

$$i{v}_{\eta }+u+(\tilde{A}v{u}^{\ast }+\tilde{B}u{v}^{\ast })u=0,$$

$$i{u}_{\xi }+(1-\gamma )v+(u{v}^{\ast }+Bv{u}^{\ast })v=0,$$

Full text

1 Introduction

Rooted in the nonhermitian quantum mechanics [1], the $\mathcal{PT}$-symmetric extensions of conservative systems are becoming increasingly relevant in applied disciplines [2]. In optics, the $\mathcal{PT}$-symmetric structures are brought about by the balanced application of gain and loss. Behaviours afforded by the $\mathcal{PT}$-symmetric arrangements and unattainable in standard set-ups, include the unconventional beam refraction [3, 4], loss-induced transparency [5], and nonreciprocal light propagation [6]. $\mathcal{PT}$-symmetric systems are expected to promote an efficient control of light, including all-optical low-threshold switching [8, 7] and unidirectional invisibility [8, 4, 9]. There is also a growing interest in the context of plasmonics [10], optomechanical systems [11] and metamaterials [12].

Recent theoretical analyses of the distributed $\mathcal{PT}$-symmetric systems were exploiting various forms of the nonlinear Schrödinger equation, continuous [13, 14] or discrete [15]. The fundamental object under scrutiny was a particle-like bunch of energy — a soliton, breather or localised internal mode [16]. The present study is concerned with a $\mathcal{PT}$-symmetric extension of another workhorse of the wave theory, namely the nonlinear Dirac equation.

A particularly simple $\mathcal{PT}$-symmetric Schrödinger system consists of two coupled modes, of which one component gains and the other one loses energy at an equal rate [14]. Like this Schrödinger dimer, the (1+1)-dimensional Dirac field consists of two symmetrically arranged components and is ideally suited for the symmetric application of gain and loss. The Lorentz invariance of the Dirac field is an additional built-in symmetry which can be preserved by the gain and loss terms. The two components of the field transform as a Lorentz spinor.

As the authors of the Schrödinger-based studies, we will be focussing on the localised solutions of the nonlinear Dirac equation. That is to say, our interest lies in the $\mathcal{PT}$-symmetric spinor solitons.

The nonlinear Dirac equation with the scalar self-interaction was introduced by Ivanenko [17]; the vector self-interaction is due to Thirring [18]. In the 1950s, Ivanenko [19] and, independently, Heisenberg [20] adopted the former equation as a basis for the unified nonlinear field theory. In the 1970s, Soler tried to use the scalar self-interaction model to describe extended nucleons [21] while Gross and Neveu employed its one-dimensional version to explain the quark confinement [22]. Outside the realm of elementary particles, the Gross-Neveu theory was utilised in the study of polymers [23] while the massive Thirring model appeared in the context of optical gratings [24]. Mathematically, the Thirring model was shown to be completely integrable via the Inverse Scattering Transform [25].

The latest wave of interest in the Dirac equation in the condensed-matter context concerns the electronic structure of two-dimensional materials graphene and silicene [26] as well as the transition metal dichalcogenides [27]. A closely related topic is bosonic evolution in honeycomb lattices [28]. The recent applications of the Dirac equation in optics, are to the light propagation in honeycomb photorefractive lattices (photonic graphene) [29] and conical diffraction in such structures [30]. The spin-orbit coupled Bose-Einstein condensates is yet another area of utilisation of the Dirac-type equations [31]. We should also mention a renewed interest in the stability properties of the Dirac solitons [32, 33, 34, 35, 36] that have been commonly seen as a mystery [37] despite some early progress [38].

There can be a variety of physically meaningful $\mathcal{PT}$-symmetric perturbations of the Dirac equation — that is, perturbations by terms that break each of the parity- and time-reversal symmetries of the equation but remain invariant under the joint action of the $\mathcal{P}\;$and $\mathcal{T}\;$operators. The present study is confined to $\mathcal{PT}$-symmetric perturbations that preserve the invariance under the Lorentz rotations (velocity boosts) and the U(1) gauge transformations.

We scrutinise the general recipe of the $\mathcal{PT}$-symmetric extension of the Dirac equation, explore properties of the extended nonlinear models, construct exact soliton solutions and examine their stability. To crystallise common properties of the $\mathcal{PT}$-symmetric spinor solitons, we consider three different spinor equations together with their $\mathcal{PT}$-symmetric extensions. These include the massive Thirring and Gross-Neveu models as well as a novel spinor model that we introduce for comparison purposes. Like its two rivals, the new model is U(1)-invariant, Lagrangian and has a cubic nonlinearity.

We will show that the $\mathcal{PT}$-symmetric extensions of the three models have different numbers of conservation laws — from infinitely many to none. Despite this difference in the regularity of the dynamics, the three $\mathcal{PT}$-symmetric models will prove to be surprisingly coherent as far as their localised solutions are concerned.

The paper is organised as follows. Section 2 provides a motivation for our choice of the $\mathcal{PT}$-symmetric extension of the linear Dirac equation. Having demonstrated that the plain addition of gain to one mode and loss to the other leads to an unstable equation, we select a particular stable variant of the $\mathcal{P}\;$and $\mathcal{T}\;$-breaking perturbation. This perturbation can also be interpreted as a balanced gain and loss — yet for a different pair of modes.

In section 3, we identify a two-parameter family of cubic spinor systems that are invariant under the full Lorentz group and the U(1) phase transformations, consider its $\mathcal{PT}$-symmetric extension and choose three representatives of this family. These are the massive Thirring and Gross-Neveu models, as well as a new simple spinor equation with cubic nonlinearity.

The physical sources of the $\mathcal{PT}$-symmetric perturbations of these spinor systems are elucidated in section 4. Another aim of that section is to illustrate the occurrence of three particular types of cubic nonlinearity in simple model settings.

The $\mathcal{PT}$-symmetric Thirring model is scrutinised in section 5. We establish that this nonlinear Dirac equation is gauge-equivalent to the original (“parent”) massive Thirring model. Accordingly, the $\mathcal{PT}$-symmetric model represents a completely integrable system and has infinitely many conserved quantities. We derive, explicitly, the first three of these. An explicit expression for the $\mathcal{PT}$-symmetric soliton is also produced in that section.

The following section focusses on the $\mathcal{PT}$-symmetric Gross-Neveu equation. Here our contributions include (a) the local momentum conservation law and (b) an exact explicit soliton solution for this model.

In sections 7 and 8 we derive an exact soliton solution for our novel spinor model. Section 7 deals with the parent system (the equation with no gain and no loss). In this case the soliton is obtained in explicit form. The subsequent section considers the $\mathcal{PT}$-symmetric extension of the model; here the solution is obtained as a quadrature. The stability of the novel solitons is examined in section 9.

Finally, section 10 compares all three spinor models and draws general conclusions.

2 $\mathcal{PT}$-symmetric extension of the Dirac equation

The covariant form of the linear Dirac equation in the free space is

[TABLE]

Here $\partial_{\mu}=\partial/\partial x^{\mu}$, where $x^{\mu}$ is the space-time two-vector, with $x^{0}=t$ denoting the temporal and $x^{1}=x$ spatial coordinate. In $\gamma^{\mu}\partial_{\mu}$, the Einstein summation convention is implied: $\gamma^{\mu}\partial_{\mu}=\gamma^{0}\partial_{0}+\gamma^{1}\partial_{1}$. The $\psi$ is a Lorentz spinor,

[TABLE]

where the components $u(x,t)$ and $v(x,t)$ are complex, and $\gamma^{\mu}$ are the Dirac $\gamma$-matrices. We use the following representation for the $\gamma$-matrices:

[TABLE]

When written in components, the Dirac equation has the form of a system

[TABLE]

or simply

[TABLE]

where we have introduced the light-cone coordinates

[TABLE]

In terms of the light-cone variables, the proper Lorentz transformation has the form

[TABLE]

where $\beta$ is a real boost parameter. (The velocity of the moving reference frame is $\tanh\beta$.) The components of the spinor $\psi$ transform as

[TABLE]

As one can readily check, the system (2.4) is invariant under the spatial reflections

[TABLE]

time reversals

[TABLE]

and the proper Lorentz transformations (2.5)-(2.6). Our aim is to identify a physically meaningful and mathematically consistent $\mathcal{PT}$-symmetric extension of the Dirac equation. The extended system must be invariant under the combined $\mathcal{PT}$ transformations — but not under the $\mathcal{P}$- or $\mathcal{T}$-reflections individually. The physical requirement is that it remain invariant under the Lorentz boosts (2.5)-(2.6) and the U(1) rotations $u\to ue^{i\theta}$, $v\to ve^{i\theta}$ with real constant $\theta$. The added perturbation terms are expected to admit the gain-and-loss interpretation — in some physical contexts, at least. Mathematically, the equation with small perturbation (small gain-loss coefficient) should be stable; that is, all its solutions with initial conditions satisfying

[TABLE]

with some $C>0$, should remain bounded in some norm as $t\to\infty$.

The pair of equations (2.3) bears some similarity with a vector Schrödinger equation

[TABLE]

which is commonly used as a model of the diffractive waveguide coupler [39]. The $\mathcal{PT}$-symmetric extension of the system (2.9) describes the coupler with gain and loss [14]:

[TABLE]

Here $\gamma$ is a positive gain-loss coefficient, $0\leq\gamma\leq 1$ (not to be confused with the $\gamma$-matrices in (2.1) and (2.2)). The system (2.10) is invariant under the $\mathcal{PT}$ transformation where the $\mathcal{P}\;$-operator is given by (2.7) and $\mathcal{T}\;$has the form

[TABLE]

(This transformation is different from (2.8) but still acceptable because $u$ and $v$ are not required to transform as components of a spinor in (2.9)-(2.10).)

Modelling on the Schrödinger dimer (2.10), one could add gain and loss to the Dirac system (2.3):

[TABLE]

However, unlike the Schrödinger system (2.10), this “naive” $\mathcal{PT}$-symmetric extension turns out to have an unstable dispersion relation:

[TABLE]

The instability is caused by the disbalance between gain and loss in (2.11); indeed, (2.11) is not $\mathcal{PT}$ invariant under the spinor transformations (2.7)-(2.8).

There is a whole range of $\mathcal{PT}$-symmetric perturbations of the system (2.3) with stable dispersions — for example, the system

[TABLE]

with the dispersion relation

[TABLE]

or a pair of equations

[TABLE]

with

[TABLE]

While these $\mathcal{PT}$-symmetric systems may be of interest in some physical contexts, in the present study we are focussing on a different extension of the Dirac equation:

[TABLE]

The corresponding dispersion relation is

[TABLE]

Unlike equations (2.13) and (2.15), the relation (2.17) exhibits the symmetry breaking as $\gamma$ exceeds $\gamma_{c}=1$. This suggests that the $\gamma$-terms in (2.16) may account for the gain and loss of energy.

This conjecture turns out to be indeed correct. Defining

[TABLE]

equations (2.16) are transformed into

[TABLE]

According to the representation (2.19), the system comprises two interacting modes, where the $u_{1}$ mode is gaining and $u_{2}$ losing energy at an equal rate $\gamma$. Therefore equation (2.16) is more likely to occur in a situation of the balanced pump and dissipation than equations (2.12) or (2.14).

Another reason for favouring the extension (2.16) over (2.12) and (2.14), is that the $\mathcal{PT}$-symmetric terms in (2.16) preserve and those in (2.12), (2.14) break the Lorentz invariance of the Dirac equation (2.3). This can be readily verified by changing to the light-cone variables and using the transformation rules (2.5)-(2.6).

The model (2.16) has originally appeared in [40] as a $\mathcal{PT}$-symmetric free-fermion quantum theory, outside the gain and loss context. We should also note two recent publications [41] and [42] where a closely related $\mathcal{PT}$-symmetric extension of the Dirac equation was introduced.

For the sake of completeness, we need to mention another stable and $\mathcal{PT}$-symmetric extension of (2.3):

[TABLE]

with the dispersion $\omega^{2}=(1-\gamma)^{2}+k^{2}$. This system is reducible to (2.3) by a trivial coordinate scaling and its $\gamma$-terms do not have any gain and loss interpretation.

3 Nonlinear spinor models

Turning to the nonlinear Dirac equations, we restrict ourselves to considering the simplest, cubic, nonlinearity. The most general cubic Dirac equation that is invariant under the proper Lorentz transformations and U(1) rotations, has the form

[TABLE]

where $A,B,{\widetilde{A}},{\widetilde{B}}$ are complex parameters. If we insist that the system (3.1) be invariant under the spatial reflections, we will have to let ${\widetilde{A}}=A$ and $\widetilde{B}=B$. If the system is required to be invariant under the time reversals, the parameters should satisfy $\widetilde{A}=A^{*}$, $\widetilde{B}=B^{*}$. Therefore the most general cubic U(1)-symmetric spinor system that is invariant under the full Lorentz group, including the $\mathcal{P}\;$- and $\mathcal{T}\;$-transformations, has the form (3.1) with real $A=\widetilde{A}$ and $B=\widetilde{B}$.

Adding the $\mathcal{P}$- and $\mathcal{T}$-breaking terms as in (2.16) we arrive at the $\mathcal{PT}$-symmetric perturbation of the general model (3.1):

[TABLE]

In (3.2) we have scaled the parameter $A$ out. This can be done without loss of generality — except when $A=0$; the latter case has to be considered separately.

The above family of models and their soliton solutions is the topic of our interest in this paper. We will consider three representatives of this family, compare their properties, derive exact expressions for the solitons and examine the soliton stability.

The first representative of the family (3.2) is the $\mathcal{PT}$-symmetric extension of the massive Thirring model, selected by letting $B=0$ in (3.2):

[TABLE]

In the covariant notation, equations (3.3) have the form

[TABLE]

Here $J_{\mu}$ is a two-vector of current: $J_{\mu}={\bar{\psi}}\gamma_{\mu}\psi$; the ${\bar{\psi}}$ stands for the Dirac-conjugate spinor: ${\bar{\psi}}=\psi^{\dagger}\gamma^{0}$, and $\gamma^{5}$ is the matrix defined by

[TABLE]

(We alert the reader to a slight abuse of notation here. In the equation (3.4) — and later in (3.6) and (3.8) — we use the traditional letters $\gamma^{\mu}$ and $\gamma^{5}$ for the Dirac $\gamma$-matrices, whereas $\gamma$ without superscripts is just a scalar — again, a traditional notation for the gain-loss coefficient in literature on $\mathcal{PT}$ symmetry.)

Our second nonlinear spinor system is the $\mathcal{PT}$-symmetric extension of the single-component massive Gross-Neveu model ($B=1$ in the list (3.2)):

[TABLE]

The covariant formulation of the Gross-Neveu equation is

[TABLE]

The third $\mathcal{PT}$-symmetric spinor model we scrutinise in this paper, has the form

[TABLE]

This pair of equations results by setting $A=0$ (and scaling $B$ out) in the general system (3.1) with real $A=\widetilde{A}$ and $B=\widetilde{B}$. Alternatively, one can write $u=B^{-1/2}\widetilde{u}$, $v=B^{-1/2}\widetilde{v}$ and send $B\to\infty$ in the system (3.2). In this sense, equations (3.7) represent the $B=\infty$ special case of (3.2). When written in covariant notation, the system (3.7) reads

[TABLE]

As we have already mentioned, both the massive Thirring and Gross-Neveu models are utilised, extensively, in quantum field theory, condensed matter physics and nonlinear optics. The novel spinor system, equation (3.7), has the nonlinearity as simple as Thirring’s, and this fact suggests that it may also find physical applications. Below, we derive this system in a simple model context.

4 Spinors as amplitudes for counter-propagating waves

In this section we derive our three nonlinear spinor models as equations for the amplitudes of the back- and forward-propagating waves in a medium that supports waves travelling in both directions. The coupling of the two linear waves is achieved by inserting a grating (or two gratings) in the system. (This is not a unique way to produce the coupling; one could alternatively consider a time-periodic parameter variation.) Nonlinear effects also contribute to the coupling.

4.1 Oscillator lattice with periodic grating

As the first prototypical system, we adopt a chain of oscillators coupled, symmetrically, to their left and right nearest neighbours. In the absence of perturbations, this system only involves second-order time derivatives and allows waves propagating in either direction. We are assuming that the atoms in the chain are moving in the external periodic potential (a grating) with the period much larger than the lattice spacing. In the continuum limit, the above discrete system reduces to the Klein-Gordon equation:

[TABLE]

The second last term in (4.1) represents the grating, with the wavelength $\pi$. The first-derivative ($\phi_{t}$) term looks similar, but it has a different physical meaning. This term describes damping with a variable coefficient, changing from positive to negative, and back. Both periodic terms are considered to be small perturbations; accordingly, we have entered a small parameter $\epsilon^{2}$ in front of each of these. Finally, the spatial modulation of the cubic term in (4.1) ensures that the nonlinear terms in the resulting amplitude equations transform as the Lorentz spinors. If we do not enter the $\cos(4x)$ factor, we will end up with a system of amplitude equations where only the linear part is Lorentz-covariant. (It is fitting to note that this is not the only way to achieve the covariance. We could have employed a time-periodic variation of the cubic self-coupling instead.)

We expect the evolution to occur over a hierarchy of space and time scales. Defining

[TABLE]

we denote

[TABLE]

By the chain rule,

[TABLE]

Expanding $\phi$ in powers of the small parameter:

[TABLE]

and substituting, along with the expansions (4.2), in equation (4.1), we equate coefficients of like powers of $\epsilon$.

At the lowest, $\epsilon^{1}$-, order, we have

[TABLE]

We take the solution in the form of a superposition of two counter-propagating waves with equal wavenumbers:

[TABLE]

where $\mathrm{c.c.}$ stands for the complex conjugate of the preceding terms. In (4.3), the amplitudes $u$ and $v$ depend on $T_{1},T_{2},...,$ and $X_{1},X_{2},...$ — but not on $T_{0}$ or $X_{0}$. The factor of $1/\sqrt{3}$ is introduced for later convenience.

The order $\epsilon^{3}$ yields

[TABLE]

where

[TABLE]

Substituting for $\phi_{1}$ from (4.3), the right-hand side in (4.4) becomes a linear combination of resonant and nonresonant harmonics. Setting to zero the coefficients of the resonant harmonics gives the system (3.7),

[TABLE]

where

[TABLE]

The above analysis explains how a particular type of cubic nonlinearity in the model (3.7) can come into being. It also illustrates one possible source of the $\mathcal{PT}$-symmetric perturbation term: a balanced gain and loss of energy in the system.

A shortcoming of the above model system (4.1) is an artificial implementation of the gain and loss. It is not immediately clear how the periodically varied damping coefficient in (4.1) can be realised physically.

4.2 Diatomic chain with periodic coupling

To describe a more realistic source of the $\mathcal{PT}$-symmetric perturbation, we turn to a slightly more complex nonlinear bi-directional medium. This time, the system consists of two chains of oscillators with linear and nonlinear coupling. A common example of such a system is given by a diatomic chain [43]. Confining the consideration to the continuum limit, we write

[TABLE]

Here, the $\epsilon^{2}\cos(2x)$ term describes a weak periodic modulation of the linear inter-chain coupling. (A temporal variation of the coupling produces an equivalent set of amplitude equations.) The two nonlinear terms (proportional to $A$ and $B$, respectively) are introduced for generality. The coefficients $A$ and $B$ are real; depending on the physical setting, one may choose a particular value for each of these.

Note that the $\chi$- and $\phi$-equations in (4.7) have different linear coupling amplitudes, $(1-\gamma)$ vs $(1+\gamma)$. This dissonance breaks the reflection symmetry between the two chains. We will show, however, that the asymmetric coupling preserves the $\mathcal{PT}$-symmetry of the underlying amplitude equations.

Expanding

[TABLE]

and substituting in (4.1), we obtain, at the lowest order of $\epsilon$:

[TABLE]

In (4.8)-(4.9), we use the notation of the previous subsection. We take the solution of (4.8) describing the wave of the unit wavenumber, travelling to the left:

[TABLE]

The solution of (4.9) is taken in the form of a wave travelling to the right:

[TABLE]

The order $\epsilon^{3}$ yields

[TABLE]

where

[TABLE]

Substituting for $\phi_{1}$ and $\chi_{1}$ from (4.10)-(4.11), and setting to zero coefficients of the resonant harmonics $e^{i(\sqrt{2}T_{0}+X_{0})}$ and $e^{i(\sqrt{2}T_{0}-X_{0})}$, we obtain

[TABLE]

with the light-cone variables as defined in (4.6).

Scaling out the parameter $A$ gives the $\mathcal{PT}$-symmetric spinor model (3.2). Therefore the diatomic lattice (4.7) with $B=0$ gives rise to the $\mathcal{PT}$-symmetric Thirring model, the chain with $A=B$ brings about the $\mathcal{PT}$-symmetric Gross-Neveu, while the equations (4.7) with $A=0$ produce the new spinor model (3.7).

Note that the $\mathcal{PT}$-symmetric perturbation term in (4.14) is no longer owing to the damping sign variation. This time, the $\gamma$-term is due to the difference in the linear coupling amplitudes of the $\phi$ and $\chi$ fields. That is, the breaking of the $\mathcal{P}\;$and $\mathcal{T}\;$invariances is caused by the coupling asymmetry.

5 $\mathcal{PT}$-symmetric Thirring model

Having selected three representative $\mathcal{PT}$-symmetric spinor models and exemplified possible sources of the associated nonlinearities, we scrutinise each of the three models individually. We start with the $\mathcal{PT}$-symmetric Thirring model, equation (3.3).

Defining a new pair of the light-cone coordinates by

[TABLE]

and scaling the components of the spinor as in

[TABLE]

casts equations (3.3) in the form

[TABLE]

(where we have dropped the tildes).

The system (5.3) has been considered previously [44, 45]. In particular, the one-soliton solution of (5.3) has been obtained [45]. Using the scaling transformation (5.1)-(5.2) we can readily determine the soliton solution of the $\mathcal{PT}$-symmetric Thirring model (3.3):

[TABLE]

where

[TABLE]

and $\alpha$ is a free parameter, $0<\alpha<\frac{\pi}{2}$.

It is fitting to note that the general $N$-soliton solution of the system (5.3) is also available in literature. This solution is expressible in terms of determinants of $N\times N$ and $(N+1)\times(N+1)$ matrices [45]. Using the scaling (5.1)-(5.2) it is straightforward to obtain the corresponding explicit $N$-soliton solution of the $\mathcal{PT}$-symmetric model (3.3).

The $\mathcal{PT}$-symmetric extension (3.3) is gauge-equivalent to the original (“parent”) Thirring model. Indeed, the following local conservation law follows from (5.3):

[TABLE]

Equation (5.5) implies that there exists a potential $W(\eta,\xi)$ such that

[TABLE]

The gauge transformation [45]

[TABLE]

takes (5.3) to the “parent” Thirring model (equation (3.3) with $\gamma=0$):

[TABLE]

Therefore, the $\mathcal{PT}$-symmetric extension of the Thirring model, equation (3.3), is a completely integrable system — like the original Thirring model itself [25]. This implies, in particular, that equation (3.3) has infinitely many functionally-independent conserved quantities.

The physically-meaningful conservation laws are the electric charge conservation

[TABLE]

where

[TABLE]

energy conservation

[TABLE]

where

[TABLE]

and conservation of momentum

[TABLE]

where

[TABLE]

Note that the above conservation laws cannot be established using the Noether theorem as the $\mathcal{PT}$-symmetric Thirring model (3.3) does not admit a Lagrangian in its $u$ and $v$ variables. We have obtained (5.8), (5.10) and (5.12) by means of the gauge transformation (5.6) and scaling (5.1)-(5.2) from the corresponding conservation laws of the original Thirring model.

We close this section with a comment on stability. The soliton of the original Thirring model (5.7) is stable due to the complete integrability of that equation. (For rigorous stability analysis, see [35].) In view of the gauge equivalence of (5.3) and (5.7), the soliton (5.4) of the extended model is also linearly and nonlinearly stable — regardless of the gain-loss coefficient $0\leq\gamma<1$ and frequency $-\sqrt{1-\gamma^{2}}<\omega<\sqrt{1-\gamma^{2}}$.

6 $\mathcal{PT}$-symmetric Gross-Neveu model

The $\mathcal{PT}$-symmetric extension of the Gross-Neveu model, equation (3.5), was introduced in [41] (though in a different formulation). Like the present study, the earlier investigation focussed on solitons.

6.1 Explicit soliton solution

The soliton solution of the original ($\gamma=0$) Gross-Neveu model is known explicitly [46]. Using a Newtonian path-following algorithm, the authors of [41] continued it to nonzero $\gamma$ and established the domain of existence of the resulting numerical solution. In what follows, we obtain an exact analytical expression for that localised solution. An analytical solution has numerous advantages over its numerical approximation; in particular, the domain of existence of the $\mathcal{PT}$-symmetric Gross-Neveu soliton will be demarcated exactly and explicitly.

Our derivation of the explicit solution exploits two conservation laws of the system (3.5) The authors of Ref [41] observed that the equation conserves energy. It is not difficult to derive the associated flux in the local conservation law (5.10). We have:

[TABLE]

We also establish the conservation of the field momentum. The momentum density and flux in the local conservation law (5.12) have the form

[TABLE]

To construct the soliton, we decompose

[TABLE]

Substituting in (6.1)-(6.2) and assuming that $a(x),b(x)\to 0$ as $|x|\to\infty$, we establish two useful relations:

[TABLE]

where

[TABLE]

Taking a product of (6.4) and (6.5) we obtain

[TABLE]

This equation implies that $\gamma^{2}$ has to be smaller than $(1+\sigma/2)^{2}$ for all $x$, including $x=\pm\infty$, where $\sigma=0$. Therefore, the solution that we are going to construct, will be valid for $\gamma^{2}\leq 1$.

Substituting (6.3) in (3.5) gives an equation for $a$,

[TABLE]

a similar equation for $b$, and two equations for the phases of the fields:

[TABLE]

(In obtaining (6.8), we made use of (6.4) and (6.5).) Comparing (6.8) to (6.4) and (6.5), we observe that the function $\theta(x)$ is monotonically decreasing and $\varphi(x)$ monotonically growing.

Differentiating the relation (6.4) in $x$ and using (6.7), we arrive at

[TABLE]

This equation implies that regions of growth of $\sigma(x)$ correspond to $\omega\sin(\varphi-\theta)>0$ whereas regions of decay are those where $\omega\sin(\varphi-\theta)<0$.

The factor $\sin(\varphi-\theta)$ in the right-hand side of (6.9) is determined, up to a sign, by $\cos(\varphi-\theta)$ which, in turn, can be written as

[TABLE]

With the help of (6.6), this relation gives

[TABLE]

The numerator in (6.10) needs to be smaller than the denominator for all $x$ — in particular, for $x=\pm\infty$. Accordingly, the parameters $\omega$ and $\gamma$ have to be constrained by $\omega^{2}+\gamma^{2}\leq 1$.

Equations (6.6) and (6.10) can be used to express the right-hand side in (6.9) in terms of a single variable, $\sigma(x)$. This converts (6.9) to a pair of simple separable equations:

[TABLE]

where

[TABLE]

One of these equations is valid in the region of growth of $\sigma(x)$ and the other one is valid in the complementary region of its decay.

The compatible nonsingular solution of (6.11a) and (6.11b), approaching zero as $|x|\to\infty$, is

[TABLE]

where

[TABLE]

The absolute values of the $u$ and $v$ components are obtained from (6.4)-(6.5),

[TABLE]

From equations (6.11m) it is clear that $\omega$ has to be positive. Therefore, the admissible range of $\omega$ is $0<\omega<\sqrt{1-\gamma^{2}}$ for each $0\leq\gamma<1$.

The phase variables are obtained from (6.8), by integration:

[TABLE]

where

[TABLE]

The solution (6.11l) satisfies equation (6.11a) in the region $x<0$ and equation (6.11b) in the region $x>0$. Accordingly, when recovering $\theta(x)$ and $\varphi(x)$ from (6.8), the integration constants were chosen so that $\sin(\varphi-\theta)<0$ in the region $x<0$ and $\sin(\varphi-\theta)>0$ in the region $x>0$.

The soliton of the original Gross-Neveu model ($\gamma=0$) is known to be stable for all values of its frequency $\omega$. For the analytical proof using the Evans function, see [33]; the comprehensive numerical study is in [34, 36]. The soliton of the $\mathcal{PT}$-symmetric extension (3.5) was also found to be stable — for all $\gamma$ and $\omega$ [41]. (The analysis of [41] appealed to the stability eigenvalues of the numerically-determined soliton.) Accordingly, we conclude that our explicit soliton solution (6.3)+(6.11m)+(6.11n) is stable regardless of the gain-loss coefficient $0\leq\gamma<1$ and frequency $0<\omega<\sqrt{1-\gamma^{2}}$.

6.2 Charge nonconservation

Despite their common possession of explicit soliton solutions, there is an important difference between the $\mathcal{PT}$-symmetric Thirring and Gross-Neveu models. Whereas the Thirring model has an infinity of conservation laws, in the case of the Gross-Neveu equation we were unable to determine any other conserved quantities in addition to energy and momentum. We therefore conjecture that the system (3.5) with $\gamma\neq 0$ has only two conservation laws, equations (6.1) and (6.2).

7 Novel spinor model

The system (3.7) is the third on our list of spinor models. Since the model is new, we start with its “original”, i.e. $\gamma=0$, version:

[TABLE]

The model (6.11a) admits a Lagrangian, with the density

[TABLE]

or, in the covariant formulation,

[TABLE]

The charge, energy and momentum conservation laws are straightforward by means of the Noether theorem. The local charge is governed by equation (5.8) with

[TABLE]

The energy and momentum conservation laws have the form (5.10) and (5.12), where

[TABLE]

Proceeding to the soliton solutions of the equations (6.11a), we consider stationary (nonpropagating) solitons of the form

[TABLE]

The local charge conservation (5.8) + (6.11c), together with the vanishing boundary conditions $u,v\to 0$ as $|x|\to\infty$, requires $|f|=|g|$. Hence we write

[TABLE]

Substituting in (6.11a) we observe that

[TABLE]

so that $\theta+\varphi$ is a constant. Using the U(1) invariance of the model we can set this constant to zero. The resulting system has the form

[TABLE]

The momentum conservation (5.12)+(6.11d), together with the vanishing boundary conditions, furnishes an invariant manifold of the system (6.11g)-(6.11h):

[TABLE]

According to (6.11i), $\cos 2\theta(\pm\infty)=\omega$; hence $|\omega|\leq 1$. Substituting for $a^{2}$ in (6.11h) gives a simple separable equation

[TABLE]

with two solutions,

[TABLE]

and

[TABLE]

Here

[TABLE]

The above solutions should be filtered using the equation

[TABLE]

which is a consequence of (6.11i) and (6.11j). Starting with the monotonically growing solution (6.11k), we observe that the corresponding $\cos 4\theta$ passes through $1$ as $x$ goes through the origin. According to (6.11m), this solution produces negative $a^{2}$ and should be discarded.

On the other hand, the function (6.11l) is monotonically decreasing from $\arctan\lambda$ to $-\arctan\lambda$. The right-hand side of (6.11m) is nonnegative provided $|4\theta|\leq\pi/2$. Therefore $\arctan\lambda$ must not exceed $\pi/8$. This gives the range of admissible $\omega$:

[TABLE]

The absolute value of the $u$ and $v$ components is determined from (6.11i) and (6.11l):

[TABLE]

The function (6.11n) with $\omega$ in the range $\frac{3}{4}\leq\omega<1$ is unimodal (bell-shaped). In the remaining part of the parameter interval, $\frac{1}{\sqrt{2}}<\omega<\frac{3}{4}$, the function $a(x)$ has two humps placed at $x=\pm x_{m}$, where $x_{m}$ is the positive root of

[TABLE]

As $\omega\to 1/\sqrt{2}$, the two humps diverge to infinities: $x_{m}\to\infty$.

Equations (6.11e), where $f=ae^{i\theta}$ and $g=ae^{-i\theta}$, while $a$ and $\theta$ are as in (6.11n) and (6.11l), provide an explicit soliton solution to the new spinor model (6.11a).

8 $\mathcal{PT}$-symmetric extension of the new model: solution by quadrature

Finally, we consider the $\mathcal{PT}$-symmetric extension of the model (6.11a):

[TABLE]

Neither the local charge (5.8)+(6.11c) nor the energy-momentum conservation laws (5.10)+(6.11d) and (5.12)+(6.11d) persist as $\gamma$ is taken away from zero. In fact we were unable to establish any conservation laws for this $\mathcal{PT}$-symmetric spinor equation. We conjecture that there aren’t any.

The lack of conservation laws deprives us of prior knowledge of the invariant manifold that harbours the homoclinic trajectory in the four-dimensional phase space of the stationary system. As a result, we will be able to construct the exact soliton solution as a quadrature — but not explicitly.

8.1 Invariant manifold

Assuming stationary solutions of the form (6.11e) and letting

[TABLE]

equations (6.11a) reduce to a four-dimensional stationary system

[TABLE]

where

[TABLE]

Here we assume that $a(x),b(x)\to 0$ as $|x|\to\infty$.

The system (6.11c)-(6.11f) can be conveniently analysed using the Stokes vector ${\vec{\mathcal{R}}}=(X,Y,Z)$, where

[TABLE]

(For review and references on the Stokes coordinates see [47].) When transformed to the Stokes variables, equations (6.11c), (6.11d) and a combination of (6.11e)-(6.11f) form a self-contained dynamical system in three dimensions:

[TABLE]

Here the overdot indicates derivative w.r.t. $x$.

Another linear combination of (6.11e) and (6.11f) is a stand-alone equation for the variable $\beta=\theta+\varphi$:

[TABLE]

Once $X,Y$ and $Z$ have been determined, $\beta$ can be obtained from (6.11l) by simple integration.

The quantity $\mathcal{R}$ is governed by equation

[TABLE]

which is a consequence of (6.11i)-(6.11k). From (6.11i) and (6.11m) we obtain a simple separable equation

[TABLE]

The solution curve satisfying the initial condition $\left.X\right|_{\mathcal{R}=0}=0$, is

[TABLE]

Equation (6.11o) gives an invariant manifold of the system (6.11i)-(6.11k) containing the homoclinic trajectory that we are trying to determine.

8.2 High-frequency soliton ($\omega^{2}>1/2$)

Assuming, first, that

[TABLE]

we define $\rho>0$ such that

[TABLE]

The implicit curve (6.11o) can be parametrised by letting

[TABLE]

where $\chi=\chi(x)$ is a real parameter. Equations (6.11q) and (6.11r) constitute one of the two branches of the hyperbola (6.11o) — specifically, the branch with $\mathcal{R}<2\omega$. (The branch with $\mathcal{R}>2\omega$ is considered in subsection 8.3.)

Substituting (6.11q) and (6.11r) into (6.11m) we obtain

[TABLE]

Comparing this to (6.11k) gives

[TABLE]

where $\chi_{0}$ is the value of $\chi$ attained as $|x|\to\infty$:

[TABLE]

The first equation in (6.11u) infers that $\omega$ has to be positive. The second one tells us that $\chi_{0}>0$.

Differentiating (6.11s) with respect to $x$ and using (6.11j) gives an equation of motion of a fictitious classical particle with the coordinate $\chi$:

[TABLE]

Here the potential $U(\chi)$ can be cast in the form

[TABLE]

(We have chosen the zero of the potential to be at $\chi_{0}$: $U(\chi_{0})=0$.) The first integral of equation (6.11v) is a sum of the kinetic and potential energy of the particle:

[TABLE]

where we have taken into account that $\chi=\chi_{0}$ is an equilibrium.

The potential $U(\chi)$ has a double zero at $\chi=\chi_{0}$ and two simple zeros, at $\chi_{1}$ and $\chi_{2}$. One can readily check that $\chi_{0}$ is a point of maximum if

[TABLE]

In what follows we assume that the inequality (6.11y) is satisfied — for if it were not, the particle would oscillate about the minimum of $U(\chi)$ and the corresponding solution $\chi(x)$ would not be localised.

Since $U(-\chi_{0})>0$, the other zeros are to the left of $\chi_{0}$ in this case: $\chi_{2}<-\chi_{0}<\chi_{1}<\chi_{0}$. This implies, in particular, that

[TABLE]

and

[TABLE]

The point $\chi_{1}$ is to the left of the origin if $U(0)<0$ and to the right of the origin if $U(0)>0$. Here

[TABLE]

It is not difficult to check that when $\gamma$ is smaller than a certain $\gamma_{*}$, there are $\chi_{a}$ and $\chi_{b}$, $0<\chi_{a}<\chi_{b}$, such that the quantity (6.11ab) is negative in the interval $\chi_{a}<\chi_{0}<\chi_{b}$ and positive outside it. On the other hand, when $\gamma>\gamma_{*}$, the expression (6.11ab) is positive for all $\chi_{0}$. The critical value $\gamma_{*}$ is given by

[TABLE]

where $\mathcal{X}=1.374$ is the positive root of the equation

[TABLE]

Since $\omega=\frac{1}{\sqrt{2}}\,\mathrm{coth}\,\chi_{0}$, we conclude that for small $\gamma<\gamma_{*}$, there are $\omega_{a}=\omega_{a}(\gamma)$ and $\omega_{b}=\omega_{b}(\gamma)$, $1/\sqrt{2}<\omega_{b}<\omega_{a}$, such that the point $\chi_{1}$ lies to the left of the origin if $\omega_{b}<\omega<\omega_{a}$ and to the right of the origin if $\omega$ is outside the interval $(\omega_{b},\omega_{a})$.

There is a homoclinic trajectory connecting the saddle $\chi=\chi_{0}$, ${\dot{\chi}}=0$ to itself. The fictitious particle following this trajectory leaves the unstable equilibrium in the infinite past, reaches $\chi=\chi_{1}$ at “time” $x=0$, and returns to $\chi_{0}$ in the infinite future ($x\to\infty$).

The variable (6.11q) pertaining to this trajectory has one of two possible behaviours depending on the values of $\gamma$ and $\omega$. Frequencies $\omega$ lying outside the interval $(\omega_{b},\omega_{a})$ correspond to bell-shaped functions $\mathcal{R}(x)$ decreasing from

[TABLE]

to zero as $|x|$ changes from 0 to infinity. On the other hand, frequencies satisfying $\omega_{b}<\omega<\omega_{a}$ correspond to bimodal (double-humped) functions $\mathcal{R}(x)$. As $|x|$ grows from zero in the latter case, $\mathcal{R}$ increases from the value (6.11ae) to $\sqrt{2}\rho(\cosh\chi_{0}-1)$ and only then decays to zero. Due to (6.11aa), we have $\mathcal{R}(x)>0$ for all $x$ ($-\infty<x<\infty$); hence $\mathcal{R}(x)$ does represent the magnitude of the vector $(X,Y,Z)$. Accordingly, the homoclinic trajectory defines a localised solution of the system (6.11i)-(6.11k) for all $\frac{1}{\sqrt{2}}<\omega<\sqrt{1-\gamma^{2}}$. The domain of existence on the $(\gamma,\omega)$ plane is illustrated in Fig 1 (b).

Note that $d\chi/dx<0$ in the region $x<0$ and $d\chi/dx>0$ in $x>0$. Keeping this correspondence in mind, we integrate (6.11x) to obtain

[TABLE]

Equation (6.11af) defines the function $\chi(x)$ over the entire real line $-\infty<x<\infty$. The function $\chi(x)$ is even. As $|x|\to\infty$, we have $\chi(x)\to\chi_{0}$; at the origin, $\chi(0)=\chi_{1}$.

Having determined the function $\chi(x)$, equations (6.11q) and (6.11t) can be used to reconstruct the moduli of the $u$ and $v$ components of the spinor soliton:

[TABLE]

To reconstruct the corresponding phase variables, we need to determine their linear combinations, $\alpha$ and $\beta$. Equation (6.11r) reads $X=\rho(\sinh\chi_{0}-\sinh\chi)$; the inequality (6.11z) implies then $X(x)>0$ for all $-\infty<x<\infty$. Because of that, the angle $\alpha$ can be taken to lie between $-\pi/2$ and $\pi/2$ and we can let $\alpha=\arctan(Y/X)$. Consequently,

[TABLE]

where $U(\chi)$ is as in (6.11w). In the above expression, we have taken into account that $Y(x)=\frac{1}{\sqrt{2}}\dot{\chi}$ is positive and negative for $x>0$ and $x<0$, respectively.

The function $\beta(x)$ can be found from (6.11l) by integration:

[TABLE]

where

[TABLE]

Both $\alpha(x)$ and $\beta(x)$ are odd functions, bounded as $x\to\pm\infty$.

Once $\alpha$ and $\beta$ have been determined, the phases of $u$ and $v$ are found simply as

[TABLE]

8.3 Solitons with negative frequencies?

We return to equation (6.11o) with $2\omega^{2}-1>0$ and consider the second branch of this hyperbola. Instead of equation (6.11q), the corresponding $\mathcal{R}$-component is given by

[TABLE]

while the $X$-component is given by equation (6.11r), as before. This time, the parameter $\chi$ satisfies

[TABLE]

and $Z$ is found to be

[TABLE]

where $\chi_{0}$ is the value of $\chi$ attained as $|x|\to\infty$:

[TABLE]

The equations (6.11an) imply that $\omega$ has to be negative this time while $\chi_{0}$ remains positive.

Differentiating (6.11al) in $x$ and using (6.11j), we arrive at the same Newton’s equation (6.11v) as in the analysis of the hyperbola branch given by (6.11q)-(6.11r). The potential energy of the fictitious particle is given by the same equation (6.11w) as before. As we have established, there is a homoclinic trajectory connecting the saddle $\chi=\chi_{0},\dot{\chi}=0$ to itself. (We assume that inequality (6.11y) is in place.) The fictitious particle following this trajectory moves from its equilibrium position at $\chi_{0}>0$ to the point $\chi_{1}<\chi_{0}$ and then returns to $\chi_{0}$.

This time, the homoclinic trajectory does not furnish a localised solution of the system (6.11i)-(6.11k) though. Indeed, in view of the inequality (6.11aa), the variable

[TABLE]

remains negative along the entire trajectory and cannot represent the magnitude of the vector $(X,Y,Z)$. We conclude that the system (6.11a) does not have solitons with frequencies $\omega<-\frac{1}{\sqrt{2}}$.

8.4 $\gamma=0$ reduction of quadrature

It is instructive to follow the transformation of the quadrature (6.11af)-(6.11aj) to the explicit solution (6.11l), (6.11n) once $\gamma$ is set to zero. When $\gamma=0$, we have $Z=0$ so that $a=b$ and $\cos\alpha=X/\mathcal{R}$. With the help of (6.11u), equations (6.11q) and (6.11r) give

[TABLE]

Using (6.11ao) the integration over $\chi$ in (6.11af) can be changed to integration over $\alpha$:

[TABLE]

In transforming (6.11af) to (6.11ap), we made use of

[TABLE]

Here the top sign corresponds to the region $x>0$ (where $\sin\alpha>0$) and the bottom sign to $x<0$ (where $\sin\alpha<0$). We have also used the relation

[TABLE]

When $\gamma=0$, we have $\beta=0$ and so $\alpha=-2\theta$. With this observation, equation (6.11ap) is nothing but an integral of the separable equation (6.11j) with an explicit solution (6.11l).

8.5 Low-frequency soliton ($\omega^{2}<1/2$)

We proceed to the situation $2\omega^{2}<1$ and define $\rho>0$, where

[TABLE]

In this case the hyperbola (6.11o) admits a unique parametrisation consistent with the boundary condition $X(\pm\infty)=0$:

[TABLE]

As in subsection 8.2, $\chi$ satisfies (6.11s) and $Z$ satisfies equation (6.11t) where $\chi_{0}$ is the value of $\chi$ attained as $|x|\to\infty$:

[TABLE]

Equations (6.11at) allow both signs of $\omega$: positive $\omega$’s correspond to $\chi_{0}>0$ and negative $\omega$’s correspond to $\chi_{0}<0$.

Differentiating (6.11s) in $x$ we obtain the Newton’s equation (6.11v) where the potential

[TABLE]

The point $\chi_{0}$ is a point of maximum of $U(\chi)$ if $\omega$ and $\gamma$ satisfy the inequality (6.11y). This is the first necessary condition for the existence of the homoclinic orbit.

When $\gamma=0$, the point $\chi=\chi_{0}$ is the only maximum of the potential (6.11au). The potential decreases monotonically in either direction away from $\chi_{0}$. Assume we now keep $\rho$ unchanged and raise $\gamma$. As $\gamma$ reaches a certain $\gamma_{\mathrm{max}}>0$, the potential function develops the second maximum at the point $\chi=\chi_{\mathrm{max}}$. A simple graphical analysis of the derivative

[TABLE]

indicates that the point $\chi_{\mathrm{max}}$ is on the right of $\chi_{0}$ when $\chi_{0}<0$ and on the left of $\chi_{0}$ when $\chi_{0}>0$.

As $\gamma$ exceeds a critical value $\gamma_{c}>\gamma_{\mathrm{max}}$, the second maximum reaches above zero: $U(\chi_{\mathrm{max}})>0$. In the parameter region $\gamma>\gamma_{c}$, the function $U(\chi)$ has two simple zeroes in addition to the double zero at $\chi_{0}$: $U(\chi_{1,2})=0$. When $\chi_{0}<0$, we have $\chi_{0}<\chi_{1}<\chi_{2}$, and when $\chi_{0}>0$, the arrangement is $\chi_{2}<\chi_{1}<\chi_{0}$. It is not difficult to realise that in the region $\gamma>\gamma_{c}$, the equation (6.11v) has a homoclinic orbit. However, this does not necessarily mean that the three-dimensional dynamical system (6.11i)-(6.11k) has one.

Indeed, let $\chi_{0}<0$. As $x$ changes from minus infinity to zero, $\chi(x)$ grows from $\chi_{0}$ to $\chi_{1}$. Equation (6.11ar) implies then that the corresponding $\mathcal{R}(x)$ is negative; this disqualifies the choice $\chi_{0}<0$. As a result, there are no solitons with $\omega<0$.

In contrast, solitons with $\omega>0$ do exist. In this case $\chi_{0}>0$; as $x$ varies from $-\infty$ to [math], the parameter $\chi$ decreases from $\chi_{0}$ to $\chi_{1}$. According to (6.11ar), the corresponding $\mathcal{R}(x)$ grows, monotonically, from [math] to its maximum value

[TABLE]

Therefore $\mathcal{R}$ does give the length of the vector $(X,Y,Z)$ and the homoclinic orbit of (6.11v) defines a localised solution of (6.11i)-(6.11k).

Note that, unlike the high-frequency solitons considered in section 8.2, the solitons in the range $0<\omega<\frac{1}{\sqrt{2}}$ all have a unimodal function $\mathcal{R}(x)$.

For each $\chi_{0}>0$, the critical value $\gamma_{c}=\gamma_{c}(\omega)$ can be determined as a root of the system

[TABLE]

Substituting from (6.11au) this system acquires the form

[TABLE]

where $\xi=(\chi-\chi_{0})/2$. Eliminating $\gamma$ between (6.11av) and (6.11aw), we arrive at a simple transcendental equation

[TABLE]

where

[TABLE]

The subscript $\omega$ in $G_{\omega}$ serves to remind that there is a one-to-one correspondence between $\omega$ and $\chi_{0}$: $\tanh\chi_{0}=\sqrt{2}\omega$.

The even function $F(\xi)$ has a single extremum (a minimum) at $\xi=0$ and grows to infinity as $|\xi|\to\infty$. The function $G_{\omega}(\xi)$ also has a single extremum (a maximum) at some $\xi<0$, and approaches $1/2$ as $|\xi|\to\infty$. The functions $F$ and $G_{\omega}$ intersect at $\xi=0$ where $F(0)=G_{\omega}(0)=1$. At the point of intersection, $G_{\omega}(\xi)$ has a negative slope,

[TABLE]

whereas $\left.dF/d\xi\right|_{\xi=0}=0$. Therefore for each $\chi_{0}>0$, there is one more intersection, at $\xi=\xi_{c}<0$. See Fig 1 (a). The point $\xi_{c}$ depends on $\omega$.

Having computed $\xi_{c}$ for a sample of $\omega$ values by means of the standard Newtonian iteration, we use (6.11av) to determine the corresponding $\gamma_{c}(\omega)$:

[TABLE]

The inverse function, $\omega=\omega_{c}(\gamma)$, gives the lower boundary of the soliton’s domain of existence; see Fig 1(b).

Note that equation (6.11az) implies that $\gamma_{c}>\rho$. This furnishes a simple lower bound on the function $\omega_{c}(\gamma)$:

[TABLE]

Thus we have established the domain of existence of solitons with frequencies $\omega<\frac{1}{\sqrt{2}}$. Specifically, the system (6.11i)-(6.11k) has a localised solution provided $\omega$ lies between $\omega_{c}(\gamma)$ and $\sqrt{1-\gamma^{2}}$. The domain of existence of the low-frequency solitons seamlessly adjoins the domain of solitons with $\omega>\frac{1}{\sqrt{2}}$; see Fig 1 (b).

The soliton is expressible in terms of the function $\chi(x)$ given by the quadrature (6.11af), with $U(\chi)$ as in (6.11au). Once the function $\chi(x)$ has been determined, we recover $a(x)$ and $b(x)$ using (6.11ar) and (6.11t):

[TABLE]

The function $\alpha(x)$ is determined from

[TABLE]

where $U$ is as in (6.11au). The function $\beta(x)$ can be found from (6.11l) by integration:

[TABLE]

where

[TABLE]

Both $\alpha(x)$ and $\beta(x)$ are odd functions, bounded as $x\to\pm\infty$. Once $\alpha(x)$ and $\beta(x)$ have been constructed, the phases $\theta(x)$ and $\varphi(x)$ are determined from (6.11aj).

9 Stability of solitons in the new model

To classify the stability of the soliton of the new spinor model and its $\mathcal{PT}$-symmetric extension, we linearise equations (6.11a) and (6.11a) about the solution (6.11e). Choosing perturbations of the form

[TABLE]

gives an eigenvalue problem

[TABLE]

Here ${\bf z}=(z_{1},z_{2},z_{3},z_{4})^{T}$; the operator $H$ is defined by

[TABLE]

and $J$ is a diagonal constant matrix

[TABLE]

In (6.11k)-(6.11l), $\sigma_{1,2,3}$ are the Pauli matrices, and $\sigma_{0}$ is the $2\times 2$ identity matrix.

The continuous spectrum of $\lambda$ lies on the imaginary axis and consists of two branches. The first branch has a narrow gap, $|\mathrm{Im}\,\lambda_{1}|\geq\sqrt{1-\gamma^{2}}-\omega$. The second branch’s gap is wider: $|\mathrm{Im}\,\lambda_{2}|\geq\sqrt{1-\gamma^{2}}+\omega$.

We approximate the boundary conditions in (6.11a) by ${\bf z}(\pm L)=0$; the bulk of our calculations was done with $L=40$. The Chebyshev differentiation on a nonuniform mesh with $N=1200$ nodes converts (6.11a) to an eigenvalue problem for a $4N\times 4N$ matrix. The matrix eigenvalues are then computed using a standard numerical routine.

The soliton solution (6.11ag)-(6.11aj) was examined on a grid of $\gamma$ and $\omega$ parameters covering the two-dimensional domain $0<\gamma<1$, $\frac{1}{\sqrt{2}}<\omega<\sqrt{1-\gamma^{2}}$. The soliton (6.11ba)-(6.11bc) was similarly studied by sampling on $0<\gamma<1$ and $\omega_{c}(\gamma)<\omega<\frac{1}{\sqrt{2}}$. We paid special attention to the explicit solution (6.11l)+(6.11n), with $\omega$ changing from 1 to $1/\sqrt{2}$. The solution with a particular choice of $\gamma$ and $\omega$ was classified as unstable if the spectrum contained eigenvalues with $\mathrm{Re}\,\lambda>0.003$.

The soliton’s stability properties were found to be qualitatively similar for all $0\leq\gamma<1$. Assume the soliton frequency $\omega$ is decreased while the gain-loss coefficient $\gamma$ is kept constant. As $\omega$ passes through a certain $\omega_{\mathrm{inst}}=\omega_{\mathrm{inst}}(\gamma)$, a complex quadruplet $\pm\lambda,\pm\lambda^{*}$ bifurcates from $\pm i(\sqrt{1-\gamma^{2}}+\omega)$, the edges of the wider gap in the continuous spectrum. As $\omega$ is further decreased towards $\omega_{c}(\gamma)$, the lower boundary of the soliton’s existence domain, the real parts of $\lambda$ grow in absolute value (but never exceed $0.3$). The instability band $\omega_{c}(\gamma)<\omega<\omega_{\mathrm{inst}}(\gamma)$ is quite narrow; see Fig.1(b) where it is tinted brown.

The conclusions of the spectral analysis were found to be consistent with the direct numerical simulations of the equations (6.11a) and (6.11a). In either case the initial condition was chosen in the form of the corresponding soliton solution. Figure 2 illustrates the growth of the instability of the soliton with $\omega$ near the boundary of its existence domain, with $\gamma=0$ and $\gamma\neq 0$. The simulations were performed using Lakoba’s method of characteristics [36].

10 Concluding remarks

The key results of our paper can be summarised as follows.

(1) We have formulated a recipe for the $\mathcal{PT}$-symmetric extension of spinor models that is consistent with the relativistic invariance and gain-loss interpretation of the $\mathcal{P}\;$- and $\mathcal{T}\;$-breaking terms.

(2) The $\mathcal{PT}$-symmetric extension of the massive Thirring model was shown to be gauge equivalent to the original Thirring model. This proves that the $\mathcal{PT}$-symmetric model is a completely integrable system and has infinitely many conserved quantities. We have derived, explicitly, the first three of these and produced an explicit expression for the soliton solution.

(3) We have established the local momentum conservation law for the $\mathcal{PT}$-symmetric Gross-Neveu equation and obtained an exact explicit soliton solution for that model.

(4) A novel nonlinear Dirac equation was introduced, along with its $\mathcal{PT}$-symmetric extension. We have determined an exact soliton solution of the new model. In the $\gamma=0$ case this solution is explicit while in the $\gamma\neq 0$ sector it is obtained as a quadrature. The soliton was found to be stable in most of its domain of existence in the $(\gamma,\omega)$ parameter plane. The instability is only present in a narrow strip along the boundary of this domain.

The central message of this study is that of a remarkable ubiquity as well as structural and dynamical stability of spinor solitons. Solitons are supported by the Lorentz-invariant Dirac equations with a broad range of cubic nonlinearities. They persist under the addition of the $\mathcal{P}\;$- and $\mathcal{T}\;$-breaking $\mathcal{PT}$-symmetric terms. No matter how orderly or disorderly the $\mathcal{PT}$-symmetric extension is — whether it has infinitely many conservation laws or none — the solitons are expressible in exact analytic form. Finally, the spinor solitons are stable in the Lyapunov sense — either in the entirety or in the vast majority of their parameter domain, both in the original model and in its $\mathcal{PT}$-symmetric extension.

We close this section with two remarks.

The first one concerns travelling solitons. In this paper, we have restricted ourselves to considering the relativistically invariant nonlinear Dirac equations. Accordingly, each of our stationary localised solutions represents a one-parameter family of travelling solitons which can be retrieved by the Lorentz boost (2.5)-(2.6). The stability properties of the moving solitons do not depend on their velocity. In contrast, obtaining travelling solitons in the Dirac equations with non-invariant nonlinearities would be a nontrivial affair. Stability would also have to be examined for each value of the velocity individually.

The second remark is on the interpretation of the $\mathcal{P}$- and $\mathcal{T}$-breaking terms. In the Schrödinger equations governing the amplitudes of optical beams, these terms are associated with gain and loss. The particular type of the $\mathcal{PT}$-symmetric Dirac perturbations that we scrutinised in this study, have a similar nature. However the Dirac modes that gain and lose energy, are not the components of the spinor but their linear superpositions. (See equation (2.19)). In the underlying physical system (e.g. two coupled oscillator chains), the symmetry-breaking terms may arise due to the coupling asymmetry rather than plain gain and loss. See section 4.2 above.

Acknowledgments

The authors gratefully acknowledge useful discussions with Georgy Alfimov, Abdul Kara and Boris Malomed. We thank Dmitry Pelinovsky and Taras Lakoba for reading the paper and giving their comments. This work was supported in part by the US Department of Energy. NA and IB would like to thank Center for Nonlinear Studies, Los Alamos National Laboratory, for warm hospitality during their stay. NA and IB were also supported by the National Research Foundation of South Africa (grants 105835, 85751 and 466082) and the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie Grant Agreement No. 691011. Computations were performed at the UCT HPC Cluster.