String-breaking dynamics in a quantum simulator

MainAtomic nuclei are described at a fundamental level by quantum chromodynamics of quarks and gluons. Quarks and gluons, which carry a net colour charge, are not found in isolation, but only in colourless agglomerates known as hadrons (for example, protons and neutrons), tied together in an effect referred to as colour-charge confinement1,2. This phenomenon is commonly conceptualized via a paradigmatic thought experiment, whereby a pair of probe charges are pulled apart. In nature, pulling apart an electron and a positron to an infinite distance requires a finite electrostatic energy. By contrast, in quantum chromodynamics, the potential energy between a quark and an antiquark increases indefinitely with distance, stretching a gluon string that binds them together. As the separation increases, the energy eventually becomes larger than twice the quark mass, driving the creation of new quark–antiquark pairs, thus breaking the string. String-fragmentation dynamics are conjectured to govern the hadronization processes in relativistic heavy-ion collisions3,4 and in the cooling of the expanding universe after the Big Bang5,6. However, such phenomena remain poorly understood at a fundamental level. While classical simulations based on sampling methods have successfully computed the static string breaking in quantum chromodynamics7, they are severely inefficient in simulating string dynamics8. Other classical methods, such as tensor networks, may be more suited to study this problem9,10,11,12,13,14,15,16. However, such methods have yet to be developed for, and applied to, the theory of the strong force17,18,19, and they will probably be inefficient for high-energy processes with abundant entanglement generation.By using quantum simulators—in analogue, digital or hybrid modes—the same dynamical processes can be simulated efficiently, that is, their simulation requires computational resources that scale only polynomially with system size and evolution time. As a result, quantum technologies may present a promising path towards simulating quantum chromodynamics17,20,21,22,23,24,25,26,27,28. While a quantum-chromodynamics simulator remains beyond current capabilities, existing platforms have started to investigate charge confinement in simplified lower-dimensional models. These experiments have been realized using trapped ions29,30,31,32,33,34, neutral-atom tweezer arrays35,36, optical lattices37,38,39 and superconducting-qubit processors40,41,42,43. However, experimental quantum simulations of lattice gauge theories beyond proof-of-principle demonstrations have so far been limited to probing the real-time response to homogeneous quantum quenches from a vacuum state. Investigations of more directly relevant phenomena such as hadron collisions44,45,46,47 or string-fragmentation dynamics demand a higher level of control precision and programmability than that achieved so far.In this work, we develop, and experimentally implement, an efficient strategy to study the dynamics of confined charges and strings using an analogue trapped-ion quantum simulator. Precise control over individual ions allows us to emulate the effects of virtual static environments, which we use to induce a net physical charge, or a physical string stretched between heavy probe charges. A rapid change of Hamiltonian parameters then drives non-equilibrium charge or string dynamics, which we observe with full spatiotemporal resolution. Our results push the boundaries of quantum simulation, providing new insights into quantum many-body dynamics, with potential relevance to nuclear and high-energy physics.Confinement, a hallmark feature of quantum chromodynamics in (3 + 1) dimensions, can also be studied using (1 + 1)-dimensional gauge theories48,49,50. In this work, we study dynamics of a \({{\mathbb{Z}}}_{2}\) lattice gauge theory consisting of fermionic matter residing on the sites of a one-dimensional lattice, coupled to \({{\mathbb{Z}}}_{2}\) gauge fields residing on the bonds (links) between lattice sites. The fermionic charge at each site is constrained by the electric fields at adjacent lattice bonds via Gauss’s law. Alternatively, Gauss’s law can be used to eliminate the fermions from the description of dynamics51. The resulting fully bosonic dual theory, as shown in the Methods and illustrated in Extended Data Fig. 1, can be mapped to a one-dimensional quantum Ising chain with the Hamiltonian52,53,54,55$$H=-\sum _{i 0 is the Ising-interaction coupling between spins i and j.A charge in the original formulation is encoded as a kink in the dual formulation, where a kink is a pair of anti-aligned spins ↑↓ or ↓↑. The charge density is quantified by \({q}_{i}=\langle 1-{\sigma }_{i-1}^{z}{\sigma }_{i}^{z}\rangle /2\). Meanwhile, the local longitudinal spin polarization encodes the electric field, denoted \({\epsilon }_{i}=\langle {\sigma }_{i}^{z}\rangle\). An electric-field string is associated with a down-spin domain ↓⋯↓. Figure 1a depicts degrees of freedom in the original lattice gauge theory and the dual spin theory. The couplings Ji,j control the particle mass and short-range interactions. The longitudinal field h generates a string tension, that is, the energy stored in the string per unit length. The transverse field introduces quantum fluctuations that couple charge and string dynamics, with coupling strength g. The tendency of strings to break dynamically and form new charge pairs is enhanced through both an increased coupling strength and an increased string tension.Fig. 1: Charge and string dynamics with a trapped-ion quantum simulator.Full size imagea, Mapping between the degrees of freedom of a (1 + 1)-dimensional \({{\mathbb{Z}}}_{2}\) lattice gauge theory and an Ising spin model implemented by our simulator. Charges in the gauge theory are represented by kinks in the spin model, and strings are represented by down-spin domains. b, The initial state prepared for the simulation protocol of Fig. 2, implementing the evolution of an isolated charge. The virtual chains (grey-shaded regions) realize a semi-infinite static string stretching to the right and a semi-infinite static vacuum on the left, enforcing a net dynamical charge in the physical region. c, The initial state prepared for the simulation protocol of Fig. 3. The static-spin configurations at the edges of the two semi-infinite virtual chains realize a pair of static charges (physically corresponding to heavy charges). d, Sketch of a spin configuration representing a state with dynamically generated charge pairs. Such states can arise during the simulation protocol of Fig. 3, which implements far-from-equilibrium string-breaking dynamics. e, Schematic of the trapped-ion apparatus, showcasing two arrays of tightly focused laser beams that can be independently controlled. One of the arrays (purple) generates long-range spin–spin Ising interactions. The second array of beams (red) implements a site-dependent magnetic field. Each of the beam amplitudes (different shades) is controlled to generate a uniform interaction profile across the ion chain as well as a site-dependent local magnetic field (Extended Data Fig. 2). Frequency and phase control of each beam (wiggly lines inside the beam) enables full tunability of the magnetic-field direction along any axis of the Bloch sphere x, y, z (dashed arrows on top).Using a quantum simulator to study the charge and string dynamics induced by the Hamiltonian in equation (1) poses two main challenges. First, experiments can only simulate finite systems. To make the best use of quantum resources, we analyse a simpler system composed of infinitely long arrays of static spins to the left and right of a finite domain of dynamical spins14. While only the dynamical spins are encoded in the quantum simulator, an additional site-dependent longitudinal field, Δhi, can be applied to the dynamical spins to emulate interactions with the static spins (Methods). This scheme enables the realization of external static strings or charges, as depicted in Fig. 1b,c, and alleviates the undesired effects of hard boundaries. It, nonetheless, requires site-dependent control of h in the experiment. This requirement is related to the second experimental challenge, which is the need for simultaneous implementation of all the Hamiltonian terms with programmable strengths.We address these challenges by using a trapped-ion quantum simulator with enhanced control. We encode the spins into two internal levels of 171Yb+ ions and implement the transverse-field Ising model (the first and third terms in equation (1)) with a linear array of L = 13 spins addressed via bichromatic laser beams. These beams off-resonantly couple to the collective normal modes of ion motion56,57,58. Simultaneous control of the amplitude, phase, and frequency of each beam in the array allows us to realize the uniform field g, as well as to program the interaction profile Ji,j to be translationally invariant and exponentially decaying with distance, Ji,j = Je−β(∣i−j∣−1) with β ≈ 0.78 (Methods). To realize the site-dependent longitudinal field (the second term in equation (1) with h → hi = h + Δhi), we apply a second array of tightly focused beams that drive carrier Raman transitions, with independent control over the amplitude, phase and frequency of each beam, as shown in Fig. 1e. This set-up enables simultaneous and independent control of the local magnetic-field vector at each spin site, in addition to programmable long-range spin–spin interactions. This new Hamiltonian-engineering capability enables the emulation of virtual degrees of freedom, such as those in Fig. 1b,c, opening the door to controlled investigations of string-endpoint and bulk-fragmentation dynamics.In preparation for studying string dynamics, we first focus on the evolution of an isolated charge. While an isolated quark cannot be observed in particle colliders, in our quantum simulator an isolated charge is realized by emulating the effect of a semi-infinite static string extending to the right and a semi-infinite vacuum extending to the left of the simulated physical region, as illustrated in Fig. 1b. The expression of the corresponding site-dependent field, \(\Delta {h}_{i}^{\,\text{(charge)}\,}\), is provided in the Methods. We study the spreading dynamics of the string endpoint, which may precede and affect bulk string breaking. The interplay between these two dynamical processes was first theoretically addressed in refs. 14,16 for the specific case of resonant string breaking.In Fig. 2a–d, we present the evolution of the charge distribution governed by several Hamiltonians with non-vanishing interactions, starting from the classical (that is, fluctuation-free) initial state shown in Fig. 1b, with charge localized at the centre of the simulated region. The corresponding electric-field spatiotemporal dynamics are shown in Fig. 2e–h. In the absence of string tension (h = 0), charge spreads ballistically across the physical lattice until it hits the static boundaries, as shown in Fig. 2a,c,e,g. The measured velocity matches well with the value \({v}_{\max }=2g\), calculated by retaining charge-hopping processes and neglecting charge-pair creation processes, which is justified for small g/J. As the string tension is increased, we observe a Wannier-Stark localization phenomenon59, shown in Fig. 2b,d,f,h, which can be understood within the same single-charge approximation60. The string tension imparts a constant acceleration to the charge. In continuum space, a string can pull a charge indefinitely far, as its momentum can grow unbounded. However, on a lattice, the momentum is limited to a finite range. A constant acceleration, therefore, generates a periodic variation of momentum in time. This gives rise to coherent oscillations of the charge around its initial position, akin to Bloch oscillations, originally predicted for an electron moving in a defect-free crystal lattice under a constant electric field61. This single-charge description of dynamics becomes quantitatively accurate for small g/J and h/J, where the spatial amplitude and temporal period of the oscillations (in units of inverse J) are predicted to be 2g/h and πJ/h, respectively60,62. We find good agreement between these theoretical predictions and the data, even for non-perturbatively large g/J and h/J, and between the experiment and numerical simulations, as shown in Extended Data Figs. 3 and 4. Overall, Fig. 2 demonstrates that the main qualitative effect of charge confinement is to halt charge spreading dynamics, thereby localizing string endpoints14,54,62.Fig. 2: Non-equilibrium charge dynamics.Full size imageTime evolution of a perfectly localized charge at the centre of the physical-spin chain (Fig. 1b) in the presence of a range of coupling-strength and string-tension values. a–d, The evolution of the charge density \({q}_{i}=\langle 1-{\sigma }_{i-1}^{z}{\sigma }_{i}^{z}\rangle /2\) for different values of string tension (h/J) and coupling strength (g/J). e–h, The evolution of the electric field \({\epsilon }_{i}=\langle {\sigma }_{i}^{z}\rangle\) for the same values of string tension and coupling strength. Note that, initially, the charge distribution is qi = δi,0, the region to the left of the charge is the classical vacuum state with the electric field \({\epsilon }_{i}=\langle {\sigma }_{i}^{z}\rangle =1\), and the region to the right of the charge is a (semi-infinite) classical string state with \(\langle {\sigma }_{i}^{z}\rangle =-1\). The superimposed lines and arrows correspond to a single-charge approximation of dynamics: lines \(i=\pm {v}_{\max }t\) in a, c, e and g, with \({v}_{\max }=2g\), highlight the maximum speed of propagation of the charge, and arrows in b, d, f and h highlight amplitude (2g/h) and period (\(\pi{J/{\bar{h}}}\)) of the coherent oscillations. Dashed vertical lines mark the boundaries between the simulated physical region and the virtual static environments. Further exploration of the parameter space for both charge and electric-field dynamics, along with their comparisons to numerical simulations, is presented in Extended Data Figs. 3 and 4, respectively.To study the non-equilibrium evolution of strings, it is convenient to disentangle the spatial dynamics of string endpoints, observed in Fig. 2, from the genuine charge-pair-creation processes arising from string breaking. To do so, we will consider a string delimited by static charges and surrounded by a static vacuum, as illustrated in Fig. 1c. To accurately realize this configuration experimentally, we include a set of site-dependent longitudinal-field values \(\Delta {h}_{i}^{\,\text{(string)}\,}\) in the simulated physical region. This field emulates the coupling with the two semi-infinite virtual chains shown in Fig. 1c. Similar to our simulation protocol in Fig. 2, we prepare the physical region at time t = 0 in the classical (that is, fluctuation-free) string state—which is the ground state with vanishing string tension and coupling strength (g = h = 0)—and abruptly increase both g and h. The observed out-of-equilibrium string dynamics are shown in Fig. 3. We have explored a wide range of parameter values encompassing stronger abrupt changes in the coupling strength and the string tension compared to our isolated-charge experiment in Fig. 2. These stronger parameter values correspond to an enhanced probability of charge-pair creation in the ensuing time evolution.Fig. 3: Non-equilibrium string dynamics.Full size imageTime evolution of a classical string state (Fig. 1c) with a pair of external static charges. a–i, The evolution of the charge density \({q}_{i}=\langle (1-{\sigma }_{i-1}^{z}{\sigma }_{i}^{z})/2\rangle\), with coupling strength (g/J) and string tension (h/J) varying across the panels. The initial charge distribution is qi = δi,−7 + δi,8, with the yellow bars denoting the position of static charges at the boundary. The corresponding numerical results are shown in Extended Data Fig. 5. j–l, The electric field \({\epsilon }_{i}=\langle {\sigma }_{i}^{z}\rangle\), corresponding to the panels c, f and i, measured at different times (colours denote different values of Jt) for h = 0.6J (Exp.), together with the corresponding numerically simulated profiles (Sim.), and contrasted to the expected profiles in thermal equilibrium (dashed red lines). The electric-field dynamics for all the parameters and the corresponding numerical simulations are reported in Extended Data Figs. 6 and 7, respectively. All experimental measurements were averaged over 300 repetitions, and error bars in j–l represent standard deviations around the mean value.The abrupt changes in the Hamiltonian parameters inject a finite energy density into the system. The strong many-body interactions are expected to dynamically lead the system from its out-of-equilibrium classical-string initial state to a thermal-equilibrium state at long times. The relaxation dynamics toward this equilibrium state is conjectured to be initiated by a uniform, spontaneous creation of charge pairs in the bulk of the system, via a mechanism attributed to Schwinger63. A hallmark of this mechanism is the extreme sensitivity of the characteristic time scale of bulk charge-pair creation to variations in the system parameters—a manifestation of the non-perturbative origin of this effect. The observed string evolution in our experiment, nonetheless, reveals a mechanism distinct from this conventional expectation.Explicitly, we observe charge-pair formation systematically occurring at the string edges, that is, near the static charges, as schematically portrayed in Fig. 1d. This is demonstrated in the spatiotemporally resolved evolution of the charge distribution plotted in Fig. 3a–i, as well as in the corresponding electric-field evolution shown in Extended Data Fig. 6. For vanishing or weak string tension, we find that such dynamical charge pairs perform coherent oscillations confined to the edges, as shown in Fig. 3a,d,g. For larger string tension, however, these charge pairs propagate and spread from the edges towards the bulk, as shown in Fig. 3c,f,i. This phenomenon—confirmed by numerical simulations shown in Extended Data Figs. 5 and 7—occurs over a time s cale that does not sensitively depend on Hamiltonian parameters, all the way from the region of strong to weak quantum fluctuations. For example, the onset and the oscillation period of pair creation at the edges remain consistent over the string-tension and coupling-strength values studied. This is in stark contrast to the exponential enhancement, as a function of string tension and coupling strength, of the rate of spatially uniform charge-pair formation in the bulk associated with the Schwinger mechanism63,64. We conclude that what we observe is a new, non-conventional string-breaking mechanism.To elucidate the nature of the observed string-breaking dynamics, we develop a perturbative approach, thoroughly discussed in the Methods. This perturbative analysis demonstrates that, for small g/J and h/J, the main contribution to non-equilibrium string dynamics arises from the quantum spreading dynamics of a single pair of charges generated by our dynamical protocol, as sketched in Fig. 4. The associated two-body wave function initially peaks at the two configurations with the charge pair localized at the left or right edge of the string, respectively, where the influence of the external vacuum makes the energy cost of charge-pair creation smallest. The two-body wave function subsequently evolves in an effective potential-energy landscape shaped by the coupling to the external regions and by the value of the string tension, shown in Fig. 4a–c. The plots show that, as the string tension is increased, two equipotential trajectories originating from the two string edges open up in configuration space, demonstrating the possibility for the charge pair to quantum-mechanically propagate and spread from the edges to the bulk, as sketched in Fig. 4d and calculated numerically in Extended Data Fig. 8a–f.Fig. 4: Edge-facilitated string-breaking mechanism.Full size imagea–c, Colour maps show the effective two-body potential landscape experienced by a pair of dynamical charges, as a function of their positions −6 ≤ i1 0). This refined observable thus provides a non-perturbative metric to confirm the validity of our string-breaking picture throughout the range of dynamical control parameters. Altogether, this analysis provides strong theoretical support to our previous conclusion that the observed string-breaking process differs from the Schwinger mechanism. Such an edge-facilitated string-breaking mechanism is distinct from previously studied processes. In the Methods and in Extended Data Fig. 8g–i, we further analyse this novel mechanism in a continuum limit, highlighting its generality. Furthermore, in Extended Data Figs. 9 and 10, we report additional numerical computations of string dynamics for the same parameters as in Fig. 3 but in the presence of a dynamical fluctuating exterior vacuum. These computations confirm that the edge-facilitated string-breaking mechanism is not an artefact of the simplified set-up with fully static environments as simulated in our experiment.Fig. 5: Time-dependent domain-size probability.Full size imagea–i, The time evolution of the probability distribution Pd(s) of the largest ‘vacuum bubble’ (that is, spin-up domain) size s = 0, 1, …, L within the same instances of far-from-equilibrium string dynamics for the values of coupling strength and string tension as in the corresponding panels of Fig. 3.Lastly, to confirm that the observed transient string dynamics precede the system’s eventual thermalization, we compare the time evolution of a local observable with its thermal-equilibrium counterpart (Methods) in Fig. 3j–l. We focus on the electric field, whose sign indicates the local predominance of either a string state (as in the initial state) or a vacuum state as the system evolves in time. For weak quantum fluctuations (Fig. 3j), the string character of the state persists throughout the observed evolution, with slight variations from the motion of charges at the edges. For stronger quantum fluctuations (Fig. 3k,l), the string character is lost quicker, indicating a faster thermalization process. However, even in this regime, the system remains out of equilibrium during the entire observation time.The experiments reported in this work demonstrate that trapped-ion architectures with fully programmable dual beam arrays have the capability to perform analogue simulations of far-from-equilibrium real-time dynamics of interacting charges and strings in a confining gauge theory, with full spatiotemporal control of Hamiltonian interactions. Our strategy of emulating probe charges and virtual environments using spatiotemporal control can be adapted to other quantum-simulation platforms that feature individual qubit control, such as two-dimensional ion crystals65,66 or neutral-atom arrays67,68,69, and combined with existing proposals for simulating non-Abelian and higher-dimensional lattice gauge theories17,21,22,23,24,25,26,27,28, to enable string-breaking studies in theories of further relevance to the Standard Model.Our results reveal a new mechanism for transient charge-formation and string-breaking dynamics, distinct from the conventional Schwinger mechanism. Future experiments with a receding pair of probe charges, or fully dynamical strings and surrounding environment, can elucidate the potential relevance of the string-breaking mechanism observed in this work to high-energy colliders and cosmology. Further studies of string breaking in controlled dynamical protocols, and in particle collisions, can probe these and more features of string-breaking dynamics9,10,11,12,13,14,15,16,44,46,47,70. Such simulations are also within reach with programmable trapped-ion simulators46,71,72.In summary, our experimental realization of string breaking in a gauge theory marks substantial progress towards the ultimate goal of first-principles quantum simulations of string fragmentation and hadronization in high-energy collisions and cosmological settings.Some of the results, including clear experimental signatures of string breaking, were reported by Or Katz at a KITP conference in November 2023 (https://doi.org/10.26081/K66Q5V). All results were presented by Arinjoy De in ‘Observation of string-breaking and charge confinement in a trapped-ion quantum simulator‘ at DAMOP 2024, Abstract C06.00001.After completion of our experiments, ref. 73 reported observation of dynamics of charges and strings in a two-dimensional array of superconducting qubits. During the preparation of this Article, we were also made aware of a manuscript on the observation of string breaking in a (2 + 1)-dimensional lattice gauge theory implemented in a Rydberg-atom array by QuEra Computing Inc. and collaborators74.MethodsTheoretical model and analysisIn this first part of the Methods, we detail the (1 + 1)-dimensional \({{\mathbb{Z}}}_{2}\) lattice gauge theory that we consider and its connection with the experimentally simulated quantum spin-chain dynamics according to equation (1). We then discuss non-equilibrium string dynamics at the lowest order in the inverse particle mass and provide a theoretical understanding of the experimentally observed charge-pair formation and spatiotemporal evolution. Lastly, we describe our calculation of thermal expectation values of the electric field, which are compared with the experimental and numerical time-evolving electric field profiles during string-breaking dynamics.Mapping a \({{\mathbb{Z}}}_{2}\) lattice gauge theory to the Ising spin modelHere, we introduce a \({{\mathbb{Z}}}_{2}\) lattice gauge theory and illustrate how its dynamics in the gauge-invariant sector are equivalent to those of the quantum Ising chain in equation (1). Our description adapts the derivation of refs. 44,54 to general variable-range Ising interactions. The Hamiltonian of the (1 + 1)-dimensional \({{\mathbb{Z}}}_{2}\) lattice gauge theory we consider reads$$\begin{array}{ll}{H}_{\mathrm{LGT}}= & -g\mathop{\sum}\limits _{l}\left({c}_{l}^{\dagger }({b}_{l}+{b}_{l}^{\dagger }){c}_{l+1}+{c}_{l}^{\dagger }({b}_{l}+{b}_{l}^{\dagger }){c}_{l+1}^{\dagger }+\mathrm{h.c.}\right)\\ & +m\mathop{\sum}\limits _{l}{c}_{l}^{\dagger }{c}_{l}+\kappa \mathop{\sum}\limits _{l}{n}_{l}-\mathop{\sum}\limits _{l}\mathop{\sum}\limits _{r > 1}{v}_{r}{n}_{l}{n}_{l+r},\end{array}$$ (2) where \(l\in {\mathbb{Z}}\) are sites of a one-dimensional lattice. In this equation, \({c}_{l}^{\dagger }\) and cl are creation and annihilation operators of fermionic particles on site l. We call an occupied fermionic site, that is, \({c}_{l}^{\dagger }{c}_{l}=1\), a charge at site l. Similarly, \({b}_{l}^{\dagger }\) and bl are hard-core-boson creation and annihilation operators residing on the lattice bond connecting sites l and l + 1, with \({n}_{l}={b}_{l}^{\dagger }{b}_{l}\). These fields correspond to gauge-field degrees of freedom—a \({{\mathbb{Z}}}_{2}\) electric field. The first line of equation (2) represents a minimal coupling between the matter and gauge fields, of strength g; the second line is associated with the rest mass m of the fermionic particles; the third line represents an electrostatic energy, including a uniform energy cost κ for the excited electric-field configuration, as well as a variable-range self-interaction of this field, described by the couplings vr.The Hamiltonian in equation (2) is invariant under \({{\mathbb{Z}}}_{2}\) gauge transformations generated by the local Gauss-law operators$${G}_{l}={(-1)}^{{n}_{l-1}+{n}_{l}+{c}_{l}^{\dagger }{c}_{l}},$$ (3) that is, [Gl, HLGT] = 0 for all l. We will consider the gauge-invariant sector, spanned by eigenstates of the Gauss’s-law operators with eigenvalue one, that is, the sector where Gl = 1 for all l. Due to these local constraints, a configuration of the bosonic gauge fields fixes a unique configuration of the fermionic matter (see, for example, Extended Data Fig. 1). One can thus eliminate the fermionic degrees of freedom and obtain an exact description of the model in terms of the gauge fields only.To formally obtain such a representation, first one applies the Jordan–Wigner transformation to turn fermions into spin-\(\frac{1}{2}\) operators \({\tau }_{l}^{\alpha }\), for α = +, −, z, with$${\tau }_{l}^{-}=\prod _{m 1}{v}_{r}\right)\mathop{\sum}\limits _{l}{\sigma }_{l}^{z}-\frac{1}{4}\mathop{\sum}\limits _{l}\mathop{\sum}\limits _{r > 1}{v}_{r}{\sigma }_{l}^{z}{\sigma }_{l+r}^{z},\end{array}$$ (6) with the constraint \({G}_{l}=-{\sigma }_{l-1}^{z}{\tau }_{l}^{z}{\sigma }_{l}^{z}=1\). We now introduce a unitary transformation U such that the transformed Gauss’s law \(G{\prime} =U{G}_{l}{U}^{\dagger }=1\) depends only on the matter degrees of freedom, whereas the transformed Hamiltonian \(H{\prime} =U{H}_{{\rm{LGT}}}{U}^{\dagger }\) depends only on the gauge-field degrees of freedom. This is accomplished by choosing$$U=\prod _{l}\exp \left[\frac{i\pi }{2}({\tau }_{l}^{x}-1)\frac{1-{\sigma }_{l-1}^{z}{\sigma }_{l}^{z}}{2}\right]\,.$$ (7) The transformed constraint \({G}_{l}^{{\prime} }=-{\tau }_{l}^{z}=1\) decouples the τ spins. Subsequently, the transformed Hamiltonian depends only on the σ spins:$$\begin{array}{ll}{H}_{{\rm{LGT}}}^{{\prime} }=\,-g\mathop{\sum}\limits _{l}{\sigma }_{l}^{x} \,+m\mathop{\sum}\limits _{l}\frac{1-{\sigma }_{l-1}^{z}{\sigma }_{l}^{z}}{2}\,-\frac{1}{2}\left(\kappa -\mathop{\sum}\limits _{r > 1}{v}_{r}\right)\mathop{\sum}\limits _{l}{\sigma }_{l}^{z}\\\qquad\qquad-\frac{1}{4}\mathop{\sum}\limits _{l}\mathop{\sum}\limits _{r > 1}{v}_{r}{\sigma }_{l}^{z}{\sigma }_{l+r}^{z}.\end{array}$$ (8) We now identify the bonds of the fermionic chain connecting sites l and l + 1 with sites i of the dual spin chain, and the parameters as$$m=2{J}_{1},\,\,\,\,\kappa =2h+\mathop{\sum }\limits_{r=2}^{\infty }{J}_{r},\quad {v}_{r}=4{J}_{r}.$$ (9) Here, Jr ≡ Ji,i+r is the spatial profile of translationally invariant Ising interactions. One can then readily observe that equation (8) coincides (up to an irrelevant additive constant) with the Hamiltonian in equation (1) implemented by our trapped-ion simulator.Perturbative analysis of string dynamicsHere, we provide theoretical understanding of non-equilibrium string-breaking dynamics in the limit where the particle mass is the dominant energy scale, g, h, J2, J3, J4, ⋯ ≪ J1 = J (see equation (9)). In this regime, the evolution of the classical-string initial state \(\left\vert {\varPsi }_{0}\right\rangle\) (depicted in Fig. 1c) is accurately described for times 0 < t ≪ J/g2 by a perturbative expression of the form$$\left\vert \varPsi (t)\right\rangle =\mathrm{e}^{-i{E}_{0}t}\left[\left\vert {\varPsi }_{0}\right\rangle +\left(\left\vert {\varPsi }_{1}(t)\right\rangle -\left\vert {\varPsi }_{1}(0)\right\rangle \right)+{\mathcal{O}}\left({g}^{2}t/J\right)\right],$$ (10) where \({E}_{0}=\left\langle {\varPsi }_{0}\right\vert H\left\vert {\varPsi }_{0}\right\rangle\) is the unperturbed string energy, and the leading non-trivial time dependence, \(\left\vert {\varPsi }_{1}(t)\right\rangle\), describes the spreading dynamics of a conserved number of charge pairs generated by the abrupt perturbation at t = 0.To calculate \(\left\vert {\varPsi }_{1}(t)\right\rangle\), let us split the infinite-spin-chain Hamiltonian as H = H0 + Vdiag + Voffdiag, where the dominant part \({H}_{0}=-J\sum _{i}{\sigma }_{i}^{z}{\sigma }_{i+1}^{z}\) defines subspaces with different numbers of charges; the remainder is split into a block-diagonal component \({V}_{{\rm{diag}}}=\)\(-\mathop{\sum}\nolimits _{i}\mathop{\sum}\nolimits _{r\ge 2}\,{J}_{r}{\sigma }_{i}^{z}{\sigma }_{i+r}^{z}\)\(-h\mathop{\sum}\nolimits _{i}{\sigma }_{i}^{z}-g\mathop{\sum}\nolimits _{{i}_{0}\le i\le {i}_{0}+L-1}\)\(\left({P}_{i-1}^{\uparrow }{\sigma }_{i}^{x}{P}_{i+1}^{\downarrow }+{P}_{i-1}^{\downarrow }{\sigma }_{i}^{x}{P}_{i+1}^{\uparrow }\right)\), which conserves the number of charges (that is, [H0, Vdiag] = 0), and an off-block-diagonal component \({V}_{{\rm{offdiag}}}=\)\(-g\mathop{\sum}\nolimits _{{i}_{0}\le i\le {i}_{0}+L-1}\)\(\left({P}_{i-1}^{\uparrow }{\sigma }_{i}^{x}{P}_{i+1}^{\uparrow}+\right.\) \(\left.{P}_{i-1}^{\downarrow }{\sigma }_{i}^{x}{P}_{i+1}^{\downarrow }\right)\), which creates or destroys pairs of charges (here, \({P}_{i}^{\uparrow }\) (\({P}_{i}^{\downarrow }\)) projects spin i on state \(\left\vert \uparrow \right\rangle\) (\(\left\vert \downarrow \right\rangle\))). Notice that the transverse-field perturbation acts only on the spins in the simulated physical region [i0, i0 + L − 1]; in our experiments, we have L = 13 and choose a labelling of spin sites such that i0 = −(L − 1)/2 = −6.We seek a unitary that transforms the Hamiltonian into a block-diagonal form, order by order in perturbation theory. At the lowest order, a transformation of the form eiS is sought such that \({{\rm{e}}}^{iS}H{{\rm{e}}}^{-iS}={E}_{0}+{H}_{\mathrm{eff}}+{\mathscr{ \mathcal O }}({g}^{2}/J)\), where [Heff, H0] = 0. This requires choosing a block-off-diagonal generator S satisfying i[S, H0 + Vdiag] = −Voffdiag. By computing the time evolution of the transformed initial state \(\mathrm{e}^{iS}\left\vert {\varPsi }_{0}\right\rangle\) under the transformed Hamiltonian, we obtain a perturbative expansion of the time-evolving state,$$\left\vert \varPsi (t)\right\rangle =\mathrm{e}^{-it{E}_{0}}\left[{\mathbb{1}}+i\left(\mathrm{e}^{-it{H}_{{\rm{eff}}}}-{\mathbb{1}}\right)S+{\mathcal{O}}\left({g}^{2}t/J\right)\right]\left\vert {\varPsi }_{0}\right\rangle .$$ (11) From equation (11), we read off$$\left\vert {\varPsi }_{1}(t)\right\rangle =i\mathrm{e}^{-it{H}_{{\rm{eff}}}}S\left\vert {\varPsi }_{0}\right\rangle .$$ (12) This procedure can be iterated to arbitrarily high orders in perturbation theory, to find an accurate description of time evolution over an increasingly long time window. We will show that the lowest-order analysis is sufficient to grasp the qualitative features observed in our experiment.The picture of string dynamics resulting from this time-dependent perturbation theory has a transparent physical interpretation: (1) The abrupt perturbation at time t = 0 generates a quantum superposition of charge pairs at different locations, locally breaking the string; (2) subsequently, these charges undergo quantum-mechanical motion across the physical region during time evolution. Processes 1 and 2 are encoded in the leading time dependence \(\left\vert {\varPsi }_{1}(t)\right\rangle\) in equation (12). Specifically, process 1 is described by the action of the operator S—which is by construction an off-block-diagonal operator—on the classical string state \(\left\vert {\varPsi }_{0}\right\rangle\). At lowest order, S flips a single spin, that is, it creates a single pair of charges adjacent to each other. The amplitude of this process is affected by the position-dependent energy cost of creating a pair, resulting in a spatially inhomogeneous initial wave function of the charge pair (see below for its expression). Process 2 is described by the time-evolution operator \(\mathrm{e}^{-it{H}_{{\rm{eff}}}}\) associated with the number-conserving effective Hamiltonian Heff. These dynamics are affected by energy of the two-charge configurations, which depend on the two positions, resulting in a quantum diffusion through an inhomogeneous potential landscape (see below for its expression).At the lowest order in perturbation theory, string dynamics thus reduce to a two-body problem (for the two charges created), which can be straightforwardly analysed. Its solution provides key insights into string-breaking dynamics and determines the time-dependent probability \({\mathcal{P}}(t)\) of finding the system in a broken-string configuration,$${\mathscr{ \mathcal P }}(t)=2\left(\langle {\varPsi }_{1}(0)| {\varPsi }_{1}(0)\rangle -\,\mathrm{Re}\,\langle {\varPsi }_{1}(0)| {\varPsi }_{1}(t)\rangle \right)+{\mathscr{ \mathcal O }}\left({({g}^{2}t/J)}^{2}\right)$$ (13) (notice that the norm of \(\left|{\varPsi }_{1}(t)\right\rangle\) is perturbatively small).Let us now analyse the two-body problem. We denote by \(\left\vert {l}_{1},{l}_{2}\right\rangle\) the configuration with the two charges sitting at lattice sites l1 and l2, with i0 ≤ l1 < l2 ≤ i0 + L (corresponding to the L + 1 spin bonds involving at least one dynamical spin). Considering an exponentially decaying interaction profile Jr = Je−β(r−1), the energy of the charge-pair configuration relative to the intact classical string, that is, \({V}_{{l}_{1},{l}_{2}}=\langle {l}_{1},{l}_{2}| H| {l}_{1},{l}_{2}\rangle -{E}_{0}\), is given by$$\begin{array}{l}{V}_{{l}_{1},{l}_{2}}=\frac{4J}{{(1-\mathrm{e}^{-\beta })}^{2}}\left(1+\mathrm{e}^{-\beta {l}_{2}}-\mathrm{e}^{-\beta {l}_{1}}-\mathrm{e}^{-\beta ({l}_{2}-{l}_{1})}\right.\\ \left.\qquad+\mathrm{e}^{-\beta (L+2-{l}_{1})}-\mathrm{e}^{-\beta (L+2-{l}_{2})}\right)\\\qquad -2h({l}_{2}-{l}_{1}).\end{array}$$ (14) The initial state of the two charges can then be expressed as$$\left\vert {\varPsi }_{1}(0)\right\rangle =iS\left\vert {\varPsi }_{0}\right\rangle =i\mathop{\sum }\limits_{l={i}_{0}}^{{i}_{0}+L-1}\frac{g}{{V}_{l,l+1}}\left\vert l,l+1\right\rangle ,$$ (15) and the effective Hamiltonian reads$$\begin{array}{l}{H}_{{\rm{eff}}}=\mathop{\sum}\limits _{{l}_{1} < {l}_{2}}{V}_{{l}_{1},{l}_{2}}\,\left\vert {l}_{1},{l}_{2}\right\rangle \left\langle {l}_{1},{l}_{2}\right\vert -g\mathop{\sum}\limits _{{l}_{1} < {l}_{2}}\left(\left\vert {l}_{1}-1,{l}_{2}\right\rangle \left\langle {l}_{1},{l}_{2}\right\vert\right.\\\left.\quad\quad\;+\left\vert {l}_{1},{l}_{2}+1\right\rangle \left\langle {l}_{1},{l}_{2}\right\vert +{\rm{h.c.}}\right).\end{array}$$ (16) Dynamics governed by this Hamiltonian can be viewed as a single particle hopping in a two-dimensional lattice of triangular shape with uniform amplitude g and hard-wall boundaries, under the influence of an inhomogeneous potential given by equation (14). Higher-order terms contain corrections to the potential and longer-range hopping processes.Figure 4a–c illustrates this description of charge dynamics, with a colour map of the potential \({V}_{{l}_{1},{l}_{2}}\). The appearance of negative values of the potential in the top-right corner for sufficiently large h reflects static string breaking, that is, a change from a string-like to a broken-string-like ground state. The initial state is represented by a wave function supported along the diagonal side of the domain (shaded dots), with amplitude peaks at the two edges. From these initial positions, the particle may then quantum-mechanically propagate through the landscape. Superimposed arrows in Fig. 4a–c highlight approximate equipotential lines, which appear as the string tension h is increased, thereby opening up a channel for charge pairs generated at the edges to spread into the bulk. This mechanism underlies the qualitative features of the spatiotemporal string-breaking patterns observed in our experiment (Fig. 3).Instances of spatiotemporal charge and electric-field evolution as approximated by the two-body time-dependent Schrödinger equation are shown in Extended Data Fig. 8a–f for several values of h/J and the lowest value g/J = 0.75 for the coupling strength in our experimental simulations. These parameters are far from the theoretical perturbative limit, where the approximation would be quantitatively accurate. Nonetheless, the dynamical behaviour exhibits qualitative similarity with the exact string dynamics, approximately capturing, in particular, the observed dynamical string-breaking mechanism.Long-time string dynamics, far beyond the transient scale 0 < t ≪ J/g2, are expected to be described by the Schwinger mechanism—a mechanism associated with the spontaneous materialization of resonant charge pairs. The spatiotemporal rate of this non-perturbative bulk process is suppressed faster than any power of g (refs. 9,54,63), precluding observation over the accessible timescales in our experiment. Within the picture in Fig. 4a–c, resonant charge-pair configurations satisfy \({V}_{{l}_{1},{l}_{2}}\approx 0\) and are indicated in the plots by black stars.Continuum limitThere are two possible continuum descriptions of the novel edge-facilitated string-breaking mechanism: (1) the ‘standard’ continuum-field-theory limit as the system approaches its quantum critical point [g = gc, h = 0], where the correlation length diverges, and a continuum relativistic quantum field theory of Majorana fermions emerges (see, for example, the review in ref. 75); (2) the limit where the spin–spin interactions Ji,i+r = Je−β(r−1) have a long range compared with the lattice spacing, β−1 ≫ 1, and simultaneously the overall strength J becomes weak, in such a way that the total ferromagnetic interaction energy is kept constant, that is, \(J/{\beta }^{2}\equiv \,\text{const}\,\equiv \bar{J}\): this is analogous to a Kac limit in statistical mechanics. Loosely speaking, the former approach takes a continuum limit of the charges’ kinetic energy, whereas the latter approach takes a continuum limit of the charges’ potential energy. We note that the latter route has immediate relevance to our quantum-simulation experiment58.Here, we discuss continuum limit 2 within the effective two-body problem analysed above. On taking β → 0 and L → ∞ with ℓ ≡ βL fixed, and defining continuum spatial variables x1,2 = βl1,2 ∈ [ − ℓ/2, ℓ/2], one obtains the following two-body Hamiltonian:$$H({x}_{1},{x}_{2},{p}_{1},{p}_{2})\,\to \,-2g\cos {p}_{1}-2g\cos {p}_{2}+V({x}_{1},{x}_{2}),$$ (17) with the potential energy$$\begin{array}{l}V({x}_{1},{x}_{2})=4\overline{J}\left({{\rm{e}}}^{-{x}_{2}}-{{\rm{e}}}^{-{x}_{1}}+{{\rm{e}}}^{-(\ell -{x}_{1})}-{{\rm{e}}}^{-(\ell -{x}_{2})}\right.\\ \left.\qquad\qquad+1-{{\rm{e}}}^{-({x}_{2}-{x}_{1})}\right)-2\overline{h}({x}_{2}-{x}_{1}),\end{array}$$ (18) where the parameters, defined as \(\bar{J}\equiv {\beta }^{-2}J\), \(\bar{h}\equiv {\beta }^{-1}h\), take finite values as β → 0. Note that [x1, p1] = [x2, p2] = iβ, so wave-packet dynamics are semiclassical in this limit.Extended Data Fig. 8g–i shows contour plots of the potential energy, together with the equipotential trajectories originating from the initial conditions localized at x1, x2 ≈ − ℓ/2 or at x1, x2 ≈ ℓ/2, corresponding to a pair of charges created near the left or right string edge, respectively. The highlighted equipotential trajectories make the appearance of an edge-facilitated breaking mechanism particularly clear in this continuum limit.Testing the predicted string-breaking scenario in the scaling region close to the quantum critical point, that is, in the continuum limit 1 above, would probably require adiabatic state preparation of a near-critical system with a large number of spins, so that the hierarchy of length scales 1 ≪ ξ ~ β−1 ≪ L be in place, where \(\xi =1/(\beta \mu ) \sim {({g}_{{\rm{c}}}/g-1)}^{-1}\) is the correlation length in lattice spacing units. This regime is unfortunately beyond the scope of our current experiment.Thermal expectation valuesHere, we outline the calculation of thermal expectation values of the electric field in the main text (Fig. 3). The state of a quantum many-body system far from equilibrium is described by a superposition of a large number of high-energy eigenstates of the system’s Hamiltonian. In generic interacting systems, local observables are expected to relax to stationary values compatible with a thermal-ensemble average. The corresponding thermal-equilibrium state is given by the canonical Gibbs ensemble, \({\rho }_{{\rm{Gibbs}}}(T)=\mathrm{e}^{-H/({k}_{B}T)}/{\rm{Tr}}\left[\mathrm{e}^{-H/({k}_{B}T)}\right]\), with Boltzmann factor kB and with the temperature T fixed by the total energy. For each of the systems analysed in Fig. 3, the temperature is determined by solving \(\left\langle {\varPsi }_{0}\right\vert H\left\vert {\varPsi }_{0}\right\rangle ={\rm{Tr}}\left[H{\rho }_{{\rm{Gibbs}}}(T)\right]\), where \(\left\vert {\varPsi }_{0}\right\rangle\) is the classical-string state depicted in Fig. 1c. The thermal-equilibrium profiles for the electric field displayed in Fig. 3j–l are computed via \({\langle {\sigma }_{z}^{i}\rangle }_{{\rm{thermal}}}={\rm{Tr}}[{\sigma }_{z}^{i}{\rho }_{{\rm{Gibbs}}}(T)]\).Experimental apparatus and Hamiltonian engineeringIn this section, we first introduce the basic experimental features of the trapped-ion apparatus used in this work. We then describe the methods used to generate the programmable spin–spin interaction and the three-dimensional (3D) magnetic field used to engineer the Hamiltonian in equation (1).Trapped-ion apparatusThe trapped-ion apparatus used in this work features a linear array of L = 15 171Yb+ ions trapped in a microfabricated surface ion trap 56,57,58. We encode the spin states in the two clock levels of the 2S1/2 ground state, where \({\left\vert \uparrow \right\rangle }_{x}\equiv \left\vert F=1,{m}_{F}=0\right\rangle\) and \({\left\vert \downarrow \right\rangle }_{x}\equiv \left\vert F=0,{m}_{F}=0\right\rangle\) with an energy separation of ω0 = 2π × 12.6 GHz. At the beginning of each experimental cycle, the ions are cooled using Doppler cooling. The transverse motional modes used to generate long-range interactions are further cooled down to their motional ground state using resolved-sideband cooling. The initial state is prepared by optically pumping the ions to the \({\left\vert \downarrow \cdots \downarrow \right\rangle }_{x}\) state, and rotating them to the z-basis with Ry(−π/2) gates, that is, single-qubit rotations around the y axis of the Bloch sphere by angle −π/2. Such gates are implemented by tuning the Raman beatnote to drive the carrier transition. After the state preparation, a spin-spin interaction is generated to realize the Hamiltonian in equation (1) (see the following subsections), and the spins evolve under this Hamiltonian for a variable duration. We measure the final spin-state population along the z basis by applying a Ry(−π/2) gate, followed by spin-dependent fluorescence detection.The optical set-up consists of two global beams and a dual array of tightly focused laser beams, each with independently controllable optical frequency, phase, and amplitude. All beams originate from a single optical frequency comb pulsed laser at 355 nm. The global beams are generated by the passage of one beam through an acousto-optic modulator (AOM), which splits the beam into two spatially and frequency-distinct components by shifting the radio-frequency (RF) beatnotes by \({\omega }_{0}^{\,\mathrm{glob}}\pm ({\omega }_{L}+\mu )\). ωL is the frequency of the lowest-frequency radial phonon mode, the zig-zag mode, and μ is the detuning of the Raman beatnote relative to this phonon mode. These two components are combined and directed onto the ions via a telecentric optical path, ensuring uniform illumination across the ion chain.The dual beam arrays originate from a single beam that passes through a diffractive optical element, initially generating an array of more than 30 tightly focused laser beams. Each beam then passes through an independent channel of a multichannel AOM, which applies two RF frequencies, generating two first-order diffracted beams that are spatially separated, with controllable frequencies, amplitudes and phases. The optical set-up is precisely designed and aligned to overlap the two distinct beam arrays onto the ion plane, ensuring that each pair of tightly focused laser beams addresses the ions perpendicular to the combined global beam. The system is stabilized using an active vibration control platform to minimize the sensitivity of the optical phase to acoustic vibrations. We routinely calibrated the relative phase at the ion positions to account for any changes or drifts in the optical phase between the pairs of beams in the dual array, due to differences in their optical paths. In addition, by controlling the radial trapping potential, we aligned the wavevector difference between the (combined) global beam and the dual beam array so that only a single set out of L phonon modes, associated with a specific direction of motion, is addressed by the beams.Generation of spin–spin interactionThe spin–spin Ising interaction in the trapped-ion quantum simulator is generated by applying spin-dependent forces, achieved by simultaneously driving the red and blue motional sidebands. This is facilitated by the two components of the global beams and one array of tightly focused laser beams. The second array is used to generate effective magnetic fields, which will be described in the following subsections. The multichannel AOM shifts the frequency of each beam in the first array by \({\omega }_{0}^{\,\text{ind}\,}\), which satisfies the condition \(| {\omega }_{0}^{\,\text{glob}}-{\omega }_{0}^{\text{ind}\,}| ={\omega }_{0}\). Consequently, the Raman beatnotes ω0 ± (ωL + μ) of the global beam, along with the first array, simultaneously drive a blue and red motional sideband transition, generating an effective spin–spin Ising interaction \({H}_{ZZ}=\sum _{i,j}{J}_{i,j}{\sigma }_{i}^{z}{\sigma }_{j}^{z}\) via the Mølmer–Sørensen protocol76. The experimental Ji,j matrix for L spins has the form56,76$${J}_{i,j}=\sum _{k=1}^{L}\frac{{\eta }_{i,k}{\eta }_{j,k}{\varOmega }_{i}{\varOmega }_{j}}{{\omega }_{L}+\mu -{\omega }_{k}},$$ (19) where ηi,k = 0.08bi,k are the Lamb–Dicke parameters describing the coupling between spin i and motional mode k, with bi,k being the mode-participation matrix elements. Ωi is the carrier Rabi frequency at ion i, and ωk is the radial motional frequency of mode k, labelled in decreasing order with 1 ≤ k ≤ L.By precisely tuning the axial trap potential56,58, we achieve a nearly uniform ion spacing of 3.75 μm for the 13 central ions in the 15-ion chain. In this experiment, we control the phase and amplitude of the 13 individual beams of the first array via the RF signals driving the corresponding AOM channels, while keeping the remaining beams off. Therefore, the two edge ions assist in trapping and participate in the motion, but they are not driven by the Raman interaction, and their spins do not contribute to the Hamiltonian in equation (1). This capability allows for varying the system size without changing the trapping potential56, which might otherwise alter the interaction matrix77.One can numerically calculate the interaction matrix elements Ji,j in equation (19) using the axial trapping parameters. We find that these elements approximately decay with a functional form best described by Jr = Je−β(r−1)r−α, where r = ∣i − j∣, and the average interaction strength is \({J}_{r}=\frac{1}{L-r}\sum _{i}\left\vert \,{J}_{i,i+r}\right\vert\). The range of interaction can be adjusted by changing the symmetric detuning ∣μ∣ of the global beam relative to the red and blue sideband—for large detuning, the exponential term dominates at short distances, and for small detuning, the interaction becomes longer range56. For our chosen trap parameters, the zig-zag mode frequency is ωL = 2π × 2.78 MHz, and the detuning is μ = −2π × 35 kHz. For this configuration, we find that the numerically calculated interaction matrix follows the decay pattern of Jr with α = 0 and β = 0.78.The participation matrix elements (bi,L) of the zig-zag mode have alternating signs for even and odd sites, which causes the Ji,j matrix elements to be staggered, that is \(\,\text{sgn}\,\left({J}_{r}\right)=-\,\text{sgn}\,\left({J}_{r+1}\right)\). To correct this form of the interaction and ensure that all the Ji,j elements have the same sign, we shift the optical phase of every other beam in the first array by π while generating the spin–spin interaction.The amplitudes of the mode-participation matrix of the zig-zag mode are non-uniform (Extended Data Fig. 2a, black triangles), leading to a non-uniform interaction matrix across the chain. To correct for this non-uniformity, we adjust the amplitudes of the individual beams by controlling the RF amplitude driving each AOM channel. For a chain of L spins, we experimentally measure L − 1 nearest-neighbour (NN) and L − 2 next-nearest-neighbour (NNN) couplings. By tuning the individual beam amplitudes, we achieve approximately uniform NN couplings and minimize fluctuations in the NNN couplings. The Rabi frequencies of the individual beams are represented by purple dots in Extended Data Fig. 2a. A comparison of the mode-participation factors with the beam-amplitude profiles indicates that the amplitude profile is the inverse of the participation factors. After applying this amplitude profile, up to a scaling factor, we experimentally observe near-uniform NN and NNN couplings as shown in Extended Data Fig. 2b. Across the chain, the average NN coupling reads J = J1 = 2π × 0.34(1) kHz, and NNN coupling reads J2 = 2π × 0.10(2) kHz, where the error indicates the standard deviation across the chain.When the spin–spin interaction is quickly turned on, an off-resonant coupling to the carrier transition can be detected. To avoid this undesired coupling, the individual beam amplitudes are linearly ramped on with a 5 μs duration. This ramp time is substantially longer than the off-resonant carrier oscillation time scale \(1/\sqrt{{\varOmega }_{i}^{2}+{({\omega }_{L}+\mu )}^{2}}\approx 1/({\omega }_{L}+\mu )\), but remains much shorter than 1/J1. As a result, the spin–spin interaction is still effectively turned on instantaneously.Generation of transverse magnetic fieldIn addition to the Ising interaction, an effective local transverse field is applied by adjusting the frequencies of individual beams. Specifically, if the optical beatnote frequencies addressing ion i are shifted by gi relative to the qubit resonance, this detuning induces an effective static local transverse field term \({g}_{i}{\sigma }_{i}^{x}\) in the Hamiltonian, under the condition gi ≲ J (refs. 56,78). We realize this term by shifting the RF frequencies of the dual-array beams. The use of inhomogeneous beam amplitudes, as discussed earlier, inevitably results in an inhomogeneous a.c. Stark shift across the ion chain. These shifts manifest as unwanted \({\epsilon }_{i}{\sigma }_{i}^{x}\) terms in the interaction Hamiltonian. We experimentally measure these shifts and apply an equal but opposite frequency shift for each tightly focused beam in the dual arrays, effectively generating Hamiltonian terms \(-{\epsilon }_{i}{\sigma }_{i}^{x}\), which exactly cancel the Stark shift effects and ensure a uniform transverse field across the chain.Due to the slight misalignment of the spatial overlap between the two frequency components of the global beam, and their finite sizes, the power imbalance between these components results in a large Stark shift error on the edge ions. To mitigate this, we implemented dynamical-decoupling pulses on the edge spins in the 13-spin chain, minimizing their response to shot-to-shot fluctuations in the Stark shifts79.Generation of longitudinal magnetic fieldThe dual-array set-up enables the implementation of fully programmable local longitudinal fields, applied simultaneously with the transverse-field Ising Hamiltonian terms. To achieve this, we introduce additional frequency shifts of \({\omega }_{0}^{\,\text{ind}\,,i}+{\omega }_{L}+\mu\) in the multichannel AOM, generating a second array of tightly focused beams. Through a telecentric beam path, they are overlapped with the beams generating the spin-spin interactions, as shown in Fig. 1e. The second array of beams, along with the \({\omega }_{0}^{\,\text{glob}\,}+{\omega }_{L}+\mu\) frequency component of the global beam, produces a beatnote that drives the qubit transition on resonance at frequency ω0. This is equivalent to a site-dependent magnetic field \({h}_{i}^{\phi }(t){\sigma }_{i}^{\phi }\), where \({\sigma }_{i}^{\phi }=\cos \phi \,{\sigma }_{i}^{z}+\sin \phi \,{\sigma }_{i}^{y}\), enabling control of all three components of the local magnetic field of the ith spin.For each individual beam, we calibrate the Rabi frequency versus the RF amplitude on the AOM, and additionally shift the phase of the RF signal to align the field in the z axis (owing to a difference in the optical path relative to the first array) such that \({\sigma }_{i}^{\phi }={\sigma }_{i}^{z}\). This allows us to apply the tunable inhomogeneous longitudinal-field profiles necessary for the experiments.In this work, we use two different configurations of virtual static spins for studying isolated-charge dynamics (Fig. 1b) and dynamical string breaking (Fig. 1c,d). The inhomogeneous fields corresponding to these configurations are, respectively, given by$$\Delta {h}_{i}^{({\rm{charge}})}=\mathop{\sum }\limits_{r=i}^{\infty }{J}_{r}-\mathop{\sum }\limits_{r=L+1-i}^{\infty }{J}_{r},$$ (20) $$\Delta {h}_{i}^{({\rm{string}})}=-{J}_{i}+\mathop{\sum }\limits_{r=i+1}^{\infty }{J}_{r}-{J}_{L+1-i}+\mathop{\sum }\limits_{r=L+2-i}^{\infty }{J}_{r},$$ (21) with Jr = Je−0.78(r−1) in our specific set-up. The numerical and measured values associated with the experimental parameters are reported in Extended Data Fig. 2c,d.