Bond Order via Light-Induced Synthetic Many-body Interactions of Ultracold Atoms in Optical Lattices

We show how bond order emerges due to light mediated synthetic interactions in ultracold atoms in optical lattices in an optical cavity. This is a consequence of the competition between both short- and long-range interactions designed by choosing the optical geometry. Light induces effective many-body interactions that modify the landscape of quantum phases supported by the typical Bose-Hubbard model. Using exact diagonalization of small system sizes in one dimension, we present the many-body quantum phases the system can support via the interplay between the density and bond (or matter-wave coherence) interactions. We find numerical evidence to support that dimer phases due to bond order are analogous to valence bond states. Different possibilities of light-induced atomic interactions are considered that go beyond the typical atomic system with dipolar and other intrinsic interactions. This will broaden the Hamiltonian toolbox available for quantum simulation of condensed matter physics via atomic systems.


Introduction
Ultracold gases loaded in optical lattices are an ideal tool for studying competing phases of quantum matter. Engineering the effective potential seen by the atoms using light beams allows to realize with optical lattices simple models of condensed matter, particle physics and even biological systems [1]. Moreover, the realization of these models in experiments would aid in the development of applications towards quantum information processing (QIP) and the development of novel quantum materials via quantum simulation [2]. Typically, one can realize effective Hamiltonians which contain short-range physical processes such as tunneling between neighbor lattice sites and on-site interactions according to a prescribed lattice potential engineered by classical light fields, such as the Bose-Hubbard (BH) model. Long range interactions are experimentally challenging but accessible in principle via polar molecules [3,4] or Rydberg atoms [5,6,7]. However, the nature of the interaction is fixed by the characteristics of its constituents. In addition, finite range interactions by other approaches of light-matter interactions [8,9,10,11], and extended Bose-Hubbard models via dipolar interactions [12,13] are also possible.
In contrast to the above, loading an optical lattice inside a cavity allows to engineer effective synthetic many-body interactions between light induced atomic modes with an arbitrary spatial profile [14,15,16]. These interactions are mediated by the light field and do not depend on the nature of the atoms considered, making them extremely tunable and suitable for realizing quantum simulations of effective many-body longrange Hamiltonians. In principle fermionic, bosonic, molecular systems, etc. can be studied. This allows to explore the interplay between additional non-conventional quantum many-body phases, besides from the typical superfluid (SF) and Mott-insulator (MI). It is now experimentally possible to access the regime where light-matter coupling is strong enough with cavity decay rates of MHz [18,17] and kHz [20,19], in the range to compete with typical short-range processes (tunneling and on-site interactions). Moreover, bosonic ultracold atoms loaded in an optical lattice inside an optical cavity have been recently realized [19,17], opening a new venue to analyse the interplay between competing orders of quantum matter by design. The light inside the cavity can be used to control the formation of many-body phases of matter even in a single cavity mode by properly choosing the arrangement of the cavity, optical lattice and light pumped into the system [21,22,14,15]. Additional freedom can be achieved by multimode cavities or multiple cavities extending the possibilities to condense into exotic quantum phases even further [24,23,25,16]. Recent advances [26,27,28], will enable the experimental realization of synthetic interactions by design with additional freedom in the near future. This can be realized as the cavity parameters (decay rates) and detunnings with respect the cavity modes can be externally modified with respect to the atomic system. Moreover, the spatial profile of the cavity modes can be designed depending on the geometry of the coherent light beams pumped into the high-Q cavity [16]. In this article we present how different arrangements involving multiple probes and/or multiple light modes configurations, lead to the competition of atypical quantum many-body phases via synthetic interactions. In the past, single cavity competition between typical SF and MI phases [29,30], density wave orders [31,32,14,15,16,33,34], and disorder [35] has been addressed. In contrast to other works, here we consider the interplay with what we call "bond order" and other orders in the system. Bond order is a form of self-organization [36] of matter-wave coherences or "bonds" due to cavity backaction to compensate for the phase difference imposed by the pattern of light pumped into the cavity and scattered by the atoms [14]. In the effective Hamiltonian, as the number of cavity or pump modes increases different physics are possible due to the light induced atomic mode structure. Thus, competition of quantum many-body phases triggers due to the induced atomic mode structure (breaking symmetries, e.g. time reversal and translational) and the regular BH Hamiltonian processes ( homogenous tunnelling and on-site interactions). Particularly difficult is the regime where strong correlations will be present in addition to the well known physics of the BH Hamiltonian. This occurs when there is a large effective light-matter coupling relative to on-site interactions. In this limit, standard mean-field theory becomes unreliable as strong imbalanced configurations can occur and large on-site fluctuations take place due to the broken symmetry of the ground state [16]. As the interaction has a global character with non-trivial structure, and in fact is of infinite range, symmetries are broken even in 1D. The ground state can acquire states with different competing orders, Fig.1. Thus, quantum phase transitions different from the usual type in 1D systems, of the Berezinsky-Kosterlitz-Thouless (BKT) type, can occur. In what follows, we analyse the system by using exact diagonalization for small number of sites in 1D. Our simulations are an indicative picture of the expected behaviour in a larger system. Additionally, we will study the behaviour related with the competition between between bond-order and other orders present in the system. In particular, we find numerical evidence to support the analogy between valence bond states (VBS) [37,38] and dimerised states that arise in the cavity system due to bond ordering via different mechanisms. Our results could be used as basis for the quantum simulation of analogous dimer states important in quantum magnetism [39]. Using classical optical lattices the AKLT Hamiltonian can be implemented in principle as the large interaction limit of the Bose-Hubbard Hamiltonian [40]. Non-trivial entanglement properties have been found [41,42], while bulk-boundary duality with entangled pair states occurs [43,44] and spin glasses are also possible [45]. The states are potentially useful for measurement based quantum computation [46]. Other ultracold systems where the possibility of bond ordered states has been explored include dipolar gases with nearest neighbour and truncated finite range density interactions [47,48], in the framework of the so-called extended Bose-Hubbard models [4,13]. Bond order can also occur via density dependent frustration of the hopping amplitudes with Raman assisted tunelling [49]. Moreover in condensed matter systems, bond order was first explored in low dimensions via extended Hubbard models with two species fermions [50,51,52,53].
In our treatment, we provide an alternative route to explore physics of this kind using cavity fields that relaxes the constraint on very strong on-site interaction (U ), as cavity coupling can be tailored to simulate effective Hamiltonians with similar properties. In ultracold atoms in the fermionic version of our system, a closer analogy to investigate and simulate resonance valence bond (RVB) states important in high-Tc superconductivity [54], might be possible. Moreover, the competition between superconductivity and density wave orders is actively studied [55,56], and light induced superconductivity is being researched [57,58,59]. We analyse the emergence and competition between superfluid, supersolid, insulating and dimerised quantum manybody phases of matter by means of the behaviour in their order parameters.
Our findings, will foster the study of competing orders in multicomponent optomechanical systems [60]. Moreover, the interplay of the quantum phases we study and their generalization, may appear in hybrid system networks [61,62,63,64]. In connection to our work, non interacting fermions in cavity systems have been studied [65,66,67] and even chiral states have been found [68]. Towards quantum state engineering via measurement back-action, competition with other correlated quantum many-body states [ [82,83,84,85] can occur. Moreover, feedback control [86,87,88,89,90,91,92] can also be explored in relation to the dynamical stabilization of quantum many-body phases. As such, the behaviour of the emergent phases in the cavity system we will show, might aid towards the design of novel quantum materials with analogous properties. It follows, that the use of the mechanisms described here could be incorporated in the future development of real materials and composite devices in hybrid solid state systems [93], where both light and matter are in the quantum limit and quantum coherence can be exploited.
The article has the structure that follows. We introduce the general model of ultracold atoms in high-Q cavity(ies) where the atoms are in the regime of quantum degeneracy. We continue by stating the effective models we will consider due to synthetic interactions between light induced atomic spatial modes. Then, we define the different competing orders that can arise in the system. Next, we present our results via the phase diagrams of competing phases. Finally, we conclude our manuscript by summarising our findings.

The model
The system consists of atoms trapped in an OL inside a cavity (single/multi mode) or many cavities with the cavity mode frequency(ies) ω c and decay rate(s) κ c [22,74,94,14,15,16,95,96]. The atomic system is subject to additional light beam(s) pumped into the system in off-resonant light scattering. The off-resonant light scattering condition means that Γ |∆ pa |, where Γ is the spontaneous emission rate of the atoms, where ∆ pa = ω p − ω a is the detuning between the light mode(s) frequency(ies) ω p and the atomic resonance frequency ω a . The scattered light from the ultracold atoms in the OL is selected and coupling is enhanced by the optical cavity(ies), generating a quantum potential. The light pumped into the system has amplitude(s) Ω p ∈ C (in units of the Rabi frequency). The pump-cavity detunnings are ∆ pc = ω p − ω c . The light is pumped from the side of the main axis of the high-Q cavity(ies), at an angle not necessarily at 90 • which allows for arbitrary control of the overlap of light induced spatial modes [16]. The cavity modes couple with the atoms via the effective coupling strengths g p = g c Ω p /(2∆ pa ), with g c the light-matter coupling coefficient of the cavity. The lightmatter Hamiltonian describing the system after the light-field has been adiabatically eliminated [21,15] in the good cavity limit (κ c ∆ pc ) is: where H b is the regular Bose-Hubbard (BH) Hamiltonian [101,102], with t 0 the nearest neighbour tunneling amplitude, U the on-site interaction and µ the chemical potential. The operators b † i (b i ) create (annihilate) bosonic atoms at site i, the number operator of atoms per site is given byn i =b † ibi . The on-site interaction and hopping amplitude terms are short-range local processes. The BH Hamiltonian contains the effective parameters forming the classical optical lattice [102]. The emergent effective light-induced interaction is [15,16], withg pc = g c Ω p /(2∆ pa ) where {η, ν} ∈ {D, B}. The sum over "p" ans "q" go over the number of pumps and "c" goes over the cavity modes (for a multi-mode cavity/several cavities). The couplings J pc η,ϕ ∈ C, correspond to the possible values of J pc ij (Wannier overlap integrals) [16,15,14,74,96,95,94] for each mode of the cavity system through the inter-site amplitudes, labeled B, or through the site density, labeled D. These can either be for a single mode cavity with one pump and one cavity or a multi-mode cavity, and even multiple cavities and multiple pumps. These coupling constants are given by, where "i" and "j" can be the same site for density coupling or be nearest neighbours for bond coupling (inter-site densities), where u c,p (x) are the cavity(ies) and pump(s) mode functions, typically travelling or standing waves. The w(x) are the Wannier functions given by the classical optical lattice in the lowest band. The light induced "density"N ϕ and "bond"Ŝ ϕ mode operators are such that: The sums go over illuminated sites N s and nearest neighbour pairs i, j that belong to the light-induced atomic spatial mode ϕ. As it has been shown [16] the coupling constants can be designed with great freedom by choosing the angle of incident light with respect to the classical optical lattice plane and the cavity axis. The spatial structure of light is useful as a natural basis to define these atomic modes, as the coupling coefficients J pc ij can periodically repeat in space [14,15,16,74,75,96,94]. The atoms that belong to a particular light-induced atomic mode scatter light with the same phase. Thus, one can use the distribution of values of J pc ij , to define the light induced spatial atomic modes. As the pump and cavity modes are external to the internal structure of the system (the BH model), they provide a large set of independently tuneable parameters. This allows to tailor the effective light-induced atomic mode interaction with an arbitrary spatial profile. By addressing the density via the couplings J pc D,ϕ one can generate multicomponent density orders. Density wave orders correspond to different groups of atoms for each light induced atomic mode. In the case of J pc B,ϕ one can generate dimer, trimers, tetramers, etc. which will form as a consequence of the pattern induced to the matterwave coherences or "bonds". These two kinds of orders will compete in addition with the superfluid order in the system and the Mott insulating phase of the BH model as we will see. The current effective model disregards additional density dependent Wannier functions modified dynamically by light, which are difficult to calculate selfconsistently. However the proper redefinition and self-consistent determination of these functions won't alter the essential structure of the effective Hamiltonian. This will only renormalise the effective coupling strengths and parameters of the Bose Hubbard model. Thus our results are applicable in a frame of reference with this renormalized parameters. In addition, coupling between cavity modes has not been included, as these processes have much smaller amplitudes compared to the pump modes.

Effective Hamiltonians
In contrast to previous works, here we focus in the large effective light matter interaction where quantum fluctuations cannot be accounted for in mean-field theory regarding bond order. Moreover, we will consider the interplay between density coupling and bond order in the strong-coupling limit. To do this we will analyse the following Hamiltonian corresponding to a single cavity and a single pump, where the incident light illuminating from the side has been designed to scatter through the bonds and densities as a staggered field (at 90 • with respect to the cavity axis [14,16,94]) with components effectively tuned by the couplings of the bond J B and densities J D , the effective interaction strength is g eff /N s ∼γ = ∆ c |g| 2 /(∆ 2 c + κ 2 c ), which depends on amplitude of the light pumped into the system and the light detunnings, Eq.(3). The bondB − and density operatorsD − are: For ultracold atoms in an optical cavity in the adiabatic limit, cavity decay rates are of the order of MHz. The effective interaction strength, g eff , can be typically be made of the same order of magnitude or larger than on-site interactions, [17], with E R the recoil energy. Note that the ratio t 0 /U can be tuned via the classical optical lattice depth and/or Feshbach resonances [102]. Essentially the sign of the light induced interaction can be chosen via the cavity-pump detunning ∆ c and the amplitude by the pump strength Ω p [14,15,16]. In addition, without loss of generality, {J B , J D } ∈ [0, 1]. Depending on the lattice depth of the classical optical lattice, e.g. the Bose-Hubbard Wannier functions and the choice of illumination, the magnitude of the J B,D coupling constants can be tuned using real Wannier functions [94]. Beyond a gaussian ansatz this gives J B = 0 depending on the lattice depth of the classical optical lattice. Typically for a lattice depth of 5E R (where the single band approximation is valid) then ∆φ (Σφ) is the difference (sum) of phases between two crossed standing waves (with phases φ 0,1 ) pumped from the side at 90 • with respect to the classical optical lattice potential. The beams are arranged such that k 0x = 0 and k x,1 = π/a, with a the lattice spacing (typically a = λ/2 for a standing wave in the classical optical lattice). Therefore, the ratio between the value of the two contributions can be adjusted arbitrarily, e.g for J B /J D = 0.25 we have ∆φ ≈ 0.844π with Σφ = 0 for simplicity. Thus, any ratio between the coefficients J D and J B is possible and can be modified in addition by changing the classical optical lattice depth, maximal bond coupling is achieved by ∆φ = π/2 while maximal density coupling is achieved by ∆φ = 0 with Σφ = 0. By increasing the depth of the classical optical lattice the coefficient J B becomes smaller and it is basically negligible for a lattice depth of 15 E R . Moreover, we consider the previously unexplored scenario where via multiple cavities or multiple pumps one can perform the quantum simulation of the following light induced atomic mode interactions [16]: • Bessel type potential, V ϕm = j 0 (π(x m − 1)), see Fig.2(a).
• Morse type potential, The function j 0 is the zero order spherical Bessel function of the first kind and ϕ m with m ∈ [1, R] with R light-induced spatial atomic modes. x m − 1 are the locations of the maxima and minima of j 0 (y) with y ∈ [0, 3]. The maximum amplitude of the interactions has been chosen such that, max(V ϕm ) ∼ 1, for simplicity. Morse type potentials are a typical phenomenological tool to model effective molecular systems. The Bessel type potential we consider, shows an example of the flexibility of the construction with respect to the degree of control that can be achieved via the synthetic light-induced atomic mode interactions. The relationship with the pump and cavity coupling via inverse discrete Fourier transforms can be found in general in [16]. Certainly other types of potentials can be tailored with great flexibility depending on the quantum many-body system we would like to simulate. For many cavity modes (multiple cavities/ multimode cavity) we have, where the interaction depends on the mode distance |ϕ−ϕ |. This kind of effective manybody interaction is physically motivated to account for finite range effective interacting potentials. For many pumps in a single mode cavity, we have, where the interaction depends on the position between light-induced atomic modes. This is potentially useful for the simulation of biological systems and other hybrid networks [61,62,63,64]. Here in contrast to the many cavity mode case, the interaction is position dependent and corresponds to the interaction between different branches, channels or nodes in the network. Without loss of generality, we will consider the case of R = 4 light-induced atomic modes, such that ϕ m with m ∈ [1, 4] for simplicity. The correspondence rules between light induced modes ϕ m and lattice sites i is shown in Fig.2(c) for 1D lattice with N s = 8. In principle, the number of pump modes can be arbitrarily increased by shining the light at different angles with respect to the cavity axis in combination with beam splitters. We call our interactions synthetic, as they are artificially designed by the choice of the spatial profile of the light pumped into the system and the cavity modes [16]. The properties of the light pumped into the system and the cavity parameters are external to the intrinsic properties of the atoms (t 0 and U ) and easily tuned in the range of the atomic processes of the order of the recoil energy.

Order Parameters
Bond order occurs whenever dimerized structures appear in the ground state of the Hamiltonian. These bosonic dimerised structures, akin to valence bond states (VBS) [37], appear in the particular case where the structure of light-matter coupling alternates sign in the inter-site amplitudes or bond between two neighbouring sites [14]. Concretely, the ground state of the system is such that in order to maximise light scattering the inter-site coherences self-organise to minimise the energy. This can be extracted from the ground state configuration via a function of the operatorB − . In exact diagonalisation, if B − = 0, it inherently implies a degenerate quantum superposition (Schrödinger "cat-state"), if B 2 − = 0. Thus, a useful order parameter regarding bond-order can be defined as, akin to a staggered magnetisation, a bond order structure factor [49]. In the case when O B = 0 there is imbalance between mater wave coherences in the ground state. This is a manifestation of a broken time reversal symmetry in the ground state. This is not necessarily coexisting with broken translational invariance, e. g. a ground state with density wave order. It is worth noting that this order parameter will signal bondorder whenever we are not in a MI [94]. Deep in the MI with exact diagonalization, we have: O 2 B | MI = 2(ρ + 1)ρ/N s with ρ ∈ Z + , which in large N s limit vanishes. For a density wave insulator with maximal imbalance we have: O 2 B | DW = 2ρ/N s . In order to tell the difference between a MI and a bond ordered state we will use the fact that the on-site fluctuations ∆(n i ) 2 = n 2 i − ρ 2 are zero in the MI. Thus, we have bond order when ∆(n i ) = and O B = 0. When bond order emerges matter-wave coherence patterns can be from a slight imbalance between matter-wave coherences, . Maximal matterwave coherence amplitude difference is stablished when b † nbn+1 + H.c. = 0 and b † n+1bn+2 + H.c. = 0. The typical matter-wave (MW) coherence patterns found with bond order in the onebody reduced density matrix of sites (i, j) are: • Partial MW amplitude imbalance: ρ 0c cc cc cc c ρ 1c cc cc c cc ρ 2c cc cc c cc ρ 3c cc c cc cc ρ 4c cc c cc cc ρ 5c c cc cc cc ρ 6c c cc cc cc ρ 7 forc and c positive real constants, where each entry in the matrix corresponds to b † ibj for {i, j} ∈ 0, . . . , N s − 1. Thus, distant matter-wave amplitudes are correlated. In a perfect SF (U = 0, e.g. g eff = 0), c =c = ρ i = ρ. Deep in the MI (t 0 = 0, e.g. g eff = 0) ρ i = ρ ∈ Z + , c =c = 0. A pictorial representation is given in Fig.1. On the other hand the matrix representing, the product of nearest neighbour coherences can be constructed, a typical structure for bond-ordered states with maximal phase difference between MW is the following: where α,β,λ,λ are positive real constants. Here each entry corresponds to the product of elements ŝ nŝm withŝ m = (b † mbm+1 + H.c.). The alternating character of the sign of its elements is characteristic of bond ordered states, e.g. for SF and SS all elements are positive. Note that O 2 B = (1/N 2 s ) n,m (−1) n+m ŝ nŝm (where we have used periodic boundary conditions, e.g.b Ns =b 0 ). In the large N s limit we have for a bond ordered state: In terms of the above matrix elements, deep in the MI or DW insulators we have λ =λ = β ≈ 0, thus O 2 B | MI/DW = α/N s for N s 1. Moreover, the above order will compete and coexist with density wave order, typically given by the structure factor [97], DW order breaks translation invariance in the ground state and signals a Z 2 symmetry between odd and even sites (Schrödinger "cat state" in the density configurations).
Deep in a density wave insulator the bond order parameter is: We use the condensate fraction as an estimator of the SF fraction in the system, Alternatively one could use the difference in energy with respect to a phase twist [100].
For an ideal SF, f SF = ρ. At commensurate fillings, in addition to a MI the system can present hidden string order, the string order parameter is given by, with δn k =n k − ρ and ρ the average density per-site (the filling factor). In order to distinguish the MI and string ordered states, it is necessary to define the parity order parameter, In combination with the other order parameters, the string and parity order parameters allow to distinguish the emergence of a Haldane insulator (HI) [38,98,97]. While O DW = 0, O P = 0, and f SF = 0, if O S > 0, with a gapped spectrum, the system is a HI. If the spectrum is gapless (f SF = 0) or gapped (insulator, f SF = 0) and O S > 0, O P = 0, then we have a type of VBS, a dimerised phase. For filling ρ = 1, θ = π, otherwise θ needs to be determined with the help of additional methods [99].
If the system is in the SF state with coexisting bond order O B = 0, it is in the superfluid dimer phase (SFD). The system has matter wave coherence patterns but is homogenous in the density. The typical difference between the order in the ground state for either O B = |O 2 B | or O DW = |O 2 DW | is shown in Fig.1. Note that whenever O DW > O B the system will be in a DW phase either an insulator if the SF component is zero or a supersolid phase (SS) if f SF = 0 . If O DW < O B and with SF component different from zero the system will be in a supersolid dimer phase (SSD). In the SSD, the system has density variation and matter wave coherence pattern with finite superfluid fraction. Whenever the system is in the SS, SSD or SFD phases, the spectrum is gapless as there is a finite SF component in the system. In addition, it can occur that the system is in a bond insulator (BI) phase, where O DW = 0 and O B > 0 and the SF component is zero at incommensurate fillings. Here the system, has phase pattern but there is no SF fraction or DW order and it is not a MI. In BI, dimerised structures form the ground state and it is homogenous in space. It can also happen that, we have a coexistence insulating phase where, 0 < O DW ≤ O B , while the superfluid fraction is completely suppressed, this is a BI+DW insulating phase. This phase is a dimer insulator akin to VBS with density imbalance between components. For commensurate filling ρ = 1 the BI+DW phase presents O S = 0 and O P = 0. Thus bond ordered phases can be gapped (BI) and gapless (SFD, SSD). SSD and SFD phases are bear similarities of RVB states [54], being bosonic gapless ground states with dimerised structures. On the other hand, BI states are similar to VBS, being gapped.
As finite size effects are considerable for small number of sites in the order parameters, to circumvent this problem, we have used the fact that in the large t 0 /U limit the system will tend to be a perfect SF. Thus, all other order parameters besides f SF should approach zero. Therefore, the finite size spurious contribution in other order parameters is eliminated by renormalising with respect to the perfect SF value (f SF in the limit t 0 /U 0). We subtract the SF fraction profile multiplied by the large t 0 /U limit off-set due to finite size. Besides from this, some intermediate phases found will be harder to observe as the number of sites increases, concretely: SFD, SSD and SS which appear as the system moves from insulating states when t 0 /U = 0 to the ideal SF in the limit t 0 /U 1. In table 1, we summarise the quantum many-body phases of the system and the relation with the order parameters defined.
In what follows, we will constrain our discussion on the half-filled ρ = 1/2, and integer filling ρ = 1 cases, while considering simulations for N s = 8 and renormalized order parameters as previously explained.

Bond order vs Density wave order
In this section we will analyse the results from simulations performed using the effective Hamiltonian (6).
At integer filling ρ = 1, the simplest case to understand is when there is only density coupling (J D = 0 and J B = 0), as the density wave instability forms for negative g eff . The system in addition to SF and MI states is able to support the emergence of DW insulator and SS phases, see Fig.3 (a) and (d). The SF and MI exist for g eff J 2 D > −U while the transition point shifts to higher values of t 0 /U as g eff J 2 D > 0 and smaller . values for g eff J 2 D < 0, compared to the system without cavity light. This occurs as on-site fluctuations are enhanced because light scatters minimally being a quantum optical lattice effect, not recoverable by simple mean-field analysis as corrections must be included [14] .Once g eff J 2 D < −U , the system has a discontinuous transition to the DW insulator state. While increasing the effective tunneling, the system goes from DW through SS to SF smoothly. This can be seen in the behaviour of f SF and O DW , Fig.3(g). The insulating character of the MI is confirmed by the absence of string order parameter O S = 0 away from the DW order phases, while having O P = 0, Fig.3(d), while the onsite fluctuations are also minimal Fig.3(j).
When only dimer coupling occurs (J D = 0, J B = 0), the system for strong on-site interactions t 0 /U 1 supports DW for g eff J 2 B 0 and for g eff J 2 B < 0 BI+DW states. For g eff J 2 B > 0, the system evolves from DW to SF via an intermediate SFD phase as tunneling increases. As g eff J 2 B 0 the system goes from the BI+DW state the SF phase as t 0 /U increases rather sharply, Fig.3(h). Complementarily, Fig.3(e), O S is different from zero as SFD and BI+DW phases emerge. We have that O B > O DW = 0 and f SF = 0, thus a BI+DW, a coexisting bond Insulating with density wave insulating phase. This state is different from a SSD state as there is no SF component (the state is gapped), and it is neither a MI, neither a HI. Summarising, the transitions from the insulating phases towards the SF state are sharp for g eff J 2 B > 0 and continuous for g eff J 2 B < 0. The O P smoothly decreases to zero as DW is approached for g eff > 0, while there is sharp change as g eff < 0 in the BI+DW phase, where O P = 0. We consider this as numerical evidence to support that the BI+DW phase is analogous to a gapped VBS

state.
The analogy between dimer states in the cavity system and VBS states can be traced back to the relationship between spin operators and bosonic operators via the Schwinger mapping [103]. Thus, the bond operators with alternating sign coupling in the bosonic system induce an analogous staggered field interaction. However, typically in spin systems interactions are of local character. In contrast to this, our interactions are global but not trivial, as they are structured. This produces a similar mechanism for the formation of an anti-ferromagnetic like state. However, as we are considering soft-core bosons the analogy is not complete to spins. The dimers in the system are similar to the typical spin singlets of the original VBS [38]. The typical hardcore bosonic representation of the Bose Hubbard Hamiltonian (limit U → ∞) via Matsubara-Matsuda mapping [104] which is an anisotropic Heisenberg model [103] while Now making the above contributions isotropic, only keeping the J B terms (J D = 0) and at fixed particle number we have: where S = (S x , S y , S z ). On the other hand the AKTL Hamiltonian [37] is: Therefore, one can see (17) is a global relative to (18) from which our Hamiltonian is an anisotropic relative without the hard core constraint formally. As such, some similarities might be expected between their ground states. When density coupling and bond coupling act simultaneously (J B = 0 and J D = 0, with J B /J D = 0.25), the situation interpolates in between the above two limits. However the additional bond-density coupling terms have strong effect even for small J B . Bond ordering can take over the system behaviour instead of DW order, see Fig.3(c), (f), (i) and (l). Interestingly, for large on-site interactions (small t 0 /U ) and g eff J 2 D < 0, we find that a state with both DW and bond order occurs while being insulating, the DW does not destroy bond order. When we increase the effective tunneling for g eff J 2 D 0, the system smoothly transitions from the BI+DW state to the SF via a mixture of SSD and SS phases. In contrast, when g eff > 0 and t 0 /U = 0 the system is a DW insulating state. The system transitions smoothly form this state one increases t 0 /U to the SF state via an intermediate SFD phase. In general, bond order takes over and competes with DW order as the ratio J B /J D increases for strong on-site interactions while smoothly transitioning to the SF state as t 0 /U increases via SSD and SS (g eff J 2 D < 0) and superfluid dimer phases SFD (g eff J 2 D > 0). In the current parameter range explored there is no indication of a HI phase for ρ = 1. The phase diagram of the above cases is shown in Fig.4.
It follows to consider the emergent phases at the half filled case (ρ = 1/2). In the system without cavity light, we only have SF phase as there is no gap in the excitation spectra due to incommensuration for the homogenous system. However, even when J B = 0 and J D = 0, the induced symmetry breaking by light will foster the formation of insulators with broken translational and time reversal symmetry. As function of the effective light matter strength for t 0 /U 1 the system has a sharp transition from a DW insulator (g eff J 2 D < 0) to a BI (g eff J 2 D > 0). In the limit of g eff J 2 D 0, the system smoothly goes from a DW insulator to the SS phase, Fig.5(a). In the opposite limit (g eff 0), the system for strong on-site interactions is a BI. As t 0 /U increases, the system smoothly reaches the SF state via an intermediate SFD phase, Fig.5(d). Surprisingly, even with only density coupling, bond ordered phases arise in the large U limit. This can be traced back to the fact that the one-body reduced density matrix has the structure where maximal amplitude MW coherence occurs in this case. This is a consequence of minimising the DW order and the fact that there is incommensuration.
Bond coupling (J B = 0 and J D = 0) at ρ = 1/2 has the effect of stabilising DW ordered phases (DW and SS) when g eff J 2 B > 0. In contrast, the system supports BI+DW . phases for g eff J 2 B < 0 and strong on-site interactions. Even at incommensurate fillings, and addressing through the bonds, the effect of cavity light is a suppression effect upon the SF component. This leads to have a sharp transition from the BI+DW phase to the SF state for g eff J 2 B 0, Fig.5(b). For g eff J 2 B 0, the transition from the DW insulator to the SF state is smoothly connected via SS and later a SFD phase, Fig.5(e).
Similar to the case with only bond coupling, bond ordered phases take over in the case of simultaneous addressing (J B /J D = 0.25). However, density coupling stabilises DW order, while BI+DW phases disappear for g eff J 2 D 0.The DW insulator takes over. The large t 0 /U limit phase for g eff < 0 is a SS and not a homogenous SF, Fig.5(c).The opposite effect occurs in the limit of g eff J 2 D 0, where instead DW ordered phases are strongly suppressed. The state of the system changes from BI+DW to SF via an intermediate SFD phase, Fig.5(f). The phase diagram of the cases considered at halffilling is shown in Fig.6.
It is expected that SS phases will shrink as the number of sites increases for Fig.6 (a) and (c). In general effects in the phase diagram when g eff J 2 D,B > 0 are a manifestation of the quantum optical lattice potential induced by light to the atoms, as light is scattered minimally [14]. The insulating phases (DW, MI) and intermediate phases that appear g eff J 2 D,B > 0 are driven by quantum fluctuations.

Bond order and Synthetic Potentials
In this section we will study the results from simulations of effective Hamiltonians (8) and (9). In the case where synthetic interactions via density coupling are considered, the situation is qualitatively the same for filling ρ = 1 and either for many pumps (8) or many cavities (9) and either Bessel or Morse like potentials. There are slight changes in critical points but overall the phase diagram can be summarised qualitatively in the phase diagram for Bessel-like interactions with many cavities, Fig.7. For light-induced many-body effective interaction strengths such that g eff > g c , with g c ≈ −2.25U the system is in the MI for small t 0 /U . The critical effective tunneling t 0 /U of the SF-MI transition shifts to lower values with respect to the system without cavity light for g eff < 0, while for g eff > 0 the MI is stabilised to larger critical t 0 /U . The situation is very similar to the density coupling case in diffraction minima (J B = 0, J D = 0), with a relevant shift in g c from g c ≈ −U to g c ≈ −2.25U . This can be attributed to the additional mode dependency of the synthetic interaction with R = 4 modes. The additional density modes suppress further the on-site fluctuations stabilising in a larger region of phase space the MI state. For g eff < g c , the density wave instability sets in and competes with dimer order close to g c , O B ≈ O DW . As the interaction strength is further decreased, the system condenses into a DW insulator for t 0 /U 1, where O DW > O B . When g eff 0 the system transitions smoothly to the SF phase via an intermediate SS phase. Bond order is relevant just near the transition to DW ordered states. The intermediate SS phase is suppressed faster as the number of sites increases with respect to the case of Hamiltonian (6).
In contrast to unit filling, the half filled case (ρ = 1/2), offers intriguing behaviour. In the case of light-induced synthetic interactions via cavity modes, we have that the Bessel potential and the Morse potential have different insulating phases for g eff 0 and t 0 /U 1, Fig.8(a) and (b). The Bessel potential supports a DW insulator that smoothly transitions via a SS to a SF state. This is similar to the behaviour seen in the diffraction minima coupling Fig. 6 (a). However, here the system reaches the SF state and is not SS for large t 0 /U . In contrast, the Morse potential supports a BI+DW phase, even though that here there are only density-density interactions. Different from the above, the situation for g eff 0 is qualitatively similar in both potentials, Fig.8(a) and (b). Here the system for t 0 /U = 0 is in a BI and a continuous transition occurs via a very narrow intermediate SFD (not shown) eventually reaching the SF for t 0 /U 0. When light-induced interactions are constructed via pump modes Fig.9, a qualitative similarity between phase diagrams for both potentials occurs with respect to Fig.8. In general, we find suppression of bond order with respect to the cavity mode case, the transition to the SF occurs at lower critical t 0 /U . In contrast to the cavity mode case, the pump case is easier to implement experimentally requiring many pump beams compared to a multi-mode cavity or many cavities and the behaviour is qualitatively similar.
The difference in behaviour between the potentials can be traced to the staggered like nature of the Bessel like potential and the smoother pattern given by the Morse like potential. This combines with the position dependent interaction nature of the many pump/cavity configuration, leading to the effects seen. In the case of Bessel like potential, even though simulating a finite range, the effect in the system is dominated by the staggered field like nature of its form, alternating coupling between density modes. Thus, the results are similar to density coupling from the single cavity case in the diffraction minima, a staggered density global interaction. In contrast to this, the effect of the Morse like potentials is manifest by stabilising bond ordered phases. Importantly, in the Morse cases Fig.8(b) and 9(b) one can have BI+DW insulators with density-density interactions. The BI phases for g eff > 0 in all the different cases in this section are driven by the quantum fluctuations induced by quantum optical lattice potential, see the previous section.

Conclusions
In this paper we have shown the interplay between bond ordered states (bond insulators, supersolid dimers and superfluid dimers), density wave ordered states (supersolid and density wave insulators), superfluid and Mott insulators due to strong cavity induced interactions in 1D, via exact diagonalization. We have shown numerical evidence to support the identification of bond ordered states with valence bond states, via the calculation of string and parity order parameters in states with bond addressing. We have investigated the suppression and stabilisation of density wave ordered phases and their competition with bond ordered phases due to different choices of synthetic lightinduced many-body matter interactions. We have found that using multiple cavity modes and multiple pump modes, one can modify the behaviour of the supported quantum many-body phases in the system. We have found that one can induce bond ordering even by density addressing. In general, the interplay of the BH model with the cavity induced interaction can change the nature of the quantum phase transitions that appear in the system. These can be either discontinuous or continuous depending on the design of the spatial profile of the interactions, modifying the typical scenario of the SF-MI transition (BKT type in 1D). We have shown the typical quantum manybody phases the system can support and how do these compete as relevant experimental parameters are changed and different geometries of light are chosen: tunneling, on-site interaction and the effective light induced interaction strength. The interplay between different order parameters demonstrates the connection between the designed light induced interactions and the supported quantum many body-phases in the system. This provides a rich landscape of phases to explore experimentally with intriguing properties.
Our results support the possibility to use synthetic many-body matter interactions for the quantum simulation of analogous strongly correlated states related to quantum magnetism from condensed matter. Our findings also suggest, the possibility to use them in the study of the fermionic variant of our system to perform quantum simulations of other interesting states such as, RVB states, which are relevant in mechanisms related to the on set of Hi-temperature superconductivity in real materials. Beyond ultracold atoms in cavities, our results can extended to other arrays of naturally occurring, synthetic, hybrid systems or solid state devices with quantum degrees of freedom with strong light-matter coupling in low dimensions. The results in principle can be applied to systems of fermions, spins, molecules (including biological ones) [105], atoms in multiple cavities [106], ions [107] and semiconductor [108] or superconducting qubits [109]. The setup we study might aid in the design of novel quantum materials, where the concepts that we describe can be exploited for the design of properties of real materials and composite devices in solid state systems with strong light-matter coupling. This opens a new chapter on what can be achieved by quantum simulation via atomic systems.