# Nonlocal order parameter of pair superfluids

Nitya Cuzzuol1, Luca Barbiero1 and Arianna Montorsi1 1Institute for Condensed
Matter Physics and Complex Systems, DISAT, Politecnico di Torino, I-10129
Torino, Italy<EMAIL_ADDRESS><EMAIL_ADDRESS>,

###### Abstract

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.

##  1 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
[1, 2]. Notable examples are topological phases [3, 4, 5, 6, 7], as well as
insulators with continuum symmetries, i.e. low dimensional Mott insulators
[8]. 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 [9]. Such intrinsic nonlocality allows for
insulators with continuous symmetries both in 1D and 2D to be accurately
captured by nonlocal parity operators [10, 11, 12, 13, 14], while nonlocal
string order parameters [15, 16, 17] result the fundamental quantity to
capture one dimensional topological insulators. Relevantly, from SSB [18, 19,
20, 21, 22, 23, 24] to low dimensional Mott [25, 26, 27, 28] and topological
[29, 30] insulators, experimental setups working with ultracold atoms [31, 32]
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 [33, 34] or quantum gas microscopy [35]
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 [12, 36, 37, 38]. In this
regard, nonlocal parity [11, 12] and string [17] operators relative to the
spin degrees of freedom have been efficiently employed to characterize Luther
Emery liquid [11] and superconducting phases [36, 38, 37] as well as
topological conducting regimes [36, 38] 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 [39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50,
51, 52]. 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 [39, 40]. Previous theoretical
works [39, 41] 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 [39].
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 [41, 43, 44, 46], we extend the definition of odd parity by means of the
brane generalization [53, 37] 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 [54, 55, 56, 57, 58, 59].

##  2 The model

Figure 1: Cartoon of a phase with finite odd parity order (see (8) and (10))
for a 1D (a) and 2D (b) lattice respectively . It amounts to a background of
holons (empty sites) and doublons (doubly occupied sites) in which density
fluctuations (single bosons) occur only in pairs with finite correlation
length, highlighted by red shaded regions. Panels (c) and (d) present the
resulting phase diagrams obtained through iDMRG simulations for the 1D system
(with $\chi_{\rm max}=200$) and a lattice system with $M=6$ (with $\chi_{\rm
max}=1500$), respectively. The SF region (depicted in blue) corresponds to
zero expectation value of odd parity, while the PSF region (depicted in red)
is characterized by a non-zero value.

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

$H=-t\sum_{\langle
R,R^{\prime}\rangle}(b_{R}^{\dagger}b_{R^{\prime}}+\textit{\rm
h.c.})+\frac{U}{2}\sum_{R}n_{R}(n_{R}-1)+\mu\sum_{R}n_{R}$ (1)

where $R=(x,y)$ is a generic site of a lattice composed by $L$ rungs and $M$
legs ($1\leq x\leq L$; $1\leq y\leq M$, $M=1$ in 1D) schematized in figure
(1)(a,b). Here, $\langle R,R^{\prime}\rangle$ indicates nearest-neighbor sites
and $b_{R}$, $b_{R}^{\dagger}$ are the bosonic creation and annihilation
operators on the site $R$. In addition, $\mu$ 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}^{\dagger})^{3}=0$ which
allows to prevent the system collapse. As shown both in the 1D and 2D case in
figure (1)(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)=\langle b_{x,y}^{\dagger}b_{x+r,y}\rangle$ (2)
$D(r)=\langle(b^{\dagger})^{2}_{x,y}(b)^{2}_{x+r,y}\rangle,$ (3)

respectively 111In (2) and (3), 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 [39].
These intuitive aspects reflect the specific behavior of two different energy
gaps defined as:

$\Delta_{\alpha=1,2}(L\times M)=E(N+\alpha;L\times M)+E(N-\alpha;L\times
M)-2E(N;L\times M)$ (4)

where $E(N;L\times M)$ is the ground state energy of a system of $N$ particles
on $L\times M$ sites. Specifically, one can find that, for a finite size
system, $\Delta_{2}>\Delta_{1}$ in SF while $\Delta_{2}<\Delta_{1}$ in PSF
[42]. On the other hand, in the termodynamic limit one finds that
$\Delta_{2}=0$ for any value of $U$, while $\Delta_{1}>0$ when entering the
PSF regime and thus explaining the exponential decay of $S(r)$.

###  2.1 Bosonization analysis

The 3-body constraint Bose-Hubbard model at filling one is directly connected
to spin-1 XXZ chain with single ion anisotropy [60, 61, 10, 42, 62]. This
aspect is at the core of the analytical treatment of (1) in 1D by means of a
bosonization approach [63]. Specifically, by following [64, 10], the
Hamiltonian (1) 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 $\phi_{1,2}(x)$, and their dual fields
$\theta_{1,2}(x)$ such that
$[\phi_{\alpha}(x),\theta_{\beta}(y)]=i\pi\delta_{\alpha\beta}\Theta(x-y)$,
with $\alpha,\beta=1,2$ and $\Theta(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 $\phi_{1,2}(x)$, namely
$\phi_{\pm}(x)=\phi_{1}(x)\pm\phi_{2}(x)$ and
$\theta_{\pm}(x)=(\theta_{1}(x)\pm\theta_{2}(x))/2$. Explicitly

$H\rightarrow H_{+}+H_{-}\quad,$ (5)

with

$H_{+}\equiv\frac{u_{+}}{2\pi}\int
dx\left[K_{+}(\partial_{x}\theta_{+})^{2}+\frac{1}{K_{+}}(\partial_{x}\phi_{+})^{2}\right]+\frac{g}{(\pi
a)^{2}}\int dx\cos(2\phi_{+})$ (6)
$\begin{split}H_{-}\equiv&\frac{u_{-}}{2\pi}\int
dx\left[K_{-}(\partial_{x}\theta_{-})^{2}+\frac{1}{K_{-}}(\partial_{x}\phi_{-})^{2}\right]-\frac{g}{(\pi
a)^{2}}\int dx\cos(2\phi_{-})\\\ &-\frac{t}{\pi a}\int
dx\cos(2\theta_{-}).\end{split}$ (7)

Here $a$ is the lattice spacing, $g=-aU/2$ is the coupling constant, and
$K_{\pm}=2(\sqrt{1\pm U/t\pi})^{-1}$ are the Luttinger parameters, while
$u_{\pm}=2ta/K_{\pm}$ are the velocities in the $\pm$-channels. We recognize
in $H_{\pm}$ two sine-Gordon models where the condition $U<0$ implies
$K_{+}>2$ and $K_{-}<2$. A renormalization group analysis[10] suggests that in
this regime the field $\phi_{+}$ never pins, whereas $\phi_{-}$ (for
$U<-2t\pi$) 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 [11, 12, 17]. In particular, it was shown that the MI
phase of the 1D (1) can be characterized by the non vanishing of the
expectation value of the parity operator [10], namely $O_{\rm
P}(j)\equiv\prod_{x=1}^{j-1}\exp[i\pi(n_{x}-1)]$, which in the continuous
limit becomes $O_{\rm P}(x)\sim\cos(2\phi_{+}(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 $\phi_{-}$ pins to the value $0$ a
different parity operator $O_{\rm P}^{\rm(o)}$ should become finite, defined
as

$O_{\rm P}^{\rm(o)}(j)=\prod_{x=1}^{j-1}\exp\left[i\pi n_{x}\right]\quad.$ (8)

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
$\phi_{-}$, i.e. the antisymmetric combination of fermionic fields
$\phi_{1,2}$. Explicitly

$O_{\rm P}^{\rm(o)}(x)\sim\cos(2\phi_{-}(x))\quad,$ (9)

which expectation value is finite when $\phi_{-}(x)=0$. Thus the latter should
be finite in the whole PSF phase, and vanishing as soon as $\phi_{-}$ unpins.
Figure (1)(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 [13, 53, 14] 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_{\rm
P}^{\rm(o)}(j)=\prod_{x=1}^{j-1}\exp\left[i\pi\sum_{y=1}^{M}n_{x,y}\right]\quad.$
(10)

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

##  3 Numerical results

Figure 2: Analysis of a 1D system ($M=1$) at $\bar{n}=1$ by performing iDMRG
simulations with a maximum bond dimension of $\chi_{\rm max}=200$. The panels
(a) and (b) present the analysis of the correlators in (2),(3) respectively.
In the panel (c), the odd parity $O_{\rm P}^{\rm(o)}$ is evaluated.

We present in this section the results we obtained by performing iDMRG [65,
66, 67] 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 $\bar{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 (2, 3). Then we present our resulting phase
diagrams in $(\mu,U)$ plane, shown in figure (1)(c,d). Next, for both the 1D
and 2D cases, we further investigate the filling dependence of the SF-PSF
phase transition by imposing $\bar{n}\neq 1$ throughout the transition.

###  3.1 Results in 1D

Figure 3:  Analysis of the gaps defined in (4) by increasing the size of the
chain $L=60,80,120,200,300$. One can see in (a) the single particle gap while
in (b) the double particle gap. Only the first one remains finite by
extrapolating the thermodynamic limit (TDL) by considering a second order
polynomial in $1/L$. We keep as maximum bond dimension $\chi_{\rm max}=200$.

Our investigation starts with the study of local correlators in (2) and (3) by
fixing $\bar{n}=1$ and enforcing particle number conservation. As illustrated
in figure (2)(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 (2)(c). Here, $O_{\rm P}^{\rm(o)}$ shows a sharp
and well defined transition point. Specifically, the odd parity becomes finite
at a specific $U_{\rm 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 $\Delta_{1,2}$. The
thermodynamic extrapolation in figure (3) shows indeed that in the SF
$\Delta_{1}=\Delta_{2}=0$ while, as long as PSF takes place, $\Delta_{2}=0$
but $\Delta_{1}$ becomes finite. Notably, the condition $\Delta_{1}>0$ turns
out to be fulfilled exactly when $O_{\rm P}^{\rm(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 $\mu$. The resulting phase diagram is reported in figure
(1)(c), based on the analysis of $O_{\rm P}^{\rm(o)}$. It confirms the
presence of the PSF phase (in red), in which the odd parity order remains
robust irrespectively to the particle density.

###  3.2 Results in 2D

Figure 4: Analysis for a 2D system. The first two panels (a) and (b) present
the the correlators in (2),(3) respectively, by increasing the onsite
attractive interaction $U<0$ for $M=6$. Panel (c) presents the plot of the SF-
PSF transition for 2D systems as a function of the number of legs $M$. The
inset shows the extrapolation for a square 2D lattice of the critical
interaction $U_{\rm c}({\rm 2D})=-17.65\pm 0.05$ imposing $\chi_{\rm
max}=1500$. Similar to the 1D case, a comparison of the three plots reveals
that, coinciding with the change in decay observed in $S(r)$ the odd parity
becomes finite.

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
[68, 67]. 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 (2) and (3) at $\bar{n}=1$. As in the previous
case, figure (4)(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 (4)(c) this transition to coincide with the
critical $U_{\rm 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 222We 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 (4)(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
$(\mu,U)$ plane. The resulting phase diagram, which should be representative
of the 2D limit, is shown in figure (1)(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 $\bar{n}=0$ and $\bar{n}=2$
cases.

###  3.3 Arbitrary Filling

Figure 5:  Evaluation of the odd parity $O_{\rm P}^{\rm(o)}$ for a fixed
$\bar{n}\neq 1$ for the 1D system in (a) and 2D in (b), both by performing
iDMRG simulations for $\chi_{\rm max}=200,1500$ respectively. The $\bar{n}=1$
case is included for comparison with previous results. As evident from
figures, for the considered densities there is no change in the nature of the
transition when moving away from filling one.

The previously discussed results for 1D and 2D systems focused first on the
canonical ensemble, where we enforced a filling of $\bar{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
$(\mu,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
$\bar{n}\neq 1$.
The results are shown in figure (5) for a 1D (a) and a 2D (b) system. There
the odd parity (8) 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 cases333This is at variance with the case of fixed chemical
potential, in which traversing a transition in the ($\mu$, $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_{\rm c}$ on filling is
obtained.

###  3.4 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 [54, 55,
56, 57, 58, 59]. Although, the regime of PSF has not yet been observed, Bose
mixtures in optical lattice have been the subject of intense experimental
activity [69, 70, 71, 72, 73, 74, 75, 76]. The Hamiltonian we are interested
in reads

$H=-t\sum_{i,\sigma=\rm{A,B}}(b_{i,\sigma}^{\dagger}b_{i+1,\sigma}+\text{h.c.})+U\sum_{i}n_{i,\rm{A}}n_{i,\rm{B}}$
(11)

where $i$ is a generic site of the chain for boson species $\sigma=\rm A,B$.

Figure 6: (a) Schematic representation of the two hard-core bosons chains
model. One can see the direct connection to the previous representation in
figure (1)(a,b). Now the couple of bosons are composed by A (in blue) and B
(in green) bosons. (b) Evaluation of the odd parity in (12) for SF-PSF phase
transition. Differently from previous case, the transition occurs for small
negative interaction.

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 (8) in this case becomes

$O_{\rm P}^{\rm(o)}(j)=\prod_{i=1}^{j-1}\exp\left[\sum_{\sigma=\rm A,B}i\pi
n_{i,\sigma}\right]\>.$ (12)

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 (6)(a).
In agreement with results in [77, 59], and with bosonization analysis in the
fermionic case[11], $O_{\rm P}^{\rm(o)}$ becomes finite as soon as $U<0$, as
shown in figure (6)(b).

##  4 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)[66].

## References

## References

  * [1] Landau L D, Lifshitz E M and Pitaevskii M 1999 Statistical Physics (Butterworth-Heinemann, New York)
  * [2] Wilson K G and Kogut J 1974 Phys. Rep. 12 75–199 ISSN 0370-1573 URL https://www.sciencedirect.com/science/article/pii/0370157374900234
  * [3] Altland A and Zirnbauer M R 1997 Phys. Rev. B 55(2) 1142–1161 URL https://link.aps.org/doi/10.1103/PhysRevB.55.1142
  * [4] Kitaev A 2009 AIP Conference Proceedings 1134 22–30 ISSN 0094-243X (Preprint https://pubs.aip.org/aip/acp/article-pdf/1134/1/22/11584243/22_1_online.pdf) URL https://doi.org/10.1063/1.3149495
  * [5] Wen X G 2012 Phys. Rev. B 85(8) 085103 URL https://link.aps.org/doi/10.1103/PhysRevB.85.085103
  * [6] Pollmann F, Berg E, Turner A M and Oshikawa M 2012 Phys. Rev. B 85(7) 075125 URL https://link.aps.org/doi/10.1103/PhysRevB.85.075125
  * [7] Chen X, Gu Z C, Liu Z X and Wen X G 2013 Phys. Rev. B 87(15) 155114 URL https://link.aps.org/doi/10.1103/PhysRevB.87.155114
  * [8] Giamarchi T 2004 Quantum Physics in One Dimension (Oxford University Press) ISBN 0198525001 URL https://www.ebook.de/de/product/3263637/thierry_giamarchi_quantum_physics_in_one_dimension.html
  * [9] Mermin N D and Wagner H 1966 Phys. Rev. Lett. 17(22) 1133–1136 URL https://link.aps.org/doi/10.1103/PhysRevLett.17.1133
  * [10] Berg E, Dalla Torre E G, Giamarchi T and Altman E 2008 Phys. Rev. B 77(24) 245119 URL https://link.aps.org/doi/10.1103/PhysRevB.77.245119
  * [11] Montorsi A and Roncaglia M 2012 Phys. Rev. Lett. 109(23) 236404 URL https://link.aps.org/doi/10.1103/PhysRevLett.109.236404
  * [12] Barbiero L, Montorsi A and Roncaglia M 2013 Phys. Rev. B 88(3) 035109 URL https://link.aps.org/doi/10.1103/PhysRevB.88.035109
  * [13] Rath S P, Simeth W, Endres M and Zwerger W 2013 Annals of Physics 334 256–271 ISSN 0003-4916 URL https://www.sciencedirect.com/science/article/pii/S0003491613000717
  * [14] Fazzini S, Becca F and Montorsi A 2017 Phys. Rev. Lett. 118(15) 157602 URL https://link.aps.org/doi/10.1103/PhysRevLett.118.157602
  * [15] den Nijs M and Rommelse K 1989 Phys. Rev. B 40(7) 4709–4734 URL https://link.aps.org/doi/10.1103/PhysRevB.40.4709
  * [16] Manmana S R, Essin A M, Noack R M and Gurarie V 2012 Phys. Rev. B 86(20) 205119 URL https://link.aps.org/doi/10.1103/PhysRevB.86.205119
  * [17] Montorsi A, Dolcini F, Iotti R C and Rossi F 2017 Phys. Rev. B 95(24) 245108 URL https://link.aps.org/doi/10.1103/PhysRevB.95.245108
  * [18] Mazurenko A, Chiu C S, Ji G, Parsons M F, Kanász-Nagy M, Schmidt R, Grusdt F, Demler E, Greif D and Greiner M 2017 Nature 545 462–466 ISSN 1476-4687 URL https://doi.org/10.1038/nature22362
  * [19] Tanzi L, Roccuzzo S M, Lucioni E, Famà F, Fioretti A, Gabbanini C, Modugno G, Recati A and Stringari S 2019 Nature 574 382–385
  * [20] Guo M, Böttcher F, Hertkorn J, Schmidt J N, Wenzel M, Büchler H P, Langen T and Pfau T 2019 Nature 574 386–389
  * [21] Chomaz L, Petter D, Ilzhöfer P, Natale G, Trautmann A, Politi C, Durastante G, van Bijnen R M W, Patscheider A, Sohmen M, Mark M J and Ferlaino F 2019 Phys. Rev. X 9(2) 021012 URL https://link.aps.org/doi/10.1103/PhysRevX.9.021012
  * [22] Zahn H P, Singh V P, Kosch M N, Asteria L, Freystatzky L, Sengstock K, Mathey L and Weitenberg C 2022 Phys. Rev. X 12(2) 021014 URL https://link.aps.org/doi/10.1103/PhysRevX.12.021014
  * [23] Su L, Douglas A, Szurek M, Groth R, Ozturk S F, Krahn A, Hébert A H, Phelps G A, Ebadi S, Dickerson S, Ferlaino F, Marković O and Greiner M 2023 Nature 622 724–729 URL https://doi.org/10.1038/s41586-023-06614-3
  * [24] Shao H J, Wang Y X, Zhu D Z, Zhu Y S, Sun H N, Chen S Y, Zhang C, Fan Z J, Deng Y, Yao X C, Chen Y A and Pan J W 2024 Observation of the antiferromagnetic phase transition in the fermionic hubbard model (Preprint 2402.14605)
  * [25] Haller E, Hart R, Mark M J, Danzl J G, Reichsöllner L, Gustavsson M, Dalmonte M, Pupillo G and Nägerl H C 2010 Nature 466 597–600 URL https://doi.org/10.1038/nature09259
  * [26] Endres M, Cheneau M, Fukuhara T, Weitenberg C, Schauß P, Gross C, Mazza L, Bañuls M C, Pollet L, Bloch I and Kuhr S 2011 Science 334 200–203 (Preprint https://www.science.org/doi/pdf/10.1126/science.1209284) URL https://www.science.org/doi/abs/10.1126/science.1209284
  * [27] Wei D, Adler D, Srakaew K, Agrawal S, Weckesser P, Bloch I and Zeiher J 2023 Phys. Rev. X 13(2) 021042 URL https://link.aps.org/doi/10.1103/PhysRevX.13.021042
  * [28] Hur J, Lee W, Kwon K, Huh S, Cho G Y and Choi J y 2024 Phys. Rev. X 14(1) 011003 URL https://link.aps.org/doi/10.1103/PhysRevX.14.011003
  * [29] de Léséleuc S, Lienhard V, Scholl P, Barredo D, Weber S, Lang N, Büchler H P, Lahaye T and Browaeys A 2019 Science 365 775–780 (Preprint https://www.science.org/doi/pdf/10.1126/science.aav9105) URL https://www.science.org/doi/abs/10.1126/science.aav9105
  * [30] Sompet P, Hirthe S, Bourgund D, Chalopin T, Bibo J, Koepsell J, Bojović P, Verresen R, Pollmann F, Salomon G, Gross C, Hilker T A and Bloch I 2022 Nature 606 484–488 URL https://doi.org/10.1038/s41586-022-04688-z
  * [31] Lewenstein M, Sanpera A, Ahufinger V, Damski B, Sen(De) A and Sen U 2007 Adv. Phys. 56 243–379 URL https://doi.org/10.1080/00018730701223200
  * [32] Gross C and Bloch I 2017 Science 357 995–1001
  * [33] Altman E, Demler E and Lukin M D 2004 Phys. Rev. A 70(1) 013603 URL https://link.aps.org/doi/10.1103/PhysRevA.70.013603
  * [34] Argüello-Luengo J, Julià-Farré S, Lewenstein M, Weitenberg C and Barbiero L 2024 Probing spontaneously symmetry-broken phases with spin-charge separation through noise correlation measurements (Preprint 2404.08374)
  * [35] Gross C and Bakr W S 2021 Nat. Phys. 17 1316–1323 URL https://doi.org/10.1038/s41567-021-01370-5
  * [36] Fazzini S, Barbiero L and Montorsi A 2019 Phys. Rev. Lett. 122(10) 106402 URL https://link.aps.org/doi/10.1103/PhysRevLett.122.106402
  * [37] Tocchio L F, Becca F and Montorsi A 2019 SciPost Phys. 6 018 URL https://scipost.org/10.21468/SciPostPhys.6.2.018
  * [38] Montorsi A, Fazzini S and Barbiero L 2020 Phys. Rev. A 101(4) 043618 URL https://link.aps.org/doi/10.1103/PhysRevA.101.043618
  * [39] Daley A J, Taylor J M, Diehl S, Baranov M and Zoller P 2009 Phys. Rev. Lett. 102(4) 040402 URL https://link.aps.org/doi/10.1103/PhysRevLett.102.040402
  * [40] Mazza L, Rizzi M, Lewenstein M and Cirac J I 2010 Phys. Rev. A 82(4) 043629 URL https://link.aps.org/doi/10.1103/PhysRevA.82.043629
  * [41] Diehl S, Baranov M, Daley A J and Zoller P 2010 Phys. Rev. Lett. 104(16) 165301 URL https://link.aps.org/doi/10.1103/PhysRevLett.104.165301
  * [42] Dalmonte M, Di Dio M, Barbiero L and Ortolani F 2011 Phys. Rev. B 83(15) 155110 URL https://link.aps.org/doi/10.1103/PhysRevB.83.155110
  * [43] Lee Y W and Yang M F 2010 Phys. Rev. A 81(6) 061604 URL https://link.aps.org/doi/10.1103/PhysRevA.81.061604
  * [44] Bonnes L and Wessel S 2011 Phys. Rev. Lett. 106(18) 185302 URL https://link.aps.org/doi/10.1103/PhysRevLett.106.185302
  * [45] Perali A, Neilson D and Hamilton A R 2013 Phys. Rev. Lett. 110(14) 146803 URL https://link.aps.org/doi/10.1103/PhysRevLett.110.146803
  * [46] Singh M, Mishra T, Pai R V and Das B P 2014 Phys. Rev. A 90(1) 013625 URL https://link.aps.org/doi/10.1103/PhysRevA.90.013625
  * [47] Macia A, Astrakharchik G E, Mazzanti F, Giorgini S and Boronat J 2014 Phys. Rev. A 90(4) 043623 URL https://link.aps.org/doi/10.1103/PhysRevA.90.043623
  * [48] Astrakharchik G E, Zillich R E, Mazzanti F and Boronat J 2016 Phys. Rev. A 94(6) 063630 URL https://link.aps.org/doi/10.1103/PhysRevA.94.063630
  * [49] Bilitewski T and Cooper N R 2016 Phys. Rev. A 94(2) 023630 URL https://link.aps.org/doi/10.1103/PhysRevA.94.023630
  * [50] Hubert C, Baruchi Y, Mazuz-Harpaz Y, Cohen K, Biermann K, Lemeshko M, West K, Pfeiffer L, Rapaport R and Santos P 2019 Phys. Rev. X 9(2) 021026 URL https://link.aps.org/doi/10.1103/PhysRevX.9.021026
  * [51] Zimmerman M, Rapaport R and Gazit S 2022 Proceedings of the National Academy of Sciences 119 e2205845119 (Preprint https://www.pnas.org/doi/pdf/10.1073/pnas.2205845119) URL https://www.pnas.org/doi/abs/10.1073/pnas.2205845119
  * [52] Padhan A, Parida R, Lahiri S, Giri M K and Mishra T 2023 Phys. Rev. A 108(1) 013316 URL https://link.aps.org/doi/10.1103/PhysRevA.108.013316
  * [53] Degli Esposti Boschi C, Montorsi A and Roncaglia M 2016 Phys. Rev. B 94(8) 085119 URL https://link.aps.org/doi/10.1103/PhysRevB.94.085119
  * [54] Paredes B and Cirac J I 2003 Phys. Rev. Lett. 90(15) 150402 URL https://link.aps.org/doi/10.1103/PhysRevLett.90.150402
  * [55] Kuklov A, Prokof’ev N and Svistunov B 2004 Phys. Rev. Lett. 92(3) 030403 URL https://link.aps.org/doi/10.1103/PhysRevLett.92.030403
  * [56] Kuklov A, Prokof’ev N and Svistunov B 2004 Phys. Rev. Lett. 92(5) 050402 URL https://link.aps.org/doi/10.1103/PhysRevLett.92.050402
  * [57] Hu A, Mathey L, Danshita I, Tiesinga E, Williams C J and Clark C W 2009 Phys. Rev. A 80(2) 023619 URL https://link.aps.org/doi/10.1103/PhysRevA.80.023619
  * [58] Hu A, Mathey L, Williams C J and Clark C W 2010 Phys. Rev. A 81(6) 063602 URL https://link.aps.org/doi/10.1103/PhysRevA.81.063602
  * [59] Grémaud B and Batrouni G G 2021 Phys. Rev. Lett. 127(2) 025301 URL https://link.aps.org/doi/10.1103/PhysRevLett.127.025301
  * [60] Chen W, Hida K and Sanctuary B C 2003 Phys. Rev. B 67(10) 104401 URL https://link.aps.org/doi/10.1103/PhysRevB.67.104401
  * [61] Dalla Torre E G, Berg E and Altman E 2006 Phys. Rev. Lett. 97(26) 260401 URL https://link.aps.org/doi/10.1103/PhysRevLett.97.260401
  * [62] Ejima S and Fehske H 2015 Phys. Rev. B 91(4) 045121 URL https://link.aps.org/doi/10.1103/PhysRevB.91.045121
  * [63] Giamarchi T 2003 Quantum Physics in One Dimension (Oxford University Press) ISBN 9780198525004 URL https://doi.org/10.1093/acprof:oso/9780198525004.001.0001
  * [64] Schulz H J 1986 Phys. Rev. B 34(9) 6372–6385 URL https://link.aps.org/doi/10.1103/PhysRevB.34.6372
  * [65] McCulloch I P 2008 Infinite size density matrix renormalization group, revisited (Preprint 0804.2509)
  * [66] Hauschild J and Pollmann F 2018 SciPost Phys. Lect. Notes 5 (Preprint 1805.00055) URL https://scipost.org/10.21468/SciPostPhysLectNotes.5
  * [67] Schollwöck U 2011 Annals of Physics 326 96–192 ISSN 0003-4916 january 2011 Special Issue URL https://www.sciencedirect.com/science/article/pii/S0003491610001752
  * [68] Eisert J, Cramer M and Plenio M B 2010 Rev. Mod. Phys. 82(1) 277–306 URL https://link.aps.org/doi/10.1103/RevModPhys.82.277
  * [69] Catani J, De Sarlo L, Barontini G, Minardi F and Inguscio M 2008 Phys. Rev. A 77(1) 011603 URL https://link.aps.org/doi/10.1103/PhysRevA.77.011603
  * [70] Trotzky S, Cheinet P, Fölling S, Feld M, Schnorrberger U, Rey A M, Polkovnikov A, Demler E A, Lukin M D and Bloch I 2008 Science 319 295–299 (Preprint https://www.science.org/doi/pdf/10.1126/science.1150841) URL https://www.science.org/doi/abs/10.1126/science.1150841
  * [71] Best T, Will S, Schneider U, Hackermüller L, van Oosten D, Bloch I and Lühmann D S 2009 Phys. Rev. Lett. 102(3) 030408 URL https://link.aps.org/doi/10.1103/PhysRevLett.102.030408
  * [72] Weld D M, Miyake H, Medley P, Pritchard D E and Ketterle W 2010 Phys. Rev. A 82(5) 051603 URL https://link.aps.org/doi/10.1103/PhysRevA.82.051603
  * [73] Amato-Grill J, Jepsen N, Dimitrova I, Lunden W and Ketterle W 2019 Phys. Rev. A 99(3) 033612 URL https://link.aps.org/doi/10.1103/PhysRevA.99.033612
  * [74] Jepsen P N, Amato-Grill J, Dimitrova I, Ho W W, Demler E and Ketterle W 2020 Nature 588 403–407 URL https://doi.org/10.1038/s41586-020-3033-y
  * [75] Chung W C, de Hond J, Xiang J, Cruz-Colón E and Ketterle W 2021 Phys. Rev. Lett. 126(16) 163203 URL https://link.aps.org/doi/10.1103/PhysRevLett.126.163203
  * [76] Jepsen P N, Ho W W, Amato-Grill J, Dimitrova I, Demler E and Ketterle W 2021 Phys. Rev. X 11(4) 041054 URL https://link.aps.org/doi/10.1103/PhysRevX.11.041054
  * [77] Safavi-Naini A, Capogrosso-Sansone B and Kuklov A 2014 Phys. Rev. A 90(4) 043604 URL https://link.aps.org/doi/10.1103/PhysRevA.90.043604