The phase diagram and vortex properties of PT-symmetric non-Hermitian two-component superfluid

We discuss the phase diagram and properties of global vortices in the non-Hermitian parity-time-symmetric relativistic model possessing two interacting scalar complex fields. The phase diagram contains stable PT-symmetric regions and unstable PT-broken regions, which intertwine nontrivially with the U(1)-symmetric and U(1)-broken phases, thus forming rich patterns in the space of parameters of the model. The notion of the PT-symmetry breaking is generalized to the interacting theory. At finite quartic couplings, the non-Hermitian model possesses classical vortex solutions in the PT-symmetric regions characterized by broken U(1) symmetry. In the long-range limit of two-component Bose-Einstein condensates, the vortices from different condensates experience mutual dissipative dynamics unless their cores overlap precisely. For comparison, we also consider a close Hermitian analog of the system and demonstrate that the non-Hermitian two-component model possesses much richer dynamics than its Hermitian counterpart.


I. INTRODUCTION
Quantum-mechanical systems are traditionally described by Hermitian Hamiltonians which ensure the realvaluedness of the full energy spectrum and, therefore, the unitary evolution of the system as a whole. It turns out, however, that the Hermitian description can be extended with a large class of non-Hermitian terms which are invariant under combined parity-time (PT ) transformations. The PT -symmetric non-Hermitian systems are as meaningful as the conventional Hermitian quantum mechanics in the regions where the PT -symmetry is not broken spontaneously [1,2].
Mathematically, in the PT -symmetric systems the familiar Hermiticity condition H † = H is replaced by the requirement of the PT -symmetricity (PT )H(PT ) = H, which is equivalent to the commutation of the Hamiltonian with the combined parity P and time-reversal inversion T operation [2], [PT , H] = 0. This combined symmetry leads to real-valued energy spectrum that ensures the stability of the system. One can show that all PT -symmetric non-Hermitian Hamiltonians belong to the class of the so-called pseudo-Hermitian Hamiltonians ηHη −1 = H † where η is a Hermitian linear automorphism [3]. The pseudo-Hermiticity is a generalization of the PT -symmetry which, in turn, depends crucially on the fact that the time-reversal transformation T is an anti-linear operation.
As it was shown in Ref. [4] and most recently emphasized in Ref. [5], it is the anti-linearity property rather than Hermiticity which is important for the selfconsistent description of stable quantum-mechanical systems. The non-Hermitian PT -invariant quantum systems can be mapped to their Hermitian counterparts via a non-Unitary transformation [6,7] (the existence of the map is not guaranteed as there are known exceptions in quantum mechanics [8]). These extensions broaden the class of stable physical systems beyond the tight Her-miticity constraints and open new horizons for the research.
The non-Hermitian description has been extended to interacting relativistic field theories, including the systems of fundamental particle interactions. The PTsymmetric interactions which explicitly break the non-Hermiticity of the system can arise in fermionic theories [9], contribute to the mass gap generation in the NJL model and affect the phase structure of the model [10]. The non-Hermitian Dirac fermions allow for the realization of an anomalous equilibrium transport [11]. The ordinary Hermitian models can generate a new, non-Hermitian ground state which could potentially be formed, for example, in fireballs of quark-gluon plasma created after heavy-ion collisions [12]. In the context of the extensions of the Standard Model of particle interactions, anti-Hermitian Yukawa interactions may lead to an anomalous radiative mass-gap generation in a model of the right-handed sterile neutrinos [13,14]. The concept of non-Hermitian quantum theory allows an extension via the gauge-gravity duality well beyond the scope of the field-theoretical models [15].
The PT -symmetric non-Hermitian Hamiltonians arise in the description of various open quantum systems in optics and solid state physics where this symmetry can be interpreted as a result of a perfect balance between the gains and losses as the system interacts with the external environment [16,17]. The recent works also include the studies of the effect of non-Hermitian terms in topological superconductivity which leads to nonlocal anomalous transport effects [18] as well as in the conventional superconductivity which gives rise to the unusual first-order phase transition between the phases [19]. The possibility of non-Hermitian superfluidity with a complex-valued, non-Hermitian interaction constant naturally arises from inelastic scattering between fermions [20]. The associated non-Hermitian BCS-BEC crossover of Dirac fermions in field-theoretical models of many-body systems reveals a nontrivial phase diagram as a function of the complex coupling [21].
In our paper, we work with vortex topological defects in a bosonic non-Hermitian model which involves a pair of scalar fields associated with interacting condensates. The topological solutions in the multicomponent scalar models are interesting because they appear in the models which have applications from condensed matter to high energy physics. Some of these models can serve as viable extensions of the Standard model of fundamental particle physics [23][24][25][26]. Similarly to the Grand Unification particle models and their close counterparts, they host 't Hooft-Polyakov monopole configurations [27] along with complex skyrmions [28] and kink/anti-kink solutions [29] with real-valued energies. As in the Hermitian models, these classical solutions are associated with the saddle points of the corresponding partition functions.
At the condensed matter side, the many-condensate systems possess richer dynamics than their onecondensate counterparts. For example, the standard classification of superconductivity into types I and II fails to describe the phases of multiband condensates so that a proposal to adopt a new terminology, a type-1.5 superconductivity, appeared in the theoretical community [30]. Experimentally, the existence of the type-1.5 superconductivity has been demonstrated shortly afterwards [31]. The semi-Meissner state of a type-1.5 superconductor demonstrates non-pairwise interaction between the vortices which leads to formation of a multitude of complicated vortex states [32]. We discuss the non-Hermitian extension of the two-component model possessing a global, rather than local, continuous symmetry, appropriate for the two-component superfluidity. We concentrate on stability of the ground state, fate of the PT symmetry in the interacting model, and the properties of the vortex configurations. In a different context, the vortices in the weakly interacting superfluid Bose-Einstein condensates with complexvalued PT -symmetric potentials have been investigated in Ref. [33]. As for the realistic model, one could expect the model can apply to the non-Hermitian version of twocomponent Bose-Einstein condensates in different hyperfine spin states similarly to the Hermitian analogues [50].
Non-Hermitian two-field models appear, for example, in the description of two-component out-of-equilibrium condensates in a non-Hermitian system of electron-hole pairs and photons in a semiconductor microcavity system. This open quantum many-body system resides in steady-state regimes characterized by a nontrivial phase diagram which contains an exceptional point that marks an endpoint of the first-order phase boundary [46] and exhibits anomalous critical phenomena [47].
The plan of our paper is as follows. In Section II we briefly overview the Lagrangian and its symmetries, and discuss the ground state of the minimal non-Hermitian theory with two scalar fields. The special attention is paid to the extension of the analysis of the PT sym-metries to the case of interacting model. In particular, we discuss the fate of the PT symmetry breaking in the presence of interactions. In Section III we consider the vortex properties for condensates with frozen radial degrees of freedom for which, in the absence of vortices, the condensate density is approximately uniform and the dynamics is encoded in the phase of the field. In Section IV we describe the examples of the vortex solutions at finite quartic couplings. The last section is devoted to our conclusions.

A. Lagrangians
We consider a simplest example of a scalar non-Hermitian theory which describes a PT -symmetric dynamics of two complex scalar fields φ 1 and φ 2 conveniently grouped into the single doublet field, The Lagrangian of the theory [7], includes the classical Hermitian self-interaction potential for the scalar fields: The non-Hermiticity is encoded in the real-valued mass matrixM 2 of the Lagrangian (2): provided the off-diagonal element 1 of this matrix is a nonzero, m 2 5 = 0. To see how the non-Hermiticity enters the theory, it is instructive to write the Lagrangian in terms of the individual fields φ 1 and φ 2 : The first term of the second line in Eq. (5) takes a purely complex value: −2im 2 5 Im (φ * 1 φ 2 ) if the off-diagonal component of the mass matrix (4) is a real-valued nonzero quantity. The complex valuedness of the Lagrangian (5) is consistent with the non-Hermiticity of the mass matrix in Eq. (2):M 2, † NH =M 2 NH . The off-diagonal mass with m 2 5 = 0 sets up the non-Hermitian regime in (4), while the point m 2 5 = 0 corresponds to a Hermitian theory (withM 2, † NH =M 2 NH ) which describes two non-interacting scalar fields φ 1 and φ 2 .
The model (2) describes two relativistic superfluids which interact with each other via the off-diagonal non-Hermitian coupling. We consider the potential (3) in the form which explicitly breaks the U (2) symmetry, Φ → ΩΦ with the 2 × 2 matrix Ω ∈ U (2), down to its Cartan [U (1)] 2 subgroup since the U (2) group is explicitly broken by the mass matrix (4) anyway provided In order to highlight the features of non-Hermiticity, we briefly discuss the Hermitian version of the model with the following mass matrix in the Lagrangian (2): Notice thatM 2, † H =M 2 H as expected and the Hermitian Lagrangian is a real-valued expression for all values of the parameter m 5 : As we discuss in Section III B 3, the Hermitian counterpart (7) of the non-Hermitian system (5) has some parallels with effective models of the QCD strings in highdensity quark matter, relevant to the physics of neutron stars, and the global cosmological strings, which reveal itself in the cosmological context [38][39][40].
Below, we will consider the classical equations of motion concentrating on the non-Hermitian theory described by Lagrangian (2). While we work with the classical solutions, we would like to notice that on quantum level, the non-Hermiticity propagates into the loops of perturbation theory. The quantum corrections could, therefore, induce a complex term in the interaction potential (3) and literally complexify the phase diagram of the theory. Leaving aside the quantum corrections, which were discussed in Ref. [25], in our paper we concentrate on classical properties of the Lagrangian (2) and compare it with the classical properties of its Hermitian counterpart (7).

B. Symmetries
We consider the non-degenerate model (2) with a nonzero off-diagonal mass m 5 = 0 which possesses a couple of continuous and discrete symmetries. Both Hermitian and non-Hermitian versions of the model are invariant under the global U(1) transformation in which the single phase factor (with a real-valued parameter ω) is shared by both complex scalar fields φ 1 and φ 2 . If m 5 = 0, both components of the doublet field Φ can be transformed independently so that the symmetry is enlarged to the global [U(1)] 2 . We do not consider this trivial case. The non-Hermitian theory (4) is also invariant under a discrete transformation corresponding to the product of the parity inversion P and the time conjugation T operators, respectively: Due to the presence of the Pauli matrix σ 3 , the upper φ 1 component transforms under the parity inversion P as a genuine scalar while the lower component φ 2 transforms as a pseudoscalar field. The P and T operations whose product leaves invariant the Hermitian theory (6), H : indicate that both fields φ 1 and φ 2 behave as true scalars under the parity inversion. The model (2), along with its extensions, possesses interesting features of the Goldstone modes associated with the spontaneous breaking of the continuous symmetry (8) as well as the unusual properties of the conserved currents [7,22,23]. Below we consider the ground state and the vortex solutions of the model.

Non-Hermitian ground state
The analysis of the ground state of the two-field model (2) has already been done analytically in Ref. [22] for the special case of the potential (3) in which one of the fields was not self-interacting (λ 1 = 0 and λ 2 = 0). In our paper, we complement this work with the numerical analysis of the system in which both fields experience the self-interaction, λ 1,2 = 0. We show that the simple extension (or, better to say, completion) of the model makes the analysis of the phase diagram very complicated.
Let us start from the simplest case when the quartic interaction is absent: λ 1 = λ 2 = 0. The Hermitian (6) and non-Hermitian (4) mass-squared matrices have the following eigenvalues, respectively: The vaccua in these models are stable provided the eigenmasses have no imaginary parts. For the Hermitian model with the mass matrix (6), this requirement implies M 2 H,− 0 or, naturally, H : These equations determine the region in the parameter space were the instability due to a negative mode (or, modes) at the trivial minimum, φ 1 = φ 2 = 0, does not occur. In an interacting theory with λ 1,2 = 0, this instability corresponds to the spontaneous symmetry breaking.
In the non-Hermitian model (4), the spectrum of the free (i.e., with λ 1 = λ 2 = 0) theory does not contain complex energy eigenvalues provided NH : where we also took into account a possibility that the squared masses m 2 1 and m 2 2 can take negative values. The first two requirements in Eq. (13) correspond to the instability related to the spontaneous symmetry breaking rather than to the non-Hermiticity of the model. The last condition in Eq. (13) highlights the region of the parameter space where the PT symmetry is said to be unbroken [7,22,23]. Together, these conditions guarantee the stability of the ground state.
The classical equations of motion of the model (2) are obtained by the variation of the Lagrangian (5) with respect to the independent fields φ * 1 and φ * 2 , respectively: A variation with respect to the fields φ 1 and φ 2 gives us an inequivalent set of equations: which differs from Eq. (14) by the sign flip of the offdiagonal mass term, m 2 5 → −m 2 5 . The striking inconsistency of the two pairs of equations of motion, (14) and (15), poses the natural question of how to treat this non-Hermitian system at the classical level. One of the earlier proposals [7] suggests to use only one set of equations of motion, either (14) or (15), and omit the complementary set. Indeed, the choice of the set flips the sign of the off-diagonal mass, m 2 5 → −m 2 5 , which does not affect the physical spectrum neither in free model (11) nor in the interacting model, as we will see below. Moreover, the non-Hermitian system must be open via coupling to an external source which equilibrates non-vanishing surface terms in the variation of the action. When the off-diagonal mass vanishes, the system becomes Hermitian, and the single set of equations is enough to describe the solutions consistently.
The next proposal was to resolve the apparent problem with the inconsistency of the classical equations of motion was put forward in Ref. [22] where a similarity transformation for the action has been used to achieve harmony between the two sets of equations of motion. This elegant idea initially required an extension of the field space using two complex components (with four degrees of freedom) for every original complex field (which hosts two degrees of freedom).
The same mapping strategy has later been adapted and extended to many-field theories in Ref. [23] where two real-valued fields were used to represent one complex field. However, this procedure, used for the mapping of the non-Hermitian theory to the Hermitian one via the similarity transformation, leads to the appearance of a negative kinetic energy term for one of the fields. While it was argued that this artifact does not change the signature of the appropriate Hilbert space [22], the appearance of the negative kinetic action leads to a negative contribution to the energy-momentum tensor so that the energy of a classical configuration (which is not necessarily a classical solution) would become unbounded from below. This pseudo-Hermitian method can be adapted, in specifying the appropriate physical regions, to give the classical solutions for non-Abelian monopoles in a spontaneously broken theory [27] (we refer the Reader to Refs. [36,37] for complementary discussions of a non-Hermitian non-Abelian gauge theory).
In our paper we follow the approach of Ref. [7] where only one set of equations of motion -either (14) or (15), but not the both of them -is considered. While this approach may seem restrictive, it still gives the complete description of classical and quantized theories [24]. In addition to the invariance of the physical solutions with respect to the flip m 2 5 → −m 2 5 , the choice of the equations of motion coincides with the choice of the Hamiltonian operator in the quantum theory since bothĤ(m 2 5 ) and H † (m 2 5 ) =Ĥ(−m 2 5 ) can be used to promote the time evolution of the system. Moreover, the choice does not affect the non-Hermitian physics: the ±m 2 5 versions of the classical equations of motion give the physically equivalent classical solutions of the theory and, in parallel, both original and Hermitian-conjugated Hamiltonians determine the very same evolution of the quantum theory [24].
The classical non-Hermitian theory becomes selfconsistent if one supplements the Hermitian conjugation with the subsequent PT symmetric operation. Then the complex conjugation along with the PT flip m 2 5 → −m 2 5 leaves the classical equations intact [24]. The new combined operation is also important at the quantum level since the Hermitian-conjugated conversion between the bra and ket states should now be supplemented by the additional PT invariance operation in order to maintain the orthonormality of the eigenvectors of the non-Hermitian Hamiltonian [7]. Thus, the would-be apparent non-equivalence of the original (14) and conjugated (15) equations does not lead to an inconsistency of the non-Hermitian theory. It gives the opposite: both formulations are consistent with lath other and lead to the same result both at classical and quantum levels. At a certain stage, we use a numerical method to find the classical solutions using the criteria of the energy minimization as a selection principle of the right, "minimal energy" solution among all other available solutions. In our approach, the energy of the classical configuration in the non-Hermitian theory is bounded from below so that the numerical approach is self-consistent in finding the correct configuration. If we would otherwise employ the pseudo-Hermitian procedure of mapping the non-Hermitian theory to the Hermitian one, then the classical energy becomes unbounded due to the negative sign in the kinetic terms, and the numerical procedure fails to converge to a reasonable solution.
In the ground state the condensates are coordinateindependent quantities, and Eqs. (14) reduce to the nonlinear algebraic relations: The use of the complementary set of equations (15) instead of Eqs. (14) would lead to an equivalent physical solutions. Indeed, the swap of equations leads to the sign flip in Eq. (16) which corresponds to a simple swap of the fields φ 1 and φ 2 .
It is convenient to represent the fields φ a in the radial form φ a = v a e iθa for a = 1, 2. Equations (16) can possess nontrivial solutions provided the phases θ a satisfy one of the following relations: One can cover both these cases assuming θ 1 = θ 2 and, simultaneously, allowing the amplitudes v a to take both positive and negative values, v a ∈ R. The classical equations (16) then determine the amplitudes: The canonical energy-momentum of the model (2) can be calculated by endowing the model with the gravitational background g µν , considering the variation of the theory action, S = d 4 L with respect to the metric, and then setting the metric to its Minkowski values back again, g µν → η µν ≡ diag(+1, −1, −1, −1). The energy density of the model (5) is given by the µ = ν = 0 component of the energy-momentum tensor (19): One gets the following simple expression for the energy density of a uniform time-independent (ground) state: Both solutions (17) for the phases θ a ≡ arg φ a of the condensates φ a lead to the vanishing non-Hermitian term in the energy density (21). Thus, the inconvenient imaginary terms do not enter the ground state of the model (we will see below that any solution of the single set of the classical equations of motion (14) has a real-valued energy density). However, despite the non-Hermitian mass m 5 not explicitly appearing in the vacuum energy (21), it affects the ground state indirectly via the solutions of the equations of motion (18).

Hermitian ground state
The equations of motion of the Hermitian model (7) take the following form: which differ from Eq. (14) of the non-Hermitian model (5) as the both signs in front of the off-diagonal mass terms m 2 5 are the same in the Hermitian case (22). This small difference naturally propagates into the equations for the amplitudes in the ground state, as compared to its non-Hermitian analogue (18). However, in the striking dissimilarity with non-Hermitian case, the variations over the fields and their complex conjugates are obviously consistent with each other in the Hermitian model. Moreover, the energy density of the ground state of the Hermitian model, explicitly includes the m 2 5 coupling between the amplitudes of the different condensates. We remind that this coupling does not enter the non-Hermitian energy density (21).

Hermitian model vs. non-Hermitian model
In the Hermitian model, a negative diagonal mass can trigger a natural spontaneous symmetry breaking of the global U(1) symmetry (8) which leads to a non-zero expectation value of the doublet field, Φ 0 = 0. The runaway of the scalar field from the symmetric state is balanced by energy density condensate 1 condensate 2 quartic (self)-interactions so that the symmetry-broken ground state of the system is stable. In an interacting model, the stability of the symmetric, Φ = 0, state is determined by Eq. (12).
The interacting non-Hermitian theory possesses two types of instabilities. In addition to the mentioned spontaneous symmetry breaking stipulated by the first two relations in Eq. (13), the symmetric ground state can also experience the U(1) symmetry breaking which could be caused, in turn, by the broken discrete PT symmetry. This purely non-Hermitian effect is dictated by the third relation in Eq. (13), the violation of which leads to the complex mass spectrum and, consequently, to the instability. The PT -induced instability could be well seen when both diagonal masses are positive so that the first two relations in Eq. (13) are satisfied. As in the purely Hermitian case, the symmetry-broken ground state of the non-Hermitian system is expected to be stabilized by the quartic (self)-interactions.
In Fig. 1 we compare the U(1) symmetry breaking patterns in the non-Hermitian model with m 2 1 > 0 (the upper panel) and m 2 1 < 0 (the middle panel) as well as the Hermitian model with m 2 1 < 0 (the lower panel). We show the total energy density E determined by Eqs. (21) and (24), and the condensates v 1 and v 2 . Notice that the non-Hermitian (18) and Hermitian (23) equations on the ground state imply that the relative sign of the condensates is determined by the ground state while the overall sign is not fixed. To remove this arbitrariness, we require the condensate v 2 to be positive.
If both diagonal squared masses are positive, m 2 1 > 0 and m 2 2 > 0, and the off-diagonal mass is zero, m 5 = 0, then the ground state resides in the symmetric phase with vanishing condensates. As the off-diagonal mass increases in its absolute value, the theory should experience the PT -symmetric instability and we expect the appearance of the condensates in the (m 2 1 > 0, m 2 2 > 0) plane. However, this does not happen, as one can see in Figs. (1)(b) and (c). The presence of the "PT -induced" condensates is forbidden by the requirement of the minimization of the energy density (21) as these condensates would have made the energy density positive while the trivial ground state with v 1 = v 2 = 0 has zero energy. This conclusion appears to be the consequence of the fact that the expression for the non-Hermitian energy density (21) does not contain a "compensating" m 2 5 term which is present, on the contrary, in the expression of the Hermitian energy density (24). Therefore, the U(1) symmetric ground state with vanishing energy is preferred over the U(1) broken state with the nonvanishing condensates. We will see below that some of these symmetric states belong to the PT -broken phases of the interacting model and, therefore, are physically meaningless.
In the symmetry-broken phase with m 2 2 < 0 (while we still keep m 2 1 > 0), the conventional U(1) symmetry breaking does occur. The upper panel of Fig. 1 shows that the non-Hermitian mass influences the conventional symmetry breaking in a somewhat controversial way: at sufficiently large m 4 5 , the third requirement of Eq. (13) is violated, the PT symmetry gets broken and the U(1) symmetry gets restored because the ground-state condensates v 1 and v 2 vanish. In the symmetry-broken region, the flip of the sign in the m 2 5 mass flips the sign of the v 1 condensate (we remind that we always keep the v 2 condensate non-negative).
A rather similar picture occurs when one of the diagonal mass squared is taken to be negative, m 2 1 < 0, as shown in the middle panel of Fig. 1. The spontaneous symmetry breaking occurs at both signs of the remaining mass squared, m 2 2 , while the increase of the absolute value of the off-diagonal mass m 2 5 leads to the restoration of the U(1) symmetry. We would like to stress again that this particular picture is enforced by the requirement of minimization of the non-Hermitian energy density (21) which allows us to choose the correct ground state from the multitude of the solutions of the non-Hermitian equations of motion (18).
Finally, the ground state of the Hermitian two-scalar model with m 2 1 < 0 is shown in the lower panel of Fig. 1. The condensates appear almost at every point of the phase diagram except for the line m 5 = 0, where the condensates decouple and the v 2 condensate ceases to exist at m 2 2 > 0. At this semi-infinite line, the condensate v 1 takes its minimum. 4. Stability of the ground state and the PT symmetry in the non-Hermitian model Before proceeding to the discussion of the vortex solutions, let us address the formal stability issues of the ground state. Usually, the local stability of a classical configuration is probed by expanding the scalar fields in the vicinity of the configuration, φ a = v a +φ a with |φ a | |v a |. The configuration is unstable if the fluctuation matrix corresponding to the variation of the action with respect to the fluctuation of the fields contains negative modes.
In the Hermitian theory, the global minimum in the total energy of the solution corresponds to an absolutely stable state. In the non-Hermitian theory, this criterion may not work since even in the classical theory, the ground state is determined by a single set of classical equations of motion out of the existing two sets. Therefore, mathematically, one could expect the emergence of negative directions in the space of fields. This expectation is coherent with our physical intuition since the non-Hermitian system resides in the steady state which is generally not a thermodynamic equilibrium. In the remaining part of this section, we argue that the formal stability criteria can still be applied to the interacting non-Hermitian systems. These criteria give us the generalization of the PT -symmetric and PT -broken regions for the interacting theory.
From the classical equations of motion (16), we obtain that the field fluctuations around the condensates are governed by the following equation: where the fluctuation matrix M 2 NH of the non-Hermitian model is as follows: In the U(1) broken phase, this matrix has one zero eigenvalue which corresponds to the Goldstone mode. In the symmetric U(1) phase, all eigenvalues are generally nonzero.
For the ground state to be stable, one expects that the all eigenvalues of the fluctuation matrix are non-negative.
In the non-Hermitian model, this requirement is not always satisfied. In Fig. 2 we show the (in)stability phases for various quartic couplings λ 1 and λ 2 .
The (in)stability phase diagram has a rather curious form. One notices that the borderlines between the stable and unstable areas involve both U(1) symmetric and broken regions which possess, on the contrary, rather featureless, smooth behavior of the condensates as shown in Fig. 1. In order to understand the appearance of the negative modes, let us consider any unstable point in the symmetry-restored region with masses m 2 1 > 0, m 2 2 > 0 (these regions are marked by the white color in the upper panel of Fig. 2). Since the both condensates are zero and both diagonal masses as well as the couplings are positive, the energy density (21) takes an absolute minimum. On the other hand, in these regions, the third criterion of Eq. (13) is not satisfied, indicating that the symmetric state resides in the PT -broken region and is thus unstable. The instability is not captured by the minimization of the energy density (21) since this particular expression is derived for the classical solutions while the instability can drive the configuration out of the classical subspace of configurations. The non-classical configuration can acquire even a complex value thus invalidating the criterion of the energy minimization outside of the classical subspace of field configurations.
Thus, we arrive to the conclusion that the ground state of the non-Hermitian model is formally not stable in certain regions of the parameter space. One can easily check that these unstable regions become the PT -broken regions when the quartic interaction couplings are set to zero, λ 1,2 → 0. Therefore, similarly to the free noninteracting model, the unstable regions indicated in Fig. 2 can be interpreted as the PT -broken regions in the interacting case where the model cannot be used as a invalid prescription of any steady state in a physical system. The PT -symmetric regions, on the contrary, are valid stable zones where the steady-state physics can be realized.
The stability, now associated with the PT -symmetric regions, is thus determined by the standard requirement that the quadratic fluctuation matrix, does not possess negative eigenvalues. Here the vector χ = (φ 1 , φ 2 , φ * 1 , φ * 2 ) denotes the original fields and their conjugates.
For the sake of completeness, we show the fluctuation matrix of the Hermitian model (7), which has only non-negative eigenmodes, as expected. The Hermitian ground state does not have any ambiguities in the stability criteria.

III. VORTICES IN FIXED-LENGTH CONDENSATES
A. Vortices in one-component superfluids: brief review In this section, we consider in the bosonic condensates with the fixed radial (amplitude) part |ψ| of the scalar condensate ψ = |ψ|e iθ when the entire dynamics of the system is encoded in the condensate phase θ. We can neglect the fluctuations of the amplitude in the long-range limit for fluctuations with a wavelength longer than the appropriate healing length ξ [51].
The healing length corresponds to the wavelength at which the kinetic energy of the boson equals the chemical potential. It determines the size of vortices and plays a role similar to the one played by the correlation length in the superconductivity. On the superconductivity side, a condensate with a fixed amplitude appears naturally in the London limit which corresponds to the extreme type II superconducting regime. In the London limit, the penetration depth of the magnetic field is much longer than the correlation length so that the latter can be neglected in many physical situations. Equivalently, this limit poses a constraint on the modulus of the superconducting condensate which can be considered as a rigidly fixed, non-fluctuating quantity. The gauge invariance implies that the phase of the superconducting condensate is not constrained by this limit, thus sharing an analogy with a superfluid. However, the phase is absorbed by the electromagnetic gauge field, which becomes -consistent with the Meissner effect -a massive vector field via the Anderson-Higgs mechanism. [34] In the background of sufficiently strong magnetic field, the type-II superconductors enter the mixed (Abrikosov) phase where the magnetic flux penetrates the condensate in a form of parallel vortices which form the regular Abrikosov lattice. In the London limit, the vortex is characterized by a singularity in the phase of the condensate which reflects itself in a singular behavior of the gauge field close to the vortex core. [35] In our paper, we consider the model with the global, rather than local/gauge, U(1) symmetry. In the global case, the condensate is associated with the electrically neutral complex scalar field which corresponds to the phenomenon of superfluidity rather than superconductivity. The Anderson-Higgs mechanism is evidently absent and the massless mode in the phase of the condensate appears to be the Goldstone boson associated with the global U(1) symmetry breaking. In the case of a nonrelativistic Bose gas, one identifies the Goldstone mode with a phonon. Instead of the Abrikosov strings with magnetic flux, the superfluid condensates host the vortices characterized by singularities in the phases of the condensate. These vortices are also called the "global vortices" since they appear in theories possessing a global symmetry. Below, we will briefly review the theory of vortices for a non-relativistic one-component superfluid following Ref. [35] and then we proceed to the generalization of our approach to the long-range ("London") limit of the relativistic non-Hermitian two-field model.
Consider a superfluid condensate of a non-relativistic Bose gas of particles with the mass M . It is convenient to describe by the wave function in the radial coordinates ψ = |ψ|e iθ . The long-range limit corresponds to the fixed radial degree of freedom |ψ| which is a good approximation at very low temperature. Then the energy of the fluctuations in the superfluid can be written in the following form (we work in units c = = 1): where the radial part |ψ| of the condensate determines the superfluid density ρ s = M |ψ| 2 while its phase θ provides us with the velocity of the superfluid condensate, v s = (∇θ)/M . We neglected here qualitatively inessential contributions coming from the inhomogeneity of the radial part of the condensate. An extremization of the energy (27) gives us the Poisson's wave equation, ∆θ = 0 which corresponds to the incomprehensibility condition of the superfluid, ∇ · v s = 0.
The superfluid theory (27) also incorporates singular configurations of the vortex fluid which correspond to the superfluid currents that wind around a line (called the vortex line) in three-dimensional space. Precisely at the vortex line, the superfluid condensate vanishes, φ = 0, while the phase of the superfluid condensate is equal -in the vicinity to the vortex core -to a geometrical angle (times an integer) in the plane transverse to the vortex line. The integer is a topological quantity which characterizes the vortex winding number.
The vortices can be accounted in the model (27) despite the fact that this model describes a long-range macroscopic physics with a globally uniform, constant condensate. Consider, for simplicity, a straight static vortex line along the direction x 3 located at the center x 1 = x 2 = 0 in the transverse plane. The phase of the vortex is θ(x 1 , x 2 ) = nϕ where n is the winding number of the vortex and ϕ = arg(x 1 + ix 2 ) is the azimuthal angle in the two-dimensional (x 1 , x 2 ) plane. The velocity of the superfluid, is directed along the unit vector in the azimuthal e ϕ , so that the fluid "winds" around the vortex center at x 1 = x 2 = 0. Alternatively, the location of the vortex singularity corresponds to the point in the two-dimensional plane, where the derivatives do not commute, The total energy (27) evaluated at this static straight vortex configuration is: where L is the length of the vortex. The infrared cutoff of this integral is the size of the system R while the size ξ of the vortex core -typically, of the interatomic distanceserves as an ultraviolet cutoff. The important message of Eq. (30) is that the vortex energy (mass) is proportional to the length of the vortex L. One can show that the vortex curvature provides a subleading correction. The energy density per unit length is a finite quantity since the logarithmic divergence in Eq. (30) is very mild.
Consider now the relativistic one-component model with the action: where we adopted the long-range limit with the constant radial condensate κ = |φ|. The variation of the action with respect to the phase θ gives us the equation for the propagation of the massless Goldstone particle, θ = 0, where = ∂ 2 t − ∇ 2 is the d'Alembert operator. In the non-relativistic limit, the model (31) possesses the energy (27) for static vortex configurations. It is convenient to parameterize the coordinatesx µ =x µ ( σ) of the two-dimensional vortex singularities by the two component vector σ = (σ 1 , σ 2 ). We split the phase θ into the regular and singular parts, θ = θ r + θ s . The regular part θ r of the phase corresponds to the perturbative fluctuations while the singular part θ s encodes the position of the vortex (29). In the relativistic notations, where the singularity itself is given by the tensor current: which is expressed via the differential measure at the vortex world-sheet Σ: Here ab is the fully anti-symmetrized tensor in two dimensions, 12 = − 21 = +1 and 11 = 22 = 0. The vortex tensor (33) is a two-dimensional delta-function at the surface of the vortex, with the orientation at the vortex world-sheet. For example, for a straight vortex mentioned above, one uses the parameterizationx 0 = σ 1 , x 1 =x 2 = 0,x 3 = σ 2 , and obtains where n ∈ Z is the vorticity. Note that the regular part of the phase does not contain any singularity by definition, [∂ µ , ∂ ν ]θ r ≡ 0. Integrating out or, equivalently, solving the equations of motion for the regular component the phase, θ r , allows us to rewrite the action (31) in terms of the singular part of the phase θ s , which, in turn, depends only on the vortex world-sheet (33): This action is a nonlocal functional which features two integrals that are taken over the same vortex worldsheet. In the case of many vortices, the worldsheet Σ includes all their worldsheets: Σ µν = Σ 1,µν + Σ 2,µν + . . . .
The nonlocal action (36) represents the self-interaction of the vortex line as well as the interactions of the distinct vortex segments via propagation of a massless Goldstone particle between the vortex segments. In Eq. (36), this long-range interaction is represented by the advanced Green's function D(x) of the d'Alembert operator: In the case of a static straight vortex line of large length L, the action (36) calculated for the time interval δt gives us S = Eδt, where E is, up to parameter redefinitions, the known vortex energy (30). In order to demonstrate this fact, it is convenient to make a Wick transformation in the integral (36) to the Euclidean spacetime. In the Euclidean space, the massless propagator (37) is: where x is the 4-distance. For a single straight static vortex with the surface (35), the longitudinal (along the vortex) and temporal coordinates in the action (36) can be integrated out, and we get the following formal expression for the vortex energy: where is the two-dimensional massless propagator [a solution of Eq. (37) in two Euclidean dimensions] as the function of the two-dimensional distance ρ. The parameter ρ 0 , which has the dimension "length", is introduced for the consistency reasons. The argument "0" in Eq. (39) highlights the fact that the formal expression for the energy in the long-range limit is a logarithmically divergent quantity similarly to the non-relativistic expression (30). A more accurate derivation in a finite cylindrical box of the radius R 0 leads us to where ξ is the size of the vortex core given by the healing length which does not explicitly enter the model in the long-wavelength limit. The general expression (36) also gives us the interaction energy of the two straight static vortices with the vorticities n 1 and n 2 separated by the distance R: The like-charged vortices repel each other while the vortices with opposite vorticities attract to each other. Finishing this section, we notice that in order to get Eq. (42) it is sufficient to take in the action (36), instead of the single-vortex current (35), the following expression: B. Two-component superfluids in long-range limit

Lagrangians in the long-range limit
The long-range limit of both Hermitian (7) and non-Hermitian (5) two-condensate models can be reached by expressing the diagonal masses m 2 a = −2λ a κ 2 a via the condensate densities κ 2 a > 0 and a = 1, 2 and then taking the limit of strong quartic interaction, λ 1 = λ 2 → ∞. The parameters κ a > 0 fix the radial amplitudes for each field, φ a = κ a e iθa , while leaving the phases θ a as the only dynamical variables. In the long-range limit, the Lagrangians for the Hermitian (7) and non-Hermitian (5) models reduce, respectively, to the following expressions.
The only difference between the Hermitian and non-Hermitian cases appears in the interaction between the phases of different condensates: instead of the cosine function in the Hermitian model, its non-Hermitian model has a sine function preceded by a purely imaginary coupling. According to the third criterion of Eq. (13), the non-Hermitian theory in the long-range limit at κ 1 = κ 2 would correspond to the PT unbroken phase if the theory were non-interacting. Below, we will see that in the interacting theory this criterion, unsurprisingly, does not work. Still, Eq. (45) represents a meaningful theory even if m 5 = 0.

Hermitian two-condensate model
The system described by the Hermitian Lagrangian (44) corresponds to a long-range description of two coupled condensates within the Gross-Pitaevskii formalism [51]. The condensates correspond to two different hyperfine spin states of (for example, of 87 Rb), which are driven by the Rabi frequency Ω.
In the static limit, the energy functional of the system (44) can be written in a suggestive non-relativistic form: where ρ a ≡ 2mκ 2 a is the density of the a = 1, 2 condensate and the interaction term Ω ≡ −2m 2 5 κ 1 κ 2 has the meaning of the Rabi frequency (here we restored for a moment the powers of ). Below we continue to work in the notations of Eq. (44).
From the classical equations of motion of the Hermitian model, one immediately determines the presence of the Goldstone massless mode and the massive excitation, These degrees of freedom satisfy, respectively, the following equations: Here is the mass of the mode (50), is the angle which determines the relative strength of the condensates, and is the total condensate density given by the sum of the densities of the individual components. The mode χ corresponds to the gapless sound excitation (the density wave) while mode γ is the spin density wave which is always gapped provided the Rabi frequency is nonzero, Ω = 0 (or, in our notations, m 5 = 0) [42].
Since κ a > 0, the mass (53) can always be chosen as a positive quantity. The Hermitian Lagrangian (44) can therefore be rewritten as a sum of independent contributions coming from massless (49) and massive (50) fields: The Goldstone mode χ corresponds to the massless excitation of the one-component Bose gas superfluid (31). The massive mode γ is described by the sine-Gordon Lagrangian: where the parameterκ plays a role of the amplitude corresponding to the massive condensate. In Eq. (57), the term cos γ appears naturally instead of the usual γ 2 mass term in agreement with the γ → γ + 2π periodicity coming from the symmetry of the original fields θ 1 and θ 2 . It is instructive to look at the limit of large mass, M ξ 1, which reduces the fluctuations of the field γ and makes the last term in the model (57) quadratic. Then the integration over the field γ can be done explicitly taking into account the vortex singularities.
Following the analogy of the one-component model, we conclude the two-component model should contain two types of superfluid vortices associated with singularities (winding) in the phases of the fields θ 1 and θ 2 . Using the decompositions (49) and (50), we identity the massless Σ (0) and massive Σ (M ) combinations of the vortex worldsheets which appear as the singularities in χ and γ fields, respectively: Here the individual phase windings are defined according to Eq. (32): Therefore, the effective theory of vortices is written as follows: where D M (x) is the advanced Green's function corresponding to the propagator of the massive particle, Despite its cumbersome appearance, the physical sense of the effective action (62) is relatively transparent: the vortex worldsheets interact with each other via the massless density waves and the massive spin-wave exchanges in the combinations Σ (0) and Σ (M ) . The interaction is given by combining the attractive long-range (Coulomblike) potential and the short-range (Yukawa) interaction with the strength determined by the couplings κ andκ, and the angle β which determines the relative strength of the condensates.
Let us consider what the consequence of the action (62) for two static parallel straight vortices of the length L ξ is. The Coulomb interaction of the first term in Eq. (62) is given by the long-range logarithmic potential (42). The short-range part which appears in the second term of Eq. (62) corresponds to the interaction (we neglect the energy of the single vortex): where K 0 is the modified Bessel function of the second kind. In the limit of zero Rabi frequency, Ω = 0, the coupling between the vortices in different condensates disappears, m 5 = 0, the mass (53) vanishes, and the massive spin-wave exchange (63) takes the logarithmic form (42). However, since the action (62) is derived in the limit of the massive spin-wave, M ξ 1, the last term leads only to the short-range interactions. The logarithmic potential dominates the large-distance physics leading to the following potential between the vortices of the kinds a, b = 1, 2: where α 1 = sin 2 β and α 2 = cos 2 β gives the amplitudes of, respectively, the fields θ 1 and θ 2 in the massless density wave (49). The expression (64) is a trivial generalization of the vortex-vortex interaction in the case of a single condensate (42).
The two-component model possesses, however, a new object, the domain wall in the scalar condensates. [42,43] This nontrivial topological structure has no analogue in the one-component model. Technically, its appearance can be seen from the last (cosine) term in the two-field Lagrangian (44) which leads to a soliton in the massive field γ, Eq. (50). The domain wall appears due to the coupling m 5 , which, as we mentioned, corresponds to the Rabi frequency Ω that couples the two spin states. The field γ changes by 2π across the vortex, which guarantees the topological stability of the wall within the model (44).
The boundary of a domain wall is necessarily a vortex that carries an integer winding number. Since the domain wall possesses a nonzero tension, the vortex cannot exist as an isolated object because it would otherwise have infinite energy due to the wall attached. Therefore, the vortices from different condensates are bound by a linear potential because the domain wall possesses a constant energy per unit area. Such bounded vortex pairs are sometimes called merons [50].
The existence of the domain wall provides us with yet another exciting link with the physics of strong interactions described by Quantum Chromodynamics (QCD). If the domain wall in the mentioned vortex pair gets broken, new vortices appear from the condensate. The created vortices produce additional vortex pairs that are attached to the broken ends of the domain wall [42,44]. The process is a direct analog of the confining QCD string breaking, which extends between the (anti)quarks and gets broken if a sufficiently long distance separates the quarks. Thus, the appearance of the domain wall leads to a rather nontrivial vortex-pair dynamics [44,45].
The confinement picture works only for a relatively small inter-condensate coupling given by the Rabi frequency Ω (or, for small off-diagonal mass m 5 in our notations). At a high inter-condensate coupling, the domain wall becomes unstable, and the vortex pairs are no more confined [42]. The inter-vortex potential is then given by the long-range logarithmic potential (64) determined by the density wave exchange.

Non-Hermitian two-condensate model
The energy of the non-Hermitian two-condensate model (45) can also be rewritten as a formal expression similar to its Hermitian counterpart (46). In the non-Hermitian case, the meaning of the Rabi frequency Ω becomes rather transparent: the Rabi coupling Ω should couple to the different hyperfine spin components with a different sign. Leaving aside the physical realization of such scenario, we notice that the non-Hermitian case can be analyzed closely paralleling the Hermitian case considered above. In particular, exactly the same field combinations (49) and (50) can be used to rewrite the non-Hermitian theory (45) in terms of the massless (density wave) and massive (spin wave) phase combinations: This non-Hermitian model possesses the usual Hermitian Goldstone mode χ which leads to the long-range interactions between the combinations of the worldsheets (59). However, the would-be massive excitation γ exhibits a non-Hermitian behavior described by the Lagrangian: and therefore the interaction between the would-be massive components of the vortex sheets (60) is not evident. The model (66) is nothing but a non-Hermitian version of the Sine-Gordon model in (3+1) dimensions. As a side remark, we notice that the Lagrangian (66) appears as a bosonic dual of the non-Hermitian massive Thirring Model in (1+1) dimensions [6]. According to Ref. [6], this model with the purely imaginary coupling in front of the sine-term resides in the PT -broken domain and, therefore, should be characterized by complex energy dispersions which correspond to dissipation or instability, or the both.
The second term in the Lagrangian (66) implies that the γ = 0 point is not a local extremum of the corresponding action. The stable minima could appear around the values γ ± = ±π/2. Defining γ = ±π/2 + δγ, we get the following equations of motion: which differ from the classical equations of motion of the Hermitian model (52) by the purely complex coefficient in front of the sine term. For small fluctuations around the minimum, δγ = 0, the solutions of Eq. (67) give us the dispersions for the energy: ω k = √ k 2 ± iM 2 . We find that the particle-like, positive-energy solutions with Re ω > 0 lead to an explosive behavior near the γ = −π/2 minimum which appears to be unstable. For example, the amplitude of any zero-momentum solution (k = 0) diverges with time t as δγ ∼ exp(Γt) where In the language of a Hermitian theory, the point γ = −π/2 would correspond to an extremum. However, the minimum γ = +π/2 is stable so that all particle excitations around it behave as dissipative solutions δγ ∼ exp(−Γt) that approach the minimum point γ = +π/2 should the field γ deviate from it. Since the angle γ takes a constant value in the stable minimum, it evidently means that the field γ contains no vortex singularities in this ground state. According to Eq. (60), the vortex singularities in the both phases θ 1 and θ 2 should coincide with each other: Σ Thus, the vortices in the φ 1 and φ 2 condensate can only exist provided they coincide with each other, Σ The common vortex line is described by action of the vortex in the one-component condensate (36), where the coupling κ is given in Eq. (55). The energy per unit length of the joint vortex is thus given by Eq. (41). In this state, the common vortex segments interact with other segments via a long-range interaction mediated by massless particles. For the straight static vortices separated by the distance R, the interaction is given by Eq. (42).
As we mentioned earlier, the Hermitian counterpart (7) of the non-Hermitian system (5) appears in various physical contexts, including color superconductors in highdensity quark matter described by QCD, relevant for the description of, for example, the interior of the neutron stars [38,39], and the axion cosmic strings [40]. In the field-theoretical context, both original U(1) degrees of freedom, representing the phases of the scalar condensates, mix up to form the U(1) B ("baryon") global symmetry where both condensates transform with the same phase and U(1) A ("axial") global symmetry which rotates the phases of the fields oppositely. 2 In our mod-els, these global symmetries correspond to the Goldstone massless mode (49) and the massive excitation (50), respectively.
In particle physics, the axial symmetry is usually broken either by the axial quantum anomaly [48] or explicitly, by the (strange) quark mass, depending on the physical system described by the model. This effect is important in the cosmological context as it is responsible for the rapid decay of axion cosmic strings [40]. These topological defects, which are suggested to be formed in early moments of our Universe, gets bounded to the QCD domain walls at the QCD finite-temperature phase transition. The walls decay rapidly after the QCD phase transition thus leaving no trace of the cosmic strings in the present-day Universe [41]. We will not pursue this topic further concentrating, instead, on the properties of vortices themselves in the contexts of the models. We notice that in model (7), the "axial" invariance is explicitly violated by the nonzero coupling m 5 = 0 between the two fields in the scalar doublet which makes the overall mechanism similar to the explicit breaking of the axial symmetry in QCD.

Vortices in Hermitian and non-Hermitian models in long-range limit: a brief comparison
It is worth comparing the properties of the vortices in Hermitian and non-Hermitian versions of the twofield model in the long-range limit. We have already seen that outside the long-range limit, the phase diagram of the Hermitian model (Section II C 2) has a much simpler structure as compared to the rich and complicated diagram which characterizes the phases of the non-Hermitian model (Section II C 4). The properties of the corresponding vortices, discussed in Sections III B 2 and III B 3, respectively, differ dramatically as well.
The Hermitian version of the model possesses two types of stable vortices associated with the phases of both condensates (61). The vortices interact with each other via nonlocal action (62) which describes the massive and massless interaction of the vortex worldsheets via exchanges of the massless (density wave) and massive (spin wave) excitations. The model also possesses the domain wall, a coherent structure in both condensates. However, in the limit of the sizeable inter-condensate coupling (given by the Rabi coupling Ω or the off-diagonal mass m 5 , in our notations), the domain wall becomes unstable [42] and the long-range vortex dynamics is determined by the repulsive logarithmic potential.
The long-range limit of the non-Hermitian model resides in the PT -broken phase, which indicates, according to the general arguments discussed earlier, instability.
fermions where the vector and chiral rotations in internal space transform upper (right-handed) and lower (left-handed) Weyl components with the same (opposite) phases which corresponds to the vector (baryon) and axial rotations, respectively [49].
The instability indeed appears as a dissipation provided the vortex cores of different types are not overlapping exactly. Therefore, the separated vortices from different condensates tend to attract each other until they form a tightly bound pair that is perfectly stable.
Therefore, the behavior of the vortices in the Hermitian and non-Hermitian cases are somewhat similar: both of them tend to form pairs, albeit for different reasons. However, one should notice that the similarity is not exact: the binding of vortices in the Hermitian model appears due to the domain wall, which exists at low intercondensate coupling, while the dissipative vortex binding in the non-Hermitian case is realized at the large values of this coupling.

IV. VORTICES IN GENERAL CASE
The non-Hermitian two-field model possesses vortex solutions not only in the London (long-range) limit, where the amplitudes are fixed by the large quartic couplings λ 1 and λ 2 but also at finite values of the quartic couplings. The general is particularly interesting for the properties of the system at the scales of the order or smaller than the healing length.
In this section, we consider examples of the static straight vortex solutions of the classical equations of motion (14) assuming the standard axial ansatz for the scalar fields φ a (r, θ) = v a f a (r)e inθ , a = 1, 2, where r and θ are the radial coordinates in the (x 1 , x 2 ) plane and n ∈ Z is the vorticity of the solution. These static and straight field configurations do not depend on time and x 3 coordinates. The vacuum values of the condensates, v 1 and v 2 , are the solutions of Eqs. (18). The consistency of the coupled solutions at m 5 = 0 imply that the vortices in φ 1 and φ 2 condensates should be superimposed on each other and they should have the same winding numbers n 1 = n 2 = n in Eq. (69). Below we concentrate on the non-Hermitian model which is the subject of our paper (the analysis of the Hermitian counterpart can also be performed in the same way).
The radial profiles of the vortices are described by the functions f a with the following asymptotics: which guarantee that the total energy of the vortex solution is converging both at the spatial infinity and at the origin, respectively. The classical equations of motion (14) lead to the fol-lowing system of equations for the profile functions: which, evidently, do not possess straightforward analytical solutions.
Close to the origin, the last non-linear terms in both equations (71) can be neglected and the differential equations can be linearized. The solutions can be represented in the form of the polynomials, which involve only the positive even powers of the radius starting from the power r n determined by the vorticity number n = 1, 2, . . . . The ansatz (72) is thus consistent with the asymptotics (70). One gets for the first three coefficients: where A 1 and A 2 are free parameters of the solution which cannot be fixed at this stage. The series (73a) of the f 1 and f 2 profile functions are related to each other by the flip of the sign in front of the off-diagonal mass term m 2 5 . In the large-distance region, r → ∞, the asymptotics (70) imply f a = 1 − h a where |h a | 1 at sufficiently large distances. The linearized equations of motion, suggest that their solutions can be represented in the following form: The self-consistency of the solutions provides us with the power s = −1/2 of the algebraic prefactor r s and also imposes two simultaneous constraints: These equations give us two possible solutions for the ratio of the coefficients B (0) a from Eq. (76): and also determines the common exponent: We denoted for brevity: The next-order correction to Eq. (76) can also be easily obtained, where If the off-diagonal mass vanishes, m 2 5 = 0, the non-Hermitian two-scalar model reduces to two noninteracting scalar models L(φ 1 , φ 2 ) = L 1 (φ 1 ) + L 2 (φ 2 ) with This single-field model possesses the asymptotic solutions of the form (we omit the index a for simplicity): f (r) = A r n + m 2 4(n + 1) r n+2 + m 4 32(n + 1)(n + 2) where the asymptotic behavior is controlled by the mass of the single scalar field: In this limit, the two equations (77) decouple, and the mass parameters reduce to Eq. (85) for each field. The asymptotics of the non-Hermitian solution (72), (73), (81), and (82) are consistent with the single-field solution (84) as well.
Our numerical analysis confirms the existence of the stable vortex solutions in the regions of the phase diagram with non-zero condensates. An example of the profile functions for a set of coupling constants is shown in Fig. 3. All the radial n = 1 profiles of the vortices exhibit the same qualitative features, the linear rise close to the origin and the exponentially slow approach of the corresponding vacuum expectation values at large distances. These properties reveal the generic behavior of all solutions that we have analyzed.
The energy density of the non-Hermitian vortex (calculated per unit vortex length), +m 2 1 φ 2 1 + m 2 2 φ 2 2 + λ 1 φ 4 1 + λ 2 φ 4 2 , can be simplified with the use of the corresponding equations of motion (14). It can be expressed via the profile functions f a as follows: where the energy is normalized in such a way that E = 0 in the absence of the vortex. The very same expression (87) also gives us the energy of the vortex in the counterpart Hermitian theory (7), stability areas, with m 2 1 < 0 and m 2 2 > 0 in Fig. 4(a) and with m 2 1 , m 2 2 < 0 in Fig. 4(b) at the same values of quartic couplings. The vortex energies in the Hermitian and non-Hermitian versions trivially coincide at m 5 = 0 and then they tend to separate as the off-diagonal mass m 5 increases. One notices non-monotonic behavior of energies in the completely broken (m 2 1 , m 2 2 < 0) part of the phase diagram.

V. DISCUSSION AND CONCLUSION
Our paper discusses the phase diagram and reveals the properties of line-like topological excitations and vortices in the non-Hermitian relativistic model of two interacting scalar complex fields. We concentrate on non-Hermitian parity-time (PT ) symmetric realization of the model, which describes an open system that communicates with an external environment with precisely balanced gains and losses. The exact balance is encoded in the PT symmetry of the model Lagrangian. This symmetry ensures the stability of its steady ground state provided that the PT -invariance of the model is not broken spontaneously.
The model describes the properties of two scalar condensates that can exhibit the spontaneous breaking of the discrete PT symmetry in addition to the spontaneous breaking of the continuous global U(1) symmetry. Given the presence of the two symmetries that can be either broken or unbroken, the model thus contains four different phases, which makes the phase structure of the model very nontrivial.
The phase diagram of the non-Hermitian twocondensate model is shown in Figs. 1 and 2. We compared the non-Hermitian model with its Hermitian counterpart with identical couplings (Fig. 1), and we have observed the striking difference between the two cases. While the Hermitian model hosts the stable U(1) broken phase almost at every point of the phase diagram (except for a thin line, where the condensates decouple), the non-Hermitian diagram shows a wide variety of intertwined stable and unstable phases with both broken and respected U(1) symmetry. The richness of the phase diagram is guaranteed by the nontrivial pattern of the spontaneous PT -symmetry breaking, which is generalized, in our paper, to the interacting theory.
On the practical level, the stability of the ground state can be identified via the absence of the negative modes in fluctuation matrix M 2 NH that describes the quadratic fluctuations of the fields (25) over the ground-state condensates. On the contrary, in the spontaneously PTbroken regions, the quadratic fluctuation matrix contains at least a single negative eigenvalue.
In addition to the phase diagram, we studied the basic properties of global vortices in the non-Hermitian model in various regimes.
Firstly, we considered the model in the limit where the lengths of the condensates are frozen and the analytical analysis is feasible. To reveal the vortex properties, we used a set of exact transformations of the field variables that did not involve the explicit solution of the equations of motion. Noticing that the model resides in the PTbroken phase where the vortices are expected to be unstable, we show that the superfluid vortices can propagate non-dissipatively if and only if the vortex singularities in different condensates have the same vorticity (winding number) and, in addition, they overlap. In the case of the strong inter-condensate coupling determined by the non-Hermitian mass m 5 , the joint vortex segments interact via a long-ranged exchange of a massless excitation, similarly to the vortices in a Hermitian one-condensate model. The dissipation rate of the individual (separated) vortex segments is controlled by the off-diagonal mass m 5 , which, in turn, determines the interaction between the condensates.
The behavior of the vortices in the Hermitian and non-Hermitian cases is qualitatively similar: both of them are forming pairs. However, the similarity is not exact for two reasons: In the Hermitian model, the binding of vortices proceeds via a formation of the domain wall that emerges between the vortices. In the non-Hermitian case, the binding appears due to the dissipative dynamics of the vortices until they overlap. Moreover, the confining domain walls exist at low inter-condensate coupling values contrary to the non-Hermitian case in which the dissipative vortex binding is realized at the large intercondensate coupling.
Secondly, we also studied the classical vortex configurations at finite quartic couplings. In order to identify the classical configurations, we used a single set of classical equations of motion, which is obtained by the variation of the action with respect to the original fields. We omitted the equivalent but the incompatible, complimentary set of equations that correspond to the action variation with respect to the conjugated fields. This procedure, which follows Refs. [7,36], seems appropriate for the open systems residing in a steady-state regime, which does not necessarily coincide with the (thermal) equilibrium. Moreover, in this approach, the classical solutions possess a real-valued energy spectrum bounded below. The latter property is essential on the practical level as we search for the classical states using a (numerical) procedure based on the energy minimization as a criterion for the true (ground) state.
An alternative approach, based on the similarity transformation, does not possess the incompatibility of the two sets of equations of motion. This property makes the analytical procedure of finding the classical solutions more elegant [22,23]. However, the same approach makes the kinetic term (of, at least, one of the fields) negatively valued, leading to the emergence of a negative quadratic mode for the classical the solutions. Therefore, the mentioned class alternative solutions correspond to different, sphaleron-type saddle-point configurations, which can be significant for the thermal properties of the system.
We found that the PT -symmetric two-component model admits the vortex solutions inside and at the border of the PT -broken regions. These two-condensate vortex solutions share similar behavior with the vortices in the one-component relativistic superfluids. For consistency of the classical solution, the vortices of different condensates should have the same position in space-time and possess the same vorticity (winding number).
Finally, we proposed that a non-Hermitian twocondensate system can be realized in the Bose-Einstein condensates of atoms in two hyperfine spin states with the asymmetric Rabi coupling Ω, which couples to the different spin components with a different sign.