Vortex lattice melting in a boson-ladder in artificial gauge field

TL;DR

This paper investigates how interleg attractive interactions in a boson ladder under an artificial gauge field lead to a sequence of phase transitions, including vortex melting, characterized by dislocations, in a cold atom experimental setup.

Contribution

It reveals a novel sequence of phase transitions involving Ising and disorder points, and predicts observable signatures in cold atom experiments.

Findings

Vortex state can be melted by dislocations due to attractive interactions.

Sequence of phase transitions includes Ising transition and disorder point.

Predicted observable effects include changes in spin current and structure factor.

Abstract

We consider a two-leg boson ladder in an artificial U(1) gauge field and show that, in the presence of interleg attractive interaction, the flux induced Vortex state can be melted by dislocations. For increasing flux, instead of the Meissner to Vortex transition in the commensurate-incommensurate universality class, first an Ising transition from the Meissner state to a charge density wave takes place, then, at higher flux, the melted Vortex phase is established via a disorder point where incommensuration develops in the rung current correlation function and in momentum distribution.Finally, the quasi-long range ordered Vortex phase is recovered for sufficiently small interaction. Our predictions for the observables, such as the spin current and the static structure factor, could be tested in current experiments with cold atoms in bosonic ladders.

Equations

$$H_{\lambda }$$

$$H_c$$

$$H_s$$

$$J_s=-it\sum _{j,\sigma }\sigma (e^{i\lambda \sigma }{b}_{j,\sigma }^{†}{b}_{j+1,\sigma }-e^{-i\lambda \sigma }{b}_{j+1,\sigma }^{†}{b}_{j,\sigma }),$$

$$n(k)=\sum _{\sigma }{n}_{\sigma }(k)=\frac{1}{L}\sum _{\sigma }\sum _{i,j}^{L-1}{e}^{ik(r_i-r_j)}⟨b_{i,\sigma }^{+}{b}_{j,\sigma }⟩$$

$$S^s(k)=\frac{1}{L}\sum _{\sigma ,{\sigma }^{\prime }}\sum _{j,i=0}^{L-1}sgn(\sigma {\sigma }^{\prime })⟨n_{j,\sigma }{n}_{i,{\sigma }^{\prime }}⟩.$$

$$S^c(k)=\sum _{j,\sigma ,{\sigma }^{\prime }}⟨n_{j,\sigma }{n}_{0,{\sigma }^{\prime }}⟩e^{-ikj}$$

$$E_{\pm }(k{)}^2$$

$$\frac{\cos 2\theta }{\pi a}=-i({\psi }_R^{†}{\psi }_L-{\psi }_L^{†}{\psi }_R),\frac{\cos 2\varphi }{\pi a}=i({\psi }_R^{†}{\psi }_L^{†}-{\psi }_L{\psi }_R),-\frac{1}{\pi }{\partial }_x\theta ={\psi }_R^{†}{\psi }_R+{\psi }_L^{†}{\psi }_L,$$

$$H_s$$

$$h=-\frac{\lambda u_s{K}_s}{a},m=-2\pi \Omega A_0^{2}a,\Delta =-\frac{\pi }{2}{U}_{\perp }(B_{1}a{)}^2.$$

$$H_s$$

$$E_{\pm }(k{)}^2=(u_{s}k{)}^2+m^2+h^2+{\Delta }^2\pm 2\sqrt{h^2(u_{s}k{)}^2+h^2{m}^2+{\Delta }^2{m}^2}.$$

$$G_{\alpha ,\beta }(x,\tau )=-⟨T_{\tau }{\zeta }_{\alpha }(x,\tau ){\zeta }_{\beta }(0,0)⟩=\frac{1}{L\beta }\sum _{k,i{\nu }_n}{G}_{\alpha \beta }(q,i{\nu }_n)e^{i(qx-{\nu }_n\tau )},$$

$$G_{\alpha \beta }(q,i{\nu }_n)=(i{\nu }_n-\mathcal{H}(q){)}_{\alpha \beta }^{-1}.$$

$$G(k,z)$$

$${\int }_{-\infty }^{+\infty }\frac{dk}{2\pi }{e}^{ikx}{\int }_{-\infty }^{+\infty }\frac{d\nu }{2\pi }G(k,i\nu )=-{\int }_{-\infty }^{+\infty }\frac{dk}{2\pi }{e}^{ikx}\left[\frac{H(k)}{2[E_{+}(k)+E_{-}(k)]}+\frac{({\sigma }_1\otimes \text{openone})H(k{)}^2({\sigma }_1\otimes \text{openone})H(k)}{2E_{+}(k)E_{-}(k)[E_{+}(k)+E_{-}(k)]}\right].$$

$$I_1(x)$$

$$I_2(x)$$

$$⟨J_s⟩=\frac{1}{L}⟨\frac{\partial H}{\partial \lambda }⟩=\frac{1}{L}\frac{\partial E_{GS}}{\partial \lambda }=\frac{1}{L}\frac{\partial }{\partial \lambda }\left(-\sum _{0<k<\Lambda }{E}_{+}(k)+E_{-}(k)+L\frac{h^2}{2\pi u_s}\right),$$

$$⟨J_s⟩=\frac{h}{\pi u_s}-\sum _{r=\pm 1}{\int }_0^{\Lambda }\frac{dk}{2\pi }\frac{h+r\frac{h[(u_{s}k{)}^2+m^2]}{\sqrt{h^2(u_s^2{k}^2+m^2)+m^2{\Delta }^2}}}{\sqrt{(u_{s}k{)}^2+m^2+{\Delta }^2+h^2+2r\sqrt{h^2(u_s^2{k}^2+m^2)+m^2{\Delta }^2}}}.$$

$$⟨J_s⟩=\frac{h}{2\pi u_s}\left[3-\frac{m^2-{\Delta }^2}{2m\Delta }\ln \left|\frac{m-\Delta }{m+\Delta }\right|\right],$$

$$\frac{1}{L}{\left(\frac{\partial E_{GS}}{\partial h}\right)}_{h=m}=\frac{m}{\pi u_s}\left[1-\frac{\Delta }{m}\arctan \left(\frac{m}{\Delta }\right)\right],$$

$$⟨J_s{⟩}^{(sing.)}$$

$$\frac{2a}{u_s}\frac{\partial ⟨J_s⟩}{\partial h}$$

$$⟨\cos \sqrt{2}{\varphi }_s(x)\cos \sqrt{2}{\varphi }_s(0)⟩\propto {\left(\frac{1}{{\xi }_s}\right)}^{1/4}\frac{\cos (\pi nx)}{{\left[L\sin \left(\frac{\pi x}{L}\right)\right]}^{K_c}}$$

$$S_c(2\pi {\rho }_0)-S_c(2\pi {\rho }_0+q)$$

$$⟨\cos \sqrt{2}{\varphi }_s(x)\cos \sqrt{2}{\varphi }_s(0)⟩={\left(\frac{\alpha }{{\xi }_s}\right)}^{1/4}\Phi \left(\frac{x}{L},\frac{L}{{\xi }_s}\right).$$

$$E_{-}(k{)}^2\simeq {\Delta }^2\left(1-\frac{m^2}{h^2}\right)+\frac{v^4}{4h^2}(k^2-k_L^2{)}^2,$$

Full text

Recently, artificial gauge fieldsruseckas05_gauge; *dalibard2011gauge and artificial spin-orbit couplinglin2011_soc; *galitski2013_soc have been achieved in cold atomic gases using Raman coupling, allowing to probe the effect of external gauge fields on interacting bosons. The analog of the Meissner to vortex (M-to-V) phase transition for superconductorstinkham_book_superconductors was predicted for the bosonic two-leg ladder in Refs. kardar_josephson_ladder; *orignac01_meissner. The original proposal was made in the context of Josephson junction ladders, where dissipation spoils quantum coherence fazio_josephson_junction_review, affecting their use for superconducting qubits based circuitsdevoret. In the ultracold atomic gas a simple but already nontrivial realization is the bosonic two-leg ladder in artificial fluxatala2014, where the M-to-V transition was observed in non-interacting bosonic ladders at fixed flux $\pi/2$ per plaquette and varying interleg hopping.

From the theoretical point of viewkardar_josephson_ladder; *orignac01_meissner, for bosons with in-chain repulsive interactions, the M-to-V transition falls in the commensurate-incommensurate (C-IC) universality class japaridze_cic_transition; *pokrovsky_talapov_prl. Recently it was investigated by Density Matrix Renormalization Group (DMRG) and bosonization approach for hard-core bosons on a two-leg ladder as a function of fluxour_2015; piraud2014b showing that the region of stability of the M phase over the V one is largely enhanced with respect to the non–interacting casepiraud2014b. Moreover, besides the incommensuration already predicted for low-fluxkardar_josephson_ladder; *orignac01_meissner at fluxes of the order of $\pi n$, with $n$ the number of particles per rung, a second incommensuration (2-IC) in the correlation functions is induced by the interchain hoppingour_2015; our_2016. However, in statistical mechanics, it is known that transitions in C-IC universality class can be turned into different universality classes by various relevant perturbationsbohr1982; *bohr1982b; schulz83_cic_vortices; horowitz_renormalization_incommensurable; haldane83_cic, thus for a bosonic ladder in an artificial gauge field it remains a relevant question to investigate the robustness of the M-to-V phase transition.

In this Letter, we consider the effect of an interchain interaction ($U_{\perp}$) and show that it can spoil the M-to-V transition, leading to the appearance of an intermediate charge density wave phase (CDW) that can be interpreted as a melted vortex phase. The melting of vortices is accompanied by a disorder point that gives rise to incommensuration of correlation functions when the density of dislocations becomes large enough to permit a greater energy gain from the applied flux than from the pinning potential of the vortices. We recall that while in two dimensions dislocations appear at finite temperatureBerezinskii2; *kosterlitz_thouless; haldane83_cic, in one dimension even at zero temperature their formation can be driven by quantum fluctuations only.

We consider a two-leg hard-core boson ladder in a flux $\lambda$, with Hamiltonian:

[TABLE]

where $b_{j,\sigma}$ annihilates a boson on chain $\sigma=\pm 1/2$ at site $j$, $n_{j,\sigma}=b^{\dagger}_{j,\sigma}b_{j,\sigma}$ is the associated number operator, $te^{i\lambda\sigma}$ the hopping amplitude along the chains, $\Omega$ the rung hopping, and $U_{\perp}$ the interchain interaction. This Hamiltonian can be mapped onto a system of spin-1/2 bosons with spin-orbit coupling in a transverse magnetic field with each spinor state corresponding to one leg of the ladder. In the rest of this Letter, we will consider the attractive case ($U_{\perp}<0$).

When $\Omega$ and $U_{\perp}=0$, the two chains are decoupled and their low energy properties are described by the Tomonaga-Luttinger liquid modelefetov_coupled_bosons; *haldane_bosons. For weak interchain couplings a low-energy description can be derivedkardar_josephson_ladder; *orignac01_meissner leading to a Hamiltonian $H_{\lambda}=H_{c}+H_{s}$ where:

[TABLE]

The gapless Hamiltonian $H_{c}$ controls the fluctuations of the total particle density (charge) while the Hamiltonian $H_{s}$ controls those of the difference of density between the chains (spin). In Eq. (2)–(Vortex lattice melting in a boson-ladder in artificial gauge field), $A_{0},B_{1}$ are dimensionless non-universal constants, $u_{c}K_{c}=u_{s}K_{s}=2t\sin(\pi n/2)$, $K_{c}=[1+U_{\perp}/(2\pi t\sin(\pi n/2))]^{-1/2}$, $K_{s}=[1-U_{\perp}/(2\pi t\sin(\pi n/2))]^{-1/2}$ and we have assumed an incommensurate filling of $n<1$ particles per rung. We perform DMRGwhite_dmrg; *schollwock2005 simulations for this system as a function of flux and interaction between the chains for selected fillings and values of interchain hopping $\Omega$. We show results for $n=0.5$ and fixed interchain hopping $\Omega/t=0.5$ that are summarized in the phase diagram of Fig. 1. Simulations are performed for sizes up to $L=64$ in Periodic Boundary Conditions (PBC) keeping up to $m=1256$ states during the renormalization procedure. The truncation error, that is the weight of the discarded states, is at most of order $10^{-6}$, while the error on the ground-state energy is of order $5\times 10^{-5}$ at most.

In absence of interchain interaction and in applied flux, for moderatepiraud2014b; our_2015 interchain hopping $\Omega$, a commensurate-incommensurate (C-IC) transitionjaparidze_cic_transition; *pokrovsky_talapov_prl; schulz_cic2d occurs as a function of $\lambda$kardar_josephson_ladder; *orignac01_meissner; tokuno2014. In the commensurate phase, the Meissner state, the current flowing along the rungs $J_{r}(l)=i\Omega(b_{l,\uparrow}^{\dagger}b_{l,\downarrow}-b_{l,\downarrow}^{\dagger}b_{l,\uparrow})$ has exponentially decaying correlations $\langle J_{r}(l)J_{r}(0)\rangle\sim e^{-|l|/\xi}$, while the expectation value of the spin current:

[TABLE]

i.e the difference between the currents flowing along the chains, increases linearly with the flux. The occurrence of the incommensurate phase, i.e. the establishment of a vortex phase with quasi-long range order (QLRO)cha2011; piraud2014b; greschner2015; greschner2016; didio2015a, is signalled by a rapid drop in $\langle J_{s}\rangle$, the simultaneous appearance of two separate peaks at $k=\pm q(\lambda)/2$ in the momentum distribution

[TABLE]

and two peaks at $k=\pm q(\lambda)$ in the Fourier Transform (FT) $C(k)=\sum_{l}e^{-ikl}\langle J_{r}(l)J_{r}(0)\rangle$ of the rung current correlation function. In the vicinity of $\lambda=\pi n$, it is possible to observe, for sufficiently large $\Omega$, the occurrence of a 2-ICour_2015; our_2016, characterized by two satellite peaks in $n(k)$ and $C(k)$, together with two peaks at $k=\pm\pi n$ in the spin static structure factor:

[TABLE]

When the interchain interaction $U_{\perp}<0$ is switched on, the last two terms in (Vortex lattice melting in a boson-ladder in artificial gauge field) describe the competition between the Meissner and the in-phase CDWlecheminant2002sdsg as a function of $U_{\perp}/\Omega$. When $|U_{\perp}/\Omega|$ is small enough, a Meissner phase at $\lambda=0$ is obtained, and the application of a flux is expected to induce a C-IC transition. However, at the C-IC transition point, the scaling dimensionschulz_cic2d; chitra_spinchains_field of the $\cos\sqrt{2}\theta_{s}$ operator is $1/(2K_{s}^{*})=1$ implying that the operator $\cos\sqrt{8}\phi_{s}$ of dimension $2K_{s}^{*}$ is relevanthorowitz_renormalization_incommensurable; haldane83_cic as well. When increasing the flux, a gapped in-phase CDW phasemathey_pra_2009; *hu_pra_2009 with $\langle e^{i\sqrt{2}\phi_{s}}\rangle\neq 0$ is formed instead of a vortex state with quasi-long range order. By duality, both the boson annihilation operators $b_{j\sigma}\sim A_{0}e^{i\frac{(\theta_{c}+2\sigma\theta_{s})}{\sqrt{2}}}$ and the rung current operator $J_{r}=\Omega\sin\sqrt{2}\theta_{s}$ present exponentially decaying correlations, while the height of $n(k)$ at $k=0$ saturates with size because of the SRO of $e^{i\theta_{s}/\sqrt{2}}$ . Since the location of the peak in $n(k)$ remains $k=0$, this indicates that CDW retains the commensuration of the Meissner state and thus we call it Meissner-CDW (M-CDW) phase. In the statistical mechanics contextschulz83_cic_vortices; haldane83_cic, the commensurate phase is our Meissner phase, the incommensurate phase is the Vortex phase and the liquid phase is the M-CDW phase. This last phase can be detected via the (in-phase) charge density structure factor

[TABLE]

which develops peaks at $k=\pm\pi n$ whose heights don’t scale with $L$, since the Luttinger exponent $K_{c}>1$ (see supplementary). In the lower panels $A$ and $B$ of Fig. 2 we follow this transition: in the Meisnner phase (panel $A$) $S^{c}(k)$ is smooth and $n(k)$ shows a power-law divergence, while in the M-CDW phase (panel $B$) $S^{c}(k)$ acquires the above mentioned peaks and $n(k)$ shows a Lorentzian-like peak at $k=0$.

Meanwhile, $C(k)$ and $S^{s}(k)$ retain a Lorentzian shape on both sides of the transition. For sufficiently large interchain interaction the Meissner phase is replaced by the M-CDW even at zero flux, as shown in Fig. 1. Let us note that the opening of a gap in the spin-sector prevents the observation of the 2-ICour_2015; our_2016 at increasing flux for attractive interchain interation.

The universality class of the flux-driven Meissner to M-CDW transition can be obtainedbohr1982; *bohr1982b by fermionizing the Hamiltonian (Vortex lattice melting in a boson-ladder in artificial gauge field) in the vicinity of the C-IC transition, as done in Refs.wang2003field; tsvelik_field; *essler04_spin1_field; *citro02_dm_ladders for anisotropic spin-1 chains and spin-1/2 ladders. The corresponding Majorana fermion Hamiltonian supplementary has dispersionbohr1982; *bohr1982b; wang2003field:

[TABLE]

where $m=-2\pi\Omega A_{0}^{2}a$, $\Delta=-\frac{\pi}{2}V_{\perp}(B_{1}a)^{2}$ and $h=-\frac{\lambda u_{s}K_{s}}{a}$. For $h_{eff}=\sqrt{h^{2}+\Delta^{2}}=|m|$, the $-$ branch of Majorana fermions dispersion becomes gapless showing that the M-CDW quantum phase transition falls in the Ising universality class tsvelik_field; mccoy_revue_qft. For $h_{eff}<|m|$ $e^{i\theta_{s}/\sqrt{2}}$ is long range ordered, while it is short range ordered for $h_{eff}>|m|$.

At the Ising transition, the derivative of the spin current $J_{s}$, $\partial_{\lambda}J_{s}\sim-\ln[\mathrm{max}(|\lambda-\lambda_{c}|,1/L)]$ diverges logarithmically like the specific heat of the two-dimensional Ising modelmccoy_revue_qft. The absence of the square root threshold singularityjaparidze_cic_transition; *pokrovsky_talapov_prl characteristic of the C-IC transition, $\langle J_{s}(\lambda)-J_{s}(\lambda_{c})\rangle\sim\sqrt{\lambda-\lambda_{c}}\theta(\lambda-\lambda_{c})$, can be noted in the upper panel in Fig. 2. By plotting the numerical derivative of $\langle J_{s}\rangle$ a narrow peak can be spotted at $\lambda\simeq 0.15\pi$, past the maximum in the spin current (see supplemental materialsupplementary).

Another indicator of the nature of the transition is the von Neumann entropy given by $S_{vN}=c/3\ln\left[\frac{L}{\pi}\sin\left(\pi\frac{x}{L}\right)\right]+g$ for a system with PBCcalabrese04_entanglement where $c=c_{c}+c_{s}$ is the sum of the central chargesdifrancesco_book_conformal of the respectively charge and spin gapless modes and $g$ is a non-universal constant. In the Fig. 3 we show the extrapolated $c$ from fit to the numerical data of $S_{vN}$. Despite the size effects, we can observe a bell shaped curve centered around the critical value $\lambda_{cdw}$ the width of which gets smaller with increasing system size. The height of this peak extrapolates to $c=c_{c}+c_{s}=3/2$ as size increases, indicating that the critical point belongs to the Ising universality class (see supplementarysupplementary), while far from $\lambda_{cdw}$, $c=c_{c}$ extrapolates to unity. Finally, the Ising nature of the transition can also be spotted looking at the value of the peaks in the $S_{c}(k=\pi/2)$ that not too close to transition should be proportional to $(\lambda-\lambda_{cdw})^{1/4}$. We have verified this behavior for the case reported in Fig. 3 (see supplemental materialsupplementary).

Besides the Ising transition, the fermionized Hamiltonian also predictsschulz83_cic_vortices the existence of a disorder pointstephenson1970a in the crossover to the melted-Vortex phase. Indeed, beyond a critical value $\lambda=\lambda_{d}$, real space correlations function $\langle J_{r}(l)J_{r}(0)\rangle$ and $\langle b^{\dagger}_{l\sigma}b_{0\sigma}\rangle$ both acquire a periodic modulationwang2003field. Since the wavevector of the modulation is varying with $\lambda$, the disorder point is of the first kind.stephenson1970b. In reciprocal space, the modulation gives rise to a superposition of two Lorentzian-like peaks in $C(k)$ and $n(k)$ that are remnants of the divergent peakscha2011; our_2015 previously obtained in the QLRO vortex state when $U_{\perp}=0$ (see panel $C$ in Fig.2 and panels $A$ and $B$ in Fig. 4). The values of $\lambda$ for which the two peaks in $C(k)$ can be resolved are at $\lambda>\lambda_{L,C}$ with $\lambda_{L,C}$ the Lifshitz point. One has $\lambda_{L,C}>\lambda_{d}$ since resolving the peaks requires that the distance between their maxima exceeds their width. This effect is less evident in the spin resolved $n_{\sigma}(k)$ where a single Lorentzian-like peaks located at finite $k$ develops (black solid line in Fig. 1). Similar disorder and Lifshitz points have been found in one-dimensional spin-1/2bursill1995; *deschner2013 and spin-1schollwoeck1996; *pixley2014; *chepiga2016 chains as well as frustrated Ising chains in transverse fieldbeccaria2006. In this phase the spin and charge response functions retain respectively the Lorentzian shape centered around $k=0$ and the peaks at $k=\pm\pi n$. Finally, for even higher flux and sufficiently small interaction, $K_{s}^{*}$ becomes greater than $1$ and the operator $\cos\sqrt{8}\phi_{s}$ becomes irrelevant, allowing a vortex phasekardar_josephson_ladder; *orignac01_meissner with QLRO. This takes place through a Berezinskii-Kosterlitz-Thouless (BKT) transitionBerezinskii2; *kosterlitz_thouless at the point $K_{s}^{*}(\lambda_{BKT})=1$. In Fig. 4 we show the recover of the QLRO vortex phase (panel $C$ from Fig.4) from the melted one by decreasing the strenght of interchain interaction $U_{\perp}/t$ for fixed applied flux $\lambda=0.625\pi$.

To conclude, using bosonization and DMRG we have found (see Fig. 1) that, with an interchain attractive interaction, the commensurate Meissner phase and the incommensurate QLRO vortex phase leave space to a Meissner-CDW and to a melted vortex phase with SRO. Instead of having a single flux-driven M-to-V transition we are left with an Ising transition to the commensurate Meissner-CDW. On increasing the flux an exponentially damped sinusoidal modulation, incommensurate with the ladder, develops in the momentum distribution. At the Lifshitz point $\lambda=\lambda_{L,C}$ double peaks appear in the rung current structure factor. This indicates the existence of a disorder pointstephenson1970a; stephenson1970b where the bosonic Green’s functions and the rung current correlation function develop exponentially damped oscillations in real space. The Meissner-CDW is then crossing over into a melted vortex phase where a SRO with proliferation of dislocations takes over. At higher flux, a BKT transition takes place, and the quasi-long range vortex lattice is recovered.

Our predictions on the melting of vortices in Bose-Einstein condensates in optical lattices can be traced in current experiments, where static structure factorsstenger99_bragg_bec and momentum distributions can be measured, together with the spin currentatala2014; livi. Using dipolar interactions, tunable by orienting the dipoles with a fieldkollath07_dipolar; *rydberg_atoms, or Feshbach resonances, the interaction $U_{\perp}$ can be rendered attractive.

The detection of a melting transition as well as the non-trivial effects due to interactions can be relevant for atomtronic ring ladders which have been proposed for readout and gate implementation in quantum technologiesAmico, analogously to the superconducting qubits in multi-junction circuits.

I Majorana Fermions representation

Fermionizationbohr1982; *bohr1982b leads to a a detailed picture of the transition between the Meissner state and the density wave states. Rescaling $\theta_{s}=\sqrt{2}\theta$ and $\phi_{s}=\phi/\sqrt{2}$, the Hamiltonian ([Insert Eq. number from manuscript]) is fermionized using the identities:

[TABLE]

yielding:

[TABLE]

where:

[TABLE]

Once we introduce the Majorana fermion operators $\zeta_{\nu,\sigma}(x)=\frac{1}{\sqrt{2}}(\psi_{\nu}+\sigma\psi^{\dagger}_{\nu})(x)$ with $\nu=R,L$, $\sigma=\pm$, the Hamiltonian (S1) is rewritten:

[TABLE]

Hamiltonians of the form (S3) have previously been studied in the context of spin-1 chains in magnetic fieldtsvelik_field; *wang2003field; *essler04_spin1_field or spin-1/2 laddersshelton_spin_ladders; *nersesyan_biquad with anisotropic interactionscitro02_dm_ladders. The quadratic Hamiltonian (S3) can be diagonalized in momentum space to obtain the following energy eigenvalues:

[TABLE]

Since observables such as the rung current are bilinear in the Majorana fermions, obtaining their correlation functions requires knowledge of the Matsubara Green’s function for the Majorana fermions defined asabrikosov_book:

[TABLE]

where $\alpha=(\nu,j)$ with $\nu=R,L$ and $j=1,2$. In Eq.(S5), the Fourier space Green’s function is the $4\times 4$ matrix:

[TABLE]

An explicit expression in terms of Pauli matrices is:

[TABLE]

where $E_{\pm}(k)$ have been defined above in (S4).

The real space Green’s function, at equal time and zero temperature, is defined by:

[TABLE]

Thus, it can be obtained from just two integrals:

[TABLE]

by taking the appropriate number of derivatives with respect to $x$.

I.1 Ising transition

For $m=\sqrt{h^{2}+\Delta^{2}}$, $E_{-}(k)=u\frac{\Delta}{m}|k|+O(k^{2})$ a single Majorana fermion mode becomes massless at the transitionbohr1982; *bohr1982b between the Meissner and the density wave state. This transition belongs to the Isingmccoy_revue_qft universality class.

As a consequence, at the transition, the Von Neumann entanglement entropy $S_{vN}=\frac{1}{3}c\ln L=\frac{1}{3}(c_{c}+c_{s})\ln L=\frac{1}{3}(1+1/2)\ln L$, where $c$ is the central charge of the gapless modes, while away from the transition it is $S_{vN}=\frac{1}{3}c_{c}\ln L=\frac{1}{3}\ln L$ since the total density mode $\phi_{c}$ is gapless. In Fig. S1 we show size extrapolation for the central charge of two cases, at $n=0.5$, $U_{\perp}/t=-1.5$, $\Omega/t=0.5$, namely at $\lambda=0.15\pi$ where we expect the Ising transition for which we get $c=1.6(1)$ and at $\lambda=0.75\pi$ in the melted Vortex phase where $c=1.0(1)$. Since bosonization predicts that the charge mode is gapless and it is described at low energy by a free boson with $c_{c}=1$, the charge mode exhausts the central charge when $\lambda\neq\lambda_{cdw}$. At $\lambda_{cdw}$, the presence of the peak indicates the appearance of an extra gapless field with $c_{s}<1$, i. e. a critical point. Conformal Field Theories with $c_{s}<1$ form a discrete seriesdifrancesco_book_conformal with $c_{s}=1-\frac{6}{m(m+1)}$ with the integer $m\geq 3$. The size-extrapolated value of central charge of the critical theory appears to be under $7/10$, leaving us with $m=3$ and $c_{s}=1/2$ as the only possible value, indicating that the critical point belongs to the Ising universality class.

In Fig. S1 we show size extrapolation for the central charge of two cases, at $n=0.5$, $U_{\perp}/t=-1.5$, $\Omega/t=0.5$, namely at $\lambda=0.15\pi$ where we expect the Ising transition for which we get $c=1.6(1)$ and at $\lambda=0.75\pi$ in the melted Vortex phase where $c=1.0(1)$.

I.2 Current

A signature of the Ising transition can be observed also in the spin current, which is defined by

[TABLE]

from which the $\langle J_{s}\rangle$ in units of $\frac{u_{s}}{2a}$ is found:

[TABLE]

The integral (S11) is convergent in the limit $\Lambda\to+\infty$.

In the limit $\lambda\ll 1$, or $h\ll 1$,

[TABLE]

As $\Delta/m$ increases, the proportionality constant between the flux and the current decreases but remains positive, indicating that the Meissner effect is reduced by interchain repulsion. Right at $h=m$, we can obtain the exact expression of $\partial_{h}E_{GS}$ as:

[TABLE]

which shows that the current at $h=m$, i. e. where the commensurate-incommensurate transition would take place in the absence of repulsion, is always reduced by interchain interaction.

In the vicinity of $h=\sqrt{m^{2}-\Delta^{2}}$, the dominant singularity in $\langle J_{s}\rangle$ is:

[TABLE]

Thus, while $\langle J_{s}\rangle$ is continuous for $h\to\sqrt{m^{2}-\Delta^{2}}$, its derivative instead:

[TABLE]

diverges as $\frac{h}{4\pi\Delta a}\ln|h-h_{c}|$ for $h^{2}=m^{2}+\Delta^{2}$. $\langle J_{s}\rangle$ is a decreasing function of $h$, and its plot as a function of $h$ presents a vertical tangent at $h_{c}$. As $\langle J_{s}\rangle$ is increasing at small $h$ and decreasing for $h=\sqrt{m^{2}-\Delta^{2}}$, a maximum of the current must exist in the range $0<h_{\mathrm{max.}}<\sqrt{m^{2}-\Delta^{2}}$ i. e. inside the Meissner phase.

I.3 Charge density wave order parameter

The CDW order is related to the $2k_{F}$ wave-vector component of the density operator that, expressed using bosonization, is $O_{CDW}(x)=Ae^{i\sqrt{2}\phi_{c}}\cos\sqrt{2}\phi_{s}$ where $A$ is a non-universal constant dependent on the microscopic Hamiltonian. On the Meissner side, $\theta_{s}$ is long range ordered, so the CDW correlations decay exponentially while, on the density wave side and for $U<0$, $\langle\phi_{s}\rangle=0$ and $O_{CDW}(x)\simeq\bar{A}\langle\cos\sqrt{2}\phi_{s}\rangle e^{i\sqrt{2}\phi_{x}}$ where $\langle\cos\sqrt{2}\phi_{s}\rangle\sim(a/\xi_{s})^{1/8}$ with $\xi_{s}\sim|\lambda-\lambda_{cdw}|$ the critical correlation length of the Meissner-CDW transition. The exponent $1/8$ is the result of $\cos\sqrt{2}\phi_{s}$ being proportional to the order parameter of the Ising transition. For finite size and periodic boundary conditions the CDW correlations takes the form:

[TABLE]

and right at the Ising transition it becomes $\propto\left(\frac{\pi\alpha}{L\sin\left(\frac{\pi x}{L}\right)}\right)^{1/4}$ and hence that $K_{c}\to K_{c}+1/4$ in Eq. (S13).

Taking the Fourier transform with $2>K_{c}>1$ at $k=2\pi\rho_{0}$ and at finite $q$ nearby, we get:

[TABLE]

When $qL\gg 1$, one has $S_{c}(2\pi\rho_{0})-S_{c}(2\pi\rho_{0}+q)\sim|q\alpha|^{K_{c}-1}$. In the limit $L\to+\infty$, if $K_{c}>1$, $S_{c}(k)$ is finite for $k\to 2\pi\rho_{0}$ but presents a cusp at that point. This is consistent with the low energy predictions derived from bosonization where $K_{c}=[1+U_{\perp}/(2\pi t\sin(\pi n/2))]^{-1/2}$.

In the transition region between the CDW and the Ising critical point, where the correlation length is finite but comparable to $L$, we can write:

[TABLE]

For given $L$ and $\xi_{s}$, we have to distinguish two regimes in $q$ when we are close to the transition. For $q\xi_{s}\ll 1$, we recover the regime (S14). For $q\xi_{s}\gg 1$, the integral in (S15) becomes sensitive to the short distance physics, and we expect: $S_{c}(2\pi\rho_{0})-S_{c}(q+2\pi\rho_{0})\sim(q\alpha)^{K_{c}-3/4}$. To sum up, far from the transition, we can treat the Ising order parameter as a constant and obtain a scaling function that depends only on $qL$ and $K_{c}$ while near the transition (with $\xi_{s}/L$ non-negligible) the scaling function depends on $K_{c}$ and both $qL$ and $q\xi_{s}$.

In Fig. S2 we show $S_{c}(2\pi\rho_{0})$ as a function of the applied flux for the case $U_{perp}/t=-1.5$ at $n=0.5$ and $\Omega/t=0.5$, far from transition we can detect a region on both sides of the Ising transition where this quantity is $\propto\xi_{s}^{-1/4}$. A fit using $f(\lambda)=a/(\lambda-\lambda_{cdw})^{1/4}+b$ on both sides of the transition gives us as results for $\lambda_{cdw}=0.150(2)$ and $0.16(2)$, in agreement with the estimated extracted from the central central charge or from the derivative of the spin current.

I.4 Lifshitz and disorder point

From the spectrum of the Majorana fermion representation (S4) we can deduce the existence of a disorder and a Lifshitz point in some correlation functions. At a disorder point, a real space correlation function acquires a periodic modulation as a function of $x$. In reciprocal space, its Fourier transform is a the sum of two Lorentzian-shaped peaks. The value of $\lambda$ where the two peaks are resolved, i. e. where the modulation wavevector is of the order of the peak width, is the Lifshitz point. It does not coincide with the disorder point because of the short range order. This definition of disorder and Lifshitz points applies to all correlation functions, but the existence of these points can be inferred by considering the single particle Green’s function of the Majorana Fermions.

The energy $E_{+}(k)$ is always an increasing function of $k$, while for $2h^{2}>2h_{L}^{2}=m^{2}+\sqrt{m^{4}+4m^{2}\Delta^{2}}$, $E_{-}(k)$ has two degenerate minimawang2003field for $k=\pm k_{L}$ with $k_{L}^{(\pm)}=\pm\frac{\sqrt{h^{2}-m^{2}-\Delta^{2}m^{2}/h^{2}}}{u_{s}}$. For $h>h_{L}$, $E_{-}(k)$ can be Taylor expanded near $k=\pm k_{L}$ as:

[TABLE]

and inserting in the expression of $I_{1}(x)$, a sinusoidal modulation of the Green’s function of the Majorana fermions is expected wang2003field at least when $h>h_{L}$. A more detailed calculation would show that the sinusoidal modulation of $I_{1}(x)$ appears when $E_{-}(k)$ develops two disconnected branch cuts symmetric with respect to the imaginary axis in complex $k$ plane. Now, the correlation function $\langle J_{\perp}(x)J_{\perp}(0)\rangle$ is the trace of a product of two Majorana Fermion Green’s function and Pauli matrices, and thus also exhibits a sinusoidal modulation of wavevector $2k_{L}$. So a disorder point is present in the correlator $\langle J_{\perp}(x)J_{\perp}(0)\rangle$, and a Lifshitz point is expected in the Fourier transform $C(k)$. In the case of the real space Green’s function of the original bosons, it is known that it depends on a factor coming from the charge modes and a factor coming from Green’s functions of Ising order and disorder operatorswang2003field associated with the Majorana fermions. The latter factorwu_ising_correlations; wang2003field can be expressed in terms of block Toeplitz determinants the elements of which are Majorana fermion Green’s functions. Numerical calculationswang2003field show that a sinusoisal modulation appears also in the Ising order and disorder parameter correlations. This indicates that a disorder point is also present in the correlators $\langle b_{j\sigma}b^{\dagger}_{0\sigma}\rangle$. Considering the Fourier transform, a Lifshitz point is expected in $n(k)$. Since the disorder point is associated with the appearance of a modulation in the Green’s function ogf the Majorana fermions, it is expected that$\langle b_{j\sigma}b^{\dagger}_{0\sigma}\rangle$ and $\langle J_{\perp}(x)J_{\perp}(0)\rangle$ have a disorder point at the same value of $\lambda$. However, in general, their Lifshitz points are not expected to coincide since $\langle b_{j\sigma}b^{\dagger}_{0\sigma}\rangle$ depends on both “charge” and “spin” modes and since the correlation lengths of the two “spin” parts can differ.