Photon-Assisted-Tunneling Toolbox for Quantum Simulations in Ion Traps

We describe a versatile toolbox for the quantum simulation of many-body lattice models, capable of exploring the combined effects of background Abelian and non-Abelian gauge fields, bond and site disorder, and strong on-site interactions. We show how to control the quantum dynamics of particles trapped in lattice potentials by the photon-assisted tunneling induced by periodic drivings. This scheme is general enough to be applied to either bosons or fermions with the additional advantage of being non-perturbative. It finds an ideal application in microfabricated ion trap arrays, where the quantized vibrational modes of the ions can be described by a quantum lattice model. We present a detailed theoretical proposal for a quantum simulator in that experimental setup, and show that it is possible to explore phases of matter that range from the fractional quantum Hall effect, to exotic strongly-correlated glasses, or flux-lattice models decorated with arbitrary patterns of localized defects.

In general, quantum many-body systems cannot be understood by extrapolating the properties of their individual constituents [1], but rather by considering their rich collective behavior. A prototypical example is that of sound waves in solids [2], where the atoms do not move individually but vibrate collectively and give rise to propagating quasiparticles, the so-called phonons. This particular many-body problem is exceptional since the properties of the phonons can be calculated exactly, something that rarely occurs in quantum manybody systems. The usual paradigm is that the complexity of the model grows very fast with the number of particles, making both analytical and numerical methods more involved and less efficient. A radically different approach is based on the so-called experimental quantum simulations, which were envisaged several decades ago by R. P. Feynman [3], and have now evolved into a discipline that merges concepts from atomic physics, quantum optics, quantum-information science and condensed-matter physics. The central idea of a quantum simulation is to manipulate the microscopic properties of a particular experimental setup in a way that it reproduces faithfully a quantum many-body model under study. In this way, nature itself computes the properties of the model, and our measurements yield the answer to questions such as the nature of the ground-state and collective excitations, determining whether a dedicated model is sufficient to describe the relevant properties of a particular system. Accordingly, quantum simulations have the potential of solving fundamental open problems in physics, ranging from high-temperature superconductivity, to the thermalization of closed quantum systems, or the properties of spin glasses.
Imagine for a moment that the phonons in a solid were controllable to such an extent that their Hamiltonian could be engineered to target a many-body model of interest. For instance, by shaping the crystal anharmonicities, one could tune the phonon-phonon interactions and reproduce the physics of strongly-correlated models. Unfortunately, it is hard to develop techniques that allow such microscopic control in solidstate materials, and it might be beneficial to search for different experimental platforms. This exotic idea may become realized in experiments with Coulomb crystals of cold atomic ions in radio-frequency traps [4]. Current experimental tools allow for a promising control of the vibrational excitations at the quantum level, such that phonon-based quantum simulations of many-body physics may become a reality [5,6].
So far, the most exploited property of the vibrational excitations of trapped ions is their ability to transfer information between distant ions, which encode quantum bits (qubits) in their electronic states. This provides a mechanism to perform quantum logic operations between distant qubits, a fundamental building block of quantum-information protocols [7]. This type of two-qubit couplings can be interpreted as a spinspin interaction where the qubit plays the role of a pseudospin (henceforth referred as spin) [9] and can be exploited for quantum-simulation purposes. From this perspective, the phonons act as mediators of the interaction and give rise to a wide range of spin models [10], which are familiar in the field of quantum magnetism. The experimental success of these spin-based quantum simulations, either in an analog [11] or digital version [12], has motivated a variety of proposals that range from neural networks [13], to three-body interactions [14], frustrated magnetism [15], mixed-spin models [16], or topologically-ordered spin models [17].
Instead of using the phonons as a gadget to obtain the desired quantum-spin simulator, one can substantially enhance the capabilities by reclaiming phonons as the building blocks for the quantum simulation of many-body models [5]. The local vibrational excitations of each atomic ion, considered to be trapped individually, correspond to bosonic quasiparticles, and the Coulomb interaction is responsible for the interchange of these bosons between distant ions. To build interesting bosonic quantum simulators, one must complement the aforementioned scheme with additional ingredients, such as strong trapping non-linearities that lead to Mott insulating [5] and frustrated phases [18], or incommensurate trapping potentials leading to different quantum phase transitions [19]. Besides, the spins of the ions can be used as a tool to widen the applicability of this quantum simulator, which can potentially target the spin-boson model [20], the lattice Jaynes-Cummings model [21], or the phenomenon of Anderson localization due to disorder [22]. In addition to these many-body quantum simulations, there has also been a recent activity in building analogues of single-particle phenomena with a special emphasis on relativistic effects (see e.g. [23]). In this case, even if the models are tractable on a classical computer, the predicted effects are hard to access in the original experiments, and thus justify the effort in building a trapped-ion analogue.
In this article, we introduce a versatile toolbox for the quantum simulation of many-body lattice models. This toolbox is based on the phenomenon of photon-assisted tunneling (PAT) of phonons in ion traps [24], but can also be applied to different systems of bosons, and even fermions, in a lattice. As shown in this manuscript, the PAT effect has a a variety of facets that can be exploited to provide new paradigms for the aforementioned many-body quantum simulation. The idea underlying the PAT is that the tunneling of particles between the wells of a periodic lattice can be assisted by inducing resonances that correspond to the absorption/emission of photons out of an EM-field providing a periodic driving force. Even though this idea was initially introduced for condensed-matter systems [25] (see [26] and references therein), it has also been applied in the field of ultracold atoms in optical lattices [27]. Here, the analogy is straightforward since the neutral atoms can tunnel between the adjacent wells of a driven periodic potential created by light. In this work, we show how to exploit the PAT in order to induce synthetic gauge fields in the aforementioned many-body lattice models.
In the context of trapped ions, the possibility of controlling the tunneling of the vibrational excitations between ions stored in two separate potential minima, aligned within a linear radio-frequency trap, has been recently demonstrated [28]. The exchange of phonons between two ions can be controlled and optimized by matching their individual trapping frequencies. Building on these seminal experiments, we have demonstrated theoretically that driving the trapping frequencies in/out of resonance realizes a novel version of the PAT paradigm [24]. Rather than relying on a periodic force [25,27] scaled to a larger number of ions, which may be quite demanding in trapped-ion experiments, our scheme relies on the periodic modulation of the trapping frequencies. In particular, it requires an additional relative phase between the individual modulations that can be experimentally tuned. This allows for the full control of both the amplitude and the phase of the tunneling of phonons, a result that becomes specially interesting for two-dimensional arrays of micro-fabricated traps [29]. Here, the technology of ion-trap micro-fabrication provides a means of assembling arrays of ion traps in any desired geometry [30]. The combination of the capabilities for the design and fabrication of arrays of microtraps with the tool of PAT [24] opens a vast amount of possibilities for a quantum simulator of many-body physics. These range from the realization of bosonic models subjected to Abelian and non-Abelian synthetic gauge fields, to bond and site disorder leading to glassy phases, or to decorated flux lattices which can be related to anyonic excitations. In this article, we provide a detailed analysis of such PAT in arrays of microtraps.
This article is organized as follows. The main results are summarized in Sec. IV, so that this section can be consulted for a general overview. The general scheme of PAT for any lattice model is introduced in Sec. II, where we demonstrate how the synthetic gauge field arises from the assisted tunneling. The application as a toolbox for quantum simulations with trapped ions is contained in Sec. III. Here, we present a thorough description of the different many-body models that can be explored by exploiting the peculiarities of trapped ions, which widen the versatility of our quantum simulator.
In this Section, we describe the concept of PAT for a general system, and discuss how it can be exploited for the quantum simulation of lattice models incorporating synthetic gauge fields. Once the main ingredients are identified, we focus on the case of ions in micro-fabricated traps in Sec. III, and show how to implement the PAT for the vibrational excitations. Let us remark that the generality of this section may also find applications in different platforms, such as photons in arrays of cavities in circuit QED [31], or ultracold atoms in optical lattices [32]. We also note that the PAT scheme is not restricted to bosons, but works equally well for fermions. This may alleviate some of the difficulties that arise in the simulation of synthetic gauge fields via Raman-assisted tunneling of fermionic atoms in spin-dependent optical lattices [33]. Figure 1: Photon-assisted tunneling scheme: (a) Scheme for the driven one-dimensional tight-binding model in Eqs. (1) and (4). The tunneling of the particles between neighboring sites r 0 j → r 0 i , which is initially suppressed by the large gradient ∆ω σ J σ t;ij , becomes assisted by a resonant periodic modulation of the on-site energies with amplitude η d ω d . (b) Two-dimensional scheme giving rise to a synthetic gauge field such that the path around an elementary plaquette W ∝ e iφ mimics the Aharonov-Bohm phase, which is picked up by a charged particle following the path around a plaquette pierced by an external magnetic field orthogonal to the lattice.
We consider a tight-binding model describing the tunneling of particles, either bosons or fermions, between the sites of an underlying two-dimensional lattice. The lattice sites are characterized by a vector of integer numbers i = (i 1 , i 2 ), such that where {e α } are the unit vectors spanning the lattice, and {d α } are the corresponding lattice constants. The particles are represented by bosonic or fermionic creation-annihilation operators a † σ ,i , a σ ,i where σ labels some additional degrees of freedom. The dynamics of the system is described by the following Hamiltonian where ω σ ,i stands for the on-site energy (h = 1), and J σ t;ij is the tunneling amplitude of the particles between different lattice sites j → i, which usually depends on the distance between sites such that J σ t;ij = J σ t (|r 0 i − r 0 j |). The tunneling constraint i > j refers to an ordering of the sites, such that i 1 > j 1 , or i 2 > j 2 if i 1 = j 1 . Two additional ingredients are required: i) Gradient of the on-site energies: The on-site energies have the following expression ω σ ,i = ω σ + ∆ω σ i 1 . Here, ω σ is a constant energy offset, and ∆ω σ results due to a gradient along one of the lattice principal axes satisfying ii) Periodic modulation of the on-site energies: The on-site energies must be supplemented by the periodic modulation where ω d is the driving frequency, and η d ω d the driving amplitude. Note that the periodic driving incorporates a site-dependent phase φ i that shall play a crucial role for the PAT, and is described by These two ingredients, schematically represented in Fig. 1(a), are incorporated in the above description by modifying the Hamiltonian H 0 → H 0 (t) as follows Let us emphasize again that the standard formulation of the PAT relies on a periodic force acting on the particles residing on the lattice sites [25,27]. In contrast, our scheme considers a periodic modulation of the on-site energies, which turns out to be better suited to taylor the tunneling amplitudes in analogy with a background synthetic gauge field. In subsection II A, we present an analytical model, which is confronted to numerical simulations in II B.
A. Analytical description

Dressed tunneling
In this subsection, we derive a compact expression for the PAT strength, and study its dependence on the periodicdriving parameters. Let us express the tunneling Hamiltonian H t in the interaction picture with respect to H 0 (t), namely In this picture, the annihilation operators fulfill which leads to the following relation between both pictures At this point, we note that the tunneling Hamiltonian H t is invariant under U(1) gauge transformations, and thus, the last term in the above expression can be trivially gauged away a σ ,i → e −iη d sin(φ i ) a σ ,i . Accordingly, the Hamiltonian becomes H t (t) = ∑ i>j J σ d;ij (t)a † σ ,i a σ ,j + H.c., where J σ d;ij (t) = J σ t;ij Θ(t), and the time-dependence is encoded in the function In this expression, one readily observes that the tunneling amplitude becomes dressed by the photons of the periodic driving J σ t;ij → J σ d;ij (t). Besides, the fundamental role of the phase of the periodic driving also becomes apparent: only when φ i = φ j , the tunneling becomes assisted J σ d;ij (t) = J σ t;ij . In order to proceed with this analytical treatment, we use the identity where J s (η d ) are Bessel functions of the first class [34]. Hence, the dressed tunneling becomes a sum of terms with different time dependences. Let us first focus on the tunneling along the gradient, i 1 > j 1 , which can be expressed as follows (9) By tuning the driving frequency to ω d = ∆ω σ /r, where r is some positive integer, the above expression contains two different types of terms. There are resonant terms fulfilling whereas the remaining terms are far off-resonance. By applying a rotating-wave approximation (RWA) for J σ t;ij ∆ω σ (2), we neglect the rapidly-oscillating terms yielding According to this expression, both the amplitude and the phase of the tunneling are controlled by the periodic driving parameters. As shown in Sec. III, the periodic driving can be achieved by means of optical dipole forces on the ions, and the tunneling along the e 1 -axis is thus assisted by the absorption/emission of r(i 1 − j 1 ) photons from the EM-field providing the driving force, hence the name PAT. A similar analysis for the tunneling orthogonal to the gradient, i 1 = j 1 , shows that the resonant terms satisfy s = s. The complete assisted-tunneling Hamiltonian becomes with the dressed couplings expressed as follows where the function F is responsible for the modulation of the tunneling amplitude and we have defined the phase difference ∆φ ij = φ i − φ j , and the function f (i, j) = r(i 1 − j 1 ). Let us emphasize one of the important properties of this photon-assisted scheme, namely, its non-perturbative character. So far, we have not used any restriction on the values of η d or φ i , so that the dressed tunneling J d may be of the same order as the bare one J t . In Fig. 2, we represent the modulation function F f (i,j) (η d , η d , ∆φ ij ) for different periodic-driving parameters, and different ranges of the tunneling. In particular, Figs. 2(a)-(c) represent the tunneling along the direction of the gradient for r = 1, whereas Fig. 2(d) stands for the tunneling orthogonal to the gradient. As derived from Fig. 2(a), for η d ≈ 1, ∆φ = π, one gets the maximal amplitude for the assisted tunneling between nearest neighbors. In particular, we find |F 1 | ≈ 0.6 for these optimal values, which thus supports our previous claim about the non-perturbative nature of the scheme (i.e. |J σ d;i+e 1 ,i | ≈ 0.6|J σ t;i+e 1 ,i |). It is also interesting  Figure 2: Photon-assisted modulation of the tunneling amplitude: Contour plot of the modulation amplitude |F r(i 1 − j 1 ) | for the tunneling between sites i → j, such that i 1 > j 1 , as a function of the driving parameters η d , ∆φ ij for r = 1. (a) First-neighbor assisted tunneling (i → i + e 1 ), (b) Second-neighbor assisted tunneling (i → i + 2e 1 ), (c) Third-neighbor assisted tunneling (i → i + 3e 1 ), (d) First neighbor assisted tunneling orthogonal to the gradient i → i + e 2 .
to note that, as we consider longer range terms, their amplitude gets gradually diminished, as shown in Figs. 2(b)-(c). In Fig. 2(d), we represent the modulation amplitude in the direction perpendicular to the gradient, which taking the same parameters as above leads to |J σ d;i+e 2 ,i | ≈ 0.2|J σ t;i+e 2 ,i |. Before concluding this subsection, we briefly comment on another interesting property of this scheme. Let us consider that the original Hamiltonian (1) also contains on-site particleparticle interactions H +V , where and the dependence of the interaction strengths U σ σ on the internal indices depends on the bosonic/fermionic nature of the particles. It is straightforward to see that the gradient and periodic driving do not modify V . Therefore, the PAT also works in the presence of on-site interactions. It would be interesting to study how the scheme can be used to modify long-range interactions, but this lies beyond the scope of the present work.

Synthetic gauge fields
In this subsection, we demonstrate that it is possible to interpret the phase of the dressed tunneling (13) as if it was originated by a background gauge field. In particular, we show that whenever particles tunnel along a closed path in the lattice, they pick up a non-vanishing phase analogous to the celebrated Aharonov-Bohm phase for charged particles in electromagnetic fields [35]. The consecutive tunneling of a particle around a unit plaquette of the lattice i → i+e 1 → i+e 1 +e 2 → i + e 2 → i [ Fig. 1(b)] is formally given by Using Eqs. (7) and (13), it can be expressed as which leads to W (1) = |W (1) |e iφ , where φ = rφ 2 only depends on the component of the periodic-driving phase that is orthogonal to the gradient. This quantity is a gauge-invariant observable proportional to the so-called Wilson loop in lattice gauge theories [36], which can be expressed as follows Here, we have introduced an effective charge e * independent of the charged/neutral character of the particles, the synthetic gauge potential A s , and the synthetic gauge field B s . Independently on how we rearrange the phases of the tunneling strengths locally, the enclosed phase φ = e * B s · dS shall always be left invariant, and can be thus interpreted as the magnetic flux of a synthetic magnetic field B s that pierces the lattice. In order for this analogy to be complete, one should consider carefully the long-range character of the tunneling, according to which particles can follow different closed paths. We shall explore the two possible paths around the second smallest plaquettes, namely Repeating the above calculations, we find that Hence, the Aharonov-Bohm phase is doubled with respect to the unit plaquette, which is consistent with the fact that the enclosed area is also doubled for these paths. The same occurs for any other closed path, and thus the analogy to a background gauge field is consistent with the possible long-range of the tunnelings. At this point, it is worth emphasizing the generality of the scheme hereby proposed. It works for both bosons and fermions, for any range of the tunneling J σ t;ij = J σ t (|r 0 i − r 0 j |), it can incorporate local particle-particle interactions, and it is non-perturbative. In the second part of this manuscript, we shall specify this scheme to trapped-ion experiments, which provide an ideal platform where to realize a bosonic PAT (i.e. vibrational excitations). Before moving onto the numerical verification of these results, let us introduce an alternative and more compact formulation of the synthetic gauge fields. It is possible to rewrite the PAT Hamiltonian in analogy to the so-called Peierls substitution [37], which yields where we have introduced the dressed-tunneling amplitudẽ , and the synthetic gauge potential A s (x) = −B 0 ye 1 , where B 0 = rφ 2 /e * d 1 d 2 . This gauge potential, which corresponds to the famous Landau gauge for electrons in a constant magnetic field, is obtained after the following U(1) gauge transformation a σ ,i → a σ ,i e −iχ i , where  Figure 3: Photon-assisted tunneling along a link: (a) Bare tunneling for an initial state |Ψ 0 = a † 1 |vac with a single particle occupying the first site. The numerical expectation values n 1 (t) (blue circle), and n 2 (t) (orange squares) are compared to the effective description n 1 (t) eff (blue lines), and n 2 (t) eff (orange lines) described in Eq. (21). (b) Suppressed tunneling due to the gradient, and in the absence of periodic driving. (c) Assisted tunneling in the presence of a periodic driving for ∆φ 12 = π. (d) Coherent destruction of tunneling (black arrows) due to a periodic driving in the absence of the gradient. The yellow circles represent the maximum population transferred to site 2 n * following parameters of the driven Hamiltonian in Eqs. (1) and (4): the bare tunneling is J t;12 = 10 −2 ω, the gradient ∆ω = 0.5ω, and the periodic driving parameters correspond to η d = 1, and ω d = ∆ω. Note that we do not consider additional degrees of freedom, avoiding thus the index σ in the following. In Fig. 3, we represent the expectation values n i (t) = Ψ(t)|a † i a i |Ψ(t) that result from the numerical integration of the Schrödinger equation id|Ψ(t) /dt = (H 0 (t) + H t )|Ψ(t) , and compare them with the effective analytical description |Ψ(t) eff = e −iH eff t |Ψ 0 . We consider an initial state with a single bosonic particle in the first site |Ψ 0 = a † 1 |vac , where |vac stands for the vacuum. In this case, the effective description (12) can be integrated exactly, yielding the following periodic oscillations of the particle between the two lattice sites where ω eff = 2|J d;12 | is the frequency of oscillations. In the absence of any gradient or periodic driving, ω eff = 2|J t;12 | = 0.02ω, and thus the particle exchange undergoes oscillations with a period of T t;12 = π/|J t;12 | = 100π/ω, which coincides exactly with the scale shown in Fig. 3(a). In this figure, we reveal that the effective description (circles and squares) agrees with the numerical results (solid lines). Note that in order to carry out the numerical integration, we truncate the bosonic Hilbert space to n trunc = 4. We have confirmed that the truncation error after changing n trunc →ñ trunc = n trunc + 1 lies below | n i (t) n trunc − n i (t) ñ trunc | < 10 −4 . The same occurs in Fig. 3(b), where we switch on the gradient. As a consequence of the mismatch between the on-site energies, the tunneling is completely inhibited. In order to assist it, we switch on the periodic driving in Fig. 3(c) for φ 1 = φ 2 + π, which activates again the periodic oscillations but with a different period T t;12 = π/|J d;12 | = 100π/|F 1 (η d , η d , π)|ω. For the parameters used, we obtain |F 1 (η d , η d , π)| ≈ 0.6 ( Fig. 2(a)), which explains the slightly longer oscillation period. However, we remark that both dynamics occur on the same time-scale as a consequence of the non-perturbative character of the scheme. In Fig. 3(d), we study the maximal population transfer to site 2 due to the periodic driving η d = 0, but in the absence of the gradient ∆ω = 0. We set the phases to φ 1 = 0, φ 2 = π, and study the maximal transfer n * 2 = max{ n 2 (t) : t ≤ t * }, such that t * = 100π/|F 0 (η d , η d , π)|ω for a range of driving strengths η d ∈ [0, 4]. As shown in this figure, there are certain values of the driving strength, marked by black arrows, where the tunneling gets completely suppressed, a phenomenon known as coherent destruction of tunneling [25], or dynamic localization. Let us remark that this scheme leads to a perfect localization of the particles also in a longer onedimensional chain, since the dressed tunneling cancels for all pairs of nearest-neighboring lattice sites simultaneously. Once we set φ i 1 = πi 1 /2, such that ∆φ i 1 i 1 +1 = π, it is possible to find a zero of the modulating function F 0 (η 0 d , η 0 d , ±π) = 0 which is related to the particular values of the Bessel functions. In this limit, the particles shall not diffuse through the lattice.
In Fig. 4(a), we study the photon-assisted process for different driving phases ∆φ 12 ∈ [0, 2π], considering the same pa-rameters and initial state as before. We represent the maximum population transferred to site 2 n * 2 = max{ n 2 (t) : t ≤ t * }, such that t * = 50π/|F 1 (η d , η d , π)|ω corresponds to the optimal PAT with ∆φ 12 = π. The yellow dots correspond to the numerical integration of the truncated Hamiltonian (1) and (4), whereas the solid line corresponds to the analytical prediction evaluated at the related exchange period n 2 = n 2 (t ) eff for the different phases, namely t = 50π/|F 1 (η d , η d , ∆φ 12 )|ω. In this figure, we reveal the agreement between both descriptions for any of the driving phases. Note that this agreement will not be compromised by going to larger arrays. As announced previously, we can interpolate between the optimally-assisted and the completely-suppressed tunneling regimes by modifying the driving phase.
It is also interesting to consider an initial state that does not correspond to a single particle, but rather to a thermal ensemble. We focus now on bosons with independent mean number of particlesn 1 ,n 2 . This state corresponds to ρ 0 = ρ 1 ⊗ ρ 2 with such that the parameters β i are implicitly defined through the mean number of particles following a Bose-Einstein distribu- . In this case, we must integrate numerically the Liouville-Von Neumann equation , and compare the result to the effective description ρ eff (t) = e −iH eff t ρ 0 e +iH eff t . In particular, we consider the same parameters as before, and setn 1 = 0.5,n 2 = 0.25 after truncating the particle Hilbert space to n trunc = 7.
Note that due to the thermal effects, the truncation parameter has been increased with respect to the previous simulations. Again, we have checked the convergence of the results by increasing n trunc . In Fig. 4(b), we represent the expectation values n * 2 , n 2 introduced above, together with those of site 1, n * 1 = min{ n 1 (t) : t ≤ t * }, and n 1 = n 1 (t ) eff for t = 50π/|F 1 (η d , η d , ∆φ 12 )|ω. We conclude from the perfect agreement shown in the figure that the scheme also works for thermal states. In fact, in the absence of interactions, the equations of motion do not depend on whether the state is pure or a mixed thermal state. For the latter, there is a background over which the PAT phenomena will occur. The role of the additional interactions, and its interplay with the thermal states, is an extremely interesting topic that deserves a separate study.
The suitability of the PAT scheme for thermal states will turn out to be essential in Sec. III, where we discuss a realistic implementation of the quantum simulator of gauge fields with phonons in microtrap arrays. In this case, it means that cooling to the vibrational ground-state is not necessary to observe the non-trivial effects of the PAT. Another interesting topic is the occurrence of motional heating in the microtraps. In case this motional heating is global, it should not interfere with the PAT effects. However, if this heating has a local nature, it may lead to very interesting effects that mimic the role played by a reservoir providing electrons to local regions of a metallic conductor. This effect can be complemented by the controlled engineering of dephasing in the phonon dynamics by introducing noisy potentials in the trap electrodes [38]. Hence, the applicability of the PAT quantum simulator could  Figure 4: Phase-dependence and thermal photon-assisted tunneling: (a) Dependence of the PAT on the relative driving phase ∆φ 12 . The yellow circles represent the maximum population transfer to site 2, n * 2 , for a particle initially located at site 1, |Ψ 0 = a † 1 |vac , as obtained from the numerical integration of the full time-dependent Hamiltonian H 0 (t) + H t . The red line corresponds to the same quantity evaluated for the effective description H eff , referred as n 2 . Emphasized by the doted frame is the region of optimal assisted tunneling. (b) Phase-dependence of the PAT for an initial thermal state with n 1 = 0.5,n 2 = 0.25. We also represent the particle transfer to site 1.
be widened, and used to study how transport phenomena is affected by an environment that induces decoherence effects such as motional heating or dephasing.

Photon-assisted tunneling around a plaquette
Let us now apply our formalism to a two-dimensional setup where to test the application of PAT for the quantum simulation of synthetic gauge fields. We consider a model of spinless bosons in a square plaquette, and set the parameters of the total Hamiltonian in Eqs. (1) and (4) to J t;12 = 10 −2 ω, ∆ω = 0.5ω, η d = 1, r = 1, and ω d = ∆ω. Our aim is to test the accuracy of the effective description (20), trying to find clear-cut evidence of the underlying synthetic gauge field. In order to do so, we study a discrete version of the celebrated Aharonov-Bohm effect [35].
Let us briefly recall the underlying interference effect [ Fig. 5(a)]. A particle initially localized at site 1 can tunnel to site 3 following two different paths, namely γ 1 : 1 → 2 → 3, or γ 2 : 1 → 4 → 3. Due to the gauge field, each path leads to a different phase dr·A s a † 3 |vac , such that the probability of performing such trajectory is (18). From this expression, one finds a perfect destructive interference that forbids the particle to tunnel to site 3 when the enclosed phase φ = π. This Aharonov-Bohm interference will serve us as a testbed for the validity of the gauge-field analogy (20).
We now compare the dynamics obtained from the numerical integration of the driven Hamiltonian (1) and (4), to the effective description in Eq. (20). In Figs. 5(b)-(c), we consider |Ψ 0 = a † 1 |vac , and study the propagation of the particle along the square plaquette for different driving phases. In Fig. 5(b), we observe that for the phases φ 1 = π, φ 2 = 0, the particle is allowed to tunnel to every lattice site. This is consistent with the fact that for φ = φ 2 = 0, there is no interference. Conversely, we set φ 1 = π, φ 2 = π in Fig. 5(c), where one readily observes that the tunneling to site 3 is forbidden by the aforementioned Aharonov-Bohm interference. We note that for the trapped-ion case (see the sections below), the diagonal path going from sites 1 → 3 and 2 → 4 must also be accounted since the tunneling amplitude need not be small. In [24], we discussed how to cancel it by playing with the phases, so that the analogy with the Aharonov-Bohm effect is valid. In Fig. 5(e), we corroborate that the effect also holds for thermal states. Finally, we check in Fig. 5(d) that the perfect destructive interference only occurs for the so-called π-flux phase. In this figure, we represent the maximal transferred population to site 3, n * 3 = max{ n 3 (t), 0 < t < t * }, where t * = 100π/|F 1 (η d , η d , π)|ω, and show that the perfect interference only occurs when the enclosed phase φ = π. This allows us to rule out other possible sources of interference, and conclude that it is only due to the synthetic gauge field.
Let us close this subsection by highlighting the accuracy of our analytical treatment, which has been confronted to a numerical survey of different setups with a wide range of parameters. The rest of this manuscript will be devoted to analyze the experimental setups in the field of quantum optics, ions in microtrap arrays, where the necessary ingredients for this PAT are in reach with state-of-the-art technologies. Besides, we shall also introduce novel methods to to exploit additional aspects of PAT phenomenology.

III. PHOTON-ASSISTED-TUNNELING TOOLBOX FOR TRAPPED IONS
In this section, we apply the above scheme to arrays of micro-fabricated ion traps, and show that one can achieve PAT of the vibrational excitations between distant microtraps. In order to set the notation, we start by showing in Sec. III A that the vibrations of an ensemble of ions in a microtrap array can be described in terms of tunneling phonons [5]. Then, we describe in Sec. III B how the gradient and periodic driving of the trapping frequencies can be achieved using the tools of state-of-the-art experiments in ion traps (see the reviews [46]). In particular, the gradient can be implemented by the local control of DC-voltages applied to the trap electrodes, whereas the periodic driving stems from an optical dipole force. We derive a set of constraints that this dipole force must fulfill, and discuss how these requirements can be met for realistic experimental parameters. In Sec. III C, we incorporate phonon-phonon interactions in the phonon-based quantum simulator of bosonic particles exposed to synthetic gauge fields. Finally, in Secs. III D, III E, and III F, we describe how to extend the PAT scheme beyond the applications discussed so far, reaching more involved models where the spin of the ion can be exploited as an additional tool.

A. Tunneling of phonons in microtrap arrays
Let us consider an ensemble of N ions with charge e and mass m, labelled by integer numbers i = (i x , i y ). Each ion is subjected to the electric potential energy created by a planar micro-fabricated trap V t , and the Coulomb energy potential V c causing the repulsion between the remaining ions, such that the total Hamiltonian is Here, the trap is designed in such a way that the ion equilibrium positions, which are given by ∇ ∇ ∇(V t + V c )| r 0 i = 0, form a regular lattice r 0 i = i 1 d 1 e 1 + i 2 d 2 e 2 , where e 1 , e 2 are the unit vectors, and d 1 , d 2 the lattice spacings. At sufficiently small temperatures, the excursions of the ions from these equilibrium positions r i = r 0 i +δ r i are small with respect to the lattice spacing. Thus, for our purpose, we are allowed to consider the above Hamiltonian up to second order only [2,47]. This approximation yields a system of coupled harmonic oscillators with the following couplings where r 0 i−l = r 0 i − r 0 l , δ il stands for the Kronecker delta, and α, γ = x, y, z refer to the main axes of the trapping potential. Note that the micro-fabricated trap gives rise to a confinement that can be considered to be harmonic for the motional amplitudes considered here, and is therefore characterized by the frequencies ω α,i that may depend on the axis and the ion position within the lattice. The Hamiltonian (24) can be expressed in the basis of local quantized vibrations where b † α,i , b α,i are the bosonic creation-annihilation operators, and we seth = 1. In this local basis, the vibrational Hamiltonian becomes H = H 0 + H c , In this work, we consider that the microtraps fulfill which allow us to neglect the rapidly oscillating terms of the above Hamiltonian using a standard rotating-wave approximation (RWA). The first condition allows us to consider independent vibrations along each of the trapping axes, whereas the second one allows us to neglect the terms that do not conserve the number of vibrational excitations. Accordingly, where the trapping frequency of each ion is slightly modified by the electrostatic interaction with its neighboring ions ω α,i = (ω 2 α,i +V αα c;ii /2m) 1/2 , and we have defined J α c;ij = 2J αα c;ij . This Hamiltonian can be interpreted as a tight-binding model for the local phonons, which tunnel between distant microtraps according to the dipolar tunneling strengths J α c;ij . By a direct comparison to the general tight-binding model in Eq. (1), we can identify the index σ with the vibrational axis, the tunneling strength J σ t;ij with the aforementioned dipolar Coulomb couplings, and the on-site energies ω σ ,i with the microtrap trapping frequencies. In the following subsection, we describe how to obtain the gradient and the periodic driving using stateof-the-art tools in trapped-ion experiments.

B. Gradient and periodic driving of the trapping frequencies
In microtrap arrays, it is possible to design the individual trapping frequencies {ω α,i } by the control of DC-voltages applied to local micro-fabricated electrodes. In fact, the experiments [28] have made explicit use of this property in a conventionally segmented linear rf-trap, in order to switch on/off the phonon tunneling by tuning the trapping frequencies of two neighboring traps in/out of common resonance [29]. We chose the trapping frequencies distributed according to a gradient ω α,i = ω α + ∆ω α i 1 . Considering mutual trap and mutual ion distances > 10µm due to current constraints in fabrication. The correction to the trapping frequency due to the electrostatic interaction can be neglected, and we get the desired gradient of the on-site energies ω α,i ≈ ω α,i = ω α + ∆ω α i 1 which must be incorporated to the tight-binding Hamiltonian (30).
We now turn to the periodic driving of the trapping frequencies, which shall be implemented by an optical dipole force [46]. We consider a pair of laser beams with frequencies ω 1 , ω 2 , and wavevectors k 1 , k 2 , which couple to the electronic states | ↓ i , | ↑ i , |a i (see the energy-level scheme of Fig. 6 for the particular case of 25 Mg + ions). Encoding the spin degree of freedom into a pair of electronic ground states by exploiting their hyperfine or Zeeman splitting, the energy difference ω 0 /2π between the states | ↓ i , | ↑ i lies in the microwave range 1 GHz, whereas the energy splitting ω a /2π to an auxiliary state |a i is much larger ω a ω 0 , and typically lies in the optical domain 100 − 1000 THz. By tuning the beatnote of the phase-locked laser beams to a particular value, it is possible to choose out of a variety of couplings. For instance, when ω 1 −ω 2 ≈ ω 0 , such that ω 1 , ω 2 ω a , one speaks about a stimulated two-photon Raman transition between the ground and excited states [ Fig. 6(b)]. Conversely, when ω 1 − ω 2 ≈ ω α ω 0 , one obtains a running-wave realization of the so-called spin-dependent dipole forces. In this work, we are interested in yet a different regime ω 1 − ω 2 ω α [ Fig. 6(c)], where we get the desired periodic driving by the laser dipole force. We  The ground-state manifold 2 S 1/2 is split into the hyperfine levels F = 2, 3, such that an additional magnetic field allows us to isolate two magnetic sublevels | ↑ i = |2, 2 i , | ↓ i = |3, 3 i with an energy difference of ω 0 /2π ≈ 1.8 GHz. Note also that these levels are dressed by a vibrational ladder with equidistant spacing ω α /2π ≈ 1 − 10MHz. The excited manifold 2 P 3/2 contains an auxiliary level to implement a stimulated Raman scheme with two laser beams and a beatnote ω L = ω 1 − ω 2 . (b) Raman scheme for two laser beams with Rabi frequencies Ω 1 , Ω 2 in the so-called red-sideband regime ω L ≈ ω 0 − ω α . In this regime, excitations are exchanged between the electronic and vibrational states of the ions. (c) Raman scheme for the periodic modulation of the trapping frequencies ω L ω α . In this limit, only negligible coupling between the internal or vibrational states take place during the transition, but rather a periodic modulation of the trapping frequencies is realized.
consider the ion-laser interaction in the dipolar approximation where Ω Since we are assuming that the laser frequencies are far offresonant with respect to the auxiliary state, it is possible to perform an adiabatic elimination of this state to simplify the dynamics between | ↑ i , | ↓ i . Furthermore, for sufficiently small Rabi frequencies and considering that ω 1 − ω 2 ω 0 , it is possible to obtain where we have defined the beatnote ω L = ω 1 − ω 2 , ∆k = k 1 − k 2 , and ∆ = ω a − ω 0 − ω 1 such that |∆| ω 0 . Note that we have neglected the AC-Stark shifts arising from each single laser beam, since they do not play any role for the periodic driving of the trapping frequencies. By a proper choice of the laser intensities, detunings, and polarizations, one finds a regime where the two-photon Rabi frequency is a↑ /2∆, and thus the laser-ion interaction becomes ω α /2π ∆ω α /2π ω L /2π J α c;12 /2π Ω L /2π η α 1-10 MHz 50-500 kHz 50-500 kHz 1-10 kHz 0.1-1 MHz 0. 1-0.4 This expression corresponds to a Stark shift that acts equally on both electronic states, and is caused by the crossed laser beams of different frequencies.
To obtain the desired periodic driving from the ion-laser Hamiltonian (33), we express the ion position in terms of the local phonon operators r i = r 0 i + ∑ α e α (b † α,i + b α,i )/ √ 2mω α . By introducing the so-called Lamb-Dicke parameter η α = e α · ∆k/ √ 2mω α 1, we can Taylor expand the Hamiltonian Note that the first term corresponds to an irrelevant c-number that does not modify the phonon dynamics. In order to find the relevant contribution of the two remaining terms H L , we switch to the interaction picture with respect to H 0 = ∑ α,i ω α,i b † α,i b α,i . If we consider that both, the frequency gradient and the laser frequency, fulfill ∆ω α , ω L ω α , |ω α − ω γ | α =γ , a rotating wave approximation for η α Ω L ω α allow us to neglect the rapidly-oscillating terms and leads to the desired periodic driving of the trapping frequencies By direct comparison with Eq. (4), we find the following, individually tunable, mapping between the parameters of the periodic driving and the laser beams: i . Let us also note that in order to neglect the higher order terms that rotate with the laser frequency ω L , we impose that Ω L η 4 α ω L . It is important to remark that even if we are considering small Lamb-Dicke parameters η α 1, the periodic driving |η d | = Ω L η 2 α /ω L needs not to be small by selecting the appropriate laser and Rabi frequencies, so that the non-perturbative character of the PAT can be exploited.
The remaining task to demonstrate that the scheme of PAT works for phonons in microtrap arrays is to consider the typical values for the experimental parameters, and discuss whether they satisfy the constraints imposed during the above derivation. The orders of magnitude of all the relevant parameters are listed in Table I, which satisfy the different constraints made along this derivation, namely For the last inequality, it suffices to focus on the phonon modes transverse to the microtrap plane, α = z, such that |ω z − ω γ | γ =z /2π ≈ 1MHz. If the in-plane vibrational modes are to be used, the anisotropy of the trapping frequencies ω x = ω y should be considered in the microtrap design. Finally, table I allows us to estimate the parameters of the model Hamiltonian, and reveal whether the constraints for the PAT in Eqs. (2) and (3) are fulfilled. We find that J α c;12 /2π ≈ 1-10kHz ∆ω α /2π ≈ 50-500kHz, which thus satisfies the constraint in Eq. (2). Besides, we find that φ i = φ 1 i 1 + φ 2 i 2 , where φ α = −(∆k · e α )d α , and thus the constraint (3) is also satisfied. Therefore, we can conclude that the required ingredients for the PAT introduced in Sec. II can already be met with the current technology of microtrap arrays if the tunneling rates outrun the decoherence rates. Note that if motional heating turns out to degrade the results severely, we could mitigate this problem by the using a cryogenic setup.
Let us finally comment on two additional constraints different from (2) and (3), which are particular to the trapped-ion setup. Both the gradient and the the periodic driving should not modify the stability of the ion crystal, and thus must be smaller than the trapping frequencies This is also fulfilled for the parameters in Table I. Let us note that this condition sets a limit to the scalability of our proposal. For gradients ∆ω α /2π ≈ 50kHz, and trapping frequencies ω α /2π ≈1MHz, an array of 10×10 microtraps is still consistent with the maximum attainable trapping frequency. We note that the gradient can be reduced further (10kHz) without compromising the efficiency of the PAT scheme, and leading to larger systems with N =2500 microtraps. Finally, we would like to remark that the scheme could be scaled even further by considering a local gradient that only affects a few sites, and repeats periodically along one axes of the microtrap array. In order to assist the tunneling between the sites where the gradient is changed, while maintaining a homogenous synthetic flux, an additional periodic driving with a suitable frequency and phase must also be introduced.

C. Synthetic gauge fields and phonon-phonon interactions
The discussion of the PAT of phonons has focused on the quadratic Hamiltonians corresponding to the Coulombinduced tunneling (30), and the periodic driving of the trap frequencies (35). However, as discussed in Sec. II, nonlinearities corresponding to the on-site interactions in Eq. (15) can be incorporated without modifying the assisted-tunneling scheme. In fact, this broadens the applicability of our manybody quantum simulator, since it also targets models of strongly-interacting particles. To obtain such phonon-phonon interactions [5], strong non-linearities in the trapping potential are required. Remarkably, one can exploit a different realization of the above dipole forces (33) in order to give rise to such non-linearities. We denote the parameters of these new dipole forces by a wiggled bar in order to differentiate them form the dipole forces leading to the periodic driving.
We consider that two additional laser beams leading to a dipole force have the same frequency, in a way that the beatnoteω L = 0, and Eq. (33) corresponds to a standing wave. In analogy to the previous expansion (34), we consider that Ω Lηα ω α so that all the terms that do not conserve the number of phonons can be neglected. However, sinceω L = 0, the quadratic terms now correspond to a small shift of the trapping frequencies which has no effect sinceΩ Lη 2 α ∆ω σ ω α . The most relevant contribution is due to the quartic terms which corresponds exactly to the on-site interactions introduced in Eq. (15). In Table II, we summarize the mapping of the phonon Hamiltonian onto the original periodically driven tunneling Hamiltonian of Sec. II. As discussed there, since this interaction is purely local, it shall not be modified by the scheme of PAT. Therefore, the effective phonon Hamiltonian becomes H eff = K eff +V eff , where where the tunneling isJ α d;ij = J α t;ij F f (i,j) (η d , η d , ∆φ ij ), and the gauge field can be expressed as A s (x) = −B 0 ye 1 , such that B 0 = (rφ 2 /e * d 1 d 2 ). Equation (39) incorporates the central result of this section, which shows that the PAT in microtrap arrays leads to a quantum simulator of a long-range Bose-Hubbard model [42] under additional gauge fields [43].
There are several interesting regimes for this quantum simulator. In the non-interacting limit for a square microtrap array, it yields a bosonic dipolar version of the so-called Azbel-Harper-Hofstadter model [39]. The nearest-neighbor model has been studied thoroughly during the last decades, and contains several interesting properties that range from its relation to topological numbers, to the fractal and self-similar properties of its energy spectrum and wavefunctions, the existence of gapless edge excitations, or the so-called π-flux phases [40]. The addition of the long-range dipolar tunnelings introduces a new feature in the model that, to the best of our knowledge, has not been studied previously and may modify the above phenomena. Besides, most of these effects rely on a magnetic flux per plaquette on the order of the flux quantum, which cannot be achieved in solid-state materials assuming realistic magnetic fields. On the contrary, our proposal has the potential to reach these regimes since it is non-perturbative and the flux can attain arbitrary values φ ∈ [0, 2π].
Since it is possible to build any desired microtrap geometry [30], our quantum simulator can explore the physics of bosonic ladders subjected to synthetic gauge fluxes. Besides, the capability of tuning the fluxes, together with the independent control of the tunneling strength along/across the ladder rungs, dives into the phenomena of flat-band physics and edge states [41], which is typical for fermionic topological insulators that break the time-reversal symmetry [45].
Another interesting regime corresponds to |J α d;ij | ≈Ũ i,αγ , where the interactions compete with the kinetic energy and induce strong correlations in the Bose-Hubbard model [42]. With respect to the neutral-atom realizations [44], the phonon model includes the effects of longer-range tunnelings, and the possibility to address site-dependent interactions. Besides, thanks to the tunability of the functions F f (i,j) [Fig. 2], the strongly-correlated regime can be reached, in principle, regardless of how small the on-site interactions are. Let us note, however, that there is a fundamental limit to this approach, which is imposed by external sources of decoherence and heating. Accordingly, the dynamics must always be faster than the time-scale imposed by these sources of noise.
For vibrational modes transverse to the microtrap array α = z, the conditionΩ Lηα ω α can be relaxed for ions in the node of the standing wave. It then suffices to setΩ Lη 2 α ω α , which allows us to reach interaction strengths in theŨ ∼1-10 kHz-regime, which are directly on the order of the bare tunnelings. Finally, including the gauge fields in this stronglycorrelated regime further enhances the versatility of our quantum simulator. Even if the particles are bosonic and the interactions are local, one can target fractional quantum Hall states, and composite-fermion fluids [43].

D. Non-Abelian synthetic gauge fields
The synthetic gauge fields discussed so far (39) correspond to standard electromagnetism. In this theory, they are formally introduced to restore the invariance with respect to local unitary transformations in the gauge group U(1). A natural question that arises is whether these synthetic fields can be generalized to different gauge groups, possibly non-Abelian ones [48]. In this case, the local unitary also acts on some additional degree of freedom, which we shall refer to as the flavor. Here, we show that it is possible to realize such scenarios by exploiting two or three orthogonal directions of vibration as the different flavors of the non-Abelian theory.
We first exploit two main axes of vibration α = x, y within the plane defined by the microtrap array, the corresponding trapping frequencies being different ω x = ω y . In the regime where the bare tunneling strength is smaller than the frequency difference, these two directions are uncoupled (see Eq. (30)). The main idea is to use two gradients of opposite sign for each degree of freedom, and two separate periodic drivings of the same frequency but of different amplitudes According to Eq. (35), this modulation of the trap frequencies can be obtained from a single Raman-beam scheme, such that the corresponding wavevector has a component along both axes η dx = −Ω L η 2 x , η dy = −Ω L η 2 y . By repeating the analysis of Sec. II, and considering the same type of conditions, we reveal that the dressed tunneling strengths (13) must be modified for each of the vibrational axes j 1 ), and one must introduce the following axis-dependent parameter r x = r, r y = −r. It is then straightforward to rewrite the effective Hamiltonian according to a generalized Peierls substitution where the synthetic gauge potential now depends on the related vibrational axis, and is equivalent to A x s (r) = −A y s (r) = −B 0 ye 1 , such that B 0 = rφ 2 /e * d 1 d 2 . Let us introduce now a bosonic spinor field operator Ψ i = (b x,i , b y,i ) t to describe the phonon fields corresponding to vibrations in each direction. The kinetic part can be written where K d;ij is a matrix that describes the vibrational couplings for each spinor component, and A na s = −B 0 yτ z e 1 is a SU(2) non-Abelian gauge field, that we write in terms of a Pauli matrix, τ z , acting on the vibrational index (i.e. flavor index). The associated magnetic field corresponds to B na s = B 0 τ z e z , so that each flavor is subjected to an opposite flux piercing the lattice. For fermions, these types of gauge fields give rise to the so-called quantum spin Hall effect [49], which is the prototype of time-reversal preserving topological insulators in two-dimensions [45]. Let us note that Eq. (43), together with the arguments presented in Eq. (19) imply that each flavor around a plaquette accumulates a non-Abelian Aharanov-Bohm phase that is governed solely by the SU(2) gauge field A na s . However, in addition to that, each flavor x, y, is subjected to different tunnelings,J x d;i,j , J y d;i,j , such that an spin-orbit coupling is superimposed to the non-Abelian gauge. The latter effect may enrich the dynamics with respect to the usual situation in SU(2) gauge theories. We remark that the above synthetic gauge field (43) is only a particular type of non-Abelian SU(2) gauge fields. In order to consider more general fields, it is possible to exploit the quadratic terms of the Coulomb interaction (28) that mix the vibrational modes along x, y. By setting the frequency of the periodic driving to account for both the frequency difference (ω x = ω y ) and the particular gradient, the dressed tunnelings would also involve a change of the flavor index, yileding thus more general gauge fields in the group SU(2).

E. Spin-mediated disordered Hamiltonians
The properties of solids usually differ from those of perfectly periodic crystals. In realistic samples, there is a certain amount of disorder in the form of impurities, dislocations, or vacancies, which may alter dramatically the properties of the solid. The study of disorder in solids is an active and mature field of condensed matter [50], where the system Hamiltonians are usually modeled as stochastic operators. Here, the randomness is due to a statistical description of the disordered degrees of freedom. To incorporate such a randomness in a quantum simulator, which by definition should be an extremely clean and controllable setup, one can exploit the quantum parallelism by an auxiliary degree of freedom [51].
In principle, some randomness could be introduced by randomly varying the lattice constant and the trapping frequencies within the microtrap array. In the limit of large arrays, these would lead to an off-diagonal bond disorder and diagonal site disorder, respectively. In this subsection, we elaborate on two alternative directions to widen the applicability of the phonon-based quantum simulator by introducing disorder regardless of the size of the microtrap array. In both cases, we shall make use of the electronic energy levels of the ion, {| ↑ i , | ↓ i }, to introduce randomness in the phonon Hamiltonian and mimic the effects of disorder. We emphasize that this setup shall allow us to control the two usual types of disorder independently, namely, diagonal and off-diagonal disorder.

Off-diagonal bond disorder
We discuss how to induce randomness on the phonon tunneling by a slight modification of Eq. (35), so that the periodic driving will be responsible for both the disorder and the synthetic gauge fields. In the derivation of the driving, we assumed that a proper choice of the laser intensities, detunings, and polarizations would lead us from Eq. (32) to the desired expression (33). This condition can be modified so that the periodic driving becomes spin dependent. In fact, by setting ae /2∆, one obtains where we have introduced This expression corresponds to a differential Stark shift between the two electronic levels. In the regime where ω L ω α , this term generalizes Eq. (35) to the following spin-dependent driving of the trapping frequencies where the same constraints (36) over the system parameters must be fulfilled. Since no other spin operators are involved in the Hamiltonian, one can treat σ z i as c-numbers σ i ∈ {−1, 1}, and carry out the same analysis made in Sec. II. In fact, one should simply modify Hamiltonian (39) (47) and the tunneling amplitudes account for the different spin configurationsJ α d;ij (σ i , σ j ) = J α t;ij F f (i,j) (η d σ i , η d σ j , ∆φ ij ). In Fig. 7, we represent the dressed tunneling amplitudes between nearest-neighbors along the direction of the gradient. It can be observed that depending on the spin state, one obtains different values. In particular, for η d ≈ 1 and ∆φ ≈ 3π/2, the tunneling amplitude can attain four possible values If we now consider the following initial state ρ 0 = |Ψ s Ψ s | ⊗ ρ 0 ph , where |Ψ s = ∑ {σ σ σ } c {σ σ σ } |{σ σ σ } and ρ 0 ph are arbitrary spin and phonon states, its time-evolution is (49) Note that due to the superposition principle of quantum mechanics, the phonon state explores simultaneously the different tunneling paths i → J(σ i σ j ) → j with a probability that depends on the initial spin configuration p {σ σ σ } = |c {σ σ σ } | 2 . In fact, the measurement of a phonon observable O ph corresponds directly to the statistical average over all such tunneling paths which can be understood as the average over all possible realizations of the bond disorder. Therefore, our phonon-based simulator can explore the physics of interacting disordered bosons in a lattice pierced by an external magnetic field.

Diagonal site disorder
We now address a scheme to introduce randomness in the on-site energies [22]. The origin of such terms is independent of the periodic driving, and requires a pair of additional laser beams in the Raman configuration (31). In order to distinguish them from the previous laser beams, we will be denote them with an overbar. The main difference with respect to the previous case is that the beatnote is tuned close to the transition between the internal statesω L =ω 1 −ω 2 ≈ ω 0 . In this case, the laser-ion Hamiltonian becomes where we have introduced σ + i = | ↑ i ↓ i |. In analogy to the derivation of the periodic driving, we express the position in terms of the local phonon operators, and expand the Hamiltonian for a small Lamb-Dicke parameter. If the laser beatnote  is tuned as followsω L ≈ ω 0 − ω α , one obtains the so-called red-sideband excitation where the bare detuning isδ L = ω 0 − ω α −ω L ω α , and we assume thatΩ L ω α in order to neglect the remaining terms of the Taylor expansion. In the regime where the laser beams are weak enoughΩ L η α δ L , it is possible to find the following laser-ion Hamiltonian in perturbation theorȳ This term arises due to a second-order process where a phonon is virtually excited and reabsorbed by a single ion, and leads to a differential Stark shift of the atomic levels that depends on the number of phonons. When incorporated to the effective description of the PAT Hamiltonian (39), it modifies the Kinetic energy term. If one is interested in the time span t ≈ 1/J c;ij , the dynamics of the spins can be safely ignored, and σ z i can be treated once more as c-numbers σ i ∈ {−1, 1}. We note that the typical time-scales for the spin flip-flop dynamics would be J s;ij ≈10-100 Hz, whereas thevibrational couplings lie in the J c;ij ≈1-10 kHz. Accordingly, the Kinetic energy term of Eq. (39) must be modified to In this situation, it is not the tunneling strengths which depend upon the spin state (47), but rather the on-site energies. Notice that the strength of these on-site energies must be much smaller that the trapping gradient ε α ∆ω α in order not to affect the PAT scheme. According to Eq. (50), the system explores simultaneously all possible values of the on-site energies with probabilities that depend on the spin configurations of the initial state. Therefore, the measurement of phonon observables yields directly the statistical average over all possible realization of ε j ({σ σ σ }) ∈ {−ε α , ε α } with a probability distribution p {σ σ σ } = |c {σ σ σ } | 2 .
Let us now comment on the extended possibilities of our quantum simulator due to the engineered disorder in Eqs. (47) and (54). The diagonal site disorder leads to the well-known Anderson localization in the non-interacting limit [53]. This phenomenon is due to the interference of the different paths associated to the scattering of the particles from the random on-site fluctuations, and gives rise to exponentially-localized wavefunctions and absence of diffusion. The combination of Anderson Localization with strong interactions, which is also a well-studied problem [42], leads to interesting insulating, yet gapless, phases such as the Bose glass. In the case of strong bond disorder, a different gapless insulator known as the Mott glass arises, which consists of disconnected superfluid regions of random size [52]. Besides, our quantum simulator has the potential of combining both bond and site disorder, and tuning them independently, which may pave the way towards other exotic insulating phases. In addition to the aforementioned Bose and Mott glasses, the simulator can explore the random-singlet glass where the bosons form delocalized random pairs [54]. The possibility to explore higherdimensions, longer-range tunnelings, and the effect of synthetic gauge fields makes our scheme a very versatile tool.

F. Decorated synthetic gauge flux lattices
In this section, we describe an additional feature of our quantum simulator: the possibility to decorate the lattice with any desired pattern of synthetic fluxes. It is thus possible to engineer highly inhomogeneous synthetic gauge fields, even reaching inhomogeneities at the unit-cell limit. The idea is to use the spins to decorate the array with different fluxes by exploiting the differential phonon-dependent Stark shift (53).
We consider a situation where the gradient of the trapping frequencies vanishes ∆ω α = 0, so that the regime is different from that of the site-disorder case considered above ε α ∆ω α . When the periodic driving frequency (35) and the Stark shift (53) fulfill the resonance condition ω L = 2ε α /r, where r is the integer representing the umber of photons involved in the PAT, the assisted tunneling will give rise to a different phase depending on the spin states of the two neighboring ions. In order to find the correct expression for the tunneling, we readdress the derivation of Sec. II for this particular spin-dependent situation. Since we are interested in the phonon dynamics, the spins are effectively frozen, and we can consider σ z i as c-numbers σ i ∈ {−1, 1}. Hence, the dressed tunneling (13) must be modified as follows where we have introduced f (σ i , σ j ) = r(σ i − σ j )/2. By considering the possible spin configurations, we find (56) Interestingly, the phase of the tunneling between sites j → i depends on their internal spin state σ j → σ i . If their spins are parallel σ i = σ j = ±1, the phase vanishes and thus the tunneling does not contribute to the synthetic gauge field. Conversely, when the spins are anti-parallel, the phases contribute to the gauge fluxes with a sign that depends on the particular spin ordering. We have calculated the consecutive phonon tunneling around a square plaquette r 0 for all the possible spin configurations All the possible encircled fluxes W (1) = |W (1) |e iφ ({σ σ σ }) for r = 1 have been represented in Fig. 8(a), where we observe that there are 9 possible fluxes out of the 2 4 = 16 possible spin configurations. These fluxes correspond to φ ∈ {0, ±φ 1 , ±φ 2 , ±φ + , ±φ − }, where φ i = ∆k i d i , and φ ± = (φ 1 ± φ 2 )/2. Accordingly, we have 9 different tiles that can be used to decorate the underlying microtrap array with a particular distribution of fluxes. Let us emphasize that the particular distribution of tiles is completely determined by the spin state |Ψ s = |σ 1 , · · · , σ N , which can be initialized at will in trapped-ion experiments. Let us also note that we could also exploit this result to introduce randomness in the gauge fields by considering a linear superposition of the spin states. In Fig. 8(b-c), we represent two-color flux lattices that correspond to a staggered magnetic field along both principal axes. In Fig. 8(d-e), we represent three-color flux lattices where a staggered flux alternates with a vanishing flux. In Fig. 8(f), we represent a checkerboard flux lattice, and finally in Fig. 8(g), we represent a limiting case of a six-color flux lattice, which is very interesting from a physical point of view. By setting φ 1 = π, and φ 2 = π, the two fluxes ±φ 1 = ±π = πmod2π are equivalent, and lead to a homogeneous π-flux model in a square lattice. Additionally, the four remaining fluxes vanish for this choice ±φ + = ±φ − = 0mod2π. Therefore they contribute with a local defect over the π-flux lattice, which could bind excitations with anyonic statistics when the longer range tunnelings are taken into account [55].
For each of the decorated flux lattices, it is possible to find a particular inhomogeneous synthetic gauge field A s (σ i , σ j ) so that the effective phonon Hamiltonian is rewritten in a standard Peierls form The particular pattern of the spins will determine the inhomogeneous gauge field, and the way the lattice is decorated with fluxes. Let us also note that the complexity of Fig. (8) will increase when the larger plaquettes due to long-range tunnelings are also taken into account.

IV. SUMMARY AND OUTLOOK
We have introduced the two key ingredients to realize PAT experiments in micro-fabricated ion traps. The first ingredient is a gradient of the trapping frequencies achieved by the local control of the trap electrodes. The second corresponds to a Raman-beam configuration, which presents different regimes that provide (i) the periodic driving of the trapping frequencies, (ii) the on-site phonon-phonon interactions, (iii) the bond/site disorder, and also (iv) an exotic flux decoration of the microtrap array. We believe that such ingredients are at reach of current microtrap technology, and their correct combination will give raise to a very versatile quantum simulator for many-body bosonic models. Such photon-assisted-tunneling toolbox for quantum simulations can be summarized in the following general Hamiltonian H eff = ∑ {σ σ σ } (K eff ({σ σ σ }) +V eff )|{σ σ σ } {σ σ σ }|, where where the particular expression for the dressed tunneling and synthetic gauge field will depend on the configuration of the frequency gradient and the periodic driving. Let us list the possibilities that have been explored in this work: i) Dynamic localization: The tunneling amplitudeJ α d;ij = J α t;ij F 0 (η d , η d , ∆φ ij ) does not depend on the spin state, and the synthetic gauge field vanishes A α s = 0. This is achieved in the regime of vanishing gradient ∆ω α = 0, and setting the beatnote of the Raman beams ω L ω α . By tuning η d , one can find a value where the tunneling strength vanishes and thus the particles are dynamically localized, a phenomenon also known as coherent destruction of tunneling [25]. ii) Synthetic Abelian gauge fields: The spin-independent tunneling amplitude isJ α d;ij = J α t;ij F f (i,j) (η d , η d , ∆φ ij ), where f (i, j) = r(i 1 − j 1 ) depends on an integer r. The synthetic gauge potential is A s (r) = −B 0 ye 1 , and follows from the regime with a finite gradient ∆ω α J α t;ij , such that the beatnote of the Raman laser beams fulfills the resonance condition ω L = ∆ω α /r. In this case, phonons behave as charged particles that move in a two-dimensional plane pierced by an orthogonal magnetic field whose flux φ 2 = ∆k 2 d 2 can be modified by varying the Raman wavevector. Phenomena typical of integer quantum Hall samples [39,40], or bosonic flux ladders [41], can also be observed in this platform. Besides, in combination with strong phonon-phonon interactions, one can find bosonic versions of the fractional quantum Hall states [43].
iii) Synthetic non-Abelian gauge fields: The above scheme can be generalized to the non-Abelian gauge group SU(2), such that both in-plane vibrational modes play the role of a flavor component. We have described in detail a particular SU(2) gauge field, which requires the vibrations along each direction to be subjected to an opposite frequency gradient. Using the same assumptions as in the Abelian case, the tunneling amplitude becomesJ α d;ij = J α t;ij F f α (i,j) (η dα , η dα , ∆φ ij ), with f α (i, j) = r α (i 1 − j 1 ) such that r x = −r y = r. Hence, the synthetic gauge field A na s = −B 0 yτ z e 1 becomes a SU(2) operator acting in the flavor space. This scheme opens a route towards a bosonic counterpart of the quantum spin Hall effect [49].
iv) Bond and site disorder: The dressed-tunneling ampli-tudeJ α d;ij (σ i , σ j ) = J α t;ij F f (i,j) (η d σ i , η d σ j , ∆φ ij ), and the on-site energies ε α σ i , take on different values depending on the spin configuration. If the initial state is a linear superposition of different spin configurations, the phonon dynamics is determined by a random Hamiltonian with bond and site disorder. This regime is achieved for a gradient ∆ω α J α t;ij , and laser beatnote ω L = ∆ω α /r giving rise to a state-dependent periodic driving. Besides, an additional Raman beatnote tuned close to the atomic transitionω L ≈ ω 0 − ω α gives raise to a phonondependent Stark shift in the limit of large detuning. In the noninteracting regime, this tool allows us to explore the physics of Anderson localization [53]. By adding strong interactions, it yields gapless insulating phases such as the Bose glass [42], the Mott glass [52], and the random-singlet glass [54].
vi) Decorated flux lattices: In the absence of the frequency gradient, one can tune the Raman lasers beatnote in resonance to the above phonon-dependent Stark shift ω L = 2ε α /r. Once again, the dressed-tunneling amplitude becomes spin-dependentJ α d;ij (σ i , σ j ) = J α t;ij F f (σ i ,σ j ) (η d , η d , ∆φ ij ) where f (σ i , σ j ) = r(σ i − σ j )/2. Moreover, the synthetic gauge field also depends on the spin configuration and we can make decorated flux lattices as those shown in Fig. 8 by selecting a particular spin state. Phenomena related to charged particles under inhomogeneous magnetic fields can be explored, such as staggered fields or π-flux lattices with defects.