Continuous matrix product states solution for the mixing/demixing transition in one-dimensional quantum fields

We solve the mixing-demixing transition in repulsive one-dimensional bose-bose mixtures. This is done numerically by means of the continuous matrix product states variational ansatz. We show that the effective low-energy bosonization theory is able to detect the transition whenever the Luttinger parameters are exactly computed. We further characterize the transition by calculating the ground-state energy density, the field-field fluctuations and the density correlations.


I. INTRODUCTION
Understanding the physics of strongly correlated many-body systems is a formidable task, both in the lattice and in the continuum [1]. There is a fruitful synergy between condensed matter, high energy physics or quantum chemistry and the quantum information community. Ideas as tensor network states [2] or quantum simulations [3] are pursuing the goal of understanding phases and dynamics beyond the paradigm of perturbative theories.
One-dimensional many-body systems are a good example of this cooperation. Well established theoretical techniques as bosonization [4] are complemented with the density matrix renormalization group (DMRG) and matrix product states (MPS) [5,6] in the lattice and more recently, continuous matrix product states (cMPS) in the continuum [7]. Ultracold gases are a paradigmatic example of experiments realizing one dimensional quantum fields [8,9,19]. Experiments and simulations in one dimension are perfect testbeds since each of them can be used for benchmarking the other [10][11][12][13].
Of special relevance for this work is the cMPS formalism. Introduced by Verstraete and Cirac [7], these states constitute a variational class for the efficient simulation of quantum field theories that does not rely on a space discretization. The ansatz has proven to be efficient for computing the ground state, dispersion relation [14] and quantum evolution [15] of nonrelativistic theories. In addition, introducing a suitable regularization prescription, it has also been applied to the study of certain relativistic phenomena [16]. Finally, the cMPS ansatz has been already tested for bose-bose [17] and fermi-fermi [18] mixtures, both in nonrelativistic setups.
Among the different scenarios covered in experiments, we are interested here in bose-bose mixtures [20,21] . In this work, we aim to characterize the mixing/demixing phase transition occurring in repulsive bosonic mixtures, via cMPS. Here, the competition between the intra and interspecies couplings leads to the formation of two different phases. The miscible, were the two gases coexist, and the immiscible, where the two gases separate from each other [22]. This transition has been studied analytically within different approaches [23][24][25] but its numerical simulation has been elusive.
We will review the phase transition, clarify some issues on the instability occurring already in the Luttinger liquid description, and fully characterize it via the ground state energy, fluctuations and correlation functions.

II. THE PHASE TRANSITION
Two 1D bosonic gases interacting via a quartic contact potential ( = 2m = 1) are described in second quantization via the Hamiltonian whereψ † α (x) (ψ α (x)) are the bosonic field operators which create (annihilate) bosonic particles of species α at the position x ∈ [0, L]. They satisfy the commutation relations: In this work, we want to characterize numerically the mixing-demixing transition occurring in mixtures whenever two different repulsive bosonic species are trapped together. We will consider the case where the participating species are nonconvertible, i.e. the individual particle densities ρ 0α of each bosonic species are conserved separately [ρ 0α = ψ † α (x)ψ α (x) , means averages over the ground state of (1)]. We will also restrict to the symmetric case ρ 01 = ρ 02 = ρ 0 , g 11 = g 22 = c > 0 and g 12 = g 21 = g/2 > 0.
The mixing/demixing transition has been broadly studied analytically, see e.g. Refs. 23-28. The phase separation, which lies on the competition between the repulsion strengths c and g, can be understood on several grounds. The simplest approach considers a mean-field arXiv:1507.03613v1 [quant-ph] 13 Jul 2015 treatment. Here, the interaction term can be seen as a quadratic form of the densities. The latter is positive defined as long as g < 2c. When the positivity condition is violated (g ≥ 2c) an instability occurs. In more than one dimension, both species must separate in order to make the overlap integral zero, i.e., minimizing the repulsive interaction which yields the instability. In this phase (g ≥ 2c) both species are immiscible as observed by resolving the spatial density profiles in the trap [22]. In 1D, there is no possibility of spatial separation. The in- (1) are zero in the ground state.
One-dimensional systems are somehow special. Their confinement provides an enhancement of the collective behaviour, leading to a universality class in the lowenergy or long wavelength sector. The latter is known as Luttinger Liquid [4,29,30]. This regime is described by introducing the bosonic operatorsφ α andθ α in terms of which we rewrite the field operatorsψ α ( . This is nothing but the harmonic fluid approach treatment best known in the literature as bosonization [31]. Note that for high enough values of p, the exponential terms oscillate very fast and rapidly average to zero. Therefore, in order to obtain the low-energy effective Hamiltonian, we should only keep a few relevant terms. This leads to This long wavelength description is fully characterized by the dimensionless parameters K α , the velocities v α and the coupling strengths g x and g c (Luttinger parameters). For the symmetric case considered here, we have that v 1 = v 2 = v and K 1 = K 2 = K. This model can be easily decoupled by introducing the normal modesφ ± = 1/ √ 2(φ 1 ±φ 2 ) andθ ± = 1/ √ 2(θ 1 ±θ 2 ). In terms of them, the low-energy Hamiltonian reads As pointed out by Cazalilla and Ho in Ref. 24, the coupled system (2) is unstable when v 2 − becomes negative. In other words, the action is not anymore definite positive, pretty much like in the mean-field argument sketched before. This will happen whenever Kg x > v. Thus, to compute the transition point, we just need to find the Luttinger parameters from the original Hamiltonian (1). In Ref. 24, K, g x and v were approximated via expressions valid in the weak interspecies coupling g regime. In the quasi-condensate regime γ = c/ρ 1 [31], the instability is estimated to happen at g * = 2c(1 − √ γ/2π). We stress that this result deviates from the mean field value g * = 2c.
The phase separation has also been studied analytically beyond perturbation theory by Kolezhuk [25]. He found that for one and two-dimensional gases, the transition point, in the symmetric case, does not depend on the particle densities. Surprisingly enough, the nonperturbative result coincides with the mean-field description. That is, the two species demix when g ≥ g * = 2c. Following this result, one might be tempted to think that the bosonization framework is not able to predict correctly the transition point. It could be argued that the Luttinger liquid paradigm breaks down at intermediates values of g, below the critical value g * = 2c. Here, we will show that this is not the case. We demonstrate that the bosonization predicts the transition correctly when the Luttinger parameters are computed exactly instead of using approximations.

III. CMPS SOLUTION
A translational invariant cMPS of N species (N = 2 in this work) of bosonic particles is defined by the state vector [16]: whereψ α (x) are the bosonic field operators, Q and R α are a set of complex, D × D matrices acting on an auxiliary D-dimensional space and |Ω is the free vacuum state vector (ψ α (x)|Ω = 0). P denotes a path-ordering prescription and the partial trace, Tr aux , is taken over the auxiliary space. This way of writing field states is the continuous limit of a MPS [7]. Here, the dimension D of the auxiliary matrices corresponds to the so-called bond dimension, an upper bound to the entanglement entropy. Typically, the low-energy states of local Hamiltonians should possess a low amount of entanglement, consequently D is a small number. Being the bond dimension small, the state (4) represents an efficient trial for finding the ground state of one-dimensional field theories numerically.
In a previous work [17] the authors showed how to construct a two-species cMPS starting from two decoupled single species solutions. In brief, for coupled fields we considered coupled auxiliary spaces (one per bosonic field). The total auxiliary Hamiltonian was extended to where K α is the auxiliary Hamiltonian associated to bosonic species α. The parameter P accounts for the number of pairs of coupling matrices entering in the cMPS state. Consequently, the matrices R α , belonging to the auxiliary space of field α, were extended into the total product space: Denoting D the dimension of the matrices R 1 and R 2 the bond dimension is then D = D 2 . The total number of variatonal parameters is D 2 (4 + 2P ). Details can be found in Ref. 17.
In the thermodynamic limit (L → ∞), the fluctuations and correlation functions can be computed from without loss of generality, we have assumed that x > y.
Remind that, throughout this work means averages over the ground state of (1). The transfer operator T is defined as: Finally, note that the fluctuations are calculated by making x = y in (6).

IV. RESULTS
As already anticipated, our goal is to characterize the mixing/demixing transition numerically. We do it in two ways. First, we study the instability in the low-energy regime described by the effective Hamiltonian (2). The second strategy is to look directly at the ground state of (1) and compute the fluctuations and correlation functions, Cf. Eq. (6).

A. Bosonization instability
In the harmonic fluid approach the normal modes for the fields decouple (see the discussion below Eq. (2)). Each of these modes propagate with different velocities, v ± . Within the bosonization framework, these velocities can be related to the ground state energy density (e 0 ). The explicit expressions for the velocities are [31,32] with ρ ± = ρ 1 ± ρ 2 . Analytical estimations for these velocities follow from Eq. (3). In the weak-coupling regime (g c), it is safe to assume that v and K correspond to the solutions for a single bosonic field [31]. In turn, g x is approximated by g x g/π (see reference 24).
In the inset of figure 1, it can be seen that already at intermediate values of g (well below the critical value g * ) the predicted velocities v ± using weak-coupling analytical expressions deviate from the numerically computed ones [17,32]. A consequence of this deviation is the failure on the estimation of the point where v 2 − becomes negative, which in turns marks the critical value g * . In the main figure 1 we have zoomed the v 2 − around the transition point for different values of γ = c/ρ. As it has been already pointed out, within the weakcoupling treatment, the transition is estimated to happen at g * = 2c(1 − √ γ/2π). As γ = c/ρ, the latter result makes the transition point dependent on both, the intraspecies coupling c and the particle density ρ. On the other hand, once v − is exactly derived from the groundstate energy density by using relation (7), we see how the transition point becomes independent of γ. In fact, we see that the mode propagating with velocity v − becomes ill-defined at g * /c = 2 [33], in agreement with the meanfield and Kolezhuk results [25]. Therefore, once the Luttinger parameters are exactly computed, the bosonization predicts correctly the transition.

B. Characterization beyond bosonization
Having a full knowledge of the ground state of (1), we proceed now to characterize the phase transition beyond the bosonization formalism. In figure 2, we see the behaviour of the ground state energy density as a function of the interspecies coupling. It is direct to realize that after g * /c = 2 the energy remains constant. In this region, the ground state is such that the last term of (1), i.e., the one accounting for the interaction among different fields, has zero average. In other words, after the transition we have that: C 1,2 (0) = ψ † 1 (x)ψ 1 (x)ψ † 2 (x)ψ 2 (x) = 0, which is explicitly represented in the inset of figure 2 (open squares). This confirms our previous exposition for the phase transition: in one dimension, phase separation implies zero interspecies fluctuations.
Apart from the transition point estimation and the zero field-field overlapping nature for the demixed phase, we can go further in characterizing the properties of the ground state before and after the transition. Let us start with the mixed phase. By looking at the inset of figure 2 we see that the fluctuations C 1,2 (0) do not remain constant as soon as the interaction is switched on. The latter behaviour reflects a sublinear growth of e 0 in terms of g. This means that a simple mean field theory ψ † 1 (x)ψ 1 (x)ψ † 2 (x)ψ 2 (x) ∼ = ρ 1 ρ 2 is not sufficient for describing this phase.
We will discuss now the demixed phase. As explained above, after the transition C 1,2 (0) = 0. It is straightforward to see that a ground state of the form fulfils this condition (suffix dm stands for demixing). Besides, |X dm must satisfy the particle density conservation for each bosonic species: X dm |ψ † α (x)ψ α (x)|X dm = ρ, which in turn imposes that: χ 2ρ |ψ † α (x)ψ α (x)|χ 2ρ = 2ρ. Indeed, this is confirmed in Fig. 2 via the ground state energy density. We see that after the transition, e 0 is the energy of a single bosonic gas (Lieb-Liniger model) with self-interaction c but double particle density 2ρ (dashed line) [33]. Finally, by looking at the fluctuations C αα (0), we check that they coincide with those of a single bosonic gas with self-interaction c and particle density 2ρ, divided by a factor of two due to normalization in (8). The fluctuations of a single gas are shown in the inset of figure 2 with a dashed line.
We finish our phase characterization by studying the correlation functions, C αβ (x). The results are plotted in figure 3. By definition, the correlations at zero distance match the fluctuations. On the other hand, in the limit x → ∞, the correlations factorize yielding In the mixed phase the correlation length is of the order of x ∼ = 5, pretty much the same than for a single bosonic species with self-interaction c and particle density ρ.
More structure for C α,β (x) appears in the demixed phase. The interspecies correlation function C 12 (x), obviously starting at zero, has a large correlation length ∼ 10 4 (notice the logarithmic scale). To understand this large correlation length we recall that after the transition the fields are infinitely repelled. Our interpretation is reinforced by looking at C αα (x). In the range 0 < x < 10 the correlations build up to 2ρ 2 wich means that they can be approximated by C αα (x) ∼ = 1/2 χ 2ρ |ρ α |χ 2ρ 2 = 2ρ 2 . Therefore, the coherence has been lost at the single field level. However, the fully uncorrelated state will involve the full state |X dm and pretty much like for the C 12 (x) correlations, the demixed phase is equivalent to an infinite repulsive phase, explaining again the large coherence lenght to reach the asymptotic limit C αα (x → ∞) = ρ 2 .

V. CONCLUSIONS
Summarizing, by means of cMPS we have computed numerically the ground state of two repulsive onedimensional bosonic nonconvertible fields. This kind of systems exhibits the so-called mixing/demixing phase transition. We have validated previous analytical results for the transition point. Furthermore, we have demonstrated that this point can be resolved within the Luttinger liquid formalism whenever the effective parameters of the theory are calculated exactly. All this marks a step forward for the cMPS method, here, resolving a phase transition in a non-trivial quantum field theory.