Topological superfluidity with repulsive alkaline-earth atoms in optical lattices

Topological superfluids are of technological relevance since they are believed to host Majorana bound states, a powerful resource for quantum computation and memory. Here we propose to realize topological superfluidity with fermionic atoms in an optical lattice. We consider a situation where atoms in two internal states experience different lattice potentials: one species is localized and the other itinerant, and show how quantum fluctuations of the localized fermions give rise to an attraction and strong spin-orbit coupling in the itinerant band. At low temperature, these effects stabilize a topological superfluid of mobile atoms even if their bare interactions are repulsive. This emergent state can be engineered with ${}^{87}$Sr atoms in a superlattice with a dimerized unit cell. To probe its unique properties we describe protocols that use high spectral resolution and controllability of the Sr clock transition, such as momentum-resolved spectroscopy and supercurrent response to a synthetic (laser-induced) magnetic field.

Introduction.-Our understanding of many-body systems traditionally relies on the Landau classification of ordered states of matter based on global symmetries spontaneously broken within a given phase. This symmetry breaking is accompanied by emergence of an orderparameter (OP), i.e. non-zero expectation value of a local physical observable that uniquely characterizes the phase. For instance a hallmark signature of a fermionic superfluid (SF) is breaking of the particle number conservation [U (1)] symmetry which occurs as a result of Cooper pairing. The corresponding OP plays the role of a Cooper pair wavefunction and defines an energy gap in the excitation spectrum, allowing dissipationless particle currents [1]. However, many phases of matter defy the Landau paradigm. An important class of such systems are topological superfluids (TSFs) [2,3], i.e. phases that in addition to U (1), break a residual Z 2 symmetry. The latter symmetry breaking is a global phenomenon that occurs in the absence of local OP and only for appropriate boundary conditions [4][5][6].
Despite all efforts dedicated to the search for TSFs, they remain elusive with the only confirmed realization being liquid 3 He-A [7]. One reason for such scarcity is that TSFs require a very particular orbital structure of Cooper pairs, at least p-wave, which originates either from strongly spin-dependent interactions or a large spinorbit coupling (SOC) that couples particle's motion to its spin. Apparently, the coexistence of a sizeable SOC and attractive interactions (leading to Cooper pairing) is quite rare in nature [8], fundamentally because SOC and fermion pairing have very different physical origins.
In the present work, we propose a pathway towards topological superfluidity, which completely overcomes the above limitation by using the same ingredients to engineer attractive interactions and SOC. We study a model of repulsive fermions in two bands: one localized and another itinerant, and show that inhomogeneities span-ning few lattice sites (e.g. dimerization) in the localized band lead to two profound phenomena. First, they induce short-range attractive interactions among the itinerant species, by virtue of local quantum fluctuations. Second, they enlarge the unit cell in accordance with the extent of localized wavefunctions. As a result, itinerant states reside in several Bloch bands, whose index plays the role of a spin degree of freedom. These pseudospins flip whenever a fermion tunnels between unit cells, thus coupling to the orbital motion. The strength of this emergent SOC is determined by the tunneling amplitude and is of the order of the itinerant bandwidth. We show that a combination of this ultra-strong SOC and attractive interactions gives rise to a robust p-wave TSF in quasi-one dimension (1D) and a chiral p x + ip y SF in 2D.
Our TSF state can be observed in ultracold nuclearspin polarized fermionic alkaline-earth atoms (AEAs) [9], e.g. 87 Sr [10] or 173 Yb [11,12], in an optical superlattice with a few-site unit cell [13][14][15][16]. The localized (itinerant) states can be implemented with atoms in an excited 3 P 0 (ground state 1 S 0 ) clock state (respectively, e-and g-states), with a single e-atom per unit cell. We propose several experimental probes for characterizing the TSFs, including momentum-resolved spectroscopy [17,18] and generation of a particle supercurrent with a laser-induced synthetic magnetic field [19][20][21]. Our approach avoids many known experimental issues: (i) the only relevant interactions occur through the a − eg channel, and therefore the system is not affected by inelastic e-e losses [22,23] or strong scattering in the a + eg channel [24,25]; (ii) p-wave interactions in our case emerge as a result of quantum fluctuations as opposed to a p-wave Feshbach resonance, and our setup is free from the three-body losses reported in experiments [26][27][28][29]; (iii) the SOC in our system is generated as a result of the lattice structure and hence avoids heating, inherent to earlier proposals to create SOC using near-resonant Raman lasers [30][31][32][33]. Our proposal (a) 2J e λ = +1  1) can be implemented by tightly confining in an array of 1D tubes an ultra-cold gas of nuclear-spin polarized fermionic alkalineearth atoms prepared in the clock states g (in blue) and e (red color). Along the tubes, e-atoms experience a superlattice that consists of weakly-coupled double-wells (dimers) with large intra-dimer tunneling Je. The g-atoms are itinerant and experience a weaker lattice potential along the tube direction with a nearest-neighbor hopping Jg. Within each tube, a unit cell (dashed rectangle) at a position xi = i includes two lattice sites labeled with a = 1, 2 (4 wells overall). There is a e-g repulsive interaction U − eg > 0, assumed small compared to Je: U − eg Je. (b) Symmetric (λ = +) and antisymmetric (λ = −) e-atom kinetic-energy eigenstates within a dimer. (c) When e-atom dimers are prepared in the antisymmetric mode, virtual transitions to the symmetric state, caused by the e-g interaction, induce an effective attraction ugg = (U − eg ) 2 /4Je between two g-fermions within a dimer. These processes are captured by the effective model (2). is much simpler than previous works that involve more complicated laser arrays, RF pulses or additional molecular states [34][35][36][37][38][39]. Finally, our cold-atom system provides a long-sought-after realization of a pairing mechanism in repulsive fermions, which emerges because of nanoscale inhomogeneities [40][41][42][43].
p-wave SF in a 1D superlattice.-Key aspects of the emergent Cooper pairing and SOC leading to our proposed TSF state can be seen by studying a dimerized 1D optical lattice, shown in Fig. 1(a) and described by the model Hamiltonian: where i = x i = 0, . . . , N d − 1 and a = 1, 2 labels dimers and sites within a dimer, respectively. The operatorê † ia (ĝ † ia ) creates a nuclear-spin polarized e (g) atom at site a within a dimer i (n e ia =ê † iaê ia and similarly forn g ia ). The e-atoms occupy a dimerized lattice with a large intradimer hopping J e > 0 and one atom per dimer (we assume that dimers are decoupled). The g-atoms propagate in a simple (non-dimerized) lattice with a nearestneighbor tunneling J g . The second term in (1) contains a local e-g repulsion of strength U − eg > 0. We focus on the regime J e U − eg and J g when interactions and g-atom kinetic energy in (1) can be considered a perturbation to the e-atom kinetic energy. For i-th dimer, the latter has eigenstates |λ i = 1 and |vac is the vacuum state without atoms) with energies −λJ e [ Fig. 1(b)]. States of the entire e-g system can be approximately written as |Ψ eg = i |λ i ⊗ |Ψ g (|Ψ g is a state of only g-atoms), thanks to the single-dimer gap 2J e . We next assume that e-subsystem is prepared in the excited state i |λ = −1 i . This configuration is stable because of the large energy penalty 2J e that suppresses decay of individual dimers to their ground state (GS) with λ = +1 in the absence of decoherence sources (this requirement is well satisfied in cold-atom systems), for instance due to e-g scattering. The weak interactions U − eg only induce e-atom virtual transitions to dimer states with λ = +1, which we take into account via 2nd order perturbation theory (the kinetic energy of g-atoms amounts to a 1st order correction because it operates within the degenerate subspace {|Ψ eg }). These virtual processes, shown in Fig. 1(c), give rise to an effective Hamiltonian for the g-subsystem (see Methods) , σ are Pauli matrices, momentum k ∈ [−π, π] (in units of inverse lattice spacing 1/a 0 ) is defined in a dimer Brillouin zone (BZ) with N d states. u gg = (U − eg ) 2 /4J e is the strength of intra-dimer g-atom attraction mediated by quantum fluctuations of localized e-atoms. If we associate the site index a = 1, 2 inside a dimer with a spin-1 2 degree of freedom,Ĥ g 0 contains kinetic energy with a SOC that arises because any tunneling event "flips" pseudospin a.
The physical origin of the p-wave TSF is especially transparent at weak coupling u gg J g and low filling n g 1, when kinetic energy in (2) dominates and is diagonalized by statesf kτ = 1 k2 with energies k = 2τ J g cos k 2 (τ = ±1). Because relevant momenta are small, |k| π, it is allowed to keep only thef k,−1 ≡f k mode. As a result, interactions in (2) become manifestly p-wave:  (a) Zero-temperature phase diagram of Eq. (2) computed using an unconstrained Hartree-Fock-Bogoliubov mean-field theory in a system with N d = 100 dimers and periodic boundary conditions. µ is the g-atom chemical potential. Thick [thin] lines indicate 1st order transitions between topological superfluid (TSF) and insulator states [2nd order transition inside the insulating region]. In the chargedensity wave (CDW) state the unit cell has two dimers with an average density n g = 1 atom per dimer. For small ugg the CDW undergoes a transition to a band insulator with a single-dimer unit cell. Blue square (red star) corresponds to ugg/Jg = 2.7 (2.3). (b) SF gap ∆ and average density n g = 1 N d i (n g i1 + n g i2 ) plotted along the arrows shown in (a). Within the Bogoliubov mean-field theory [44] one introduces a pairing OP ∆ = − ugg 4N d q e −i q 2 f −qfq which parameterizes the single-particle excitation spectrum E k = ( k − µ) 2 + |D k | 2 with a p-wave gap D k ≈ i k∆. µ is the g-atom chemical potential (set by the Fermi energy). ∆ is a solution of a self-consistency equation Stability and topological nature of the SF state.-To assess the stability of the p-wave SF state beyond the weak coupling limit, and uncover its topological properties, we compute the phase diagram ofĤ ef within a fully unconstrained Hartree-Fock-Bogoliubov meanfield approach in real space (explained in Methods). This variational technique minimizes the grand potential Ĥ ef − µ iaĝ † iaĝ ia (µ is the g-atom chemical potential) w.r.t. local OPs ∆ i = −u gg ĝ i2ĝi1 , ξ i = ĝ † i2ĝ i1 and n g ia = ĝ † iaĝ ia , and includes the competition between SF phases with a finite gap ∆ i and various in-homogeneous states, e.g. charge-density waves (CDWs), characterized by a site-dependent n g ia . The minimization is performed at zero temperature T = 0 in a system with periodic boundary conditions (BCs). Once the GS is self-consistently determined, we open the chain and diagonalize the Bogoliubov-de Gennes (BdG) mean-field Hamiltonian with fixed OPs to determine edge modes: If a SF phase displays zero-energy (Majorana) modes, we call it topological [2]. Fig. 2(a) shows the phase diagram of model (2) as a function of chemical potential µ and interaction u gg . In agreement with the previous section, the SF phase is stable at weak coupling u gg < J g and low density, and is characterized by the mixing of singlet ĝ −k,2ĝk1 and triplet ĝ −k,aĝka Cooper pair amplitudes. To better understand this effect, let us consider properties of the model (2) under space inversion I: x → −x. In the dimerized lattice, I = σ x ⊗ I d , where I d : k → −k is an inversion acting on the dimer center-of-mass and σ x appears because I must interchange dimer sites. The SOCĤ g 0 is invariant under I but manifestly breaks I d due to the odd-momentum terms. In a SF state, Cooper pair wavefunctions inherit this feature and the system exhibits p-wave pairing between same-flavor gatoms ĝ −k,aĝka ∼ k n (n is odd) despite the s-wave nature of interactions in Eq. (2). This situation is similar to singlet-triplet mixing in non-centrosymmetric superconductors with Rashba SOC [8,19,20].
When u gg or µ is increased, the system undergoes a 1st order transition to a non-SF gapped state with an average density n g = 1 [ Fig. 2(b)]. This phase is a band insulator for small u gg , and a CDW with two dimers per unit cell for strong interactions. As shown in Fig. 2(c), the SF state is topological (i.e. possesses zero-energy edge modes) for all u gg and µ where it is stable. This happens because we used the Hartree-Fock OPs in our variational scheme: if the minimization were constrained to include only site-independent ∆ i , one would recover a well-known transition [2] from TSF to a non-topological SF state. The latter phase is unstable towards CDW formation and the transition never happens.
The phase diagram in Fig. 2(a) remains valid at finite temperature T > 0. Indeed, as demonstrated in Fig.  2(d) a typical critical temperature, above which the SF phase disappears, is T c ∼ 0.1J g (here and below we use the units with Boltzmann constant k B = 1). For T > T c , the system becomes a homogeneous Fermi liquid.
Detection of the TSF state.-The simplest way to validate our theory in cold-atom experiments, is to probe the g-atom attraction via quench dynamics in a normal state using the following protocol: (i) For times t < 0, g-atoms fill a non-interacting Fermi sea. (ii) At t = 0, J g is switched off (e.g. by increasing the lattice depth), and g-atoms are brought in contact with e-atoms, thus allowing them to experience e-g interactions described in Eq. (1). Then, one lets the system evolve for a time t 0 . As a result, basis states with doubly occupied dimers accumulate a phase −U t 0 , where U [= −u gg in Eq. (2)] is the induced g-atom interaction. (iii) At t = t 0 , the e-atoms are removed, hopping J g is restored and the system evolves with a noninteracting Hamiltonian [1st term in (2)]. The sign of U can be determined by measuring an average number of doubly occupied dimers: in the state |ψ(t) of the system at time t. For short evolution times when |U |t 0 1 and Hence, the number of double occupancies decreases (increases) for the repulsive (attractive) interaction U (see Methods for a complete derivation of this result).
The spectrum of Bogoliubov excitations can be probed by momentum-resolved spectroscopy [17,18]. Let us assume that one tube in Fig. 1(a) contains no atoms and is detuned relative to its neighboring tubes. A laser with a wavevector along the x-axis, Rabi frequency Ω, detuned (2) with a laser-assisted intra-dimer tunneling ±i by simulating a magnetic field B = byey. A supercurrent K ∼ by is induced in the TSF phase. by δ from the atomic e-g transition, transfers g-atoms from the SF phase to e-states in the empty tube (the transfer happens along the y-axis and does not change an atom's x-position along the tube). The e-lattice depth in that tube has been reduced to make flat bands at ±J e [see Fig. 1(b)] dispersive with a band structure e kτ = τ J e 1 + η 2 + 2η cos k where τ = ±1 and ηJ e is the inter-dimer hopping. The g-atom transfer rate to an e-band e kτ , R τ (δ, k), can be written in terms of the spectral density A ab (ω, k) of the single-particle normal Green function [18]: is the Fermi function at temperature T , and i0 is an infinitesimal imaginary number (see Methods).
The representative signal R is shown in Fig. 3. Its maximum (for a fixed k) occurs when ω coincides with the highest occupied BdG energy state. Therefore, one can map out the BdG band structure and extract the excitation gap ∆ [see panels (a) and (d)]. In Fig. 3(b) and (e) we plot the same signal as a function of the bare laser detuning δ, as it would be observed in an experiment. This spectroscopy technique can also be used to probe the insulating phase in Fig. 1(a). The transfer rate inside the band insulator regime is shown in Fig. 3(c) and (f). In this case, the excitation spectrum is again gapped, but as opposed to the SF state, this gap exists because all single-particle states below the Fermi level are filled, and not due to fermion pairing. Finally, we note that the signal with τ = −1 [panels (a) -(c)] is significantly stronger than the one with τ = 1 [(d) -(f)], which highlights the validity of the effective model (3).
Due to the SOC inherent in Eq. (2), the TSF state has A protocol to prepare the model in Fig. 1(a) and Eq. (1). Blue (red) color indicates g (e) atoms. Atoms are nuclear-spin polarized and hoppings are quenched until step (5). |e(g) n means a state with one e (g) atom in the n-th spatial level (n = 0 means GS). At steps (2) and (3), red-blue (blue-red) circles [from bottom to top] correspond to |e ± g states detuned from the lowest-band energy by ±Ω, respectively. remarkable features that set it apart from a usual s-wave SF and can be used as its "fingerprint". Perhaps its most revealing property is an analog of the spin-galvanic effect, when an applied Zeeman magnetic field induces a bulk supercurrent [20,21]. Since the pseudospin degrees of freedom inĤ ef correspond to a site index inside the unit cell, this "field" must couple to the motion of g-atoms and can be implemented as a laser-assisted tunneling within dimers [30,31]. This hopping can be set to have an arbitrary phase θ, but we focus on the case θ = π 2 , i.e. consider a perturbation: δĤ ef = −b y i, ab σ y abĝ † iaĝ ib and compute supercurrent response K = κ b y [ Fig. 4(a)]. We can anticipate this response based on pure symmetry arguments: as explained previously in this section, the SOC, being odd in momentum, breaks the space inversion symmetry in the dimer lattice. On the other hand, the synthetic Zeeman term δH ef with b y = 0 violates timereversal symmetry (see discussion in Methods). Breaking of these two symmetries is a necessary condition to stabilize a state with a non-zero current in the system.
The operatorK is a single-particle mass current, obtained by varying the HamiltonianĤ b =Ĥ g 0 + δĤ ef [cf Eq. (2)] w. r. t. the flux ϕ piercing the ring: K = δĤ b /δϕ ϕ=0 . This flux enters via a phase factor e iϕ on each physical link in the lattice with onesite unit cell, but in the dimerized lattice one must In Fig. 4(b) we show the magneto-electric coefficient κ as a function of the chemical potential for several interaction strengths. Remarkably, κ exists only inside the SF phase and vanishes across the TSF-insulator transition.
Due to the p-wave nature of the SF phase, the synthetic field b y induces a pseudospin polarization Ŝ y = χb y , The susceptibility χ and magneto-electric coefficient κ can be related in the weakcoupling dilute limit u gg J g , n g 1 at T = 0 [19]. Indeed, as explained in Methods, similar calculations that led to Eq. Fig.  4(b) shows χ and κ as functions of the chemical potential µ across the TSF-insulator transition. In cold-atom experiments it is possible to measure χ, for example by a Ramsey-type protocol [45]: Assuming that the system is in its GS |ψ 0 (with b y = 0), at t = 0 we quench the Hamiltonian from (2) toĤ x = J gŜx e.g. by making the intra-dimer g-atom tunneling dominant and let the system evolve for a time t 0 = π 2Jg . As a result, the state becomes |ψ = e −i π 2Ŝ x |ψ 0 . Now we measure the difference in populations on two sites of a dimer, i.e. ψ|Ŝ z |ψ . Because e i π 2Ŝ xŜ z e −i π 2Ŝ x =Ŝ y , the above protocol yields ψ 0 |Ŝ y |ψ 0 and can be used to obtain χ.
Another physical effect induced by δĤ ef is an asymmetry of the g-atom momentum distribution which can be detected in time-of-flight experiments [46]. Because these measurements involve crystal momentum p in the BZ of a single-site unit cell, we need to compute Fig.  4(c) shows ν p computed in the TSF phase of Fig. 2(a). For comparison, in a non-SF system, ν p ∼ δ k,k F − δ k,−k F (k F is the Fermi momentum). SF correlations destroy Fermi points and lead to a finite asymmetry even away from k = ±k F .
Preparation of the 87 Sr lattice clock.-The system in Fig. 1(a) and Eq. (1) can be realized with AEAs, such as fermionic 87 Sr, using steps shown in Fig. 5: Step (0) We start with a nuclear-spin polarized g-atom band insulator in a deep magic-wave lattice (in which e and g atoms experience equal light shifts and therefore same trapping potential [47]) with suppressed tunneling.
Step (1) The system is irradiated by a laser that adiabatically applies a staggered synthetic gauge field [48,49] with wavevector k = π, Rabi frequency Ω and a detuning from the e-g transition δ.
In the e-g basis, the single-atom Hamiltonian at the lattice site j isĤ j = 1 2 (−1) j Ωσ x − δσ z . Here the Pauli matrices σ act on the local g-e basis. As the detuning is decreased to zero, the Rabi frequency is simultaneously ramped up, thus adiabatically preparing atoms in their local GS 1 √ 2 |e ± g j .
Step (2) The laser wavevector is quenched from k = π to 2π, making half of the local e-g mixtures excited states.
Step (3) δ is adiabatically increased, while Ω is decreased. As a result, excited (GS) coherent e-g superpositions are transferred to e (g) states.
Step (4) A laser is used to directly transfer GS e-atoms [|e 0 ] to g-atoms [|g 0 ] and |g 0 to excited e-atoms [|e 1 ], where |e(g) n indicates an e (g) atom in n-th lattice band.
Step (5) We adiabatically enable hoppings by decreasing the magic-wavelength lattice depth. Simultaneously, we create a dimerized e-superlattice by ramping up a potential experienced only by e-atoms [50] at twice the periodicity of the magic lattice, and transfer the states |e 1 to excited antisymmetric motional states in each double-well, in order to satisfy the requirement |J e | > |J g |.
After the last step, g-atoms form a band-insulator with n g = 1 (two atoms per dimer). Their filling can be controlled spectroscopically by removing atoms from kstates near band edges with a laser which drives a narrow transition whose detuning is adiabatically changed to scan the conduction band and access atoms deeper in the Fermi sea. This can be visualized as an adiabatic injection of holes in the presence of e-atom background. The newly added holes form Cooper pairs, thus building up a SF state.
Discussion.-Topological superfluidity in Fermi liquids is intimately related to the coupling between particles' spin and orbital motion. Unfortunately, in most systems this crucial ingredient is absent or too weak to yield a measurable topological structure of SF phases. In the present work we discussed a mechanism that bypasses this "rule" and allows a coexistence of a strong spin-orbit coupling and pairing correlations. The key ingredient in our theory is the lattice modulations that host localized degrees of freedom and play a dual role. On the one hand, quantum fluctuations of localized fermions stabilize SF states in the itinerant channel, even when the bare interactions are repulsive. On the other hand, the modulations enlarge the lattice unit cell and lead to an emergent odd in momentum SOC in the conduction band. A combination of these effects always results in a topologically non-trivial SF state in a number-conserving system with potential emergence of Majorana modes. We illustrated the above mechanism by studying a system of spinless fermions in a quasi-1D lattice with a dimerized structure and showed how one can observe this physics in a quantum simulator with AEAs with a variety of probes, including momentum-resolved spectroscopy and an analog of the spin-galvanic effect, i.e. a magneto-electric phenomenon that can be used to detect a SF phase with broken inversion symmetry.
The analysis presented above can be easily extended beyond 1D. In particular, in a 2D system where g-atoms propagate in a square lattice, and e-atoms are localized inside square plaquettes, similar arguments show that quantum fluctuations of the e-atoms stabilize a p x + i p y SF state of the g-species [44].
Although our presentation illustrates main ideas behind this emergent phenomenon using mean-field approximations, the topological nature of the SF state remains intact beyond mean-field. We confirmed this by performing exact diagonalization in a single tube and verifying that the GS realizes a fermionic parity switch [4] for all fermion fillings [44].

I. METHODS
Derivation of the g-atom effective Hamiltonian (2) We focus on the regime J e J g , U − eg when the interaction termĤ int = U − eg ian g ian e ia and g-atom kinetic energyĤ g 0 in Eq. (1) can be treated as a small correction to the kinetic energy of e-atomsĤ e 0 = −J e i (ê † i1ê i2 + h.c.), and use perturbation theory to compute the effective low-energy model (2). The zero-order subspace is spanned by the states |Ψ eg = i |λ i ⊗ |Ψ g . Because J g J e , we consider this subspace degenerate for all |Ψ g . To determine 1st (2nd) order corrections to the wavefunctions (energy eigenvalues), we employ the Schrieffer-Wolff transformation (SWT) method which systematically removes off-diagonal matrix elements of H int between states |Ψ eg to each order in U − eg 1 :Ĥ = H e 0 +Ĥ int → eŜĤe −Ŝ withŜ =Ŝ 1 +Ŝ 2 + . . . being an anti-hermitian generator of the SWT andŜ n ∼ (U − eg ) n . Notice thatĤ g 0 acts only on the g-atom states |Ψ g and contains only diagonal matrix elements, thus representing a 1st order correction. On the other hand,Ĥ int has both off-diagonal and diagonal elements. The latter are proportional to ian g ia , i.e. they only correct the g-atom chemical potential and can be omitted.
To the first order in U − eg ,Ŝ is given by: whereP Ee is the projector onto the state i |λ i of the esubsystem with an energy E e , and n e ia is an off-diagonal part of the matrix i λ |n e ia |λ i = 1 2 δ λ λ + 1 2 σ x λ λ (δ a1 − δ a2 ) (δ ab is a Kronecker delta-symbol).

Unconstrained Hartree-Fock-Bogoliubov theory in position space
In the present section, we summarize the fully unconstrained Hartree-Fock-Bogoliubov (HFB) mean-field approach 2 which we used to compute the phase diagram in Fig. 2. This variational technique treats on an equal footing SF and normal phases, including various insulating states that break translational symmetry, and hence provides an additional check of robustness of the TSF phase. Our method amounts to the following linearization of the four-fermion interaction: whereĥ H,F,B refer to Hartree, Fock and Bogoliubov contributions with (in general complex) OPs n g ia = ĝ † iaĝ ia , ξ i = − ĝ † i2ĝ i1 , and ∆ i = −u gg ĝ i2ĝi1 . The full HFB mean-field Hamiltonian includesĤ int and the g-atom kinetic energyĤ g 0 [see Eq.

Momentum-resolved spectroscopy signal
Here we derive the expression for R τ (δ, k) used in the main text. Our approach almost exactly follows Ref. 3. We use units such that = k B = 1.
The transfer of g-atoms from the SF phase to e-states in an empty tube is governed by the operator: where Ω is the Rabi frequency andq † ia creates an eatom in well a of the i-th dimer in the auxiliary (empty) tube. The initial state of the system contains N g gatoms, |ψ in = |ψ i,g Ng ⊗ |vac (|vac is an empty state of the auxiliary tube). Its energy is E in = E g i + µN g + ω 0 , where ω 0 is the laser frequency (note that we work in the grand-canonical ensemble, so E g i already includes the chemical potential offset). In the final state, there is a single e-atom with momentum k and band index τ in the auxiliary tube: |ψ fi = |ψ f,g Ng−1 ⊗f † kτ |vac (|ψ f,g Ng−1 is an intermediate state of g-atoms with N g − 1 particles). This state has an energy E fi = E g f + µ(N g − 1) + ν eg + e kτ , where ν eg is the e-g transition frequency (∼ 10 14 Hz for 87 Sr). Because ν eg is a very large energy compared to SF energy scales, we introduce a laser detuning δ = ω 0 − ν eg . The operatorsf kτ are similar to the f -modes defined in the text: . η parameterizes the e-atom band structure in the empty tube via e kτ = τ r k = τ 1 + η 2 + 2η cos k.
The transfer rate R can be computed with the aid of the Fermi golden rule 3 : where summation is extended over all initial and final states, and Z is the grand partition function. The matrix elements ofV are given by where f |ĝ ka |i = ψ f,g Ng−1 |ĝ ka |ψ i,g Ng . Using the single-particle spectral density 4 , we obtain after a straightforward calculation: with ω = e kτ − µ − δ. Noting that A 21 = A * 12 , we arrive at the expression in the main text.

Magneto-electric phenomena
In this section we discuss novel magneto-electric effects that occur as a result of the odd in momentum SOC in the effective model (2) when it is subjected to a weak laserinduced synthetic magnetic field. We assume this field to be homogeneous: δĤ ef = −b · k σ abĝ † kaĝ kb . In particular, we shall derive weak-coupling expressions for the susceptibility χ and magneto-electric coefficient κ quoted in the main text.
It is well-known that a physical external magnetic field breaks time-reversal symmetry T . However, because b is artificial, it does not necessarily break T . Indeed, under time-reversal, the second-quantized operators and cnumbers transform asĝ ka →ĝ −k,a (because a is not a real spin but the position index inside a dimer) and c → c * . Therefore, only the b y -term in δĤ ef breaks T and is capable of generating a mass current. Below we focus on the case with b x = b z = 0 and b y = 0 [see Fig. 4(a)].
It is instructive to study the non-interacting system described by the HamiltonianĤ g Assuming that T = 0 and f -particles fill the Fermi sea |FS with a chemical potential µ, i.e. FS|f † kτ f kτ |FS = δ τ τ θ(µ − kτ ), the average mass current is (the operator K was defined in the main text) This result is valid even in interacting non-SF systems. Hence a magnetic field cannot generate a mass current in the absence of SF correlations 5 . This situation changes dramatically when pairing correlations are taken into account, technically because the average current can no longer be written as an integral of the quasiparticle group velocity. Consider the weak-coupling dilute limit u gg J g and n g 1. In this regime, we can project interactions in Eq. (2) onto the lowest band τ = −1. At small fields b y J g , ∆, the bare fermion operators becomeĝ k1 ≈ 1 √ 2f k and g k2 ≈ 1 √ 2 e ik/2 (1+ib y /2J g )f k , wheref k,−1 ≡f k . It is easy to check that b y enters the projected interaction term in Eq. (2) only via quadratic corrections ∼ b 2 y . Therefore, to linear order in b y the projected Hamiltonian is the same as Eq.

the BdG Hamiltonian is
respectively. However, energies of the quasiparticlesγ k andΓ k are now split by a field correction and are given by E k,± = E (0) k ∓ b y sin k 2 , respectively. If b y is small, we can assume that sign(E k,± ) = ±, and in the BCS GS γ † kγ k = ν(E k,+ ) and Γ † kΓ k = ν(E k,− ) with the Fermi function ν(x) = e x/T + 1 −1 .
The projected current and y-component of the pseu-dospinŜ y have the form: In the BCS GS the first term in both expressions vanishes at T = 0 because When T = 0, it is easy to show that κ is given by the expression in the main text. Fig. 4(b) shows κ and χ numerically computed beyond the weak-coupling dilute limit. In the simulation we used the same strategy as above: solve the BdG equations in the field, compute K and Ŝ y as functions of b y , and extract the corresponding linear coefficients.

BCS theory in a translationally invariant 1D system
The unconstrained mean-field approach presented in Methods can be handled only numerically. However, if we restrict the variational space to translationallyinvariant states in a chain with periodic boundary conditions (one can use the momentum k-space representation) and omit the Hartree-Fock terms inĤ HFB , the resulting BCS theory can be treated analytically and a number of insights can be gained on the structure of the topological SF phase. In particular, in this section we analyze the singlet-triplet mixing of Cooper pair states due to the SOC in Eq. (2).
At the BCS level 6 , the linearized interactionĤ int in Eq. (M1) contains only theĥ B term: is an s-wave OP. Because the above interaction couples fermions with opposite momenta, to avoid overcounting of k-states we will restrict momentum summations to a "positive" half of the BZ with k > 0 (indicated by a prime), and denoteĝ −k,a ≡Ĝ ka (−k belongs to the other half of the BZ): Similarly, the kinetic energy in Eq. (2) can be written aŝ Up to constant terms, the full BCS Hamiltonian iŝ The BdG Hamiltonian H BdG can be easily diagonalized when the OP ∆ is purely imaginary. This scenario is still quite general, and below we will focus on the case The eigenstates ψ ντ are labeled by two indices τ and ν = ±1, and have the form with h k = (k) − µ + i∆ 0 σ y and the spinor The corresponding energies are E kντ = νT kτ with T kτ = µ 2 + ∆ 2 0 + 2J 2 g (1 + cos k) + 2τ J g R k , R k = D 2 k + 2µ 2 (1 + cos k) and D k = ∆ 0 (1 + cos k). The parameter ∆ 0 is determined from the BCS selfconsistency condition: with the inverse temperature β = 1/T . To solve this equation, let us assume that the Fermi level (corresponding to momentum k F ) lies below zero, i.e. µ = −J g 2(1 + cos k F ) < 0. Near the Fermi level δk = |k| − k F , T k,+1 = 2|µ| and T k, In the weak coupling and zero-temperature limits, the main contribution to the sum comes from around the Fermi surface: From here, ∆ 0 = 2ρe −2πJg/ugg √ 1−(µ/2Jg) 2 with ρ ∼ J g ∆ 0 -a characteristic energy measured from the Fermi level beyond which the above expansion becomes invalid.
The OP ∆ 0 involves g-fermions with opposite pseudospin a and describes an s-wave Cooper pairing. However, due to the SOC inherent inĤ g 0 , there is also an admixture of the p-wave pairing states of two g-atoms with the same pseudospin which, because of Fermi statistics, must be antisymmetric in momentum. To quantify these p-wave correlations, we compute a spin-averaged pairing amplitude P = 1 2∆0 a ĝ kaĜka using the BdG wavefunctions (S1): This result demonstrates that p-wave correlations in the system are enslaved to the s-wave OP. Moreover, the τdependence of the quasiparticle energy is essential: if T k,+1 = T k,−1 , as it would be in the absence of SOC, there is no p-wave pairing and P ≡ 0. Finally, we derive the BCS approximation in the weakcoupling limit described by Eq.
where D k = 2i∆ sin k 2 , −k = k , as before k . . . extends over half of the BZ, andF k =f −k . This Hamiltonian is diagonalized by a Bogoliubov transformation from ffermions to quasiparticlesγ k andΓ k 7 : Here E k = ( k − µ) 2 + |D k | 2 is the quasiparticle dispersion. The SF gap has a p-wave symmetry: D k = −D −k and its amplitude ∆ is determined from the self-consistency equation which coincides with the result obtained in the multi-band case.

Topological superfluidity in a 2D lattice
The results presented in the main text, can be extended to higher dimensional systems. Here we briefly consider a generalization involving a lattice that has the plaquette structure shown in Fig. S1(a). In each plaquette, there is one e-atom that can tunnel around the plaquette (with lattice constant a 0 ). On the other hand, gfermions move in the simpler square lattice. As before, all atoms are nuclear-spin-polarized and we use units where a 0 = 1. The Hamiltonian of this system iŝ ; where i = x i = (x i , y i ), i, j = 1, . . . , N label plaquettes in the lattice, a, b = 1, 2, 3, 4 denote wells inside a plaquette, ab indicates all sides ab of a square, and ij x,y is a link connecting two plaquettes in the x or y direction [ Fig. S1(a)]. Other notations are the same as in Eq. (1).
In the rest of this section, we shall follow an analysis similar to the one in the main text: First, we consider e-atom states of an isolated plaquette and identify a spatial mode that gives rise to the g-atom attraction within that plaquette. Then, we study a weak-coupling dilute limit u gg J g , n g 1 and demonstrate the stability of a chiral p x + ip y TSF state. The study of the phase diagram, magneto-electric effects, etc. is left for a future investigation.
a. Pairing of the g-atoms.
-An e-atom localized within a plaquette, has four states λ = 0, . . . , 3 shown in . In the 1D case we saw that an attractive interaction between g-atoms occurs when the e-subsystem is prepared in the highest excited kinetic energy state. An analogous situation happens here: e-atoms must fluctuate out of the d-wave λ = 2 state with energy 2J e . Using the Schrieffer-Wolff transformation derived in Methods and performing similar calculations, we obtain the second-order g-atom Hamiltonian ) . (S5) Here u gg = (U − eg ) 2 /32J e , and density terms in the 1st (2nd) line describe attractive interactions along sides (diagonals) of the i-th plaquette [gray ellipses in Fig. S1(a)].
b. Weak-coupling dilute regime.-The g-atom kinetic energy can be written aŝ Other notations are as in Fig. 1(a). (b) States of an eatom within a plaquette. Also shown are point symmetries of the wavefunctions X (λ) a = 1 2 e i π 2 λ(a−1) with colors indicating phases ±1 and ±i. In particular, the λ = 0 (2) state has s-(d-) wave symmetry. Attractive interactions between g-atoms occur when e-atoms are prepared in the λ = 2 state.
The low-energy effective Hamiltonian can be obtained by replacingĝ k → ψ 0 f k (f k is the lowest-band fermion quasiparticle) in Eq. (S5): This expression describes a system of spinless fermions interacting via an attractive p-wave coupling. Within the BCS approximation we have: Here we switched to a grand-canonical ensemble with a chemical potential µ and performed the standard replace-mentĤ ef →Ĥ ef − µ ian g ia . ξ k = k − µ, D k = (k · ∆), and other notations are the same as in the 1D weak- k k F kfk obeys the BCS equation where E k = ξ 2 k + |D k | 2 and in the last sum we used evenness of the integrand and extended summation over the entire BZ. Whenever the SF phase develops, the grand potential G s in that phase is reduced compared to its normal-state value G n . This shift can be computed using the Hellman-Feynman theorem 7 There are two competing SF states characterized by different symmetries of the OP ∆: (i) chiral p x ± ip y phase with ∆ p+ip = ∆ 1 (e p x ∓ ie p y ) (p = p x e p x + p y e p y ), and (ii) p x (or p y ) state with ∆ px = ∆ 2 e p x (or ∆ 2 e p y ). Below we assume real amplitudes ∆ 1,2 and demonstrate that at weak coupling, the p x + ip y state is favored. c. p x + ip y SF.-The u gg -dependence of ∆ 1 can be determined by multiplying Eq. (S6) by ∆ * [notice that |∆| 2 = 2∆ 2 1 and |k · ∆| 2 = k 2 ∆ 2 1 ]: 2 = 8π , k F is the Fermi momentum, and ρ ∼ J g ∆ 1 is a characteristic energy cutoff [cf. Methods]. We have: -Proceeding in a similar manner as above, we use the BCS equation (S6) to calculate ∆ 2 where I = −2 dθ 2π cos 2 θ ln | cos θ| ≈ 0.2. δG can be computed as before Because e 2I /2 < 1, δG 1 < δG 2 and the p x ± ip y SF state is preferred over the striped p x,y phase.
Fermion parity switch in the Hamiltonian (2) In this section, we show that the exact GS of the effective Hamiltonian (2) realizes a fermionic parity switch 8,9 , i.e. the inverse compressibility, X = 1 2N E 0 (N g + 1) + E 0 (N g − 1) − E 0 (N g ), changes sign when the boundary conditions are switched from periodic to anti-periodic, thus confirming the topological nature of our system beyond mean-field. In Fig. S2 we present results of the exact diagonalization (using the Lanczos technique of Ref. 10 ) in systems with N = 24 and 26 sites (12 and 13 dimers, respectively) and even number of g-atoms N g = 2, . . . , 10. Clearly, X has different sign for the two types of boundary conditions.