Classical, semiclassical and quantum signatures of quantum phase transitions in a (pseudo) relativistic many-body system

We identify a (pseudo) relativistic spin-dependent analogue of the celebrated quantum phase transition driven by the formation of a bright soliton in attractive one-dimensional bosonic gases. In this new scenario, due to the simultaneous existence of the linear dispersion and the bosonic nature of the system, special care must be taken with the choice of energy region where the transition takes place. Still, due to a crucial adiabatic separation of scales, and identified through extensive numerical diagonalization, a suitable effective model describing the transition is found. The corresponding mean-field analysis based on this effective model provides accurate predictions for the location of the quantum phase transition when compared against extensive numerical simulations. Furthermore, we numerically investigate the dynamical exponents characterizing the approach from its finite-size precursors to the sharp quantum phase transition in the thermodynamic limit.


Introduction
The methods and ideas of quantum chaos [1,2] have provided deep insights into the way classical information conspires withh in a subtle manner. To large extent this can be understood within a semiclassical theory, explaining genuine quantum behaviour like entanglement and coherence. In this field, Shmuel Fishman made paramount contributions ranging from the celebrated explanation of dynamical localization as a type of Anderson transition in kicked systems [3] to the resummation of periodic orbit expansions to construct semiclassical approximations for individual eigenstates in chaotic systems [4]. The present contribution aims to express our admiration of his scientific work.
During the last decade, the field of quantum chaos has experienced an influx of new ideas coming from its application to the realm of interacting many-body systems. The newly emerging field of many-body quantum chaos is based on exciting developments in our understanding of fundamental problems like equilibration of closed systems [5][6][7][8][9] and the scrambling of quantum information due to classical chaos [10][11][12][13].
It is therefore not a surprise that semiclassical methods, both at the heuristic level of quantum-classical correspondence [14][15][16] and the level of asymptotic analysis of path integrals describing coherent quantum effects [17][18][19][20][21], have been lifted from its original particle-like form into the realm of quantum fields. Among the plethora of phenomena characteristic of the rich physics of interacting many-body systems, critical phenomena have always had a special place. In this new disguise, many-body semiclassical methods are a suitable tool to understand even the most delicate quantum effects related to the emergence of criticality.
A natural arena for testing this idea is the attractive Lieb-Liniger model [22] describing one-dimensional bosons attractively interacting through short-range forces and, in particular, its low-energy effective description that has been experimentally realized [23,24]. The reason for this is that this system displays a quantum phase transition [25][26][27][28] and admits a rigorous derivation of a well-defined and controlled classical limit in the form of mean-field equations, thus allowing for direct application of semiclassical techniques [21]. The semiclassical study of this system in [21] revealed the key role played by locally unstable mean-field dynamics in the corresponding dynamical and spectral quantum mechanical features.
The extension of many-body semiclassics beyond the realm of bosonic systems is still in its infancy, but a step in this direction is to first consider how the well-established picture of [21] gets modified by two new ingredients: a relativistic dispersion and the presence of spin-like degrees of freedom. Since the very possibility of having locally unstable dynamics (as opposed to global chaos) of the attractive Lieb-Liniger model is due to the integrability of the effective Hamiltonian describing its low energy regime, a natural question concerns possible non-integrable behaviour of such models and its consequences for the existence and characteristics of the quantum phase transition. In this paper, we answer some of these questions.
The paper is organized as follows. After we introduce the model and describe its general physical properties in section 2, we present the motivation for the transformation into a special Fock basis in section 3 and how this optimal transformation adiabatically fragments the Hamiltonian in section 4. After that, in section 5 the conversion of the channel containing the ground state into its classical form is examined. The most important results presented in the section 6 are the exact calculation of the critical interaction strength and the analysis of discontinuities in the functional dependence of the energy on the interaction. Finally, the asymptotic convergence of the first excited energy level towards the ground state level leading to a degenerate ground state in the mean field is quantified in section 7.

The Hamiltonian and its symmetries
The Hamiltonian of the (modified) Lieb-Liniger model with linear dispersion and contact potential is defined asĤ describing bosons on a ring with radius R with a contact interaction that can be interpreted as a mass term: The moment two bosons are at the same point they obtain a mass through the contact potential, whereas they are massless otherwise. In the following we assume attractive interactions, e.g. α > 0, and we will choose natural variablesh = 1, L = 2πR = 2π such that the unit of energy is [E] = 4π 2h2 L 2 [28]. As appealing as it is, it is important to note that the system above appears ill-defined, as its Hamiltonian (1) is not bounded from below. Unlike in fermionic systems, in this bosonic system this issue cannot be resolved by the introduction of a Fermi sea. One way out of the problem is to interpret (1) as emerging from a local approximation of a one-dimensional condensed matter or cold atom system with two crossing bands that is perturbed by an interband interaction. This naturally introduces a regularization of the noninteracting model with a single-particle momentum cutoff defining the region where the linearization is justified. In this approach, the linear dispersion is a property of excited states and has an effect on dynamical properties of states with a certain momentum. An example of such (local) Dirac bosons in two dimensions has been found in the collective plasmon dispersion in honeycomb-lattices of metallic nanoparticles [29]. In such local approximation, one has to make sure that any prediction of the model has to be independent of the cutoff, which might be realized in a quench scenario, starting with a narrow momentum distribution.
In order to proceed within a Fock space approach, we choose the eigenbasis of the non-interacting (α = 0) Hamiltonian as the single-particle basis where as orthonormal eigenbasis for the momentum operator we use plane waves as the most obvious choice. For the quasi-spin an orthonormal eigenbasis is used consisting only of "up" and "down" generated by the third Pauli-matrixσ From these definitions the Fock space is characterized through the occupation numbers n k,σ of the several states |k, σ with creation and anihilation operators satisfying canonical commutation relations where each pair of creation/annihilation operators defines an occupation number operator for the correspondig mode. With the help of these bosonic operators this leads, after truncation of the momenta from Z to {−1, 0, 1}, to the more convenient form with the relevant Fock states labeled by six occupation numbers, |n 1,+ , n 0,+ , n −1,+ , n 1,− , n 0,− , n −1,− .
One way to consider this system is to map the truncated model to a spin-one bose gas on two quantum dots (or two sites with suppressed hopping), where the physical spin takes the role of the momentum k = −1, 0, 1 and the pseudo-spin 1/2 labels the two sites that have opposite external magnetic fiels applied to them, introducing linear Zeeman splitting and thus the "three-mode linear dispersion". The interaction processes are then taking, e.g., two particles of opposite spin on the same site and distribute them into the spin-zero modes of the two sites. There are different processes, of course, but the overall interaction effect is a spin-mediated hopping of a single particle with the total spin (of the participating particles) being preserved. The noninteracting case would decouple the two sites. However, we take a different perspective here that takes the truncated model as it is, i.e., we truncate to the three lowest momenta and then assume that the ground state of this model represents a physical ground state. This Hamiltonian has a set of symmetries that will be the key for the adiabatic separation later on. We have the total number of particleŝ and the Hilbert space can be divided into sectors with the respective quantum numbers (N, L). To simplify the task we will focus on the special case of fixed N and L = 0. Except for the derivation of the effective Hamiltionian which is done for general L. In this way, the effective number of degrees of freedom is reduced from six to four. Besides these two symmetries, the energy spectrum splits up symmetrically in the positive and negative direction, as can be seen in Figure 1. For an even particle number N this observation can be explained using the operatorξ whereŜ is the total (pseudo) spinŜ that satisfies (−1) and therefore it is easy to show that Asξ is a bijection on the set of eigenstates |ψ ofĤ with energy E = ψ|Ĥ |ψ , there always exists a state |φ =ξ |ψ that is also an eigenstate ofĤ. The energy value corresponding to this state is then given by Finally, a parity operatorP can be defined which simultanously flips all spins and momenta, given by a complex conjugation to invert the momenta in the eigenbasis of plane waves followed by a spin flip, satisfyingP 2 = 1. Also, since [P,Ĥ] = 0,P represents a discrete symmetry that splits the Hilbert space into two separate subspaces leading to a separation of the energy spectrum into two independent subspectra 1 see Figure 1.
As a final remark, we note that the existence of further symmetries is ruled out by a numerical diagonalization and the analysis of avoided crossings, as indicated for N = 20, L = 0, P = 1 in Figure 2. The absence of real crossing suggests that there are no additional symmetries to be found which could be used to further reduce the dimensions of the Hamiltonian (8) [30].

Adiabatic separation of the Hamiltonian
which corresponds to the total number of particles in the zero mode, we can rearrange the Fock basis into several blocks. Figure 3 shows the wavefunction of the ground state and the first five excited states of the system for N = 120, L = 0, αN = 0.7. The vertical grid lines indicate the borders between the different blocks of the Fock basis which are arranged in ascending values of n 0 . Within one block the states are further sorted with respect to n imb ≡ n 0,+ − n 0,− which characterizes the imbalance between the occupation of the zero modes of a Fock state. 1 In general, one does have [P,L] = 0, however, for L = 0 the two operators commute. Inspection of the wavefunctions in Figure 3 indicates a further substructure: Within each n 0 -subspace the wavefunction has a form corresponding to the ground (see panels (a)-(d) and (f)) or excited (see (e)) state of a particle in a box whereas over the whole Fock space these fine structures are enveloped by an overall oscillation. Going even further, this kind of behaviour can be compared to the excitation spectrum of a molecule in the Born-Oppenheimer approximation [31]. In this picture, the behaviour within a constant n 0 -subspace corresponds to a fast degree of freedom which separates the energy spectrum into different channels [32]. Within each channel, there are smaller excitations which are determined by the slow degree of freedom corresponding to the behaviour of the oscillations in the envelope.
Based on this physical motivation we now take a look at the matrix representation of the Hamiltonian (8). If we choose N > 0 and order the Fock basis in blocks of constant n 0 including both parities P = ±1, one obtains a tridiagonal block matrix where H n 0 is the projection of the Hamiltonian (8) into the subspace with fixed n 0 , while H n 0 ±2,n 0 couples the n 0 -block to its next neighbours. Due to the form of the interaction all other blocks vanish. The next step is to define transformations U n 0 which diagonalize H n 0 and thereby the global transformation from the Fock space into a basis that diagonalizes each projection of the Hamiltonian (8) to an n 0 -subspace. This allows us to systematically select vectors solely corresponding to the ground, first or second excited states of the channels and project out all the others. This projection then neglects all possible couplings between different channels. Note that this procedure has to be repeated for every αN as the magnitude of the interaction alters the corresponding eigenvectors of the n 0 -blocks. The resulting spectrum is shown in Figure 4. Neglecting the coupling of different channels is fully justified as seen from the excellent agreement between the exact and approximated spectrum. Here, the solid blue lines show the energy levels of the full Hamiltonian (8). In comparison, the dotted lines show the excitation spectrum if the Hamiltonian is restricted to different single channels. They correspond to restrictions to the ground state (red), first (black) and second (green) excited state within the fast degree of freedom. The excellent agreement shows that this approximation provides energy levels of the original system quantitatively to very good accuracy. Furthermore, it enables us to split the spectrum into several subspectra which can be investigated independently of each other.

Effective Hamiltonian
The subsequent derivation is carried out for choosen particle number N and total momentum L, such that the respective operators are replaced by these quantum numbers. To properly derive the block structure of the Hamiltonian (8), an operatorP n 0 can be defined that projects the Hilbert space onto its subspace with constant n 0 . By definition we have which can be used to rewrite the Hamiltonian aŝ where the second term is the coupling Hamiltonian. The part diagonal in n 0 is the effective Hamiltonian that defines the adiabatically separated channels and is the starting point of all further analysis.

Redefinition of zero modes
The effective Hamiltonian (25) has an additional constant of motion, that can be easily shown to commute also withN andL. The redefinition of the creation and annihilation operators of the zero modesẑ gives whilen 0 =ẑ † +ẑ + +ẑ † −ẑ − keeps its structure. A further definition offers a new good quantum number necessary to describe the effective system: Therefore we are able to rewriteĤ 1 (n 0 ) in diagonal form aŝ where the range of this new quantum quantum number c 0 ∈ {0, 1..., n 0 } depends on n 0 . Note that the operatorĉ 0 is deliberately choosen in a way such that the resulting eigenenergy ofĤ 1 (n 0 ) is minimal for c 0 = 0.

Redefinition of kinetic modes
Now we focus on the remaining parts of the effective Hamiltonian (25) to show how it can be rendered diagonal by a redefinition of the creation and annihilation operators. Up to nowĤ eff (n 0 ) (25) consists of two parts. While the first one (E 1 (n 0 , c 0 )) has been analyzed in subsection 4.1, the second part looks comparatively difficult: whereĥ ρ is defined asĥ It is quadratic in the creation and annihilation operatorŝ suggesting to define a vector containing all annihilation operators of the (k = ±1)-modes, that allows us to rewrite the Hamiltonian (32) asĤ where and K > 0 depends on N and n 0 via The quadratic form (36) allows us to diagonalize the Hamiltonian (32). From the blockstructure of the matrix one can already conclude that the diagonalization will only mix those operators within the same k-mode, where C ± are matrices obtained from the eigenvectors of A ± . This notation is chosen in such a way that "p" corresponds to the new operators obtained from the operators acting on "positive" k-modes and "n" from the "negative" ones. Furthermore the "+", "−" indices (not to be confused with the eigenvalues of the parity operator) of the new operators refer to the associated eigenvalues of the diagonalized matrix As this redefinition is a rotation of the old operators, the sum of their occupation numbers remains unaffected,n and the same holds true for the negative k-modeŝ Finally, in view of the transformation from subsection 4.1, we are able to fully diagonalize the effective Hamiltonian (25) This can be made explicit using the eigenbasis of the operatorŝ that commute withĤ eff (n 0 ). Using one gets the explicit expression for the eigenenergies. Note that the range of the new quantum numbers c ± ∈ 0, 1, ..., is defined by N, L and n 0 , while, in the case of L = 0, (46) simplifies to Each combination of quantum numbers (c 0 , c + , c − ) then defines a different channel within the effective Hamiltonian (43). In a last step, we assume that interactions between different channels can be neglected as motivated in section 3. Within an (c 0 , c + , c − )-channel this leaves only one possible and its Hermitian conjugate, leading to an approximated Hamiltonian with a ≡ α 4 (3N + n 0 − 2).
(50) Figure 5 presents the decoupled energy spectrum of this system resulting from the Hamiltonian (50). The contribution of c + and c − to the approximate energy E eff (n 0 , c 0 , c + , c − ) depends only on their sum c + + c − for the case L = 0, and therefore c − was chosen to be always zero. The resulting spectrum is plotted (marked by dots) against the one (solid lines) from the complete Hamiltonian (8). While the higher excitations show small deviations, the results are essentially the same as without the approximation. In particular, as the main interest is in the lowest channel which corresponds to the black dots, the new quantum numbers give rise to the ability to split the spectrum into several combinations of {c 0 , c + , c − }. Further investigations will be focused on the ground state and the lowest excitations which means that these quantum numbers are always chosen to be zero.

Classical analysis
In the following, we will analyse the critical properties of our effective model (50) by means of a semiclassical analysis. Starting with the diagonal part (43) we substitute the creation/annhihlation operators by classical phase space variableŝ and neglect all terms of order O(N 0 ) in the limit of N → ∞, to obtain where "cl" refers to the classical (mean field) limit. Since the coupling between different channels can be neglected, as shown in section 3 and section 4, the classical form of the remaining interaction then gives H coup,cl (n 0 ) = α 2 (1 + tanh(γ)) n z,− √ n p,− n n,− · cos(φ z,− − φ p,− − φ n,− ).
To get an easily solvable form we reduce the Hamiltonian (51) to its channel of minimal energy by setting n z,+ = n p,+ = n n, while we reexpress {n z,− , n p,− , n n,− } in terms of N, L and n 0 through the point transformation This finally leads to a one-dimensional description with only two (conjugate) phase-space coordinates n 0 and θ, To extract the physical properties of this mean field Hamiltonian (57), valid for lim N → ∞, we define scaled variables to get the energy per particle as We are now ready to proceed with the study of the classical phase space. Obviously, it is π-periodic in θ such that the analysis can be restricted to θ ∈ [− π 2 , π 2 ]. Figure 6 shows contour plots of the energy e cl for different values of the couplingᾱ. As clearly seen, there is a qualitative change within the phase space, when the scaled interaction is increased fromᾱ = 0 andᾱ = 1. While in the non-interacting case, the phase space allows only rotations 2 , the phase space is divided into two qualitatively different regions atᾱ = 1. The regime of the lowest energies consists of vibrations/librations, separated from the rotating orbits by a separatrix. This separatrix is created at z = φ = 0 at a critical interaction α crit . Furthermore, for weak interaction (0 ≤ᾱ ≤ᾱ crit ) the energy minimum is located at z = 0 and degenerate in θ. In contrast, at a stronger interaction (ᾱ crit <ᾱ), the energy minimum consists of only one discrete point z > 0, θ = 0.
According to its definition, z represents the ratio of particles within the zero modes n 0+ and n 0− with respect to the whole particle number N. This yields the interpretation that, for an interaction greater than α crit , the occupation of the zero modes within the ground state changes from a microscopic occupation near zero to a macroscopic one at a finite value and therefore indicates that z can be taken as an order parameter characterizing a quantum phase transition. This is in complete analogy with the spin-one Bose gas without pseudospin and quadratic Zeeman shift [33] and the truncated versions of the attractive one-dimensional Bose gas [21].

Analytic analysis of the quantum phase transition
Armed with a clear signature of a phase transition in the change of morphology of the classical (mean field) limit produced by the appearance of the separatrix, we will now study the different aspects of this 2 Using the analogy to the mathematica pendulum. critical behaviour. As discussed before in section 5, the energy minimum is always located at θ = 0 for α >ᾱ crit and is degenerate in θ forᾱ ≤ᾱ crit . Therefore, this variable can be eliminated in the following discussion by setting θ = 0. The resulting energy dependence e cl (ᾱ, θ = 0, z) on z for severalᾱ is shown in Figure 7.
The range of z was deliberately chosen as {− 1 2 , 1.}, despite the fact that negative z are unphysical according to its definition, to illustrate the behaviour of the local minimum depending on the interaction strengthᾱ. Forᾱ ≤ᾱ crit this minimum would be at z * < 0. As this is not part of the allowed phase space, the minimum will simply be located at z * = 0. If the interaction strength is increased, z * increases too until it reaches z * = 0. This is exactly the point where the quantum phase transition can be expected. To find the critical value, one cane use that the derivative of the energy e cl with respect to z should vanish when evaluated at z = 0 andᾱ =ᾱ crit which provides the critical parameter as To further prove this critical behaviour, Figure 8 shows the functional dependence of the second derivative of energy minimum with respect toᾱ, The plot consists of four curves and a dashed line indicating the exact N → ∞ values of the discontinuity. All the curves are based on values of the groundstate energy for discrete sets of points ofᾱ, with the second derivative evaluated numerically. The blue dots were calculated using the energy dependence given by the classical Hamiltonian (59), whose minimal energy was numerically determined within the phase space for different values ofᾱ. They are compared to the quantum mechanical results for the ground state at various particle numbers N given by the lowest eigenvalue of the matrix representation of (50) renormalized by 1 N . The analytical result e 0,< (ᾱ), is obtained through a simple derivative of the classical energy with respect to z at the critical point e 0,< (ᾱ crit ) = ∂ 2 e cl (ᾱ crit , θ = 0, z) Extracting the second value right behind the critical threshhold is a bit harder, as the change of the z-position depending onᾱ has to be take into account. To this end a leading-order expansion in z is necessary where we used z(ᾱ crit ) = 0. Now problem is reduced to calculating the derivative of z with respect toᾱ at the critical point. For this purpose we define the function The zero of this function for a choosenᾱ gives the z-position of the energy minimum and therefore its derivative is The last step is to insertᾱ into and to calculate the second derivative Clearly, the dependence of the ground state energy is seen to be discontinuous atᾱ =ᾱ crit with α crit determined in the previous section. With the last results we even obtained an analytic expression to quantify the magnitude of the discontinuity e 0,< − e 0,> = 81 289 569 in excellent agreement with the numerical result shown in Figure 8.

Further characterization of the critical behaviour
In this last section, we will further characterize the finite-size effects in the quantum phase transition by means of the way the critical parameters approach their sharp values in the mean field limit N → ∞. Our choice of the appropriate observables comes from the behaviour of the spectrum when we approach the critical region. As seen in Figure 10, and in accordance with what happens in the attractive Lieb-Liniger model [21], one observes a strong accumulation of excited states around criticality, a phenomenon that can be related to an excited-state quantum phase transition [34].
The structure of the spectrum in Figure 9 and the dependence shown on Figure 10 suggest that the approach to criticality is well captured by two parameters, namely the minimal gap and interaction value describing its position, in the form of a power laws where β, γ > 0 [35]. To this end the gap is numerically calculated in a small region between specifically chosenᾱ min ,ᾱ max for a given particle number N. Afterwards, an interpolation function is calculated within this region and the minimum of it is numerically determined. This procedure is repeated for several N. Because of its special behaviour at the phase transition, the necessary numerical effort can be reduced drastically [36]. By means of this numerical approach, we are able to present results with particle numbers between twenty and five million. The results are shown in Figure 11 using a double logarithmic scale.  Figure 11. Asymptotic behaviour of the gap in interactionᾱ (a) and energy (b) depending on the particle number N in a double logarithmic plot with linear fits.
To extract the power law the particle numbers with N ≥ 20 000 are fitted linearly. Smaller particle numbers are taken out of the fit because this power-law is found to be valid only for large particle numbers. The obtained relations are where the powers seem to coincide with the values − 1 3 and − 2 3 within small tolerance. Therefore, apart from the quantum phase transition defined in the limit N → ∞, also finite N effects occurring in the regime of large N are properly accounted for in our analytical approach.

Summary and conclusion
In this article, we have explored a (pseudo) relativistic extension of the attractive Lieb-Liniger model, by considering both particles with linear dispersion and spin degree of freedom. Our objective was to check the existence of a relativistic analogue of the well-known quantum phase transition [26] displayed by the original non-relativistic model, where the attractive potential drives a transition of the ground state from a homogeneous state into an inhomogeneous one due to the critical appearance of a bright soliton, as thoroughly study by means of semiclassical methods in [21].
As a main result we find numerically and explain analytically that the relativistic extension indeed shows clear signatures of critical behaviour and a quantum phase transition where the macroscopic occupation of the side modes (|k| = 1), characterized by the vanishing order parameter given by the occupation of the homogeneous zero modes, is destroyed by quantum fluctuations giving rise to macroscopic occupation of the zero modes, indicating a sudden broadening of the particle distribution and an increase in the interaction energy.
Given the fact that the existence of the phase transition in the non-relativistic case is essentially due to the quantum integrability of the model, the fact that the same effect can be seen in the present non-integrable system points towards universal aspects of this transition.
In order to get an analytical understanding of this transition and its connection to the integrability of the non-relativistic case, we followed a combined approach. First, extensive numerical simulations show an adiabatic separation that mimics integrability in the low-energy region. Second, a classical analysis based on this approximate separability of the model allows for understanding the critical behaviour as a consequence of the appearance of separatrix motion in the mean field limit. This combination enabled us to provide analytical results for the location and characteristics of the quantum phase transition in excellent agreement with exact diagonalization results.