Bose-Einstein condensation processes with nontrivial geometric multiplicites realized via ${\cal PT}-$symmetric and exactly solvable linear-Bose-Hubbard building blocks

It is well known that using the conventional non-Hermitian but ${\cal PT}-$symmetric Bose-Hubbard Hamiltonian with real spectrum one can realize the Bose-Einstein condensation (BEC) process in an exceptional-point limit of order $N$. Such an exactly solvable simulation of the BEC-type phase transition is, unfortunately, incomplete because the standard version of the model only offers an extreme form of the limit characterized by a minimal geometric multiplicity $K=1$. In our paper we describe a rescaled and partitioned direct-sum modification of the linear version of the Bose-Hubbard model which remains exactly solvable while admitting any value of $K\geq 1$. It offers a complete menu of benchmark models numbered by a specific combinatorial scheme. In this manner, an exhaustive classification of the general BEC patterns with any geometric multiplicity is obtained and realized in terms of an exactly solvable generalized Bose-Hubbard model.


Introduction
The conceptual appeal of PT −symmetry (i.e., of the parity times time reversal symmetry) currently influences several different areas of theoretical and/or experimental physics [1,2] as well as of the related mathematical physics [3,4]. The birth of the concept dates back to the mathematics of perturbation theory [5,6,7,8] but the real start of popularity was only inspired by the Bender's and Boettcher's 1998 conjecture [9] that PT −symmetry of a Hamiltonian H could play a key role in a "non-Hermitian" formulation of quantum mechanics of bound states (cf., e.g., reviews [10,11,12]).
One of the basic innovations connected with the use of non-Hermitian but PT −symmetric Hamiltonians with real bound-state spectra has been found to lie in a remarkable and highly exciting possibility of a direct and experimentally controllable access to the dynamical regime of a spontaneous breakdown of the symmetry. Bender with Boettcher [9] were probably the first who managed to simulate this process (leading to a quantum phase transition, i.e., to an abrupt loss of the observability of the energy) using various elementary one-dimensional single-particle local potentials. This proved inspiring and influenced the model building efforts in multiple areas of realistic phenomenological considerations. Among others, several research teams turned also attention to a family of multiparticle PT −symmetric benchmark Hamiltonians H called Bose-Hubbard models (see, e.g., [13,14]).
The source of the appeal of the latter models lied in the possibility of a sufficiently realistic simulation of a specific quantum phase transition called Bose-Einstein condensation (BEC). From our present point of view the most inspiring analysis has been performed in Ref. [15] where a thorough mathematical description of the BEC phenomenon has been performed using both the linear (i.e., solvable, mathematically less complicated) and non-linear (i.e., fully general while just perturbative or purely numerical) versions of the underlying standard PT −symmetric Bose-Hubbard Hamiltonian.
In both of these model-building arrangements the process of the condensation has been attributed, in the Kato's language [16], to the presence of the higher-order exceptional points. In our present paper we intend to study the former mechanism of the quantum phase transition in more detail, extending the scope of the approach to certain more general dynamical scenarios characterized by nontrivial, optional geometric multiplicities of the generic exceptional-point degeneracies.

Bose-Einstein condensation 2.1 PT −symmetric Bose-Hubbard model
For introduction let us follow the guidance provided by Graefe et al [15] who performed a detailed perturbation-approximation analysis of the PT −symmetric Bose-Hubbard (BH) system living in a small vicinity of its BEC dynamical singularity. One of the most elementary versions of their non-Hermitian and BEC-supporting family of Hamiltonians had the following one-parametric form written in terms of the conventional creation and annihilation operators, As long as such a Hamiltonian commutes with the operator of the number of bosons (with eigenvalues N B = 1, 2, . . .), the model becomes exactly solvable, with spectrum E n (γ) = n 1 − γ 2 , n ∈ S(N) , N = N B + 1 where At γ ∈ (0, 1) the most characteristic features of such a spectrum are its reality, equidistance and up-down symmetry: see, for illustration, Fig. 1 where we choose N = 8. In the picture we see how the spectrum shrinks from the γ = 0 maximum {−7, −5, −3, −1, 1, 3, 5, 7} to a complete degeneracy at γ = γ (BEC) = 1. A detailed derivation of formula (3) itself can be found in [15]. In the framework of our present project it is important that the operator version (1) of the unperturbed BH Hamiltonian may be studied, most efficiently, in its explicit infinite-dimensional block-diagonal matrix form as derived in [15], Each block H (N ) (γ) is a tridiagonal N by N matrix such that etc. This representation of the Hamiltonian reconfirms the conservation of the number of bosons N B = N − 1. At every N, the bound-state Schrödinger equation degenerates to a linear algebraic eigenvalue problem yielding the closed-form eigenvalues (7). In this representation of the system also the construction of the wave functions becomes tractable non-numerically [17].

BEC-formation patterns
The properties of spectrum (3) are reflected by the equidistance of the N−dependent set (4) of the energy-level quantum numbers. The degeneracy of all of the levels in the BEC limit γ → 1 is simultaneous, lim An analogous degeneracy also controls the behavior of all of the eigenstates |ψ n (N, γ) of the BH operator (1), lim Precisely such a degeneracy can be interpreted as one of the explicit realizations of the so called Kato's exceptional point of order N (EPN, [16]). The physical limit γ → γ (BEC) = 1 and the mathematical limit γ → γ (EP N ) = 1 appear synonymous. Unfortunately, such a correspondence between mathematics and physics is one-sided.
In the context of mathematics the convergence (8) towards a single limiting vector |χ (N ) (BEC) must be considered a simplification and an artifact of the model. In any sufficiently model-independent scenario it will be necessary to replace property (8) (where the EPN merger has the so called geometric multiplicity K equal to one) by the more general, K−centered eigenvector-degeneracy pattern with any geometric multiplicity K ≥ 1, The complete set S(N) of the energy quantum numbers must be admitted to be, in general, decomposed into a K−plet of disjoint subsets. Now, the main question to be addressed in our present paper is whether the generalized mathematical form (9) of the EPN limit with arbitrary K ≥ 1 could still be assigned a suitable physical multi-bosonic interpretation via an exactly solvable BEC-supporting Hamiltonian, say, of the BH type.
3 Nontrivial geometric multiplicites K at small N In our recent heuristic study [18] we addressed the problem on a purely methodical level. In a trial-and-error search for the low-dimensional toy-model matrix Hamiltonian with nontrivial geometric multiplicity of its degenerate EPN limit we used the brute-force numerical methods and performed the search among certain real and more or less randomly selected matrices H (6) . The inspection of the parametric-dependence of the spectra led us to the empirical conclusion and conjecture of a one-to-one correspondence between the trivial geometric multiplicity K = 1 of the exceptional points and the tridiagonality of the matrix H (6) in question. This turned out attention to non-tridiagonal toy models H (6) . The hypothesis appeared to work. Several models supporting the existence of exceptional points with K = 2 and K = 3 were identified. Unfortunately, our subsequent tentative extension of the brute-force empirical search to higher N encountered severe technical obstacles. Our conclusions formulated in [18] were, therefore, discouraging and sceptical. We argued that beyond N ≈ 6, the numerical EPN searches become ill-conditioned, yielding highly unstable results marred, in the methodical context, by the influence of the ubiquitous rounding errors (for details see, in particular, Appendix A in loc. cit.).
In our present paper we decided to accept a different strategy. It will be based on a nonnumerical, strictly analytic specification of the candidates for the suitable EPN-supporting toymodel Hamiltonians. For this purpose we will merely rescale the separate fixed−N components of the K = 1 Hamiltonian (5) yielding the trivially modified tilded forms of its submatrices, viz., etc. Obviously, any such a rescaling with real c k = c (N ) k = 0 preserves not only the basic features of physics behind the model but also, at any dimension N = N B + 1, the exact solvability of Schrödinger equation. At the same time, the new, alternative, non-tridiagonal candidates for the submatrices of the generalized K > 1 versions of Hamiltonian (5) can tentatively be constructed using suitable direct-sum combinations of the building-blocks (11).

Three bosons and K = 2
The first two simplest BH models with N B = N − 1 = 1 and N B = N − 1 = 2 are trivial, offering just the usual unique option with K = 1. Thus, our search for alternatives has to start from the next, N = 4 item taken from the K = 1 sequence of Hamiltonians (6), with the spectrum displayed in Figure 2. In the picture we marked the levels by symbols "I" and "II" in order to indicate that precisely the same γ−dependent spectrum can be also generated by an alternative Hamiltonian defined as a direct sum of the two N = 2 matrices taken out of the tilded menu (11). This yields the innovated, K = 2 Hamiltonian [K=2] (γ) (12) characterized by the choice of scalings c 1 = 1 (levels "I") and c 2 = 3 (levels "II"). Besides sharing the spectrum (but not the eigenvectors!) with its BH predecessor H (4) [K=1] (γ) of Eq. (6), the new model also illustrates several other merits of the present philosophy. The key point is that the two alternative isospectral Hamiltonians can be both made PT −symmetric, and that their imaginary parts can be kept the same (i.e., diagonal). In this sense we can treat the right-hand-side matrix in (12) as a canonical non-tridiagonal form of the K = 2 alternative to the original K = 1 BH building block at N = 4.
The existence of the two isospectral models with different eigenvectors may be interpreted as forcing us to introduce a new, ad hoc quantum number (say, "color"). Certainly, any less formal specification of such a quantum number (or numbers) would have to be left to the experimentalists. Near the BEC dynamical regime only a dedicated experiment performed over the N B −plet of bosons will be able to distinguish between the systems with different structures of wave functions. With their energy spectra fitted by the K−independent formula (3), the measurements would also have to extract information about some wave-function-dependent mean values of some other suitable observables.

Four bosons and K = 2
At N = 5 the K = 1 Hamiltonian H (5) (γ) of sequence (6) yields the "colorless" five-level spectrum as displayed in Figure 3. The elementary combinatorics reveals that besides such an option we can now introduce the two alternative labelings "I + II" and "III + IV" of the states with K = 2.
Having opted for "I + II" (see the picture) we are able to reproduce the spectrum by the direct (where c 1 = c 2 = 2) while in the second scenario we choose c 1 = 1 and c 2 = 4 and have We again decided to arrange the basis in such a manner that the main diagonals remain unchanged. This renders the resulting Hamiltonian submatrices sparse, multidiagonal, complex symmetric [19] and manifestly PT −symmetric.

Canonical representation
In the light of Eq. (9) the non-triviality property K > 1 of our Hamiltonian matrices means, at any N, that their canonical representation must have, in the EPN limit, the partitioned form of a direct sum of the standard M j by M j Jordan-block submatrices such that Thus, one of the basic characteristics of the dynamics of the system near γ (EP N ) will be an ordered version of partitioning (17). Up to N = 8, the sets of these partitionings (having P (N) elements) are sampled in Table 1.
The trivial partitionings π 0 (1, N) = [N] form the leftmost column of the Table. The related single, unpartitioned Jordan block (16) with M j replaced by N is also most often encountered in the literature. As the simplest special case of the canonical Hamiltonian it is also the only option available at N = 2 and N = 3. At N > 3 the option becomes complemented by the simplest K = 2 matrix (15) characterized by the m = 1 partition π 1 (2, N) = [N − 2, 2], etc. Thus, the complete menu of the available partitions has two elements at N = 4 and N = 5, four elements at N = 6 and N = 7, etc. The growth of the number P (N) of the partitions accelerates at the larger matrix dimensions N (see Table 2 and/or Ref. [20]).

Transition matrices
The BEC-and EPN-admitting BH-type Hamiltonians H [N ] (γ) of our present interest have to be constructed along the lines as sampled in section 3 above. Once we specify their characteristic partition (18) for the conventional Schrödinger equation. The superscript π m (K, N) has been added here to refer to Eq. (18), i.e., to the partition [M 1 , M 2 , . . . , M K ] specifying the direct sum J [N ] (η) of Jordan blocks in (15).
In what follows we will assume that at a fixed number of bosons N B = N − 1 the partitioning π m (K, N) is an inseparable part of a dynamical input information about the system even at γ = γ (EP N ) . Such an input information will be carried by the preselected N by N BH-type toy-model Hamiltonian matrix H [N ] (γ), tractable as a small perturbation of its EPN limit H [N ] (γ (EP N ) ). Subsequently, in a way explained in [17], the γ−independent transition matrices Q can be still perceived as certain formal analogues of unperturbed basis [21].
In a sufficiently small EPN vicinity with γ ≈ γ (EP N ) , the Hamiltonian matrices themselves may be assumed composed of an exactly diagonalizable unperturbed component H [N ] 0 (γ) (equal to a suitable BH-type direct sum of the tilded-matrix components (11)) and a perturbation. In the EPN limit, such a perturbation must disappear. This means that for our present considerations the detailed form of this perturbation is irrelevant. We will, therefore, ignore its influence and drop the zero subscript of the relevant BH Hamiltonian as, for our present purposes, redundant,

Geometric multiplicities K > 1: realization
Our attention will be now narrowed to the study of the BH-type systems of bosons in which the hypothetical available experimental information consists of our a priori knowledge of the conserved number of bosons N B and of the γ−dependence of the spectrum as given by formula (3). Let us emphasize that such a constraint is well supported not only by its occurrence in multiple experimental setups (where the equidistance of the spectrum usually reflects its "vibrational" character and interpretation) but also by the significance of its purely formal merits. In [22], for example, we showed that even the simplest non-equidistant choice of the square-well-type unperturbed spectrum makes the localization of the EPNs technically much more complicated.
The equidistance assumption (3) was, initially, a rigorous result of study of the BH model (1). Such a property of the system with trivial geometric multiplicity K = 1 was prescribed by formula (4) for indices. In our study of the K ≥ 1 scenarios we now intend to proceed by analogy. Having in mind the importance of the exact solvability we will complement the (entirely general) directsum decomposition (10) of the set of indices S(N) by a more specific, BH-motivated additional assumption by which all of the separate components S This means that the optional real multipliers c k will serve here as a source of an adaptive rescaling of the subspectra. The resulting enhanced flexibility of these subsets of quantum numbers will be paralleled by a rescaling (11) of the original BH sub-Hamiltonians (6). In this manner, in a way inspired by the small−N constructions of section 3, our general one-parametric N by N Hamiltonian matrices will be defined as follows, This is our ultimate ansatz. It admits, by construction, the BEC-related EPN singularities of arbitrary geometric multiplicities K ≥ 1. Up to N = 7, an exhaustive list of the mutually isospectral Hamiltonian matrices of this type is given here in Table 3.

Physics behind the generalized BH model
In the standard terminology the two annihilation operators a 1 and a 2 and the two creation operators a † 1 and a † 2 in the initial BH Hamiltonian (1) correspond to the two dominant bosonic modes which may tunnel through a hypothetical unit barrier [15]. The parameter γ is left variable in order to characterize the phenomenologically most relevant imaginary part of the on-site bosonic-energy difference (see also several other studies of such a most elementary BH bosonic system in [13,23,24]). From the point of view of physics the apparent robustness of the K = 1 property of such a family of solvable BH-type models has long been perceived as disappointing. This seemed to imply that the solvable picture of the BH-type quantum dynamics can only be based on the conventional BH Hamiltonian (1), not offering a description of an EPN collapse with a non-trivial geometric multiplicity.
We have shown that an appropriate K > 1 generalization of the model will always require the direct-sum block-diagonal structure (15) of the canonical Hamiltonian in the EPN limit. In this sense the highly desirable and still sufficiently realistic BH-type realization of the general scenario has been found to be offered by our present model (21). We managed to keep it exactly solvable, and we also expect a reopening of the questions of its predictive power and/or tests in the laboratory. In this direction a few immediate comments are worth adding.

Change of phase: two alternative physical interpretations
In the conventional quantum mechanics of unitary evolution the singular values γ (EP N ) were only a mathematical curiosity. Useful in perturbation theory [16] but playing hardly any significant role in phenomenology. Indeed, as long as the conventional quantum theory had to be formulated in a fixed, pre-selected Hilbert space (say, K), the basic postulate of the diagonalizability of any meaningful quantum Hamiltonian H(γ) appeared manifestly incompatible with any form of the EPN-related degeneracy of the eigenstates. This conclusion also sounded consistent with the obligatory account of the unitarity and Stone theorem [25].
The change of paradigm and the turn of EPNs into one of the fundamental physical concepts occurred when Bender with coauthors [8,9] initiated the study of certain non-Hermitian models with relevance in the theory of relativistic quantum fields. The obvious phenomenological appeal of several non-Hermitian but Hermitizable innovative forms of interactions opened a new and promising direction of research and discoveries in physics [11] as well as in mathematics [3]. With time, these developments resulted in an extension of the scope of the traditional realistic models as well as in multiple proposals of the new, non-Hermitian but Hermitizable Hamiltonians (see, e.g., several most recent reviews in [1,2]).
In particular, a compatibility of the new theory with the Stone theorem has been achieved via an ad hoc amendment of the inner product in K. As a consequence, the later space proved converted into a new, strictly physical Hilbert space (say, H) in which the Hamiltonian re-acquired its necessary self-adjoint status [12]. In our present paper the main consequence of the perspective of working with the two Hilbert spaces K and H is twofold. Firstly, the non-Hermiticity of H(γ) in K admitted the existence of the EPNs at a real parameter γ (EP N ) . Secondly, the mathematical consistence of the theory has been found achieved by the observation that H simply ceases to exist in the singular limit of γ → γ (EP N ) .
Formally speaking, the admissibility of the reality of the EPNs was a consequence of the Hermitizability of H(γ) or, indirectly, of its PT −symmetry [11]. Still, it was also necessary to take into account that in general, the physical phenomena covered by non-Hermitian models may be found extremely sensitive to random perturbations [21,26,27]. For the latter reason, it is necessary to distinguish between the implementations of the theory, mainly in the two entirely different experimental setups. In one (let us call it, for our present terminological purposes, "approach A"), many theoreticians are trying to predict the qualitative features of the dynamics of the system after its perturbation. Interested readers can find an extensive sample of the results of such a type in paper [15].
In that paper the authors described several alternative unfolding scenarios of the EPN degeneracy, i.e., several eligible patterns of the complexification of the energy spectrum under a "strong", unitarity-violating perturbation. In the context of our present paper we would slightly prefer an alternative philosophy (let us call it "approach B"). In contrast to approach A (which analyzes, basically, the non-unitary effective open-system dynamics), this approach is restricted to the analysis of the unitary, closed-system behavior of the quantum system of interest. Before it passes through its EPN-mediated phase transition, and before the reality (i.e., the observability) of the energy spectrum becomes lost.
One of the physics-related consequences is the emphasis put upon the conservation-law role of operator (2) which describes the number of bosons in the system. We declared such a number a strictly conserved quantity. Under this assumption we are allowed to work with a fixed and finite number of bosons N B . Also the number of the energy levels is kept finite and equal to N = N B + 1. Last but not least, it is worth adding that even the inclusion of perturbations changing the number of bosons need not necessarily make the model prohibitively complicated (in this respect, interested readers could consult, e.g., Ref. [28]).

The role of PT −symmetry
The modern applications of quantum theory range from condensed matter physics and materials science up to chemistry and engineering. In this context many theoretical methods often become popular with a certain delay. One of the good examples is the concept of parity-times-timereversal symmetry (PT −symmetry) which was initially perceived as a curiosity possessing just a few purely abstract mathematical applications (say, during the studies of some subtleties in perturbation theory -see, e.g., Refs. [5,6,29]). Needless to add, the phenomenological relevance of the concept of PT −symmetry is now widely recognized and quickly growing (see, e.g., several most recent reviews collected in books [1,2]). At present such a tendency still keeps being productive in several traditional areas of physics. It offers new insights also in some subtle mathematical aspects of the theory of quantum phase transitions.
In spite of the manifest non-Hermiticity of the PT −symmetric candidates H for Hamiltonians, these operators were shown eligible as generators of unitary evolution [10,12]. In this context, one of the basic methodical assumptions accepted in the current literature on BH models [17,18,30] was that the infinite-dimensional matrix (5) as well as all of its separate submatrices (6) had to be complex symmetric, tridiagonal and PT −symmetric, with P equal to an antidiagonal unit matrix, and with symbol T representing an antilinear operation of Hermitian conjugation (i.e., transposition plus complex conjugation). This led to the conclusion (or rather conjecture) that in the EPN limit (i.e., at the instant of the loss of diagonalizability), the canonical representation of every N by N submatrix H (N ) (γ (EP ) ) can be given the form of the N by N Jordan matrix (16) with, due to PT −symmetry, η = 0.
In mathematical terminology this seemed to imply that in the BH models, the geometric multiplicity K of the EPN in question had to be, at any preselected value of N B = N − 1, always equal to one [18]. And it is precisely this belief and scepticism which are disproved by our direct-sum ansatz (21).

Combinatorics behind the classification
The relevance of the BEC models with K > 1 is twofold. First, the differences between partitions (18) may inspire new experiments, especially in the dynamical regime close to the EPN singularity. Second, the number of non-equivalent partitionings (18) as well as the number of admissible Hamiltonians (21) are both quickly increasing functions of the matrix dimension N. Obviously, in the language of mathematics, the occurrence of non-minimal geometric multiplicities K > 1 may be expected to dominate, especially among the random Hamiltonian matrices. For both of these reasons, the currently open problem of the classification of the models at large N would certainly deserve an enhanced attention.

Classification scheme
At the not too large matrix dimensions N, our present condition of equidistance of the energy levels remains reasonably restrictive. At the smallest N one even encounters not too numerous alternative isospectral BH models with Hamiltonians (21). Once we denote the number of nonequivalent BEC-degeneracy scenarios by a dedicated symbol a(N), we may notice that the growth of this value at intermediate N still remains comparatively slow, especially at the even N (see Table 4). Still, asymptotically, the growth of the value of a(N) with the growth of N becomes exponentially quick. Even then, the sequence a(N) keeps exhibiting a parity-related irregularity (see [31] and [32] for details).  For the purposes of an explicit classification of the models, it is necessary to consider not only the (conserved) number of bosons N B = N B (N) = N − 1 and the geometric multiplicity K of the BEC degeneracy but also the selection of partition π m (K, N) and of the set of scaling parameters {c 1 , c 2 , . . . , c K }. All of these choices have to be made compatible with the directsum decomposition formula (21). The resulting exhaustive list of the Hamiltonians is sampled in Table 3. The Table offers a clarification of the relationship between the scarcity of our initial BH Hamiltonians in their matrix representations (5) or (6) (where we always had K = 1) and the abundance of their K > 1 EPN-supporting generalizations H [N ] (γ).
Although the construction of Table 3 was just an elementary combinatorics, its extension to the larger dimensions N is still an open question. The construction becomes tedious, well suited for the computer-assisted enumeration. Incidentally, the suitable algorithms proved different for the odd and even N. Their most efficient versions were proposed by Andrew Howroyd and may be found published, in the on-line encyclopedia of integer sequences, under the respective coordinates [31] and [32].
Our spectrum-equidistance constraint restricts, in Hamiltonian (21), the freedom of our choice of the partitions π m (K, N) and of the scalings c k . From a descriptive and phenomenological point of view, the family of the resulting generalized BH Hamiltonians still remains sufficiently rich. In the manner illustrated in Table 3 one only has to notice that at a fixed (i.e., unrestricted but conserved) value of N, a straightforward classification of these Hamiltonians is not provided by the multiplicity K (even at N = 5, the K = 2 option is already shared by the two different Hamiltonians) nor by the partitions (at N = 6, partition [3,3] is not realized at all) but only by their combination with the eligible multiplets of scalings {c 1 , c 2 , . . . , c K }.

An alternative notation
A rather abstract nature of the latter criterion indicates that a more physics-oriented guide to the classification could and might be sought in a return to the explicit spectrum (3). We know that the unique) K = 1 element H (N ) (γ) of the BH-matrix sequence (6) is unambiguously specified by the spectrum. In this sense the Hamiltonian matrix becomes uniquely identified by the γ = 0 spectrum or by the symbol S(N). Thus, with N = 5, for example, we have S  (14) and (13).
The latter, more intuitive and spectrum-representing version of notation admits an abbreviation which has been used in Refs.  (3q)(5q) . . .}. In this notation, therefore, the "colorless" K = 1 building-block items forming the sequence (6) The first few samples of the modified classification may be found displayed in our last Table  5. More comments on such a notation may be also found in [31] and [32]. In particular, these references offer an algorithm for the evaluation of the "count of the colors" a(N) at any N.

Discussion
In our paper we managed to clarify certain qualitative properties of the BEC degeneracies. For the sake of definiteness we had to assume that the set of indices S(N) was defined by Eq. (4). Nevertheless, only our independent requirement of the BH-type physics extended the same requirement to all of the components of the direct sums (10). Otherwise, the form of these subsets S (EP ) k would remain less constrained, numbering strictly just the eigenstates which had to coincide with the k−th ket |χ (EP N ) k in the limit. Thus, every such a subset would only be required to contain two or more quantum numbers n.

The specific features of bosons
In the light of the latter comment we believe that the price to pay for the enhanced freedom (consisting in the loss of the connection of physics with the BH model) would be too high and hardly acceptable. Indeed, the BH model mimics the experimentally realizable double-well arrangement in which, say, the cold bosonic atoms are injected in one well and, simultaneously, extracted from the other one [33,34]. Moreover, the model also reflects the key difference between fermions and bosons because the "natural" Pauli exclusion principle only applies to the the former class of the quantum particles. Bose with Einstein [35] were among the first to notice that in the systems of bosons, this opens the possibility of a specific quantum phase transition (which is, at present, widely known as the Bose-Einstein condensation [36]).
Both the theoretical and experimental aspects of the BEC idea proved particularly useful and inspiring in the condensed matter physics. The condensates were successfully described there, typically, by various nonlinear forms of Schrödinger equation [37]. In our present paper we paid attention to the less widespread simulation of the condensation processes in which the underlying dynamical equations are required PT −symmetric, i.e., invariant with respect to the simultaneous action of parity P and of the antilinear time reversal T [11]. As we already indicated above, the scope and impact of such a methodical innovation is not yet fully explored, with the existing results ranging from an amendment of our understanding of the topological phase transitions (say, in a non-Hermitian Aubry-Andre-Harper model [38] or in quasicrystals [39]) up to the new approaches to the perception of the conservation laws [40], and from the theoretical studies of the interference between channels [41] and of the mechanisms of squeezing [42] up to the detailed, experiment-oriented simulations of the properties of the specific BEC-type condensates [34]. In this framework, promising results are also being obtained in the area of related mathematics.
Along these lines we showed that although the model (1) is, in some sense, naive, its exact solvability represents an important advantage. In fact, precisely this property encouraged us to propose a generalization which enhanced the theoretical scope while still remaining realistic. We may summarize that our "benchmark" amendment of the model seems to offer a fairly satisfactory picture of the physical BEC-related reality which is obtained, in addition, by the purely nonnumerical means.

The problem of the non-uniqueness of the model
In the above-mentioned open-system-theory approach A (which may be found more thoroughly explained also in monograph [43]), the one-parametric PT −symmetric BH model (1) can be interpreted as one of the phenomenological descriptions of a realistic and potentially unstable quantum system of bosons in which their number N B itself is conserved. In such an effectivedescription approach the model simulates the behavior of an open system exposed, say, to the influence of an environment in a way characterized, globally, by the parameter. The manifest non-Hermiticity and non-unitarity of the model reflects, effectively, the randomness of the environment.
In our presently preferred approach B, the physical interpretation of the model is different, treated fully in the spirit of the original Bender's philosophy [11]. This, naturally, enhances the impact of PT −symmetry, and it changes also a part of the related mathematics. First of all, the conventional Hilbert space K must be declared unphysical because in this space the evolution (controlled by Hamiltonian H which is non-Hermitian) would be, in the light of the well known Stone theorem [25], non-unitary. Incidentally, the apparent emerging paradox has a virtually elementary resolution: Along the lines discovered and discussed already in older literature [10,44], it is sufficient to endow the same space K with an amended inner product. Interested readers should search for the necessary mathematical details elsewhere [3,10,12]. For introduction it is just sufficient to keep in mind that the amendment of the inner product converts the unphysical Hilbert space K into its physical Hilbert-space alternative H in which H becomes, by construction, Hermitian alias self-adjoint.
In approach B, the first task is to demonstrate the existence of a suitable inner product (see, e.,g., an extensive discussion of this aspect of the theory in [45]), the necessary condition of which is the reality of the spectrum. This means that, in some sense, the two approaches A and B meet at the boundary where γ (EP N ) = 1. In both of the neighboring dynamical regimes, after all, the hypothesis of the existence of the EPN boundary opens an exciting theoretical possibility of description of an EPN-mediated phase transition in the system.
In our study, we had mainly in mind approach B. We considered, exclusively, the positive and unitarity-compatible parameters γ < γ (EP N ) . In such a subinterval of parameters it is possible to guarantee that the necessary physical Hilbert space H does really exist because the spectrum of H is real [12]. In such a context there emerges one of the most interesting theoretical questions: Would the experimental confirmation of the BH-model-based predictions mean and imply that the validity of the theoretical model is confirmed? Obviously, it is not so because in an inverse BH bound-state problem, the candidate for the simulation of a given energy spectrum is not unique. Constructively we showed that a given solvable BH-type Hamiltonian can be complemented by its many isospectral alternatives, yielding even exactly the same parameter-dependence of the energy spectrum. This implies that in a complete experiment one would have to measure also some wave-function-dependent features of the system.

A remark on the theory of quantum phase transitions
Before the BEC/EPN collapse (i.e., whenever |γ| < 1), the spectrum of all of our present BH Hamiltonians remains real and observable. At the instant of collapse, one can speak about a quantum phase transition [46]. After the system has passed through such an EPN singularity, one has a choice between the above-mentioned open-system evolution scenarios of approach A (admitting the loss of the reality of the spectrum) and the closed-system scenario B in which the escape from the EPN singularity proceeds through a fine-tuned corridor of unitarity as sampled, in the BH context, in [17].
In a broader phenomenological context, one of the most typical aspects of the change of phase of any quantum system can be seen in the loss of observability of at least one of its observable characteristics. Among the most popular scenarios of phenomenological interest one encounters the loss of observability of the energy levels, i.e., in Schrödinger picture [47], of at least some of the parameter-controlled eigenvalues of the Hamiltonian.
From the point of view of mathematics, it makes sense to prefer the study of the Hamiltonians which are analytic functions of the parameter, H = H(γ). In the generic, infinite-dimensional and purely numerical models the loss of the observability can be then perceived as caused, in a partial analogy with our present solvable models, by the limiting transition γ → γ (EP N ) in which the limiting value of the parameter is the Kato's exceptional point of some finite order N. One can say that at in such a limit there also exists a finite subset S = S (EP N ) of quantum numbers n for which the energies as well as the wave functions merge. Due to our present intention of keeping the models solvable, it was only necessary to decompose the discrete set S (EP N ) in the very special, BH-related direct sum (10). Naturally, the same method and approach should also work in the descriptions of the quantum systems which are not solvable exactly.

Summary
The authors of Ref. [15] emphasized the deep methodical relevance of model (1) in the vicinity of its BEC degeneracy. They even developed a specific ad hoc perturbation theory and described some of the consequences of the inclusion of perturbations. Several alternative BEC-related perturbationinfluence results can be also found in Refs. [48]. These studies clarified multiple quantitative aspects of the condensation.
In our present project we managed to construct a BH-based simulation of the general BECs in a solvable-model realization. The resulting picture of a generic EPN degeneracy appeared to have two aspects. The positive one is that in our innovative family of the BH-type models the EPN limit has been made nontrivial, characterized by any kinematically admissible geometric multiplicity K. In our methodically oriented considerations, on the other hand, our results only involved the study of the Hamiltonian. Our description of the system of bosons did not involve any other quantities which would be measurable and which would enable us to clarify the differences between our nonequivalent BH-type models. Thus, the clarification of these differences (characterized here just by a rather formal proposal of an auxiliary quantum number) will certainly be a task for some forthcoming studies of the BH-type systems of bosons.
From the mathematical point of view the most exciting aspect of our generalized BH models of BEC was, certainly, the discovery of the feasibility of the enumerative classification of all of the admissible direct-sum Hamiltonians (21). This was an interesting combinatorial task which found the necessary background in the specialized literature. In such a formulation one cannot say that the task is already completed. Indeed, the two Howroyd's brute-force algorithms of the enumeration of all of the phenomenologically meaningful unions (10) of the specific equidistant and symmetric spectral subsets (20) of the energy quantum numbers would certainly deserve to be complemented by some less implicit (e.g., recurrent) alternatives (provided only that they do exist at all).