Emergent bosons in the fermionic two-leg flux ladder

TL;DR

This paper investigates how fermionic particles in a two-leg flux ladder system can form bosonic pairs under certain interactions, leading to a phase transition from vortex density waves to charge density waves, demonstrated through DMRG simulations.

Contribution

It reveals the emergence of bosonic behavior from fermionic systems in a flux ladder, highlighting a flux-induced quantum phase transition.

Findings

Fermions form tightly bound pairs acting as bosons under specific interactions.

A flux-induced quantum phase transition from vortex density wave to charge density wave.

Numerical evidence of bosonic-like phases in a fermionic flux ladder.

Abstract

We study the emergence of bosonic pairs in a system of two coupled one-dimensional fermionic chains subject to a gauge flux (two-leg flux ladder), with both attractive and repulsive interaction. In the presence of strong attractive nearest-neighbor interaction and repulsive next-to-nearest-neighbor interaction, the system crosses into a regime in which fermions form tightly bound pairs, which behave as bosonic entities. By means of numerical simulations based on the density-matrix-renormalization-group (DMRG) method, we show in particular that in the strongly paired regime, the gauge flux induces a quantum phase transition of the Ising type from vortex density wave (VDW) to a charge density wave (CDW), characteristic of bosonic systems.

Equations

$${\hat{H}}_0=-t\sum _{j=1}^{L-1}\sum _{m=\pm 1/2}{\hat{c}}_{j,m}^{†}{\hat{c}}_{j+1,m}+H.c.$$

$${\hat{H}}_{\perp }=t_{\perp }\sum _{j=1}^L{e}^{-i\Phi j}{\hat{c}}_{j,-\frac{1}{2}}^{†}{\hat{c}}_{j,+\frac{1}{2}}+H.c.$$

$${\hat{H}}_{\mathrm{int}}=\sum _{m=\pm 1/2}\left(V\sum _{j=1}^{L-1}{\hat{n}}_{j,m}{\hat{n}}_{j+1,m}+W\sum _{j=1}^{L-2}{\hat{n}}_{j,m}{\hat{n}}_{j+2,m}\right)$$

$${\hat{c}}_{r,\alpha ,m}\sim \sum _p{e}^{-i\sqrt{\pi }{\hat{\theta }}_{\alpha ,m}(r)}{e}^{-ip\pi (n/2)r}{e}^{ip\sqrt{\pi }{\hat{\phi }}_{\alpha ,m}(r)},$$

$${\hat{\phi }}_{\pm ,m}=\frac{{\hat{\phi }}_{\uparrow ,m}\pm {\hat{\phi }}_{\downarrow ,m}}{\sqrt{2}}{\hat{\theta }}_{\pm ,m}=\frac{{\hat{\theta }}_{\uparrow ,m}\pm {\hat{\theta }}_{\downarrow ,m}}{\sqrt{2}},$$

$${\hat{B}}_{r,m}^{†}\sim e^{i\sqrt{2\pi }{\hat{\theta }}_{+,m}}\sum _p(e^{i2p\pi (n/2)r}{e}^{-ip\sqrt{2\pi }{\hat{\phi }}_{+,m}}$$

$$+e^{-ip\sqrt{2\pi }{\hat{\phi }}_{-,m}}),$$

$$\displaystyle V\sum_{r}(\hat{n}_{r,\uparrow,{\color[rgb]{0,0,0}m}}\hat{n}_{r,\downarrow,{\color[rgb]{0,0,0}m}}+\hat{n}_{r,\downarrow,{\color[rgb]{0,0,0}m}}\hat{n}_{r+1,\uparrow,{\color[rgb]{0,0,0}m}})$$

$$\displaystyle\hskip 56.9055pt\sim\int dr\,\cos\left(2\sqrt{2\pi}\,\hat{\varphi}_{-,{\color[rgb]{0,0,0}m}}\right)\,\,.$$

$${\hat{\theta }}_{+,m}=\frac{{\hat{\theta }}_{B,m}}{\sqrt{2}}{\hat{\phi }}_{+,m}=\sqrt{2}{\hat{\phi }}_{B,m},$$

$${\hat{B}}_{r,m}^{†}\sim e^{i\sqrt{\pi }{\hat{\theta }}_{B,m}(r)}\sum _q{e}^{iq\pi n_{B}r}{e}^{-iq\sqrt{\pi }{\hat{\phi }}_{B,m}(r)},$$

$$\sum _z{e}^{-i2\Phi z}{\hat{C}}_{z,-\frac{1}{2}}^{†}{\hat{C}}_{z,+\frac{1}{2}}+H.c.$$

$$\sim \int dz\cos \left[\sqrt{\pi }\left({\hat{\theta }}_{B,+\frac{1}{2}}-{\hat{\theta }}_{B,-\frac{1}{2}}\right)+2\Phi z\right],$$

$$\sum _z{\hat{n}}_{z,-\frac{1}{2}}{\hat{n}}_{z,+\frac{1}{2}}\sim \int dz\cos \left[2\sqrt{\pi }\left({\hat{\phi }}_{B,+\frac{1}{2}}-{\hat{\phi }}_{B,-\frac{1}{2}}\right)\right].$$

$${\hat{\theta }}_{B,s/a}=\frac{{\hat{\theta }}_{B,+\frac{1}{2}}\pm {\hat{\theta }}_{B,-\frac{1}{2}}}{\sqrt{2}}{\hat{\phi }}_{B,s/a}=\frac{{\hat{\phi }}_{B,+\frac{1}{2}}\pm {\hat{\phi }}_{B,-\frac{1}{2}}}{\sqrt{2}}.$$

$$\displaystyle{\color[rgb]{0,0,0}\hat{H}^{(a)}=\hat{H}^{(a)}_{\rm LL}}+\tilde{t}_{\perp}\int dz\,\cos\left(\sqrt{2\pi}\,\hat{\theta}_{B,a}+2\Phi z\right)$$

$$+\tilde{U}\int dz\cos \left(2\sqrt{2\pi }{\hat{\phi }}_{B,a}\right).$$

$$r_B(\Phi ):=\frac{2n_B(\Phi )}{n},$$

$$\delta n_{B,j}(\Phi )=\left|⟨{\Psi }_{\mathrm{GS}}(\Phi )|\left({\hat{n}}_{B,j,-\frac{1}{2}}-{\hat{n}}_{B,j,+\frac{1}{2}}\right)|{\Psi }_{\mathrm{GS}}(\Phi )⟩\right|,$$

$$S(\ell )=a+\frac{c}{6}\ln \left[\left(\frac{2L}{\pi }\right)\sin \left(\frac{\pi \ell }{L}\right)\right],$$

$$C(x)=\left|⟨{\Psi }_{\mathrm{GS}}|{\hat{c}}_x^{†}{\hat{c}}_0|{\Psi }_{\mathrm{GS}}⟩\right|.$$

$${\hat{H}}_U=U\sum _j{\hat{n}}_{j,-\frac{1}{2}}{\hat{n}}_{j,+,\frac{1}{2}},$$

$${\hat{J}}_{\perp ,j}^{(B)}=-i{t}_{\perp }\left(e^{-i\Phi (2j+1)}{\hat{B}}_{j,-\frac{1}{2}}^{†}{\hat{B}}_{j,+\frac{1}{2}}-H.c.\right).$$

$${\left({\bar{J}}_{\perp }^{(B)}\right)}^2=\frac{1}{L-2\Delta L}\sum _{j=\Delta L}^{L-\Delta L}{\left|⟨{\Psi }_{\mathrm{GS}}|{\hat{J}}_{\perp ,j}^{(B)}|{\Psi }_{\mathrm{GS}}⟩\right|}^2,$$

Full text

I Introduction

Exotic phases of matter emerging from the interplay between strong interactions, magnetic fields and enhanced quantum fluctuations due to low dimensionality have been an active field of research in condensed-matter physics during the last decades, both for fermionic and bosonic systems. In the last years, a renewed theoretical and experimental interest in the realization and characterization of such intriguing phases has been triggered by the advances in the field of ultra-cold atomic gases in optical lattices with artificial gauge fields, the latter mimicking the effects of applied magnetic fields Bloch (2005); Lewenstein et al. (2007); Bloch et al. (2008); Dalibard et al. (2011); Boada et al. (2015); Goldman et al. (2016); Kang et al. (2018). Such techniques provide the ability of creating and manipulating matter (synthetic matter) with unprecedented precision.

In this respect, systems of many coupled one-dimensional (1D) chains immersed in a gauge field (flux ladders) represent a versatile platform in which such effects can be studied, in which dimensionality is controlled by the number of wires. Because of their 1D nature, the toolbox to theoretically analyze phases in these systems is provided by ad-hoc numerical algorithms based on the density-matrix-renormalization-group (DMRG) White (1992); Schollwöck (2005) or matrix-product-state (MPS) Schollwöck (2011) formalism, and effective field theories, such as bosonization Gogolin et al. (2004); Giamarchi (2003).

The minimal setup in which gauge-field effects can be obtained is the two-leg flux ladder, i.e., two connected chains. Several works have focused on this system, discussing interesting aspects both for bosons Orignac and Giamarchi (2001); Dhar et al. (2012); Crépin et al. (2011); Atzmon and Shimshoni (2011); Petrescu and Le Hur (2013); Dhar et al. (2013); Wei and Mueller (2014); Tokuno and Georges (2014); Di Dio et al. (2015a); Piraud et al. (2015); Di Dio et al. (2015b); Kolley et al. (2015); Natu (2015); Greschner et al. (2015, 2016); Orignac et al. (2016); Calvanese Strinati et al. (2018); Greschner and Heidrich-Meisner (2018); Loida et al. (2018); Buser et al. (2019) and fermions Narozhny et al. (2005); Carr et al. (2006); Mazza et al. (2015); Barbarino et al. (2015, 2016); Ghosh et al. (2017); Taddia et al. (2017); Lacki et al. (2016); Sun (2016); Haller et al. (2018). In particular, it has been shown that flux ladders can host phases that, at low energies, are analogous to fractional quantum Hall (FQH) phases Petrescu and Le Hur (2015); Cornfeld and Sela (2015); Calvanese Strinati et al. (2017); Petrescu et al. (2017); Calvanese Strinati et al. (2019), or manifest quantum phase transitions from superconducting (SC) to Mott insulating phases Orignac and Giamarchi (2001); Atzmon and Shimshoni (2011); Piraud et al. (2015); Orignac et al. (2016); Di Dio et al. (2015a).

The phase diagram of the fermionic and bosonic two-leg flux ladder has been discussed in details for different models of interactions Carr et al. (2006, 2013); Natu (2015), both attractive and repulsive. While it is expected that attractive on-site interactions in the fermionic ladder lead to the formation of fermionic pairs, which behave as bosonic particles Orignac and Giamarchi (2001); Giamarchi (2003), to the best of our knowledge, a detailed study of how such bosons emerge in the fermionic ladder for longer-range interactions is still missing. In this paper, we aim to bridge this gap, studying the emergence of bosonic particles in the fermionic two-leg flux ladder with attractive and repulsive finite-range interactions.

This paper is organized as follows. We introduce our microscopic model in Sec. II, and discuss its low-energy theory in the strongly-interacting regime in Sec. III. We then present our numerical results in Sec. IV. We draw our conclusions in Sec. V, and present additional numerical data in the Appendixes.

II Model

The system consists of two 1D chains immersed in a gauge flux (Fig. 1). In the following, we consider open boundary conditions (OBC) both in the longitudinal ($j$) and transverse ($m$) dimensions, and model our system by the Hamiltonian $\hat{H}_{F}=\hat{H}_{0}+\hat{H}_{\perp}+\hat{H}_{\rm int}$, where

[TABLE]

[TABLE]

[TABLE]

in which we set the lattice constant to unity. In the Hamiltonian, $\hat{c}_{j,m}$ ($\hat{c}^{\dagger}_{j,m}$) is the annihilation (creation) operator of a fermion on site $j$ and on leg $m$, and $\hat{n}_{j,m}=\hat{c}^{\dagger}_{j,m}\hat{c}_{j,m}$ is the fermionic density operator; $t$ and $t_{\perp}$ denote the longitudinal and transverse hopping parameters, respectively, and $\Phi$ is the gauge flux per plaquette. We denote by $L$ the number of rungs of the ladder, and $N$ the total number of particles in the system. We define the total particle density as $n=N/L$. The interaction Hamiltonian $\hat{H}_{\rm int}$ accounts for both intra-leg nearest-neighbour (NN) and next-to-nearest-neighbour (NNN) interaction, whose strengths are identified by $V$ and $W$, respectively. In the following, if not explicit, we use $t$ as reference energy scale. Since we aim at forming fermionic pairs, we consider $V<0$ and $W>0$. The first condition induces NN particles to bind, whereas the second one prevents clusters from forming.

III Low-energy theory for strong interactions

In this section, we derive the low-energy theory of the model introduced in Sec. II, in the strongly-interacting regime. For sufficiently large but finite $|V|$ and $W$, which is the case of interest, the resulting ground state (GS) is composed of tightly-bound fermionic pairs subject to a hard pairing gap Ruhman and Altman (2017); Borla et al. (2019) (see also Appendix A), with effective hopping parameters $\tilde{t}\sim t^{2}/|V|$ and $\tilde{t}_{\perp}\sim t_{\perp}^{2}/|V|$. In this strongly-coupled limit, we bosonize the model starting from the fermionic pair operator $\hat{B}^{\dagger}_{j,m}=\hat{c}^{\dagger}_{j,m}\hat{c}^{\dagger}_{j+1,m}$.

Since the NN and NNN interactions couple site with different and equal parity, respectively, one can interpret each chain of length $L$ (suppose $L$ even) as the composition of two sublattices of length $L/2$ each, identified by a pseudo-spin index $\alpha=\uparrow,\downarrow$ and lattice coordinate $r$, such that $\alpha=\uparrow$ includes the sites of the of the original lattice with $j$ odd, and $\alpha=\downarrow$ includes those with $j$ even [Fig. 2(a)]: $j=2r-R(\alpha)$, where $R(\uparrow)=1$ and $R(\downarrow)=0$. The fermionic lattice operator is then recast as $\hat{c}_{j,m}\rightarrow\hat{c}_{r,\alpha,m}$, whose bosonized version reads Giamarchi (2003); Cazalilla et al. (2011)

[TABLE]

where $\hat{\theta}_{\alpha,m}(r)$ and $\hat{\varphi}_{\alpha,m}(r)$ are the phase and density bosonic fields, respectively, which obey the canonical commutation relations $[\hat{\varphi}_{\alpha,m}(r),\partial_{r^{\prime}}\hat{\theta}_{\alpha^{\prime},m^{\prime}}(r^{\prime})]=\delta_{\alpha,\alpha^{\prime}}\,\delta_{m,m^{\prime}}\,\delta(r-r^{\prime})$, and $p$ is an odd integer. The pair operator, in the remapped lattice, reads either $\hat{B}^{\dagger}_{r,m}=\hat{c}^{\dagger}_{r,\uparrow,m}\hat{c}^{\dagger}_{r,\downarrow,m}$ or $\hat{B}^{\dagger}_{r,r+1,m}=\hat{c}^{\dagger}_{r+1,\uparrow,m}\hat{c}^{\dagger}_{r,\downarrow,m}$ (Fig. 2). We discuss now the bosonization form of $\hat{B}^{\dagger}_{r,m}$.

By introducing the fields

[TABLE]

the pair operator is bosonized using Eq. (4):

[TABLE]

for odd $p$. Accordingly, the lowest non-oscillating harmonic of the NN interaction term now reads

[TABLE]

When $V<0$, it pins the \hat{\varphi}_{-,{\color[rgb]{0,0,0}m}} fields Giamarchi (2003) in Eq. (6) to \hat{\varphi}_{-,{\color[rgb]{0,0,0}m}}=0, providing an effective $p=0$ harmonic. A further canonical transformation

[TABLE]

by introducing $q=2p$ and $n_{B}=n/2$, allows to recast Eq. (6) as

[TABLE]

for $q$ even, therefore recovering a bosonic operator Giamarchi (2003); Cazalilla et al. (2011). An analogous result is found for $\hat{B}^{\dagger}_{r,r+1,m}$. This result allows us to treat the pair as a single bosonic particle: $\hat{B}^{\dagger}_{j,m}=\hat{c}^{\dagger}_{j,m}\hat{c}^{\dagger}_{j+1,m}\rightarrow\hat{C}^{\dagger}_{z,m}$, localized at $z=j+1/2$, and therefore coarse-grain the system [Fig. 2(b)]. In the strongly-paired regime, a pair experiences a flux per plaquette equal to $2\Phi$.

The NNN interaction between fermions represent an intra-chain repulsive NN interaction $\tilde{W}$ between pairs. Moreover, even if the original fermions are not coupled by an inter-chain interaction, the presence of $t_{\perp}$, in addition to providing the inter-leg pair (Josephson) tunnelling

[TABLE]

can perturbatively generate all interactions processes allowed by symmetry. The minimal one that one expects is an on-site interaction between pairs on different legs

[TABLE]

Therefore, we expect the system in the strongly-paired regime to be described by $\hat{H}_{\rm eff}\simeq\hat{H}^{(s)}+\hat{H}^{(a)}$, in terms of independent symmetric and anti-symmetric fields Atzmon and Shimshoni (2011)

[TABLE]

Here, $\hat{H}^{(s)}=\hat{H}^{(s)}_{\rm LL}$ is a gapless Luttinger liquid, whereas $\hat{H}^{(a)}$ is the self-dual sine-Gordon model Lecheminant et al. (2002)

[TABLE]

This model belongs to the Ising universality class Lecheminant et al. (2002), and exhibits an Ising-type quantum phase transition. We will now use this result in order to validate the emergence of bosons in the fermionic chain.

IV Numerical results

We now present our numerical results on the emergence of bosonic pairs in the system. We simulate the Hamiltonian $\hat{H}_{F}$ by means of a DMRG algorithm that is the same used in Ref. Rossini et al. (2019). For fixed values of $L$, $N$, $t_{\perp}$, $V$, $W$, and $\Phi$, after the initial infinite-DMRG sweep, a number of sweeps $S$ of finite-DMRG are performed in order to variationally find the density matrix of the system. During the sweeps, we truncate the dimension of the density matrix keeping up to $M$ states, where $M$ is chosen such that the truncation error does not exceed $\sim 10^{-7}$ Schollwöck (2011). Because of the high numerical complexity of the problem, we can scan a limited range of parameters. Specifically, here, we present numerical data for $t_{\perp}=0.3\,t$, $n=1/4$, and keep $W=-V>0$. We use $120\leq M\leq 200$ and $3\leq S\leq 5$ depending on the observable that we measure, and on the value of $L$. We refer the interested reader to Appendix B for additional numerical data.

With the numerical algorithm that we use, we can measure only one- and two-point observables in terms of the original fermions $\hat{c}_{j,m}$, which means at most on-site observables for the emergent bosons $\hat{B}_{j,m}$. However, as we discuss below, the emergent bosonic physics can be detected already by looking at the pair density $n_{B,j,m}(\Phi)=\langle\Psi_{\rm GS}(\Phi)|\hat{n}_{B,j,m}|\Psi_{\rm GS}(\Phi)\rangle$, where $\hat{n}_{B,j,m}=\hat{B}^{\dagger}_{j,m}\hat{B}_{j,m}$, which is the focus of the rest of our work (for a discussion on the measurement of the inter-leg current, the reader is referred to Appendix C).

IV.1 Detecting fermionic pairs

A first evidence of the formation of pairs is provided by comparing the average local pair density $n_{B}(\Phi)=L^{-1}\sum_{j}n_{B,j,m}(\Phi)$ with the fermionic density $n$. Since we expect $n_{B}\ll n$ and $n_{B}=n/2$ [Eq. (9)] in the unpaired and paired regimes, respectively, monitoring how the quantity

[TABLE]

varies as $W$ is scanned from $W=0$ to large values provides information on the emergence of bosonic pairs in the system.

The result is shown in Fig. 3. We simulate $\hat{H}_{F}$ for $k=41$ values of $\Phi\in[0:2\pi]$, and show the flux average $\bar{r}_{B}=k^{-1}\sum_{\Phi}{r_{B}}(\Phi)$. We see that, for small $W$, $\bar{r}_{B}\simeq 0$, whereas it approaches $\bar{r}_{B}=1$ as $W$ is increased. Between these two regimes, there is a wide range of $W$ in which $\bar{r}_{B}$ smoothly interpolates between [math] and $1$. In such a region, fermions and bosonic pairs coexist, and the number of pairs $n_{B}(\Phi)$ is found to fluctuate with $\Phi$, quantified by the uncertainties on the data. For $W\gtrsim 5$, instead, such fluctuations are suppressed, and the system approaches the fully-paired regime. Because of the large number of flux values that we need for each $W$, we use $L=32,64$ in order to keep a reasonable computational complexity. Importantly, for the simulated values of $L$, the data are almost overlapped, and show no finite-size scaling.

IV.2 Flux-driven Ising-type transition

We now discuss the existence of an Ising-type transition driven by the gauge flux, in the paired regime. The first striking feature is that, in this regime, the system undergoes a flux-driven transition between a vortex density wave (VDW) ($\Phi<\Phi_{c}$) and relative charge density wave (CDW) ($\Phi>\Phi_{c}$), for some critical value $\Phi_{c}$ that depends on the system parameters. This manifests itself in the spatial patterns of $n_{B,j,m}$ and the local currents along the ladder, as in Fig. 4(a),(b), in which the fermionic intra- and inter-chain currents Piraud et al. (2015); Calvanese Strinati et al. (2017) (arrows) are shown together with $n_{B,j,m}$ (blue dots). For $\Phi<\Phi_{c}$, an ordered arrays of vortices appears (along with a vanishing relative density $\delta n_{B,j}=n_{B,j,-\frac{1}{2}}-n_{B,j,+\frac{1}{2}}$ for all $j$), which is compatible with the locking of the relative phase field $\hat{\theta}_{B,a}$ Calvanese Strinati et al. (2017). Instead, for $\Phi>\Phi_{c}$, the relative density becomes periodically modulated, signalling a (staggered) CDW order (locking of the relative charge field $\hat{\varphi}_{B,a}$) (see also Appendix C).

This allows to focus on the local density imbalance between the two legs as a function of $\Phi$:

[TABLE]

in order to quantify the transition. Specifically, we compute the space average $\delta n_{B}(\Phi)=L^{-1}\sum_{j}\delta n_{B,j}(\Phi)$ scanning $\Phi$ through the transition. The key result is shown in Fig. 4(c). We compare the results in the paired regime ($W=5,50$) with those in the unpaired regime ($W=2.5$). As evident, no transition occurs for $W=2.5$ ($\delta n_{B}=0$ for all $\Phi$, signalling no density imbalance), whereas an increase of $\delta n_{B}$ around $\Phi_{c}$ is found for $W=5,50$ ($\delta n_{B}=0$ for $\Phi<\Phi_{c}$ and $\delta n_{B}>0$ for $\Phi>\Phi_{c}$, signalling the transition from the VDW to the CDW phase).

We now try to illuminate the nature of this transition, and quantify the scaling of the correlation length $\xi$. This can be estimated by adding a localized impurity on one site $j_{0}$ of the $m=+1/2$ leg: $\hat{H}_{F,\mu}=\hat{H}_{F}-\mu\,\hat{n}_{j_{0},+\frac{1}{2}}$, and analyzing the response of $\delta n_{B,j}$. The impurity locally enforces a density imbalance: if the GS is the VDW configuration, $\delta n_{B,j}$ is locally perturbed from the balanced configuration $\delta n_{B,j}=0$, but such a configuration is recovered after a characteristic length $\xi$: $\delta n_{B,j}\sim\delta n_{B,j_{0}}e^{-|j-j_{0}|/\xi}$. Instead, if the GS is a CDW configuration, the local imbalance forced by the impurity is preserved through the whole system ($\xi\rightarrow\infty$). By simulating $\hat{H}_{F,\mu}$, one can fit the data for $\delta n_{B,j}$ vs. $j$ as $\Phi$ is varied across the transition and thus extract $\xi^{-1}(\Phi)\sim\Delta(\Phi)$, which is expected to exhibit scaling behaviour compatible with the Ising model in (1+1)-D Atzmon and Shimshoni (2011); Gogolin et al. (2004); Sachdev (2001): $\Delta\sim|\Phi-\Phi_{c}|$ for $\Phi<\Phi_{c}$ and $\Delta=0$ for $\Phi>\Phi_{c}$.

The result of the simulation is shown in Fig. 5. We show $\xi^{-1}(\Phi)$ for $W=5,50$ and $L=80,96,128$. We observe that $\delta n_{B,j}$ exhibits an exponential decay, which is on top of spatial fluctuations (see also Fig. 4). In order to extract $\xi$ and account for such fluctuations, as well as finite-size effects, we fit the envelope of $\delta n_{B,j}$ with the function $f(j)=f_{0}e^{-|j-j_{0}|/\xi}$ (using $f_{0}$ and $\xi$ as fit parameters) three times, for $j\in[0.2\,L:L-\Delta L]$ and $\Delta L=0.15,0.2,0.3$. The resulting values of $\xi^{-1}$ and uncertainties are given by the average value and $(\max_{\Delta L}\xi^{-1}-\min_{\Delta L}\xi^{-1})/2$, respectively. Within our numerical precision and limitations due to finite-size effects, our results are consistent with the linear closing of the gap $\Delta\sim|\Phi-\Phi_{c}|$, confirming the Ising transition.

IV.3 Measuring the central charge

A further observable to test the low-energy physics as in Eq. (13) is given by the central charge $c$ Giamarchi (2003). A way to extract $c$ is to measure the von Neumann entropy (VNE), defined as $S_{\rm VNE}(\ell)=-{\rm Tr}[\hat{\rho}_{\ell}\ln(\hat{\rho}_{\ell})]$, $\hat{\rho}_{\ell}$ being the reduced density matrix of a subpart of the system of size $\ell$, and fit it via the expression Calabrese and Cardy (2004)

[TABLE]

for OBC. Our numerical results for $W=50$ are shown in Fig. 4(d). In order to measure $c$ reliably, we use $M=200$, which significantly increases the computational time. We thus use $L=96$. Because of the fluctuating behaviour of the VNE, we extract the values of $c$ and relative uncertainties as in Refs. Calvanese Strinati et al. (2017, 2018). Away from the transition point, our data of $c$ are consistent with the value $c=1$ expected from Eq. (13) (which contains a single gapless mode in the symmetric sector).

V Conclusions

We analyzed the emergence of bosonic pairs in the fermionic two-leg flux ladder with competing attractive and repulsive interaction. We provided a phenomenological low-energy description in the strongly paired regime, which predicts the existence of an Ising-type transition between phases related by vortex-charge duality, and corroborated its validity by means of DMRG simulations. Although our numerics was limited to specific observables and values of parameters, due to the challenging numerical complexity of the problem, we observed a flux-driven Ising-type transition focusing on the divergence of the correlation length of the relative density order. Our work opens the possibility of creating interfaces in the flux-ladder between FQH and SC phases, thus opening a new intriguing path towards the possibility of hosting parafermions in flux-ladders Clarke et al. (2013); Lindner et al. (2012); Vaezi (2013). We leave these promising perspectives for future work.

Acknowledgements

We thank Daniel Podolsky and Jonathan Ruhman for fruitful discussions. We are grateful to Davide Rossini for support. We acknowledge support from the Israel Science Foundation (ISF), Grants No. 231/14 and 993/19 (E. S. and M. C. S.) and No. 1452/14 (M. C. S.), and the U.S.-Israel Binational Science Foundation (BSF) Grant No. 2016130 and No. 2018726 (E. S. and M. C. S.).

Appendix A Direct evidence of the pairing gap

In Sec. III, we justify our bosonization treatment by the fact that, in the regime of parameters that we use, the presence of $V$ induces a hard pairing gap $\Delta_{\rm pair}$ between fermions, which corresponds to the energy that the system gains when two fermions bind to form a bosonic pair. In this appendix, we explicitly provide direct evidence of the presence of such a pairing gap.

An observable that can monitor the presence of a pairing gap is the fermionic two-point correlation function

[TABLE]

The correlation function in Eq. (17) was also used in different contexts, for example in Refs. Ruhman and Altman (2017); Borla et al. (2019), in order to detect the occurrence of the pairing gap. At large $x$, $C(x)$ decays as a power law in the unpaired regime, and as an exponential in the paired regime, indeed due to the absence and presence of a pairing gap, respectively.

We find such a behaviour also from our numerical data. In particular, for the sake of clarity, we show in Fig. 6 the correlation function $C(x)$ for the numerical data for $L=64$ presented in Fig. 3. By analyzing the long-distance behaviour of $C(x)$, we indeed see that, when $|V|$ is sufficiently small, $C(x)$ decays as a power-law, as can be appreciated by the linear decay in log-log scale (Fig. 6, left panel), whereas when $|V|$ is sufficiently large, $C(x)$ decays as an exponential, $C(x)\sim e^{-\Delta_{\rm pair}x}$, as can be seen by the linear decay in log-linear scale (Fig. 6, right panel). Specifically, we find that up to $|V|=3\,t$, $C(x)$ decays clearly as a power law, whereas for $|V|\geq 4\,t$, it decays exponentially. The same results is found for the other simulations presented in this paper. Specifically, we show in Fig. 7 the same analysis for the data in Fig. 8 (Appendix B), for which the simulation parameters are $L=128$, $N=32$ ($n=1/4$), $t_{\perp}=0.3\,t$, fixed $W=7\,t$.

This result tells us that, in the regime of parameters considered in Sec. IV (i.e., $|V|\geq 5\,t$), where the bosonic physics in which we are interested is discussed (the Ising VDW-CDW transition), the system is always in the strongly-paired regime, with the presence of a hard pairing gap. This also justifies our bosonization treatment in Sec. III.

Appendix B Additional numerical data

In this appendix, we show additional numerical data for different ranges of parameters, in order to show that the phenomenology that we discuss in our paper is not a consequence of a fine tuning of the system parameters.

B.1 Varying the interaction strengths

First, we relax the condition $W=-V$ considered in Sec. IV. We repeat the simulations as in Fig. 3, but here we fix $W$, and vary $V$. In particular, we show in Fig. 8 the result of a simulation keeping $W=7\,t$ and by scanning $V$ from $V=0$ to $|V|=7\,t$. We compute $\bar{r}_{B}$ as explained in Sec. IV. The other simulation parameters are: $L=128$, $N=32$ (i.e., $n=1/4$), $t_{\perp}=0.3\,t$, $M=120$ and $S=5$ sweeps. As we see, also in this other case, the quantity $\bar{r}_{B}$ displays a smooth increase from $\bar{r}_{B}\ll 1$ (mostly unpaired fermions) for small $|V|$, to $\bar{r}_{B}\rightarrow 1$ (strongly-paired regime) as $|V|$ is increased.

We further show the data for $\delta n_{B}(\Phi)$, and the density and current configuration along the ladder in Fig. 9, as in Fig. 4. Specifically, the data are shown for $W=7\,t$ and $V=-9\,t$. In panels (a) and (b), we show the ladder configuration for $\Phi/\pi=0.15$ (in the VDW phase) and for $\Phi/\pi=0.25$ (in the CDW phase), respectively. We see that, apart from specific quantitative details, the same phenomenology discussed throughout the paper for $W=-V$ arises also for this other choice of $V$ and $W$.

The same numerical simulations as in Figs. 8 and 9 are repeated for a different value of $W$, namely, we fix $W=3\,t$ and scan $|V|$ from $|V|=0$ to $|V|=6\,t$, since we observed the formation of clusters for larger values of $|V|$. The result of the simulation is shown in Figs. 10 and 11. Again, the data of $\bar{r}_{B}$ in Fig. 10 show a smooth increase from almost zero to $\bar{r}_{B}\rightarrow 1$ for sufficiently strong $|V|$. In Fig. 11, $\delta n_{B}(\Phi)$ for $|V|=5.5\,t$ is shown, together with the density and current pattern configuration along the ladder, displaying once again the same phenomenology (VDW-CDW transition).

We now discuss the presence of an inter-leg density-density interaction of the form

[TABLE]

at the level of the original microscopic Hamiltonian. The original microscopic model (Sec. II) does not include such a term, but as we say in Sec. III, such an interaction term is generated in the effective model by $t_{\perp}$, $V$ and $W$. We show in Fig. 12 the result of a simulation using the interaction Hamiltonian $\hat{H}_{\rm int}$ as in Eq. (3), with the inclusion of the on-site interaction term as in Eq. (18), with $U=1\,t$. As we see, the presence of $U>0$ does not dramatically affect the physics that we discussed throughout the paper.

B.2 Varying the particle density

In our paper, since we discuss the emergence of a gap in the antisymmetric sector of the emergent bosonic pairs, what is important is that we choose a value of $n$ such that we are away from any relevant lattice commensurability condition (i.e., $n=1$), which would create a gap also in the symmetric sector. For numerical convenience, we choose $n=1/4$. In order to demonstrate that our qualitative results do not rely on this particular choice, we repeated the numerical calculations with $n=1/8$ (using $L=160$ and $N=20$) and show the numerical results in Fig. 13. The other numerical parameters are $W=-V=5\,t$, $t_{\perp}=0.3\,t$, $M=120$ and $S=3$ sweeps. We indeed see that the same phenomenology as in Figs. 4, 9, and 11 arises, apart from the fact that the VDW and CDW appear with higher spatial period due to the smaller value of $n$. This result strengthens the conclusion that the physics discussed in the our paper is not a consequence of the specific value of $n$ that we consider.

B.3 Varying the inter-leg hopping parameter

In the numerical simulations presented thus far, we use a single fixed value of $t_{\perp}$, i.e., $t_{\perp}=0.3\,t$. We now show the numerical data of $\bar{r}_{B}$, density and current configuration along the ladder, and density difference $\delta n_{B}(\Phi)$, for a different value of $t_{\perp}$, namely $t_{\perp}=0.5\,t$, in order to further show that the phenology discussed in this paper does not depend on the specific choice of $t_{\perp}=0.3\,t$.

The additional data for $t_{\perp}=0.5\,t$ are shown in Figs. 14 and 15. We report the data for $\bar{r}_{B}$, for $L=64$, $N=12$ (i.e., $n=1/4$), in order to keep a reasonable numerical complexity, and we scan $V$ from $V=0$ to $|V|=10\,t$, keeping $V=-W$. We also show the ladder density and current configuration, and density difference $\delta n_{B}(\Phi)$, for $W=-V=7\,t$. Once again, we observe the same phenomenology as in the previous cases discussed in this paper, with a slight shift in the critical value of $\Phi_{c}$

In light of all these results, we are confident in concluding that the phenomenology discussed in our paper is not a consequence of a fine tuning of the system parameters.

Appendix C Numerical results for the inter-leg pair current

In this appendix, we provide an additional evidence of the fact that the low-energy physics of our system is indeed given by Eq. (13).

An additional observable that we can measure in order to detect the VDW-CDW phase transition is the inter-leg pair current, which on the lattice is defined by the operator

[TABLE]

Such an operator, at the low-energy level, is sensitive to the pinning of the phase field $\hat{\theta}_{B,a}$. Therefore, close to the transition, we expect the inter-leg current to be suppressed for $\Phi>\Phi_{c}$, i.e., when the relative density field $\hat{\varphi}_{B,a}$ is pinned and therefore $\hat{\theta}_{B,a}$ becomes strongly fluctuating.

We notice that, by symmetry, the space-average inter-leg current is always zero. The current between the legs can be then quantified for example by the square average current as

[TABLE]

where a few sites $\Delta L$ are removed from both chain ends in order to account for open boundary conditions and finite-size effects. We show the numerical data of the square average current in Fig. 16. The data are computed from the data in Fig. 4, panel (c). We see that the current reaches a maximum value around the critical value $\Phi=\Phi_{c}$, which agrees indeed with the same value of the flux at which $\delta n_{B}$ starts to be nonzero [see Fig. 4, panel (c)], and rapidly decreases for $\Phi>\Phi_{c}$. This result is in agreement with the presence of a relative density order for $\Phi>\Phi_{c}$, which then implies a strongly fluctuating phase order, detected by the suppression of the current.