Nonlocal order parameter of pair superfluids
Nitya Cuzzuol Luca Barbiero Arianna Montorsi


Nonlocal order parameter of pair superfluids


^1Institute for Condensed Matter Physics and Complex Systems, DISAT, Politecnico di Torino, I-10129 Torino, Italy
nitya.cuzzuol@polito.it, arianna.montorsi@polito.it








Order parameters represent a fundamental resource to characterize quantum matter. We show that pair superfluids can be rigorously defined in terms of a nonlocal order parameter, named odd parity, which derivation is experimentally accessible by local density measurements. As a case of study, we first investigate a constrained Bose-Hubbard model at different densities, both in one and two spatial dimensions. Here, our analysis finds pair superfluidity for relatively strong attractive interactions. The odd parity operator acts as the unique order parameter for such phase irrespectively to the density of the system and its dimensionality. In order to enforce our finding, we confirm the generality of our approach also on a two-component Bose-Hubbard Hamiltonian, which experimental realization represents a timely topic in ultracold atomic systems. Our results shed new light on the role of correlated density fluctuations in pair superfluids. In addition, they provide a powerful tool for the experimental detection of such exotic phases and the characterization of their transition to the normal superfluid phase.



















§ INTRODUCTION



Quantumness is fundamentally rooted in delocalization and entanglement, the latter allowing quantum matter to organize in low temperature phases which escape Landau classification of spontaneously symmetry breaking (SSB) phases <cit.>. Notable examples are topological phases <cit.>, as well as insulators with continuum symmetries, i.e. low dimensional Mott insulators <cit.>. Crucially, while SSB phases are endowed with a local order, the aforementioned examples are characterized by nonlocal orderings which do not violate the Mermin-Wagner theorem <cit.>. Such intrinsic nonlocality allows for insulators with continuous symmetries both in 1D and 2D to be accurately captured by nonlocal parity operators <cit.>, while nonlocal string order parameters <cit.> result the fundamental quantity to capture one dimensional topological insulators. Relevantly, from SSB <cit.> to low dimensional Mott <cit.> and topological <cit.> insulators, experimental setups working with ultracold atoms <cit.> have allowed for an in depth description of such states of quantum matter. Besides the impressive versatility encoded in such quantum simulators, a game changing aspect is represented by highly accurate probing schemes. By indeed making use of either noise correlators <cit.> or quantum gas microscopy <cit.> both local and nonlocal orderings have been efficiently probed.  

Although the above discussion focused on insulating states of matter, the notion of nonlocal orders can be applied also to characterize conducting phases of low-dimensional interacting fermions <cit.>. In this regard, nonlocal parity <cit.> and string <cit.> operators relative to the spin degrees of freedom have been efficiently employed to characterize Luther Emery liquid <cit.> and superconducting phases <cit.> as well as topological conducting regimes <cit.> respectively. At the same time, the order parameters of gapped conducting phases in low dimensional bosonic matter are still unexplored. In this paper we tackle this important topic.   

In particular, we explore conducting phases characterized by bosonic pairing, i. e. pair superfluids (PSFs) whose existence has been predicted in large variety of physical systems <cit.>. Importantly, PSFs are neither associated to a symmetry breaking nor to the appearance of topological order, thus representing an example of a gapped superfluid with continuous symmetry. To perform our analysis, we consider a paradigmatic model of strongly interacting bosons, namely the Bose-Hubbard model with a constraint of a maximum of two bosons per site. It was initially introduced to describe regimes characterized by a large rate of three-body losses in ultracold atomics gases <cit.>.
Previous theoretical works <cit.> have established that for weakly attractive interaction a normal superfluid (SF) phase occurs while, by increasing the attraction between bosons, the system undergoes a phase transition. Here, a different SF phase, i.e. the PSF where pairs of single bosons become the effective elementary constituents, appears. As we will discuss in details, these different superfluids have been characterized by analyzing the different long distance decay of the single particle and pair-pair correlation function <cit.>. Although certainly useful, this approach poses significant challenges. On one hand it makes difficult to properly locate the transition point. On the other hand, experimental schemes to probe the long distance behavior of correlators are still characterized by a limited efficiency. Here, we propose an alternative approach based on local density measurements, which can be performed with high accuracy through quantum gas microscopy.

Specifically, our method introduces a new nonlocal parameter, named odd parity. The prediction is first confirmed by a bosonization  and renormalization group analysis. We then prove also numerically that the odd parity is the unique order parameter of the PSF phase. As PSF can occur also in 2D <cit.>, we extend the definition of odd parity by means of the brane generalization <cit.> and we verify numerically that our approach  holds also beyond the one dimensional limit. In addition, since sample imperfections might make it challenging to accurately control the particle density, we show that our method remains highly robust in a large range of bosonic fillings. Finally, we test the validity of our predictions on an experimentally feasible two component Bose-Hubbard Hamiltonian where the SF-PSF phase transition was recently predicted <cit.>.



§ THE MODEL


We analyse the Bose-Hubbard model corresponding to the following Hamiltonian

    H = - t ∑_⟨ R, R' ⟩ (b_R^† b_R' +  h.c.)
        + U/2∑_R n_R (n_R - 1) + μ∑_R n_R

where R=(x,y) is a generic site of a lattice composed by L rungs and M legs (1 ≤ x ≤ L; 1 ≤ y ≤ M, M=1 in 1D) schematized in figure (<ref>)(a,b). Here, ⟨ R, R' ⟩ indicates nearest-neighbor sites and b_R, b_R^† are the bosonic creation and annihilation operators on the site R. In addition, μ is the system chemical potential and the hopping matrix element t=1 fixes the energy scale. In order to favor the bosonic pairing we focus on regimes with on-site attractive interaction U<0 and we impose the three body constraint (b_R^†)^3=0 which allows to prevent the system collapse. As shown both in the 1D and 2D case in figure (<ref>)(c,d), bosons remain in the SF phase up to a critical value of the attractive interaction, where PSF emerges. To locate the SF-PSF transition one can investigate the long distance behavior of the single particle and pair-pair correlation functions defined as

    S(r) = ⟨ b_x, y^† b_x+r,y⟩


    D(r) = ⟨ (b^†)^2_x,y (b)^2_x+r,y⟩,

respectively [In (<ref>) and (<ref>), the y-th leg in the lattice is constant but can assume any value. However, since the model is isotropic in both hopping and interaction, any y-th leg leads to the same effects. In one dimension, obviously y=1.]. In fact, one expects that in the SF phase both correlators exhibit algebraic decay, signaling the presence of quasi-long-range order. In contrast, in the PSF phase, one can argue that while S(r) exhibits exponential decay, indicating the suppression of single bosons motion, D(r) should maintain its polynomial decay <cit.>.

These intuitive aspects reflect the specific behavior of two different energy gaps defined as:

    Δ_α = 1,2(L× M) = E(N+α; L× M) + E(N-α;L× M) - 2E(N;L× M)

where E(N;L× M) is the ground state energy of a system of N particles on L× M sites. Specifically, one can find that, for a finite size system, Δ_2 > Δ_1 in SF while Δ_2 < Δ_1 in PSF <cit.>. On the other hand, in the termodynamic limit one finds that Δ_2=0 for any value of U, while Δ_1>0 when entering the PSF regime and thus explaining the exponential decay of S(r).



 §.§ Bosonization analysis
 
The 3-body constraint Bose-Hubbard model at filling one is directly connected to spin-1 XXZ chain with single ion anisotropy <cit.>. This aspect is at the core of the analytical treatment of (<ref>) in 1D by means of a bosonization approach <cit.>. Specifically, by following <cit.>, the Hamiltonian (<ref>) is first mapped into a spin-1 model, then each spin-1 variable is written as a sum of two spin-1/2, which in turn are mapped onto two spinless fermions via a Jordan Wigner transformation. The two corresponding fermionic ladder operators are bosonized in the standard way, introducing two bosonic fields ϕ_1,2(x), and their dual fields θ_1,2(x) such that [ϕ_α(x),θ_β(y)]=iπδ_αβΘ(x-y), with α,β=1,2 and Θ(x-y) being the Heaviside step function. Apart from a coupling term negligible in the low energy limit, the resulting Hamiltonian decouples when re-written in terms of the symmetric an anti-symmetric combinations of ϕ_1,2(x), namely ϕ_±(x) = ϕ_1(x)±ϕ_2(x) and θ_±(x) = (θ_1(x)±θ_2(x))/2. Explicitly

    H→ H_++H_-    ,

with

    H_+ ≡u_+/2π∫ dx  [K_+ (∂_x θ_+)^2+1/K_+ (∂_x ϕ_+)^2  ]+g/(π a)^2∫ dx cos (2 ϕ_+)


    H_- ≡   u_-/2π∫ dx  [K_- (∂_x θ_-)^2+1/K_- (∂_x ϕ_-)^2  ]- g/(π a)^2∫ dx cos (2 ϕ_-)
       -t/π a∫ dx cos( 2θ_-).

Here a is the lattice spacing, g=-aU/2 is the coupling constant, and  K_±= 2(√(1± U/tπ))^-1 are the Luttinger parameters, while u_±=2t a/K_± are the velocities in the ±-channels. We recognize in H_± two sine-Gordon models where the condition U<0 implies K_+>2 and K_-<2. A renormalization group analysis<cit.>  suggests that in this regime the field ϕ_+ never pins, whereas ϕ_- (for U<-2tπ) can pin to the value 0. The corresponding phase can be identified with the PSF phase.

On general grounds, the bosonization approach can also be exploited to associate to each pinning value of the bosonic fields an  appropriate nonlocal operator whose expectation value becomes different from zero in the corresponding phase <cit.>. In particular, it was shown that the MI phase of the 1D (<ref>) can be characterized by the non vanishing of the expectation value of the parity operator <cit.>, namely O_ P(j)≡∏_x=1^j-1exp [ iπ (n_x -1)], which in the continuous limit becomes O_ P(x)  ∼cos (2 ϕ_+(x)). In the thermodynamic limit, its expectation value will remain finite when most sites are singly occupied, and holon (empty sites) doublon (doubly occupies sites) fluctuations occur in correlated pairs.

Here we observe that when instead the field ϕ_- pins to the value 0 a different parity operator O_ P^ (o) should become finite, defined as  

    O_ P^ (o)(j) = ∏_x=1^j-1exp[ iπ n_x ]    .

We name it odd parity, since its expectation value will be finite when fluctuations occur in the form of correlated pairs of single bosons. In fact, in the continuum limit it turns out to be connected to the bosonic field ϕ_-, i.e. the antisymmetric combination of fermionic fields ϕ_1,2. Explicitly

    O_ P^ (o) (x) ∼cos (2 ϕ_-(x))    ,

which expectation value is finite when ϕ_-(x)=0. Thus the latter should be finite in the whole PSF phase, and vanishing as soon as ϕ_- unpins. Figure (<ref>)(a) shows a simple cartoon of the case when this quantity is nonvanishing in 1D.

When extending this analysis to the 2D limit, upon exploiting <cit.>  
one may introduce a brane odd parity to capture the PSF phase. This is the product of M odd parities for each x-value in the lattice. Explicitly

    O_ P^ (o)(j) = ∏_x=1^j-1exp[ iπ∑_y=1^M n_x,y]    .

which is nonvanishing when in the disordered background of holons and doublons density fluctuations occur only in correlated pairs, as shown schematically in figure (<ref>)(b).



§ NUMERICAL RESULTS


We present in this section the results we obtained by performing iDMRG <cit.> simulations, which allows to accurately take into account the quantumness of a given system. We analyse first the system of a single chain (M=1) and then we extend our investigation to larger systems aiming to reach the 2D limit. Within this algorithm, we can either fix a priori the particle density, working in the canonical ensemble, or tune the chemical potential, in a grand canonical framework.

Firstly, for both the 1D and 2D cases, we present a validation at fixed filling n=1 of the effectiveness of exploiting the odd parity in estimating the SF-PSF transition, in comparison to the use of correlators S(r), D(r) as defined in (<ref>, <ref>). Then we present our resulting phase diagrams in (μ, U) plane, shown in figure (<ref>)(c,d). Next, for both the 1D and 2D cases, we further investigate the filling dependence of the SF-PSF phase transition by imposing n≠ 1 throughout the transition.



 §.§ Results in 1D





Our investigation starts with the study of local correlators in (<ref>) and (<ref>) by fixing n=1 and enforcing particle number conservation. As illustrated in figure (<ref>)(a) and (b) and mentioned before, we find that strong attractions generate an exponential decay of S(r) while D(r) maintains its quasi-long range order for any U. This behavior can be compared with that of the odd parity, as shown in figure (<ref>)(c). Here, O_ P^ (o) shows a sharp and well defined transition point. Specifically, the odd parity becomes finite at a specific U_ c which clearly does correspond to the one where the long-range behavior of S(r) starts to deviate from an algebraic decay.

Our numerical approach makes it also possible to accurately complement this analysis by incorporating a study of the gaps Δ_1,2. The thermodynamic extrapolation in figure (<ref>) shows indeed that in the SF Δ_1=Δ_2=0 while, as long as PSF takes place, Δ_2=0 but Δ_1 becomes finite. Notably, the condition Δ_1>0 turns out to be fulfilled exactly when O_ P^ (o)>0 therefore unambiguously proving the validity and accurateness of our approach.

In order to further enforce our results, restricted so far to the canonical ensemble, we also inspect the odd parity in the grand canonical ensemble by tuning both U and μ. The resulting phase diagram is reported in figure (<ref>)(c), based on the analysis of O_ P^ (o). It confirms the presence of the PSF phase (in red), in which the odd parity order remains robust irrespectively to the particle density. 



 §.§ Results in 2D


The simulation of 2D systems using MPS algorithms can be challenging due to the efficiency of this method relying on covering a finite ratio of the system's entanglement entropy. However, the entanglement entropy grows extensively according to the area law, diminishing the algorithm's efficiency <cit.>.
In our iDMRG simulations, we address this challenge by considering an infinite size in the x-direction, achieved through the repetition of a finite unit cell in the y-direction. The use of periodic boundary conditions accelerates convergence to a 2D configuration. The resulting infinite cylinder system can be effectively mapped to a 1D system by following a 'snake-like' path, initially traversing in the y-direction.

We start again from the canonical ensemble by considering the different behaviours of the correlators (<ref>) and (<ref>) at n=1. As in the previous case, figure (<ref>)(a) and (b) shows that for a sufficiently negative U there is a change in the decay law for S(r) while D(r) conserves his long distance behavior. As before, this point signals the emergence of the PSF phase. Again, we show in figure (<ref>)(c) this transition to coincide with the critical U_ c where the odd parity becomes finite.

Then, by increasing M from 3 to 8, we evaluated the odd parity to extrapolate the critical interaction for SF-PSF transition in the C_4 symmetric 2D limit [We consider the case M=6 as a reliable approximation of a 2D system, given that the difference between this case and our resulting limit is less than 3%.], as illustrated in figure (<ref>)(c).

For M=6 we further release the particle number conservation and evaluate again the odd parity by varying the onsite attractive interaction U in the (μ, U) plane. The resulting phase diagram, which should be representative of the 2D limit, is shown in figure (<ref>)(d): also in this case the presence of the PSF phase (in red) is in one-to-one correspondence with the region of finite odd parity. This is confined between the SF phase (in blue), and the two trivial phases (in white) corresponding to the n=0 and n=2 cases.



 §.§ Arbitrary Filling




The previously discussed results for 1D and 2D systems focused first on the canonical ensemble, where we enforced a filling of n=1. Subsequently, such condition was released upon exploring the behavior of the system in the grand canonical ensemble, where the number of bosons is not necessarily constant throughout the SF-PSF phase transitions, and we obtained the (μ, U) 1D and 2D phase diagrams. To confirm the picture, we now inspect the SF-PSF transition again in the canonical ensemble, for fixed fillings n≠ 1.

The results are shown in figure (<ref>) for a 1D (a) and a 2D (b) system. There the odd parity (<ref>) was evaluated for fillings greater and smaller than one. The SF and PSF phases are again distinguished through the behavior of odd parity: in the first, the odd parity is zero, whereas in the other the emergence of a nonlocal order is confirmed by its non zero value. A careful analysis reveals how at all fillings the transition is second order, in both 1D and 2D cases[This is at variance with the case of fixed chemical potential, in which traversing a transition in the (μ, U) plane may result in the total number of bosons not being conserved, leading to abrupt filling changes. This discontinuity can also manifest in odd parity, transitioning from zero to a finite value.]. Moreover, in this way the actual dependence of the critical value of interaction U_ c on filling is obtained.



 §.§ Two-species system

The presence of the SF-PSF phase transition we discussed in previous sections is not limited to the 3-body constrained BHM. For instance, it is observed also in case of two components hard-core bosons in optical lattice <cit.>. Although, the regime of PSF has not yet been observed, Bose mixtures in optical lattice have been the subject of intense experimental activity <cit.>. The Hamiltonian we are interested in reads  

    H = - t ∑_i, σ = A,B (b_i,σ^† b_i+1,σ + h.c.) + U ∑_i n_i,A n_i,B

where i is a generic site of the chain for boson species σ=  A,B.

For attractive interchain on-site interaction, the system should show a tendency to create pairs of A-B bosons.

As before we can measure the emergence of the odd parity order. The previous definition (<ref>) in this case becomes

    O_ P^ (o)(j) = ∏_i=1^j-1exp[ ∑_σ=  A,B iπ n_i,σ]   .

Now the same site on the two chains is either empty or occupied in both chains, whereas the fluctuations occur in the form of correlated pairs of bosons of the different species as depicted in figure (<ref>)(a).

In agreement with results in <cit.>, and with bosonization analysis in the fermionic case<cit.>, O_ P^ (o) becomes finite as soon as U<0, as shown in figure (<ref>)(b).



§ CONCLUSIONS

We unveiled the order parameter of the exotic pair superfluid phase. The latter takes the form of a nonlocal parity operator named odd parity which can be accurately detected through local density measurements. Such quantity is finite whenever in a disordered background of holons and doublons single bosons occur in correlated pairs. The prediction is confirmed by a bosonization analysis and extensive numerical simulations.

Relevantly, we have shown  our results to hold both in one and two spatial dimensions, at any density supporting a pair superfluid phase, and in two different models where such state of matter takes place. This aspect confirms the generality and the large range of applicability of our findings. Moreover, we proved the aforementioned points through extensive matrix-product-state simulations where quantum fluctuation are accurately taken into account. This makes our results highly robust and reliable in deep quantum regimes where mean-field approaches are not reliable.

We further stress that a major advantage of nonlocal operators is that they do not break continuous Hamiltonian symmetries, thus describing orders that can persist in low dimensional systems even at non vanishing temperatures. Hence the measurement of odd parity order is even more at hand in present experimental setups, especially because it can be  obtained by means of just local density probes. In conclusion, we believe that our results represent an important step towards a more complete and accurate characterization of strongly correlated many-body quantum systems.




§ ACKNOWLEDGMENTS

N.C. and A.M. thank for the financial support from the ICSC – Centro Nazionale di Ricerca in High Performance Computing, Big Data and Quantum Computing, funded by European Union – NextGenerationEU (Grant number CN00000013). 
L. B. acknowledges financial support within the DiQut Grant No.
2022523NA7 funded by European Union – Next
Generation EU, PRIN 2022 program (D.D. 104
- 02/02/2022 Ministero dell’Università e della
Ricerca).
Computational resources were provided by HPC@POLITO (http://www.hpc.polito.it). Calculations were performed using the TeNPy Library (version 0.10.0)<cit.>.




§ REFERENCES


iopart-num.bst

