Strain and pseudo-magnetic fields in optical lattices from density-assisted tunneling

Applying time-periodic modulations is routinely used to control and design synthetic matter in quantum-engineered settings. In lattice systems, this approach is explored to engineer band structures with non-trivial topological properties, but also to generate exotic interaction processes. A prime example is density-assisted tunneling, by which the hopping amplitude of a particle between neighboring sites explicitly depends on their respective occupations. Here, we show how density-assisted tunneling can be tailored in view of simulating the effects of strain in synthetic graphene-type systems. Specifically, we consider a mixture of two atomic species on a honeycomb optical lattice: one species forms a Bose-Einstein condensate in an anisotropic harmonic trap, whose inhomogeneous density profile induces an effective uniaxial strain for the second species through density-assisted tunneling processes. In direct analogy with strained graphene, the second species experiences a pseudo-magnetic field, hence exhibiting relativistic Landau levels and the valley Hall effect. Our proposed scheme introduces a unique platform for the investigation of strain-induced gauge fields, opening the door to future studies of their possible interplay with quantum fluctuations and collective excitations. Realizing and controlling artificial gauge fields in quantum systems constitutes an intriguing route towards simulating of exotic quantum theories. Here, the authors propose a tunable strategy to engineer strain and synthetic magnetic fields in optical lattices from coupling to trapped Bose-Einstein condensates.

T he rise of cold atoms in optical lattices as a versatile platform to study quantum phases [1][2][3][4] has led to the realization and investigation of rich physical models. Building on the milestone implementation of the Bose-Hubbard model 5 , a model of bosonic particles on a lattice with onsite (contact) interactions 6 , a variety of opportunities have become available. For instance, the realization of the Fermi-Hubbard model with cold atoms [7][8][9][10][11] represents a promising route to unveil the microscopic origin of high-temperature superconductivity 12 . More recently, tunable long-range interactions have been introduced in atomic lattice systems [13][14][15] , as well as more exotic features 16 , such as SU(N)-symmetric interactions 17,18 , and density-assisted tunneling 19-26 . Gauge fields are at the core of remarkable phenomena in condensed matter, as was exemplified by the discovery of the quantum Hall effects and topological materials [27][28][29] . These exciting topics have become accessible in ultracold gases through the design of synthetic gauge fields [30][31][32][33][34] . One of the key methods to realize synthetic gauge potentials in quantum-engineered systems consists in driving the system periodically in time [35][36][37] , a general scheme also known as Floquet engineering 38 ; in this driven-lattice context, tunneling matrix elements acquire welldesigned complex phase factors, known as Peierls phase factors, hence mimicking the Aharonov-Bohm effect caused by an external magnetic field. These Floquet schemes have been applied to generate, for example, Peierls phase factors in one-and twodimensional lattice geometries 39,40 , and to realize the Harper-Hofstadter [41][42][43][44][45][46][47] and Haldane-type 48-51 models in optical lattices.
An exciting scenario, which has become more concrete and realistic over the last few years, concerns the realization of dynamical gauge fields in cold gases, namely, engineered gauge fields that experience a back-action from the matter degrees of freedom [52][53][54][55][56][57][58][59] . A principal motivation behind such developments concerns the elucidation of non-perturbative effects in lattice gauge theories (LGTs) [60][61][62][63] . First realizations of densitydependent gauge fields in cold gases (which did not satisfy the constraints of a gauge theory), were reported in refs. 64,65 . Besides, major advances in the quantum simulation of LGTs have been achieved in trapped ions 66 , and more recently, in ultracold atoms 67,68 .
It is well established that artificial gauge fields can also be engineered in the solid state, for instance, by applying strain to materials 33 . In the context of graphene [69][70][71] , strain generates an effective "magnetic" field, which strongly modifies its low-energy relativistic excitations: strain induces relativistic Landau levels in the vicinity of the Dirac points 72 . The main and crucial difference with the action of a real magnetic field is that time-reversal symmetry is preserved in strained graphene. As a result, the vector potential that emerges from the strain field has opposite signs at the two valleys, thus providing the conditions for the valley Hall effect 73 . The characteristic relativistic Landau spectrum has been successfully observed in graphene 74 and in molecular graphene 75 . In synthetic systems, lattice patterning is often an intrinsic requirement, such that an external stretching is not needed to produce the effects of strain. Instead, strain can be mimicked by displacing the lattice sites according to the most convenient profile [76][77][78][79][80][81][82] . Similar strategies have been exploited to investigate the physics of strained honeycomb lattices with arrays of optical waveguides 83 , microwave resonators 84 , excitonpolaritons 85 , acoustic metamaterials 86 , dipole emitters embedded inside a cavity waveguide 87 , and ring resonators arrays 88 . In contrast, optical-lattice potentials for ultracold atoms are typically rigid: their perfect periodicity is generally fixed by the lasers wavelength. This makes the realization of strain more challenging in cold atoms than for other synthetic-matter platforms. We note that a promising proposal, which consists in displacing one of the three laser beams generating the honeycomb-lattice potential, was described in refs. 89,90 . Very recently, effects of elasticity as described by phonon modes have been experimentally explored in optical lattices by coupling a BEC to multimode cavity photons 91 .
In this paper, we introduce a different strategy to realize and investigate the effects of strain in optical lattices, which is summarized in Fig. 1a. Our scheme builds on a mixture of two atomic species, one of which is bosonic (denoted by ↑) and forms a Bose-Einstein condensate (BEC), while the second species (denoted by ↓) can be either bosonic or fermionic. As a central ingredient, the two species are assumed to be coupled through a density-assisted tunneling term, which affects the hopping of ↓ atoms through the density of ↑ atoms. When the BEC is harmonically trapped, the density of ↑ atoms is inhomogeneous, and the correlated tunneling of ↓ atoms displays the effects of a fictitious uniaxial strain: the ↓ atoms behave as electrons moving in a strained lattice. Throughout our analysis, we make the assumption that the BEC forms a static classical background and that the back-action from the second species can be neglected. Within this framework, we discuss the scheme for two different regimes of the condensate, namely the non-interacting and the Thomas-Fermi regimes. In both cases, we show that the spectrum Fig. 1 Description of the strain scheme. a Representation of the density-assisted hopping model on the honeycomb lattice of length L x = 2x c . The parameter a denotes the lattice spacing. The shaded area depicts the inhomogeneous Bose-Einstein condensate. In the zoomed panel, we depict the hopping process of the ↓ atoms between two neighboring sites occupied by a different amount of ↑ atoms. The unit cell of the lattice where we assume translational invariance in the y direction is shown by the green dotted line with zig-zag terminations. b First Brillouin zone of the honeycomb lattice where Dirac points at K ¼ ð2π=3a; 2π=3 ffiffi ffi 3 p aÞ and K 0 ¼ ÀK are indicated. associated with ↓ atoms can display pseudo-Landau levels under proper conditions. We propose a Floquet driving protocol to engineer the required coupling between the two species and we outline probing methods to extract the spectral and topological features. We note that models with similar features have been suggested in different contexts, for instance, to design an atomic dissipative bath 92,93 , in order to simulate the Su-Schrieffer-Heeger instability 94 , to study the back-action of dipolar crystal 95 and vortex lattice fluctuations 96,97 . In the present context, fluctuations of the strain field originate both from the BEC quantum fluctuations and from the back-action of the second species. We leave to future work the task to identify the impact of these fluctuations onto the static background picture that we have employed in our analysis. These effects open the door to interesting scenarios as they would allow us to explore the interplay of correlations and dynamics involving the two atomic species and the effective gauge structure.

Results and discussion
Strain on the honeycomb lattice. The Hamiltonian of a singleparticle on the honeycomb lattice in the tight-binding approximation readŝ whereâ r ;â y r (b r ;b y r ) are respectively the annihilation and creation operators at position r ≡ (x, y) in the A (B) sublattice, the quantities t j are the nearest-neighbor hopping amplitudes and a=2Þ, as in Fig. 1a. In momentum space,Ĥ 0 can be rewritten aŝ with Assuming C 3 discrete rotational invariance, namely t j = t, the Hamiltonian around the time-reversal invariant points K and −K in the Brillouin zone shown in Fig. 1b reads where ζ = ±1, q ≡ k − ζK, v F 3ta=2ℏ is the Fermi velocity and σ x , σ y are Pauli matrices. Equation (4) describes a relativistic Dirac particle whose linear dispersion relation is given by The application of a spatial deformation that changes the distance between lattice sites, also known as strain, brings new interesting effects 69,70,77 . Within the tight-binding description and for small deformations, strain affects the tunneling amplitudes, which become spatially dependent as t j → t j (r). Different types of strain can be applied to the honeycomb lattice 98 , but here we will consider the case of uniaxial linear strain along the x direction. For an intensity of strain τ ≪ 1, we assume that the hopping coefficients read with the condition τL x /3a < 1 ensuring that the strain is sufficiently small to avoid a local Lifshitz transition to a gapped state 77 . By introducing this slow space dependence of the hopping coefficients into the Hamiltonian (3), translational invariance along y is preserved and k y remains a good quantum number. In the rest of this work, we will exploit translational invariance by solvingĤ 0 for a stripe of size N x × 2, as shown in Fig. 1a where the unit cell of the y-periodic lattice is highlighted. Uniaxial linear strain on the honeycomb lattice mathematically appears in the Dirac Hamiltonian (4) as a homogeneous magnetic field 70 hðq; ζK; AÞ ¼ ℏv F ððζq y À e Ã A y Þσ x À ðq x À e Ã A x Þσ y Þ ; ð7Þ where A has the form of a vector potential in the Landau gauge In the rest of this work, we will use units where e * = 1. Since time-reversal symmetry is preserved by strain, the corresponding magnetic field B = ∇ × A has the opposite sign for the two valleys 70,71 . As in the non-relativistic case, the spectrum of a relativistic particle in a magnetic field also displays Landau levels (LLs). A major difference with respect to their non-relativistic counterparts is that they are not equispaced in energy. The full expression, which includes a momentum dependence originating from a spatially varying Fermi velocity 77 , reads corresponding to the relativistic LL eigenvectors ψ ν;q y ¼ ðψ A ν;q y ; ψ B νÀ1;q y Þ, where ψ l ν;q y ðrÞ / e iq y y e with x 0 ðq y Þ x c À ζq y ' 2 B indicating the LL center, x c the origin of the coordinate axis and l = A, B the component of the wavefunction associated to the A or B sublattice, respectively. The function H ν is the ν th Hermite polynomial and ℓ B is the magnetic length, related to the strain parameter by In Fig. 2, we compare the numerically calculated spectrum of H 0 in the presence of uniaxial linear strain with the one in the absence of strain, near the K point. We observe that straining the lattice has generated relativistic LLs as predicted by Eq. (9). In Fig. 3a, b, we compare the eigenstates ofĤ 0 corresponding to the first and second LLs at k y a ¼ 2π=3 ffiffi ffi 3 p (or q y = 0), denoted ϕ ν;q y ¼ ðϕ A ν;q y ; ϕ B νÀ1;q y Þ, with the analytical relativistic Landau states ψ ν;q y , for ν = 1 and ν = 2 respectively. Sufficiently far from K, the (almost) flat LLs become strongly dispersive. This effect originates from the dependence of the LL wavefunction center on momentum, x 0 (q y ), which we show in Fig. 3c, d. For the values of q y corresponding to wavefunctions centered near the edge of the system, the hard wall potential lifts these states in energy thus causing a strong dispersion. The energy levels therefore cross the energy gap between two subsequent Landau levels. This is indeed what we expect from the bulk-boundary correspondence for the quantum Hall effect (QHE) that predicts the existence of robust edge modes when the Fermi level sits in the gap between two LLs. Since the dispersion shows opposite slope at the two opposite edges, these modes are obviously chiral. However, we have to recall that on the other valley, the opposite effect will take place as a result of time-reversal symmetry 98 . In the end, no net current can be observed on each edge, unless valley transport can be resolved, an effect known as the valley Hall effect.
Model. In order to generate uniaxial linear strain with cold atoms in optical lattices, we propose to employ a mixture of two atomic species, which we indicate as ↑ and ↓. The ↑ atoms are weakly interacting bosons, harmonically trapped in the x direction. The corresponding Hamiltonian readŝ ";r ðn ";r À 1Þ; wheren ";r â y ";râ";r (b y ";rb";r ) if r 2 A (2 B) and x c is the position of the system's center. The parameters J, U and V x are respectively the nearest-neighbor hopping amplitude, the onsite interaction energy and the strength of the harmonic confinement.
The ↓ atoms, whose statistics does not need to be specified, hop on the same honeycomb lattice as the ↑ atoms according to the HamiltonianĤ where t is the hopping amplitude for the ↓ atoms. The two species are coupled through the interaction term where α is a dimensionless parameter quantifying the interaction strength between the ↑ and ↓ species. Since the functions F j depend on the density of the bosons, this term describes correlated-hopping (or density-assisted) processes where the tunneling of ↓ atoms between two neighboring sites depends on the number of ↑ atoms at these two sites. For the functions F j , we consider the following expression F j ðn ";r ;n ";rþδ j Þ ¼ 1 3 γ j ðn ";rþδ j Àn ";r Þ; ð15Þ where γ 1 = 1 and γ 2 = γ 3 = − 1. As we show below, these functions will generate the artificial strain for the ↓ atoms when the density of ↑ atoms is inhomogeneous. The full model readŝ which we solve in the mean-field (MF) approximation for the ↑ atoms, described by the discrete Gross-Pitaevskii equation. Specifically, we consider a regime where a large number of ↑ atoms form a weakly interacting BEC and where quantum fluctuations can be neglected; we note that such fluctuations could potentially affect the correlated hopping in Eq. (14), but leave the study of this interesting effect for future work. Furthermore, we assume that the impurities (↓ atoms) are weakly coupled to the BEC so as to neglect back-action effects. Within these approximations, the BEC of ↑ atoms acts as a static and classical background for the ↓ atoms. We therefore write the model asĤ 'Ĥ eff # þĤ MF " , whereĤ eff # Ĥ # þĤ MF "# . The density operatorn ";r is then replaced by its mean value n " ðxÞ, where we remove the dependence on y due to the assumption of homogeneity in this direction. After calculating the BEC density profile obtained by solvingĤ MF " , we input the solution intoĤ eff # . As a result, the ↓ atoms experience spatially dependent hopping parameters t eff j that read In order to recover the linear space dependence needed for uniaxial strain, we consider a parabolic profile n " ðxÞ ¼ Àη where the constants η 0 , η 1 will be specified below and depend on the microscopic parameters of the BEC regime considered. As a result, we obtain which reproduces uniaxial linear strain with τ = 2αη 1 , as described by Eq. (6). We obtain new constant terms appearing in Eq. (19) in comparison to Eq. (6), which result in the vector potential The last term in Eq. (20) shifts the position of the Dirac points. However, such a shift is negligible as long as αη 1 ≪ 1, which is the case here. The effect of this term is indeed not observable in the numerical results presented below. Furthermore, note that the LLs energy gaps are not affected because the magnetic field is determined by the derivatives ∂ x t eff j , which yields In the rest of this work, we discuss two regimes where the density profile can be approximated by a parabolic expression Eq. (18): the non-interacting and the Thomas-Fermi (TF) regimes. The former corresponds to the condition U = 0, whereas the latter is obtained for large U such that the BEC kinetic energy becomes negligible. In order to correctly identify the TF regime, we compare the energy functionals x À x c À Á 2 jχ ";r j 2 ; jχ ";r j 2 ðjχ ";r j 2 À 1Þ: with each other, where χ ";r ¼ hâ ";r i for r 2 A and χ ";r ¼ hb ";r i for r 2 B. The system enters the TF regime when E kin ≪ E trap , E int 99 , which is reached for sufficiently large values of U, as shown in Fig. 4. These two regimes are going to be the focus of our analysis, as they will lead to inhomogeneous hopping coefficients for the ↓ atoms described by Eq. (17). The BEC density profiles are shown in Fig. 5a, b, where we also show the corresponding magnetic fields obtained from Eq. (17), calculated by neglecting the space dependence of the Fermi velocity. The origin of the inhomogeneity of the magnetic fields for both regimes is discussed in the two next paragraphs.
Thomas-Fermi regime for the ↑ atoms. We start our analysis by investigating the Thomas-Fermi (TF) regime, in which the gas of ↑ atoms enters when the repulsive interactions dominate the kinetic energy. By inspecting Fig. 4, which is obtained for V x = 10 −6 J/a 2 and a number of atoms per stripe N ↑ = 1.2 × 10 5 , we see that we can safely use the TF approximation for U ≳ 10 −5 J. The corresponding density profile reads where μ is the BEC chemical potential, which we numerically compute from the relation , where E tot = E kin + E trap + E int is the total energy of the BEC. By substituting Eq. (24) into Eq. (17), we find that the effective strain intensity can be expressed in terms of the harmonic trap parameter V x and the interaction strength U as follows, In Fig. 6, we show the spectrum of the ↓ atoms for U = 10 −4 J. The agreement between the numerical results and the analytical predictions in Eq. (9) obtained for τ = τ eff is visible for the first five levels. By looking at the BEC density in Fig. 5a, we see that the TF radius R TF ffiffiffiffiffiffiffiffiffiffiffiffiffi 2μ=V x p marks a separation between a region with strain (|x − x c | < R TF ) and a region without strain (|x − x c | > R TF ). As a consequence, the spectrum displayed in Fig. 6 will also show features of a homogeneous (i.e. unstrained) honeycomb lattice, as we can conclude by comparison with Fig. 2 (gray lines). The interface between strained and unstrained regions is not a hard wall potential, thus allowing for a penetration length of the wavefunctions from both sides. This fact will therefore induce hybridization events between LLs and planewave states that will result in avoided crossings, some of which are visible in Fig. 6. To distinguish the contribution of the two regions in Fig. 6, we superimpose the spectrum (empty circles) of a strained  honeycomb system that only extends over the size of the BEC, namely L x = 2R TF with the corresponding value of the strain intensity τ = 0.003. As expected, the LL plateaus are clearly identified together with the edge states branches on the left side of the spectrum (k y a < K y ). Notice that the TF radius sharp boundary has been suggested to host edge modes in other topological interacting models as for the case of spinful bosons, see Ref. 100 .
Deviations from the ideal strain physics appear on the right side of the spectrum for k y a > 1.25 and, as indicated by the arrow in the spectrum, a distinct branch is present. As pointed out before, the TF radius introduces a separation, or an interface between the two regions. Let us focus, for simplicity, on the left interface at x ≃ 230a. If a hard wall were present and we could therefore cut the system into two separate parts, we would have a zig-zag termination for the unstrained region, which admits E = 0 edge states for k y > K y , and LL edge states at energies above the gap for sufficiently large values of k y . However, since the interface is soft, a hybridization between these two types of states takes place, which results in (i) a gap-crossing energy branch, as indicated by the arrow in Fig. 6 and (ii) a deviation from the ideal case of strained honeycomb lattice, as manifested by the empty markers in Fig. 6 not overlapping with the lines of the spectrum. We show in Fig. 7a-d that the new branch corresponds to LL edge states in the asymptotic limit of large k y , which is also spectrally observed in Fig. 6. A more quantitative analysis of the interface problem, which is beyond the scope of this work, can be addressed by employing the formalism presented in ref. 101 , which would describe the localized solutions for a relativistic particle near the interface between a region with magnetic field and a region without. However, additional care must be taken near the TF radius, as the profile smoothens out, thus causing sharp changes in the magnetic field at the interface, which are visible in Fig. 5a.
Non-interacting regime for the ↑ atoms. For the non-interacting regime, we take U = 0 and we solve the single-particle Hamiltonian for the ↑ atoms. As we will show below, this regime is less ideal for the LL physics, and this is the reason why we present it after the TF regime analysis. The density profile corresponds to the harmonic oscillator ground state, namely a Gaussian profile, as shown in Fig. 5b, which reads where ξ is the width of the cloud, x c is the trap minimum position, N ↑ is the total number of ↑ atoms and the factor 3/4 is a consequence of the honeycomb geometry. Near the center of the system, i.e. for |x − x c | ≪ ξ, the density can be approximated by the parabola n par " ðxÞ ¼ By substituting Eq. (27) into Eq. (17) for t eff j , we obtain the strain parameter Differently from the TF regime, the analytical prediction in Eq. (9) is in good agreement with the numerical spectrum only very close to the K point, as shown in Fig. 8. Important differences with the TF regime appear in the spectrum, whose origin can be identified directly from the inspection of the two types of density profiles (see Fig. 5) and traced back to (i) the severe deviations of the non-interacting density profile from an ideal parabola and (ii) the absence of a sharp transition from a region with strain to a region without it. As a result, we obtain a nonlinear space dependence of strain, which translates into an inhomogeneous magnetic field that decreases in strength away from the center, see Fig. 5b.
In order to elucidate the consequences of the inhomogenous magnetic field, let us assume that the magnetic field changes very slowly in space. Within this picture, which we will address as a local-density approximation regime and requires |x − x c |, ℓ B ≪ ξ, we can approximate the modified magnetic field with a constant B(x) ≈ B(x 0 ), where x 0 is the wavefunction center and it depends on the momentum as x 0 ¼ x c À ζ' 2 B q y . As one moves away from  the K point, the LL wavefunctions is centered further away from x c . From Eq. (9), we know that the LL energy is proportional to ffiffiffiffiffiffiffiffiffiffi ffi Bðx 0 Þ p , which therefore translates into a decrease of the LL energies away from the K point. As a result, the impact of inhomogeneous magnetic field corresponds to a deformation of the LLs which bend down away from the K point. To be more quantitative, let us expand the Gaussian profile to the next order in (x − x c )/ξ, yielding n (4) " ðxÞ ¼ The resulting strain parameter is indeed inhomogeneous and reads τ ð4Þ eff ðxÞ ¼ τ eff 1 À where τ eff is given by Eq. (28). Using the local-density approximation, we therefore replace x ! x 0 ¼ x c À ζ' 2 B q y and obtain the modified LLs energy levels which are shown in Fig. 8. The modified dispersion captures reasonably well the behavior of the lowest LLs, despite the fact that for this choice of parameters the LL wavefunctions are too broad as compared to ξ, thus weakening the validity condition of the localdensity approximation.
Numerical validation. In this section, we further provide a numerical analysis of the results presented so far by showing a direct comparison of the eigenstates ofĤ eff # with the ones expected from the ideal linear strain regime described by the HamiltonianĤ 0 and from the LLs description. In order to establish this comparison, we compute the fidelity F ðν; q y Þ ¼ jhϕ ν;q y jψ #;q y ij 2 where ψ #;q y denotes the eigenstates ofĤ eff # and ϕ ν;q y those ofĤ 0 corresponding to the ν th LL. As mentioned in the previous section, the spectrum ofĤ eff # differs from the one of H 0 due to regions without strain or with magnetic field inhomogeneities. To minimize these effects, we focus on a certain window in momentum space q y 2 ½q min y ; q max y centered around the K point where we find high values of fidelity (F ðν; q y Þ ≥ 0:8) for bulk states in each ν th LL. The results are shown in Fig. 9a, b for ν = 1, 2, 3 in the TF and the non-interacting regimes. The parameters are chosen as in the previous section.
In both cases, we observe that the eigenstates ofĤ eff # reach a high fidelity with the ones ofĤ 0 near the center of the LL, namely near the K point. However, in the non-interacting regime (Fig. 9b), the fidelity decreases as we go away from the K point. This effect originates from the inhomogeneous magnetic field that modifies the wavefunction as compared to the expected ideal LL results. In particular, the magnetic length becomes space dependent and it increases when the magnetic field decreases,  thus enlarging the tail of the wavefunctions and therefore lowering the fidelity.
We now discuss how the LL picture correctly describes the results when we change the effective magnetic field, namely the parameter α. In particular, we focus on the highest fidelity, denoted by F M ðνÞ ¼ max q y F ðν; q y Þ, as a function of α. The results are shown in Fig. 10a for the TF regime and in Fig. 10b for the non-interacting regime. The content of these plots can be understood through a lengthscale analysis. The smallest lengthscale is the lattice spacing a, whereas the largest one (besides the size of the system L x ) is related to the BEC size, namely ξ and R TF for the non-interacting and the TF regimes, respectively. The relevant lengthscale for LLs physics is the magnetic length ℓ B . We therefore conclude that the ideal situation to observe LLs requires a ≪ ℓ B ≪ ξ, R TF .
In both regimes, we find that the fidelity F M is close to 1 for large values of α whereas it drops as α decreases, see Fig. 10a, b. The best scenario is therefore reached for sufficiently large values of α. However, when α becomes too large the LL picture breaks down since the magnetic length ℓ B becomes smaller and lattice spacing effects take place. We can already see this trend for the values of α chosen in this analysis, as shown in Fig. 10c, d. We indeed observe that the fidelity F 0 M ¼ max q y jhϕ ν;q y jψ ν;q y ij 2 between the analytical relativistic Landau levels ψ ν;q y given by Eq. (10) and ϕ ν;q y , decreases as α increases. The impact of the lattice discreteness is more effective for higher LLs, which have more nodes and thus a less smooth wavefunction, which results in a lower fidelity. A second reason for the drop in fidelity F 0 M as α increases comes from the asymmetry (or parity breaking) with respect to the center x c that the wavefunctions manifest and that can be identified by inspecting Fig. 3. There one can recognize that the left peaks of the wavefunctions have different heights with respect to the right peaks, whereas the analytical LL states do not. This feature, which was already noticed in ref. 102 and caused by terms that have been neglected in the effective Dirac description, is negligible for small strain values but becomes more and more relevant for larger ones, thus causing a distinct mechanism for a mismatch with the ideal LL wavefunctions and the drop in fidelity.
We can therefore conclude that the ideal regime requires a not so large value of α because novel effects that invalidate the LL picture take place, as the condition ℓ B ≫ a is not satisfied anymore. On the other side, when α decreases, the wavefunction broadens and the condition ℓ B ≪ ξ, R TF breaks down. We may encounter situations where the lowest LL has a very good fidelity (see Fig. 11a, b) when centered near x c whereas the highest LLs are more strongly affected given their larger size (Fig. 11c, d). In the TF regime, the wavefunctions can indeed cross the interface and hybridize with the planewave solutions of the unstrained region. This effect is shown in Fig. 11c. In the non-interacting regime, one must instead consider the intermediate region where the BEC density is non-parabolic. In this region, the magnetic field is nonuniform as we discussed before, thus implying that the magnetic length acquires a space dependence and becomes larger as we go away from the center, which in turn broadens the wavefunction, as shown in Fig. 11d.
Experimental realization and probing. In this section, we outline a method to experimentally implement the model discussed in the previous section by using a time-dependent scheme and we then discuss possible detection protocols.
Floquet scheme. In order to generate the correlated-hopping parameters given in Eq. (17), we combine a Floquet engineering method inspired from ref. 21 where interactions are modulated in time, and we combine it with the resonant driving scheme analyzed in ref. 103 for a double-well system. Before discussing the coupling between the ↑ and ↓ species, let us briefly review how a  resonant Floquet driving scheme can be implemented to engineer tunneling amplitudes in a double-well system. Let us consider a single species (↓) described by the time (τ)-dependent Hamilto-nianĤðτÞ ¼Ĥ hop þVðτÞ, wherê H hop ¼ Àtĉ y #;0ĉ#;1 À tĉ y #;1ĉ#;0 þ Δn #;1 with Ω the modulation frequency,ĉ y #;i (ĉ #;i ) the creation (annihilation) operator of the ↓ atoms at site i ∈ {0, 1} and n #;i ĉ y #;iĉ#;i . The parameter Δ describes an energy offset between the two sites. Differently from the model in ref. 103 , we have imposed a time modulation for both sites. If the resonant condition Δ = ℏΩ ≫ t is met, we can follow the standard procedure of changing basis to the rotating frame through the unitary transformationR ¼R 1R2 , wherê The resulting Hamiltonian transforms intô HðτÞ ¼RĤðτÞR y À i ℏR∂ τR We obtain an effective time-independent Hamiltonian by taking the time-average ofĤðτÞ which readŝ where t eff ¼ tJ 1 ððK 0 À K 1 Þ=ℏΩÞ, with J 1 ðxÞ being the first Bessel function of the first kind. When its argument is much smaller than 1, namely K 0 À K 1 ( ℏΩ, it can be linearized as J 1 x ð Þ % x=2 thus yielding an effective hopping amplitude In order to understand how to generate the correlated-hopping term, let us replace VðτÞ by an interaction term between the ↓ species and the ↑ species that readŝ H "# ¼ U "# ðτÞðn ";0n#;0 þn ";1n#;1 Þ; ð37Þ where U "# ðτÞ U cosðΩτÞ, as in ref. 21 . We can then identify K i ¼ Un ";i and apply the resonant driving scheme described before, which yields with α 3U=2ℏΩ. We immediately conclude that the following condition αjhfn 0 ";i gjðn ";0 Àn ";1 Þjfn ";i gij ( 1, where {n ↑,i } labels the Fock states, must be satisfied in order to linearize the Bessel function.
The scheme that we have discussed so far provides the correlated-hopping term that is central to our model of strain in the subsection "Model". However, the ↓ atoms must also possess a dominant tunneling amplitude that is unaffected by the interaction with the ↑ atoms. The missing term can be obtained by considering an additional onsite energy driving termŴðτÞ ¼ Δ 0 cosðΩτÞn #;0 that yields While the terms K i are proportional to the density operatorŝ n ";i , Δ 0 is simply a c-number. We thus obtain Under the condition K 0 À K 1 ( Δ 0 ( ℏΩ, we have therefore obtained that the hopping process is described by a dominant tunneling amplitude and by a smaller density-dependent contribution, as required for the realization of strain. Here below, we will focus our discussion on the different possibilities for the second term. The double-well scheme that we have discussed must be applied to an extended honeycomb lattice in order to reproduce uniaxial strain. We have identified two possible implementations, which we sketch in Fig. 12. The first scheme consists in applying an energy offset on each lattice site that grows along the x axis, as shown in Fig. 12a. The Floquet results presented for the doublewell case apply here to the hopping between two neighboring sites of the honeycomb lattice after identifying site 1 as the one with the highest energy offset and site 0 as the one with the lowest energy offset. The generalization of Eq. (38) to the honeycomb lattice is therefore with γ 1 = 1 and γ 2,3 = − 1, which is exactly what we studied in the previous sections. The second scheme presents an offset (Δ) only for the A sublattice and none for B sublattice, as shown in Fig. 12b. In the double-well representation, this means that the site 0 is always a B site and the site 1 is always an A site. The effective hopping amplitude in this case reads Therefore, by following the same reasoning as the one leading to Eq. (20), the vector potential A for the present configuration is given by which corresponds to an effective magnetic field three times larger than the one in Eq. (21). The corresponding LLs energies are For ν = 1, the energy gap is thus larger by a factor ffiffi ffi 3 p in comparison to the energy gaps obtained with the previous scheme presented in Fig. 12a. As this scheme provides a larger gap, it would be more suitable for the experimental detection of LL physics. Further effects and probing methods. In the discussion above, we have not considered the possible renormalization of the hopping parameter J of the ↑ atoms due to the time modulation of interactions. As presented in ref. 21 , if no energy offset is experienced by the ↑ atoms and if Ω ) J; U then the hopping amplitude for the ↑ bosonic atoms is renormalized as J ! J eff ¼ JJ 0 ðαðn #;0 Àn #;1 ÞÞ. When the argument of the Bessel function is sufficiently small, the hopping amplitude is unaffected and no back-action of the ↓ atoms takes place onto the ↑ atoms. This condition can be met when the density of ↓ atoms is small or when it is homogeneous. So far, we have not specified the statistics of ↓ atoms. If we consider a Fermi gas of ↓ atoms, we may have both conditions met at once because Pauli exclusion principle will prevent the onsite density to exceed one atom per site and it will also broaden the density distribution in the presence of a harmonic trap. The latter may result in a very flat density profile over the range occupied by the BEC of ↑ atoms. If we consider a gas of bosonic ↓ atoms, we will instead have to enforce a low density or a homogeneous distribution in order to prevent the back-action on the ↑ atoms. Despite the arguments presented here to neglect back-action effects, it would nevertheless be interesting to include those, as they will provide a distinct opportunity to enrich the strain picture discussed in this work with dynamical effects.
A separate discussion for the actual implementation of the model concerns the trapping potential. The strain model that we have analyzed requires a strongly anisotropic harmonic trap experienced by the BEC of ↑ atoms with ω x ≫ ω y , where ω x,y are the harmonic trapping frequencies in the two spatial directions. In this work, we have actually considered ω y → 0 to simplify the theoretical analysis. However, we have not included trapping effects on the ↓ atoms, assuming that one can independently control the confinement of the two atomic species. In this case, several scenarios are possible, which affect the back-action discussed before and the corresponding probing methods. When the harmonic trap is absent or negligible, we obtain the picture discussed in this work for the single-particle spectrum. However, in the presence of a strong confinement and for fermionic ↓ atoms, we will have the opportunity to reveal the presence of LL physics when the corresponding Fermi level at the center of the system is fixed between the LL gaps. In this case, LLs will manifest through jumps in the density profile that will confer the typical wedding cake structure to the fermionic gas, see ref. 90 . Near these jumps, which are going to be partly smeared out because the LLs are not perfectly flat, we also expect to find valley-dependent edge modes that represent a clear signature of the valley Hall physics.
Several techniques can be employed to directly probe the properties of LLs, as well as the nature of the pseudo-magnetic fields. One possibility would be to monitor the real space dynamics of a wave packet of ↓ atoms, initially prepared in the vicinity of a Dirac point; such a wave packet will exhibit a cyclotron motion, as illustrated in refs. 90,104 . We point out that this would be a direct probe of the valley-dependent pseudomagnetic field (and thus highlight the time-reversal invariant nature of the system), since the chirality of the cyclotron orbits is opposite at the two valleys. Similarly, acting on the wave packet with a constant force will generate a Hall drift 105 , whose direction will depend on the valley considered. Another approach would be to prepare a uniform fermionic gas of ↓ atoms at half-filling. In this case, circular lattice shaking will allow us to spectroscopically resolve the LLs by measuring the absorbed energy. Moreover, band mapping techniques make possible to identify valleydependent absorption processes, thus allowing to extract the valley Hall conductivity through the corresponding valley circular dichroism 51 .

Conclusions
In this work, we have presented a strategy to generate a strain field in optical lattices that is implemented by coupling an atomic species to a trapped BEC via well-tailored density-assisted tunneling terms. By changing the shape of the BEC profile or the type of density-assisted tunneling terms, distinct strain profiles can in principle be generated. We have focused on the implementation of uniaxial linear strain in the honeycomb lattice applied along one of the three crystalline axes of the lattice, which is realized by considering a strongly anisotropic harmonic trapping potential. We have then discussed two limits of interest, namely the noninteracting and the Thomas-Fermi limits. After investigating the spectral features, we have identified the Thomas-Fermi regime as most suitable to reproduce the ideal linear strain configuration. Indeed, this regime minimizes the effects of regions with inhomogeneous magnetic field and requires smaller atomic clouds.
The Thomas-Fermi regime may also provide an interesting scenario to study the effect of quantum fluctuations originating from the phonon modes of the BEC or the effect of exciting the BEC collective modes, which would provide time-dependence to the synthetic gauge field. Some of these effects have a correspondence in solid-state system and originate from the lattice vibrations of the crystal. Another interesting scenario, which is more specific to cold atoms, is to investigate the strongly interacting regime for the bosonic gas near the Mott insulator phase. In this case, low filling and strong quantum fluctuations would provide a very different regime as compared to what is studied in solid-state materials. Our results therefore suggest a distinct direction to investigate the interplay of dynamical gauge fields, as realized through synthetic strain fields, and quantum matter with ultracold atoms.
While our work has focused on the two-dimensional honeycomb lattice where a spin-1/2 relativistic theory represents the effective description, interesting effects originating from strain can also take place in other dimensions (for the case of Weyl semimetals see, for instance, ref. 106 and references therein) or in models represented by higher-spin effective descriptions 107 . Our scheme therefore offers the opportunity to explore more generic strain effects beyond the paradigmatic honeycomb lattice (graphene-type) case. An interesting perspective concerns the possibility of shaping the density profile of the BEC, for instance using digital micromirror devices 108 , in view of realizing tunable strain configurations in optical lattices of arbitrary geometry.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.