Toward simulating quantum field theories with controlled phonon-ion dynamics: A hybrid analog-digital approach

Quantum field theories are the cornerstones of modern physics, providing relativistic and quantum mechanical descriptions of physical systems at the most fundamental level. Simulating real-time dynamics within these theories remains elusive in classical computing. This provides a unique opportunity for quantum simulators, which hold the promise of revolutionizing our simulation capabilities. Trapped-ion systems are successful quantum-simulator platforms for quantum many-body physics and can operate in digital, or gate-based, and analog modes. Inspired by the progress in proposing and realizing quantum simulations of a number of relativistic quantum field theories using trapped-ion systems, and by the hybrid analog-digital proposals for simulating interacting boson-fermion models, we propose hybrid analog-digital quantum simulations of selected quantum field theories, taking recent developments to the next level. On one hand, the semi-digital nature of this proposal offers more flexibility in engineering generic model interactions compared with a fully-analog approach. On the other hand, encoding the bosonic fields onto the phonon degrees of freedom of the trapped-ion system allows a more efficient usage of simulator resources, and a more natural implementation of intrinsic quantum operations in such platforms. This opens up new ways for simulating complex dynamics of e.g., Abelian and non-Abelian gauge theories, by combining the benefits of digital and analog schemes.


I. INTRODUCTION
Quantum field theories (QFTs) provide the underlying quantum-mechanical descriptions of physical systems, from relativistic gauge field theories of the Standard Model of particle physics [1,2], to emergent low-energy models in condensed-matter systems [3,4], to effective field theories in hadronic and nuclear physics [5,6]. Classical simulation methods have come a long way to describe phenomena emerging from these underlying theories, with notable examples in the realm of lattice gauge theory (LGT) methods applied in stronginteraction physics [7,8]. Nonetheless, there is a need for new computational strategies to overcome the limitations of the current methods, in order to achieve real-time simulations of matter, and predictions for equilibrium and out-of-equilibrium phenomena arising from stronginteraction dynamics. Hamiltonian simulation of physical systems, naturally enabled via mapping the problem to a quantum simulator or a quantum computer, is an example of such a strategy, with major advances reported in recent years in proposing  and implementing [53][54][55][56][57][58][59][60][61][62][63][64][65][66][67][68] simulations of various field theories on the limited existing quantum hardware, and in formulating suitable Hamiltonian descriptions of QFTs and LGTs for resource-efficient and robust quantum simulations [9, 12, 13, 21, 29, 36, 42-45, 47-50, 69-74]. The question we investigate here is complementary to the latter effort: Given a Hamiltonian formulation, what are the efficient implementation strategies that maximally take advantage of available hardware capacity and capability? We focus on trapped-ion quantum simulators as one of the leading quantum-hardware technologies. Nonetheless, applications of this strategy to other platforms, such as cavity QED and superconducting circuits, can be established analogously.
QFTs with bosonic fields are extremely common, from scalar field theory descriptions of phase transitions in condensed-matter systems or of the inflationary phase of the early universe and the Higgs mechanism, to gauge field theory descriptions of the fundamental forces of nature. Due to their infinite-dimensional Hilbert space, even on a single space-time point, a truncation must be imposed on the various bosonic excitations to contain their dynamics on a finite simulating hardware. Important work has emerged in understanding, qualitatively or quantitatively, the impact of this truncation, from the simple quantum harmonic oscillator [75][76][77] to scalar field theory [13,29] and LGTs [47,50,74,78]. While quantities in the low-energy subspace of the truncated theories exhibit exponential convergence to the exact values in the full theory, as guaranteed by the Shannon-Nyquist sampling theorem [29,76] in the case of scalar field theory, observables in the high-energy sector or those obtained from long-time expectation values exhibit a slow convergence to the exact values. In particular, as shown in Ref. [74] for the case of the SU (2) LGT coupled to matter in 1+1 D, for a certain threshold on the gauge-2 field cutoff quantities enter a scaling region in which the dependence on the cutoff becomes exponential and systematically improvable. However, for high-energy spectra and long-time dynamical quantities, such a threshold is reached for a large value of the cutoff, increasing the simulation-resource requirement needed to reach a fixed accuracy.
The most common encoding of bosonic degrees of freedom to qubits for digital quantum simulation is to convert their occupation number to a binary representation, the qubit-number requirement of which is logarithmic in the cutoff on the highest excitation of the boson. Most importantly for near-term applications, the number of entangling gates required to implement the dynamics associated with boson-boson or boson-fermion interactions grows polynomially with the system size, as a large number of controlled operations must be introduced on the state of the boson register, see e.g., Refs. [39,75,77]. This problem is particularly severe in the case of bosonic field theories as, for example, a scalar field in a finite discretized space represents quantum harmonic oscillators, or bosons, associated with each momentum mode, the number of which scales linearly with the volume of the reciprocal lattice. Each harmonic oscillator then requires the allocation of dedicated qubits and the implementation of associated entangling gates. A cost analysis of this kind for a fully-digital implementation is made explicit in this paper and is contrasted with our analog-digital protocol, to be introduced next.
A trapped-ion quantum simulator [79][80][81] can operate in an analog mode, with which dynamics generated by the Ising and Heisenberg spin Hamiltonians can be studied [32,[82][83][84][85][86][87][88][89][90], or a digital mode, with which universal computations expressed in terms of single-and two-qubit gates can be performed [54,[91][92][93][94][95][96][97][98][99][100][101][102][103][104]. Typically in these systems, specific pairs of internal states encode qubits, which are manipulated via ion-laser or ion-microwave interactions. The collective excitations of the motion of ions in the trap, i.e., the phonons, are often used as mediators of the interactions among the qubits, enabling the engineering of effective spin Hamiltonians or in turn entangling operations among two or more qubits [81]. In recent years, promising progress has been reported in controlling the phonon dynamics in a trapped-ion simulator, leading to experimental demonstrations of phonon hopping, phonon interference, and a phonon quantum walk in Paul traps [105][106][107][108][109]. This opens up the possibility of taking advantage of phonons, in addition to the qubits, as dynamical degrees of freedom to encode and process information [110,111]. In fact, there exist quantum-simulation proposals for an Abelian LGT, i.e., the lattice Schwinger model, which take advantage of the phonons in a fully-analog setting to encode the gauge or fermion degrees of freedom [112], but they are not experimentally feasible yet. A more viable option for near-term applications is to combine the flexibility of gate-based digital simulations with the versatility of both the controllable spin and phonon degrees of free-dom in a trapped-ion simulator [14,[113][114][115]. This idea has led to concrete gate-based protocols for simulating interacting fermion-boson models, such as the Holstein model of electron-phonon dynamics in condensed-matter physics [115], as well as the first experimental implementation of a one-site boson-fermion dynamics [116] considering only a single excitation of the boson. In the present work, the ideas in these proposals are taken to the next level to demonstrate that quantum simulation of relativistic QFTs, including LGTs, can benefit from a similar hybrid protocol.
This work introduces a complete and enhanced set of conventional as well as phonon-based gates, including but not limited to those introduced in Ref. [115], and deploys them to construct concrete circuits for the digitized time-evolution operator within a Yukawa theory of scalar-fermion interactions and the U (1) LGT coupled to matter. The gates include single-spin (qubit) operations, and entangling spin-spin, spin-phonon, and phonon-phonon operations. In order to ensure efficient implementation of all these gates in one experiment, normal modes of motion, local modes of motion, or both types of modes should be controlled and manipulated. The simulated theories considered are in 1+1 dimensions (D) but the generalization to higher dimensions, and in particular challenges involved in the case of higherdimensional LGTs, are briefly discussed. Numerical examples are provided to supplement our proposals with concrete experimental parameters, and to demonstrate their feasibility within the current technologies. Qualitative comparisons of the simulation cost within digital and hybrid analog-digital simulations of the same theories will be provided, along with a discussion of the outlook of a hybrid approach for generalization to more complex QFTs, including non-Abelian LGTs.

II. THE HYBRID ANALOG-DIGITAL BUILDING BLOCKS
The trapped-ion quantum simulator considered in this work consists of a number of ions confined in a radiofrequency Paul trap [117]. The qubit is encoded in two stable internal levels of the ion, which are separated in energy by an angular frequency ω 0 . 1 The confining potential is sufficiently stronger along the transverse axes, x and y, of the trap so that the ions form a linear crystal in the axial, z, direction. 2 With appropriate anharmonic axial confinement forces, the ions can be arranged in 3 an equally-spaced configuration for individual laser-beam addressing [118,119]. The typical spacing between adjacent ions is a few micrometers in present-day trapped-ion simulators.
Given the long-range Coulomb force among the ions and the common trapping potential applied, the motion of the ions can be described in terms of a set of collective normal modes. The displacement of the ions around an equilibrium position can then be expressed in terms of phononic degrees of freedom, whose excitation energies are quantized in units of the normal-mode frequencies. While such a normal-mode picture is sufficient for simulating the dynamics of certain fermion-boson field theories, switching to a local-mode picture is required once the bosonic field exhibits non-trivial dynamics, as is the case in the gauge-theory example considered. Local modes can be addressed by making the trapping potential tighter along the transverse directions, or conversely by relaxing the axial confinement to increase the spacing among the ions, such that the displacement of any given ion from its equilibrium position is much smaller than the distance among the adjacent ions [120]. The local modes can be addressed as long as the laser excitation on individual ions happens on timescales faster than the separation of the corresponding normal modes, which defines the timescale for phonon hopping. 3 In such a regime, the phonon hopping probability among the ion sites is small and the quanta of motion are localized. The interactions associated with each of these scenarios, along with the expanded gate-set for the hybrid analog-digital computation of this work, will be introduced in the following.

A. Phonons as excitations of normal modes of motion
The free Hamiltonian of the trapped-ion system in the absence of laser-ion interactions is Here, σ is a Pauli operator acting on the space of the ions' quasi-spin, i.e., the qubit, with the standard commutation relations. a † m (a m ) are the bosonic creation (annihilation) operators for the normal modes of motion with the associated frequencies ω m and the commutation relations [a m , a m ] = [a † m , a † m ] = 0 and [a m , a † m ] = δ m,m . We label the transverse modes of motion along the x-axis of the trap by indices m = 1, 2, · · · , N , those along the y axis of the trap by indices m = N +1, N +2, · · · , 2N , and 3 One may also consider a micro-trap where the global trapping potential is replaced by local potentials applied at the position of each ion [121]. For applications studied in this work, a Paul trap with a tight transverse potential suffices to produce the desired dynamics, and this additional possibility will not be considered.
the axial modes, i.e., those along the z axis of the trap, by indices m = 2N + 1, 2N + 2, · · · , 3N . With the introduction of two counter-propagating laser beams with wave-vector difference ∆k j , frequency difference (beatnote) ∆ω j ≡ ω L j , and phase difference ∆φ j ≡ φ j , detuned from an excited internal level of the ion, the two-photon Raman transition among the two states of the qubit can be induced with a Rabi frequency Ω j at the location of ion j. The beams are assumed to address the ions individually, hence the subscript j. In an interaction picture that rotates with the free Hamiltonian in Eq. (1), the interacting ion-laser Hamiltonian can be written as [84] where the condition |ω L j − ω 0 | ω 0 is assumed. The prime on H is to denote that this Hamiltonian is in the interaction picture. Later on, we need to adopt a different interaction picture rotating with shifted frequencies and so it is important to bear in mind the origin of Eq. (2). Multiple beatnote frequencies, amplitudes, and phases can be applied simultaneously, hence a sum over laser parameters can be introduced, a possibility that we take advantage of later. The Lamb-Dicke parameters η m,j are defined as η m,j = (∆k) 2 2Mionωm b m,j , where b m,j are the (normalized) normal-mode eigenvector components between ion j and mode m, and M ion denotes the mass of the ion. In the Lamb-Dicke regime in which η m,j (a m + a † m ) 2 1/2 1, transitions in the space of coupled spin-phonon system take a simpler form, and can be realized through quantum gates. ∆k j is assumed to be the same for lasers addressing each ion j, i.e., ∆k j ≡ ∆k.
The carrier transition corresponds to H ion-laser | O(η 0 ) and is obtained by setting ω L j −ω 0 = 0. With this setting, the Hamiltonian corresponding to ion j becomes The superscript on H σ j denotes that this Hamiltonian only acts on the qubit space with the noted σ operators. The single-spin rotations of arbitrary angle θ j , to be denoted as R σ j (θ j , φ j ), can be obtained by applying H σ j (φ j ) for time τ σ gate , to be deduced from the relation For example, rotations around the x and y axes of the Bloch sphere of the qubit j can be realized as The superscript on H j denotes that this Hamiltonian only acts on the qubit space with the noted operators. The single-qubit rotations of arbitrary angle ✓ j , to be denoted as R j (✓ j , j ), can be obtained by applying H j ( j ) for time ⌧ gate , to be deduced from the relation For example, rotations around the x and y axes of the Bloch sphere of the qubit j can be realized as respectively. Rotations along the z axis of the Bloch sphere, R Z j (✓ j ), defined as can be implemented without the need for quantum operation using a classical shift [Norbert to verify]. The blue and red sideband transitions correspond to H ion-laser 0 |O(⌘) and are obtained by setting ! L ! 0 = ! k and ! k , respectively. This setting leads to a simple coupled spin-phonon Hamiltonian. In order to achieve a Hamiltonian proportional to j x , both sets of blue and red sideband operations can be applied simultaneously with a set of Raman beams such that the phases of the two sideband beams add up to ⇡. With this setting, the Hamiltonian corresponding to ion j becomes In order to ensure only one of the phonon modes, a k , couples to the qubit in this setup, and to minimize the o↵-resonant transitions, the Raman beams could be set such that they only couple to one set of normal modes of motion. For example, by setting k = kx, the beams only couple to transverse normal modes along the x axis. The qubit-phonon rotations of arbitrary angle ✓ k,j , to be denoted as R a k,j (✓ k,j , j ), can be obtained by applying H a k,j ( j ) for time ⌧ a gate , to be deduced from the relation In particular, setting laser beam phases to proper angles can generate transitions proportional to a k + a † k and a k a † k :

2
The superscript on H j denotes that this Hamiltonian only acts on the qubit space with the noted operators. The single-qubit rotations of arbitrary angle ✓ j , to be denoted as R j (✓ j , j ), can be obtained by applying H j ( j ) for time ⌧ gate , to be deduced from the relation For example, rotations around the x and y axes of the Bloch sphere of the qubit j can be realized as respectively. Rotations along the z axis of the Bloch sphere, R Z j (✓ j ), defined as can be implemented without the need for quantum operation using a classical shift [Norbert to verify]. The blue and red sideband transitions correspond to H ion-laser 0 |O(⌘) and are obtained by setting ! L ! 0 = ! k and ! k , respectively. This setting leads to a simple coupled spin-phonon Hamiltonian. In order to achieve a Hamiltonian proportional to j x , both sets of blue and red sideband operations can be applied simultaneously with a set of Raman beams such that the phases of the two sideband beams add up to ⇡. With this setting, the Hamiltonian corresponding to ion j becomes In order to ensure only one of the phonon modes, a k , couples to the qubit in this setup, and to minimize the o↵-resonant transitions, the Raman beams could be set such that they only couple to one set of normal modes of motion. For example, by setting k = kx, the beams only couple to transverse normal modes along the x axis. The qubit-phonon rotations of arbitrary angle ✓ k,j , to be denoted as R a k,j (✓ k,j , j ), can be obtained by applying H a k,j ( j ) for time ⌧ a gate , to be deduced from the relation In particular, setting laser beam phases to proper angles can generate transitions proportional to a k + a † k and a k a † k : Implementing simultaneous blue and red sideband transitions now detuned from the normal mode frequencies will not generate dynamical phonons as long as the in this process then e↵ectively induce a spin-spin transition, which is well known as the Molmer-Sorenson transition []. 4 The corresponding Hamiltonian between ions j and j 0 is where R ⌘ ( k) 2 /2M ion is the recoil frequency of the ion with M ion being the mass of the ion, and µ ⌘ ! L ! 0 . The qubit-phonon rotations of arbitrary angle ✓ j,j 0 , to be denoted as R j,j 0 (✓ j,j 0 ), can be obtained by applying H j,j 0 for time ⌧ gate , to be deduced from the relation In reality, designing high-fidelity two-qubit entangling gates of this type requires complex pulse shaping techniques to achieve a perfect creation and annihilation of the virtual phonon in the process, which is necessary to leave the system in the same phonon-excitation state before and after the gate operation.

B. Phonons as excitations of local modes of motion
In a linear Paul trap the Coulomb energy is comparable to the potential energy associated with the axial motion of the ions. Namely, the parameter z ⌘ where ! z is the trap frequency along the axial axis of the trap, d 0 is the average spacing among the ions, and e is the electric charge. As a result, the phonons hop over the chain and are considered as the common normal modes of the system. The same parameter is typically much smaller for the transverse motions of the ions. In particular, if the trap is made extremely binding along one or both transverse axes of the ions such that x(y) ⌧ 1, the transverse phonons are bound to a single ion, with a small strength for hopping to the neighboring ion that goes as / e 2 /(d 3 0 M ion ! x(y) ) = x(y) ! x(y) []. For the quantum simulation of quantum field theories that exhibit non-trivial boson-field dynamics, it is necessary to adopt a local mode picture. Assuming only one of the transverse directions of the trap, e.g., x, supports local phonon modes, the free Hamiltonian of the ion chain can be written as ing wave beam and, e ⌘ ⌘ kx = k corresponding Lamb-Dickie parameter. standing wave is only exciting the local tion given the choice of the wave-vector in Eq. (15) denotes either higher-order te will be neglected in the regime e ⌘ h(a x j + a are phonon non-conserving operators that least the frequency ! x and can be adia nated as long as F e ⌘ 2 /! x ⌧ 1 []. Similar to the previous gates, phonon-ph of arbitrary angle can be obtained by app time ⌧ aa gate , to be deduced from the relatio An important point to notice regarding t according to Eq. (15), the relative size of of the a x j † a x j and (a x j † a x j ) 2 terms, (1) / ( 1 + e ⌘ 2 )/e ⌘ 2 , while the overall strength can be tuned arbitrarily by changing F , condition F e ⌘ 2 /! x ⌧ 1 is not violated. T dition may seem too constraining for sim bitrary Hamiltonian containing both pho will be demonstrated in Sec. IV for simula Schwinger model, by arbitrarily and cons ing the reference frame in obtaining the ture, desired ratios of the coe cients in theory can still be engineered. Last but n simulated theory requires site-dependent the phonon-phonon couplings, individual beams can replace the global beam in Eq III. A YUKAWA THEORY: SCAL COUPLED TO (STAGGERED) FE A scalar field theory coupled to fermion coupling of the Higgs boson to fermion dard Model through Yukawa interactions sible for dynamical mass generation an mixing in nature. Its non-perturbative lattice-regularized Euclidean field theory drawn considerable interest as it reveals nections to the cuto↵ scale of the Stand questions regarding triviality []. In the mology and early universe, non-perturbat are required for studies of non-equilibrium dynamics of this theory and its beyond B. Phonons as excitations of local modes of motion In a linear Paul trap the Coulomb energy is comparable to the potential energy associated with the axial motion of the ions. Namely, the parameter z ⌘ e 2 /(d 3 0 M ion ! z2 ) & 1, where ! z is the trap frequency along the axial axis of the trap, d 0 is the average spacing among the ions, and e is the electric charge. As a result, the phonons hop over the chain and are considered as the common normal modes of the system. The same parameter is typically much smaller for the transverse motions of the ions. In particular, if the trap is made extremely binding along one or both transverse axes of the ions such that x(y) ⌧ 1, the transverse phonons are bound to a single ion, with a small strength for hopping to the neighboring ion that goes as / e 2 /( For the quantum simulation of quantum field theories that exhibit non-trivial boson-field dynamics, it is necessary to adopt a local mode picture. Assuming only one of the transverse directions of the trap, e.g., x, supports local phonon modes, the free Hamiltonian of the ion chain can be written as 4 Or a geometrical phase gate in case of spin interactions propor- tions scale similar to the hopping strength and are hence suppressed by at least a factor 10 3 compared with the trap frequency in a typical trap described above. As a result, ! x j in Eq. (14) can be replaced with ! x to a first approximation.
The spin-phonon Hamiltonians in Eq. (8) can be obtained similarly by expanding the ion-laser interaction Hamiltonian in Eq. (2) using the local modes of motion. Then both the spin j and the phonon operators (a x j and a x j † ) are representing quantum operations at location j along the chain, with the associated gate operation identified as On the other hand, in the limit of a tight transverse direction, the spin-spin Hamiltonian in Eq. (12) could be derived using the normal modes of motion along the axial (or the other less tight transverse) direction so to have su ciently strong all-to-all spin-spin couplings. As will be shown, both types of the gates, one based on local modes of motions and the other based on normal modes of motion, can be employed in the simulation, exhibiting another advantage of a hybrid digital-analog setting. In other words, by suppressing the hopping of the phonons along one axis of the trap, the qubits can still be entangled with high fidelity by taking advantage of another set of modes at the same time.
For the simulation of the lattice Schwinger model, the bosons need to interact locally. Such an interaction will need to be mapped to a local phonon-phonon interaction. To engineer this new type of interactions in the trappedion simulator, as proposed in Ref. [], the ions can be put near the minima (or maxima) of a standing wave, which induces an AC Stark shift corresponding to the j m=1 m 2 tor acting on the space of the qubit, with the standard com-(a m ) is the creation (annihilarmal modes of motion with aswith the commutation relations and [a m , a † m 0 ] = m,m 0 . We label motion along the x axis of the · · · , N, those along the y axis of N + 1, N + 2, · · · , 2N , and the ong the z axis of the trap by in-· · · , 3N . With the introduction g laser beams with wave-vector di↵erence ! ⌘ ! L , and phase tuned from an exited internal honon Raman transition among it can be induced with a Rabi ion of ion j. The beams can adhain globally or individually. In t rotates with free Hamiltonian g ion-laser Hamiltonian can be ! 0 | ⌧ ! 0 must be assumed. enote that this Hamiltonian is e. Later on, we need to adopt cture rotating with shifted frertant to bear in mind the origin icke regime in which 1, transitions in the space system take simpler forms, and applications of quantum gates. responds to H ion-laser 0 |O(⌘ 0 ) and L ! 0 = 0. With this setting, nding to ion j becomes 2 respectively. Rotations along the z axis of the Bloch sphere, R Z j (✓ j ), defined as can be implemented without the need for quantum operation using a classical shift [Norbert to verify]. The blue and red sideband transitions correspond to H ion-laser 0 |O(⌘) and are obtained by setting ! L ! 0 = ! k and ! k , respectively. This setting leads to a simple coupled spin-phonon Hamiltonian. In order to achieve a Hamiltonian proportional to j x , both sets of blue and red sideband operations can be applied simultaneously with a set of Raman beams such that the phases of the two sideband beams add up to ⇡. With this setting, the Hamiltonian corresponding to ion j becomes In order to ensure only one of the phonon modes, a k , couples to the qubit in this setup, and to minimize the o↵-resonant transitions, the Raman beams could be set such that they only couple to one set of normal modes of motion. For example, by setting k = kx, the beams only couple to transverse normal modes along the x axis. The qubit-phonon rotations of arbitrary angle ✓ k,j , to be denoted as R a k,j (✓ k,j , k,j ), can be obtained by applying H a k,j ( k,j ) for time ⌧ a gate , to be deduced from the relation In particular, setting laser beam phases to proper angles can generate transitions proportional to a k + a † k and a k a † k : Implementing simultaneous blue and red sideband transitions now detuned from the normal mode frequencies will not generate dynamical phonons as long as the respectively. Rotations along the z axis of the Bloch sphere, R z j (θ j ), defined as can be implemented without the need for a quantum operation using a classical phase shift on the controller for addressing beam j. A useful gate in the circuits constructed in the following sections is the phase gate S j , which can be implemented via S j = R z j (π/4). The blue and red sideband transitions correspond to H ion-laser | O(η) and are obtained by setting ω L j − ω 0 = ω k and −ω k , respectively. This setting leads to a coupled spin-phonon Hamiltonian. In order to achieve a Hamiltonian proportional to e.g., σ y j , beatnotes associated to blue and red sidebands can be applied simultaneously with phases that add up to zero [122]. 4 With this setting, the Hamiltonian corresponding to ion j becomes where k,j denoting the red (blue) sideband laser phase. In order to ensure only one set of the phonon modes, a k , couples to the qubit in this setup, the Raman beams can be set such that the wave-vector difference ∆k is parallel to one of the principal axes of the trap. For example, by setting ∆k = ∆kx, the beams only couple to transverse normal modes along the x axis. If this condition is not met and more than one set of modes are coupled, the different frequencies of the modes can be addressed with lower power to avoid simultaneous coupling to other sets of modes. By applying multiple beatnote frequencies associated with red and blue sideband transitions, corresponding to the set {ω L j } − ω 0 = ±{ω k } for k ∈ {1, · · · , N }, and given the flexibility to set the amplitude and phase associated with 4 To achieve a Hamiltonian of the same type that is instead proportional to σ x j , the red and blue sideband laser phases can be tuned to add up to π. Finally, phonon-ion Hamiltonian proportional to σ z j can be obtained in a circuit by applying appropriate single-qubit rotations, as is demonstrated in later sections.

each beatnote independently, the Hamiltonian
can be implemented directly. Note that Ω j is promoted to Ω k,j to reflect the freedom in the choice of mode-dependent Rabi frequencies. The spin-phonon rotations of an arbitrary set of angles {θ k,j }, to be de- In particular with a single-mode addressing, opportunely setting laser-beam phases gives rise to couplings to different phonon-operator combinations proportional to a k +a † k and a k − a † k : Implementing simultaneous blue and red sideband transitions detuned from the normal mode frequencies will not generate dynamical phonons as long as the lasers are applied for certain time durations. The spin-phonon couplings in this process then effectively induce a spinspin interaction, which leads to the well-known Mølmer-Sørensen (MS) gate [92][93][94]. 5 The corresponding Hamiltonian between ions j and j is Spin-spin rotations R σσ j,j (θ j,j ) of arbitrary angle θ j,j can be obtained by applying H σσ j,j for time τ σσ gate given by the In reality, designing high-fidelity two-qubit entangling gates of this type requires complex pulse shaping techniques to close all the mode trajectories in phase space, which is necessary to eliminate any residual spin-phonon coupling, and to leave the system in the same phonon state before and after the gate operation, see e.g., Refs. [123][124][125][126][127].

B. Phonons as excitations of local modes of motion
In a linear Paul trap, the Coulomb energy is comparable to the potential energy associated with the axial motion of the ions. Namely, the parameter where ω z is the trap frequency along the axial axis of the trap, d 0 is the average spacing among the ions, and e is the electric charge. As a result, the phonons are shared across the chain as excitations of collective normal modes. The same parameter is typically much smaller for the transverse motions of the ions. In particular, if the trap is made extremely tight along the transverse axes of the trap such that β x(y) 1, the transverse phonons describe the oscillation of a single ion, with a relatively small coupling strength for hopping to the neighboring ions that is ∝ e 2 /(d 3 0 M ion ω x(y) ) = β x(y) ω x(y) [120]. For the quantum simulation of QFTs that exhibit non-trivial boson-field dynamics, it is necessary to adopt a local-mode picture. Assuming that the transverse directions of the trap support local phonon modes, the free Hamiltonian of the ion chain can be written as where in the second and third terms, the subscript j on the phonon operators is a local index corresponding to the quanta of motion of the j th ion. For the field theories considered in this work, the simulated theory does not exhibit boson hopping and so the phonon hopping in the simulator needs to be suppressed. In a typical trapped-ion system with tens of Ytterbium ions, ω x ∼ 2π × 5 MHz and d 0 is of the order of a few microns. Then β x ∼ 10 −3 − 10 −4 and β x ω x ∼ 2π × 1 kHz, and hence the dynamics associated with the phonon hopping can be neglected compared with spin-phonon (with typical strength ∼ 100 kHz) and spin-spin dynamics (with typical strength ∼ 10 kHz). If further suppression of the phonon hopping is desired to improve the accuracy of the simulated model, one can either impose tighter trapping potential along the transverse directions or spread out the ions further by reducing the axial potential. Additionally, the nearest-neighbor hopping can be actively blocked using techniques demonstrated in Ref. [107,128] or by applying local optical tweezers to vary the local confinement [129][130][131][132]. Another option, which does not require additional blockade ions, is the use of a mixed-species ion chain with a large mass ratio, in which the modes for the different species separate by mode participation [133]. In an alternating arrangement, phonon hopping will be suppressed by the localmode frequency difference between neighboring ions of different species and the increased distance between the ions of the same species. Finally, we note that the localmode frequencies can be expressed as the common trap frequency plus additional local corrections. These corrections scale similar to the hopping strength and are hence suppressed by at least a factor of 10 −3 compared with the trap frequency in a typical trap described above. As a result, ω x j in Eq. (15) can be replaced with ω x to a first approximation.
The spin-phonon Hamiltonians in Eq. (8) can be obtained similarly by expanding the ion-laser interaction Hamiltonian in Eq. (2) using the local modes of motion. Then both the spin σ j and the phonon operators (a x j and a x j † ) are representing quantum operations at location j along the chain, with the associated gate operation identified as where the single-mode addressing limit of Eq. (10) is considered. On the other hand, in the limit of a tight transverse direction where the hopping of the phonons along transverse axes of the trap is suppressed, the spin-spin gate in Eq. (14) can be still derived with high fidelity using the normal modes of motion along the axial direction so to have sufficiently strong all-to-all spin-spin couplings, given that β z β x(y) . As will be shown, both types of gates, one based on local modes of motions and the other based on normal modes of motion, can be employed in the simulation, exhibiting another advantage of a hybrid analog-digital setting.
For the simulation of the lattice Schwinger model, the bosons need to interact locally. Such an interaction will need to be mapped to a local phonon-phonon interaction. To engineer this new type of interaction in the trappedion simulator, as proposed in Ref. [120], the ions can be placed near the minima (or maxima) of a standing wave, which induces an AC Stark shift corresponding to the Hamiltonian Here, k = kx where k is the wave-vector of the standingwave beam and η ≡ kx (0) = k 2 2Mionω x is the corresponding Lamb-Dicke parameter. Note that the standing wave is only exciting the local modes of motion along one direction based on the choice of the wave-vector. Eq. (17) does not consider higher-order terms in η 2 that can be neglected in the regime η (a x j + a x † j ) 2 1/2 1. It also neglects phonon non-conserving operators that rotate with at least the frequency ω x , and can be adiabatically eliminated as long as F η 2 /ω x 1 [120]. Analogous to the previous gates, phonon-phonon rotations of arbitrary angles can be obtained by applying H s.w. j for time τ aa gate , to be deduced from the relation An important point to notice regarding this gate is that according to Eq. (17), the relative size of the coefficients of the a x j † a x j and (a x j † a x j ) 2 terms, χ (1) /χ (2) , is fixed to (−1 + η 2 )/ η 2 , while the overall strength of these terms can be tuned arbitrarily by changing F via the standingwave intensity, as long as the condition F η 2 /ω x 1 is not violated. The ratio constraint may seem too limiting for simulating an arbitrary Hamiltonian containing both phonon terms. However, as shown in Sec. IV, the desired ratio of the coefficients in the simulation of the lattice Schwinger model can be engineered by appropriately choosing the laser frequencies, which amounts to choosing a suitable interaction picture. Last but not least, if the simulated theory requires site-dependent coefficients for the phonon-phonon couplings, individual standingwave beams could replace the global beam in Eq. (17), at the expense of added experimental complexity.

III. A YUKAWA THEORY: SCALAR FIELDS COUPLED TO (STAGGERED) FERMIONS
A scalar field theory coupled to fermions describes the coupling of the Higgs boson to fermions of the Standard Model through Yukawa interactions, and is responsible for dynamical mass generation and fermion mixing in nature. Non-perturbative studies of the Yukawa theory using lattice-regularized Euclidean field-theory simulations have drawn considerable interest [134][135][136][137][138][139][140][141] as they reveal important connections to the cutoff scale of the Standard Model and questions regarding quantum triviality [142]. In the context of cosmology and the early universe, nonperturbative simulations are required for studies of nonequilibrium and real-time dynamics of this theory and its beyond-the-Standard-Model cousins [143]. A scalar field theory coupled to (non-relativistic) fermions is also the effective field theory description of the Yukawa interactions among the nucleons and pions, and enters the quantum many-body description of nuclei [144][145][146]. It is therefore highly relevant to investigate the prospects for quantum simulation of such theories, including the suitability of the hybrid analog-digital approach. In the following, we focus on a simple case: a scalar field theory coupled to one flavor of staggered fermions in 1+1 D and without the possibility of self-interactions among the scalar fields. Later, we comment on the applicability of this proposal for self-interacting scalar fields and the higher-dimensional case, and will point out the requisite extensions.
A. The Yukawa model Consider a one-dimensional spatial lattice with N sites, lattice spacing b, and with periodic boundary conditions (PBCs) imposed on the fields. 6 The Hamiltonian of the lattice-regularized Yukawa theory to be simulated with a trapped-ion quantum simulator consists of where the purely-fermionic Hamiltonian describes the hopping term and the mass term of one flavor of staggered fermions with mass m ψ . Note that PBCs impose the identification ψ N +1 ≡ ψ 1 . The free scalar-field Hamiltonian is where Π j is the conjugate momentum corresponding to the scalar field ϕ j , i.e. [ϕ j , ϕ j ] = [Π j , Π j ] = 0 and [ϕ j , Π j ] = i(N b) −1 δ j,j . These fields can be quantized in the standard way to obtain a representation in terms of the (bosonic) harmonic-oscillator creation (d † k ) and annihilation (d k ) operators, where k labels the corresponding momentum mode p k ≡ 2πk/(N b), and ε k = ( 2πk Internal states of the ion are used to encode the dynamic of fermions. field. Inputting Eqs. (22) and (23) in Eq. (21), and using the commutation relations of the bosonic operators:

one arrives at the well-known Hamiltonian
describing the energy of N uncoupled quantum harmonic oscillators. Finally, the interacting fermion scalar-field Hamiltonian is where field ϕ must be realized as the collection of quantum harmonic oscillators through Eq. (22).

B. The mapping to an analog-digital circuit
To map the Hamiltonian in Eq. (19) to the building blocks of the analog-digital trapped-ion simulators introduced in Sec. II, one may first transform the staggered fermionic fields straightforwardly to the spin operators through a Jordan-Wigner transformation: . Additionally, the harmonic oscillators can be mapped to the phonons associated with the normal modes of the motion (in either the axial or one of the transverse directions depending on the convenience of the experimental setting). Explicitly, with the labeling rule defined after Eq. (1), the mapping to the transverse normal modes along the x direction reads: d k := a k+N/2+1 and d † k := a † k+N/2+1 for k = −N/2, −N/2+1, , · · · , N/2− 1. To keep the occupation of these modes unaffected while implementing spin-spin entangling gates, the other set of transverse modes or the axial modes can be addressed to implement the MS gates. With the degrees of freedom in the simulated theory being mapped to qubit and phonon degrees of freedom of the trapped-ion simulator, what is left to identify are the gate operations that implement e −iH Yukawa t . Given the digital setting of this proposal, the time-evolution operator can be digitized using the lowest-order Trotter-Suzuki expansion [147]. Higher-order expansions [148,149], as well other state-ofthe-art simulation algorithms [150][151][152][153][154][155][156][157][158], can be similarly investigated within the current approach, but will not be discussed further in this work.
To map the Hamiltonian in Eq. (19) to the Hamiltonians associated with the gates introduced in the previous section, a re-arrangement of the terms is performed such that H Yukawa is now broken down to (29) where the mapping to spin and phonon degrees of freedom has already been carried out. To generate the free scalar-field Hamiltonian in Eq. (24) for generic modedependent coefficients ε k → ε m=k+N/2+1 , one can resort to a change of interaction-picture Hamiltonian [115]. The gates obtained so far are in an interaction picture derived using the free Hamiltonian of the ion system in Eq. (1), containing the harmonic-oscillator energy term with mode-independent and experimentallyfixed coefficients. However, one can introduce the term m ε m (a † m a m + 1 2 ) in the interacting Hamiltonian by appropriately choosing the rotating frame. ε m is the properly rescaled ε m accounting for the ratio of the model time variable and the experimental gate time, see the description after Eq. (35). In particular, setting the free Hamiltonian to H ion − m ε m (a † m a m + 1 2 ), the Hamiltonian in the interaction picture becomes H ion-laser + m ε m (a † m a m + 1 2 ). Operationally, this implies adjusting the detuning of the red and blue sideband transitions, in other words tuning the laser frequency to ω L j −ω 0 = ±(ω k − ε k ), which would implement the desired mode-dependent term ε k (a † k a k + 1 2 ) in the evolution.
However, this is not the full picture yet, since the evolution of the Yukawa theory will be implemented in a digital manner. In particular, there will be multiple sideband operations (when implementing spin-phonon gates) during each step of the Trotter evolution, and these effectively implement the phonon-energy term as well. Therefore, one must ensure that this term will be induced with the correct coefficient. As will be seen below, implementing the time evolution of the Yukawa theory requires introducing an ancilla ion and amounts to a total of N + 1 multi-mode spin-phonon gates per Trotter step. This necessitates changing the free Hamiltonian to H ion − 1 , which leads to the interaction-picture Hamiltonian where L is introduced to denote the possibility of simultaneous application of N beatnote frequencies ω L j ≡ ω 0 ± (ω k − 1 N +1 ε k ) for all k ∈ {1, · · · , N } with their associated amplitude Ω L j ≡ Ω k,j and phase φ L j ≡ φ k,j , which is taken advantage of in the application of spin-phonon gates in the following. It should be emphasized again that experimentally, the change in the interaction picture corresponds to detuning each of the spin-phonon operations. Note that the (N + 1) th normal mode is not affected by the chosen interaction picture, as its dynamics are irrelevant in simulating the scalar fields. With these adjustments, the time evolution of the Yukawa theory for duration δt can now be correctly obtained considering and upon expressing each evolution term as with θ j = 1 2 m ψ (−1) j δt, and with θ m,j = θ m,N +1 = 8N εm δt and φ m,j = 2πj N (m − N 2 − 1), and with all values of m in the range {1, · · · , N }. R σ j (θ j , φ j ) is defined in Eq. (4) and its operation remains the same despite the change in the interactionpicture Hamiltonian. R σa {m},j ({θ m,j }, {φ m,j }) is defined in Eq. (10) but must be realized with red and blue sideband detunings ω L j −ω 0 = ±(ω k − 1 (N +1) ε k ) with ε k τ σa gate = ε k δt, where τ σa gate is determined from Eq. (10) with the θ m,j values specified above. It should be noted that expectation values of observables are invariant under the change in the interaction picture as long as the time-dependent states are transformed accordingly. For simple observables such as fermion and boson occupations, measurements in the original basis obtain the correct expectation values as the corresponding operators commute with the transformation [159,160]. Finally, R σσ j,j (θ j,j ), defined in Eq. (14), can operate using transverse normal

Ancilla ion Free fermion terms
Free scalar field and fermionscalar field interactions modes other than those used in spin-phonon gates, or axial modes. Note that an ancilla qubit, labeled as N +1, is introduced in Eq. (35) to effectively implement the interactions proportional to I j in Eq. (29) [115]. This ancilla qubit is an extra ion prepared in the spin up state and can be used in all subsequent Trotter steps. A schematic of the circuit for a single Trotter step is shown in Fig. 3. We will study the benefits of this hybrid proposal for the simulation of Yukawa theory in Sec. V A.
The algorithm above can be generalized to a Yukawa theory in higher dimensions. In d spatial dimensions, the number of sites on the lattice is N d , where N is the number of sites along each Cartesian direction. The number of ions required to fully encode the dynamics is N d , plus a single ancilla ion as introduced in the 1+1 D case. This is because the number of normal modes of motion associated with each of the transverse (or axial) directions will also grow as N d (considering the ancilla ion, as (N + 1) d ), which is more than sufficient to encode the N d momentum modes of the scalar field in the harmonic oscillator basis. Each ion needs to couple to all N d phonon modes, which polynomially increases the number of spin-phonon gates required to simulate the fermion scalar-field interaction term. The more significant overhead in terms of the entangling operations is caused by the encoding of the fermionic hopping term, which in the Jordan-Wigner transformation involves implementing a chain of Pauli operators, with the number of Pauli operators growing polynomially with N . These operations can be performed through the known decomposition into spin-spin gates, either in series or in parallel. The parallel implementation can be achieved either in one go using a global MS-operation and shelving the ions that do not participate in the coupling using individual beams [161], or by more involved optical pulse-shaping methods as demonstrated in Refs. [162][163][164].
Another important generalization of the Yukawa theory is to incorporate self-interactions of the scalar fields, for example through a quartic interaction Hamiltonian: This term represents non-local interactions among quantum harmonic oscillators in momentum space, requiring all-to-all phonon-phonon couplings to be engineered in the trapped-ion simulator. While inducing phononphonon coupling among the normal modes is possible by taking advantage of the intrinsic non-linearity of the Coulomb interaction [165] or through mediating the interaction via virtual spin degrees of freedom, such an implementation will be more challenging than the other set of gates introduced so far. The use of the local phonon modes will not be optimal either, as the phonon hopping is suppressed beyond nearest-neighbor sites, and its strength cannot be made homogeneous across the sites as required by Eq. (36). Searching for efficient and feasible implementations of the Hamiltonian in Eq. (36) that are more natural to a the trapped-ion simulator will be the subject of future work. The model considered in this work is still of phenomenological interest, for example in simulating lattice effective field theory of nucleons coupled to pions. The self-interactions of pions matter only at higher orders in the chiral effective field theory of nuclear forces [145,146] and can be neglected at low energies. Besides the trivial incorporation of multiple flavors of (non-relativistic) fermions representing spin and isospin components of a nucleon, and the (local) self-interaction of fermions representing nucleons two-and three-body contact interactions, the fermion scalar-field interacting term in Eq. (25) must be promoted to a derivative coupling in order to represent a pion-nucleon coupling. Such a term can be implemented by representing the derivative of the scalar field by a finite difference among the fields at adjacent sites. This then amounts to multiple implementations of the terms of the type in Eq. (25) that are added with appropriate signs. Further, there are three different types of pions, each requiring their own harmonic-oscillator representation (two of which being electrically charged). These can either be realized through the three sets of normal modes or by enlarging the ion chain to create more normal modes of one type for encoding of all the pionic degrees of freedom. As the building blocks of this construct are identical to the ones for the simple Yukawa model above, the explicit circuit will not be discussed further.

C. An example with realistic parameters
To demonstrate the viability of near-term experimental implementations of the dynamics in this model with a trapped-ion simulator, one can investigate the range of the gate parameters required to observe interesting phenomena in this model. While only small instances of the problem can be studied classically, as is evident from the examples considered, the quantum simulator with tens of ions can still operate in the same gate-parameter range and already push the limits of classical capabilities, particularly given the native implementation of boson dynamics in the simulator. Later on in Sec. V, we demonstrate the advantage of an analog-digital approach by making a qualitative cost comparison with a fully-digital implementation of the same model.
An interesting phenomenon in the Yukawa theory considered is the dynamical generation of mass even if the bare fermion mass is set to zero originally. Such an effective mass controls the rate of fermion generation throughout the evolution and depends on the boson accumulation on the lattice. As Fig. 4 demonstrates for small lattice sizes N and small boson occupation cutoffs Λ, such a non-trivial fermion, anti-fermion, and boson occupation evolution can be revealed in the Loschmidt echo, i.e., | ψ(0)|ψ(t) | 2 , where ψ(0) is the state with no fermion, no antifermion, and no boson. After a quench, this state evolves to a superposition of states with any number of  fermions and antifermions allowed by symmetries, and involving varying boson occupation whose average is shown in the inset of the plots. The unfilled points represent the same quantities evaluated using a Trotter expansion of the evolution operator according to the split in the Hamiltonian terms outlined in the previous section. The coarsest Trotter digitization is chosen such that quantities can still be described accurately in the time window specified. The corresponding model parameters, shown in Table I, give rise to gate rotation angles denoted in Table IV of Appendix A 1 along with associated experimental gate parameters for the R σa gate, including gate's operation time, in Table V. The gate time required is shorter (a few tens of microseconds for reasonable experimental parameters) than that of the spin-spin entangling gate (a few hundreds of microseconds), which is a great improvement over fully-digital implementation that would have required multiple entangling spin-spin gates, being an order of magnitude slower in general. We will come back to gate-number scaling of the analog-digital implementation compared with the fully-digital implementation in Sec. V. The spin-spin entangling gates required in simulating the solely fermionic dynamics can be driven with the rotation angles specified in Table IV using well-known laser pulse-shaping techniques. The trap and laser characteristics, including transverse normal modes of motion and the Lamb-Dicke parameter, are provided in Table VI. A final technical detail regarding the experimental implementation is that all normal-mode vectors must have non-zero components so that all Rabi frequencies Ω j are finite when applying the R σa gate. To ensure that each ion in the simulation couples to all of the phonon modes used, the total number of ions must be even. Therefore, studying an N -site theory of the 1+1 D Yukawa model, where the number of staggered sites N is even, generally requires N + 2 ions: one ancilla, which is operated on by the specified gates in Fig. 3, and one idle ion that is only present to ensure that for any given mode, no ion is stationary along the chain.

IV. U(1) LATTICE GAUGE THEORY COUPLED TO (STAGGERED) FERMIONS IN 1+1 D
The Schwinger model, the theory of quantum electrodynamics in 1+1 D, has long served as a testbed for simulation methods (both classical and quantum) for stronglyinteracting gauge theories. It is qualitatively similar to quantum chromodynamics, the theory of the strong force in nature, as it exhibits confinement, chiral-symmetry breaking, and a nontrivial θ vacuum [166,167]. Being a low-dimensional and Abelian theory, it is simpler than its higher-dimensional and non-Abelian counterparts and its lattice-regularized form has been the subject of valuable quantum-hardware implementation benchmarks in recent years [54-57, 65, 66]. Fully-analog proposals for Hamiltonian simulation of the lattice Schwinger model  Table I. The time t is in units of lattice spacing b, which is set to one.
have been put forward using spin degrees of freedom only [32], and both the spin and phonon degrees of freedom [112], but they remain relatively challenging for hardware implementation. Here, we propose a hybrid analog-digital simulation of the lattice Schwinger model that brings the simulation proposals involving the phonon degrees of freedom a step closer to experimental realization.

A. The Schwinger model
The well-known Hamiltonian of the lattice Schwinger model, introduced by Kogut and Susskind [166,167], can be written in terms of a one-component staggered fermion ψ j defined on site j (with odd and even sites corresponding to matter and anti-matter fields, respectively), the gauge link U j and its conjugate field, the electric field E j , both defined on the link emanating from site j, and satisfying the commutation relations [E j , E j ] = [U j , U j ] = 0 and [E j , U j ] = U j δ j,j . The Hamiltonian can be written as where the gauge-matter interacting Hamiltonian, the staggered fermion-mass Hamiltonian, and the electricfield Hamiltonian are respectively. The Hilbert space of the theory is characterized by on-site quantum numbers n (g) j and n (f ) j associated with the discrete spectrum of a quantum rotor satisfying E j |n The physical states are those annihilated by the Gauss's law In the following, it is assumed that the simulation starts in a physical state and hence will remain in the same sector, as long as the gauge-symmetry violations arising from the digitization of the time-evolution operator, or from hardware imperfections, remain small.
Unfortunately, a quantum rotor representing the gauge link in the Schwinger model does not have the same algebra as the phonon modes of the trapped-ion simulator for a direct encoding. A highly-occupied boson model (HOBM) was proposed in Ref. [112] to resolve this mismatch by considering the following mappings: for an integer M . This transformation keeps the commutation relation between E j and U j intact, but modifies that between the gauge links to [U j , U j ] = δ j,j /M . Only in the limit M N , with N being the number of sites on the lattice, the HOBM will recover the lattice Schwinger model, necessitating the simulation involving the bosonic degrees of freedom to initiate in a state with a large number of bosons. Nonetheless, the numerical simulations of Ref. [112] demonstrate that for M ∼ N , a high level of accuracy is still achieved, particularly in space-and time-averaged dynamical observables. Applying the HOBM mapping as well as the Jordan-Wigner transformation, the spin-boson Hamiltonian of the lattice Schwinger model can be written as   where H (I)

B. The mapping to an analog-digital circuit
This Hamiltonian can be directly mapped to spin and phonon degrees of freedom in the trapped-ion simulator. Noting that the electric-field Hamiltonian in the HOBM, H (III) U (1) , consists of local self-interactions of the bosonic modes, it is necessary to consider a tight trap in the transverse directions, as introduced in Sec. II B, to take advantage of the local modes of the motion. In particular, the boson encoding should be realized by the mapping d j → a x j , and the local self-interaction of the phonons can be implemented by the standing-wave gate introduced in Eq. (18) The Trotterized evolution operator is where each square bracket corresponds to the time evolution of each of the four terms in Eq. (43) and θ j =

Fermion-gauge interactions
Fermion mass term Gauge-field interactions · · · · · · FIG. 6. The schematic of the analog-digital quantum circuit associated with the time evolution of a four-site Schwinger model for a single Trotter step, as expressed in Eqs. (47)- (49). The fermion-gauge interaction block must be repeated four times with various phases for the spin-phonon gates and single-spin rotations, see Eq. (47). The gate symbols are defined in Fig. 1.
with θ j = 1 2 (−1) j mδt, and with the gates defined in Secs. II A, and with χ (1) = −g 2 bM δt and χ (2) = 1 2 g 2 bδt, and with the gates defined in Secs. II B. The schematic of the corresponding circuit is shown in Fig. 6. Note that the MS gate R σσ j,j+1 in Eq. (47) can be implemented using the normal modes of the motion along the axial direction, as the phonons involved in its implementation are virtual and not contributing to the dynamics of gauge fields in the simulated theory. Furthermore, the ratio χ (1) /χ (2) in Eq. (49) is fixed in the gate design as mentioned before, and one must ensure that χ (1) /χ (2) = (−1 + η 2 )/ η 2 = −2M , with the standing-wave Lamb-Dicke parameter η defined in Sec. II B. As typically η 2 1, the requirement M N is guaranteed, but the careful tuning of η in experiment is required as the thermodynamic limit of the HOBM (N → ∞) is taken. If the necessary variation in η is not feasible in a given experimental setup, one can adjust the detuning of the laser beatnotes, as explained in Sec. III, to effectively shift the free Hamiltonian of the trapped-ion system by a term proportional to j a x j † a x j in a suitable rotating frame. In this way, the corresponding term in the interaction-picture Hamiltonian, i.e., the first term in brackets in Eq. (46), can take an arbitrary coefficient, hence relaxing the condition on the ratio of a x j † a x j and (a x j † a x j ) 2 terms. This will amount to changing the frequency associated with red and blue sideband transitions involved in the operation of R σa j (θ j , φ j ) accordingly, as detailed before for the case of the Yukawa theory.
Simulating QED in higher dimensions requires engineering the magnetic Hamiltonian, that is the sum of the product of gauge links along the elementary plaquette in each two-dimensional plane in coordinate system. In the HOBM of QED, the gauge links on these higher-dimensional lattices can still be mapped to the local modes of motion in a linear chain of ions. However, engineering the magnetic Hamiltonian requires achieving four-body interactions among the phonons pinned to four adjacent ions, which will be challenging given the experimental setup explained in this work. This is because the nearest-neighbor phonon hopping needs to be simultaneously suppressed to simulate the electric-field Hamiltonian accurately. Analog and digital simulations with quantum simulators are proposed for QED in 2+1 D with pure gauge interactions, using both the HOBM or other models [49,[168][169][170]. Future work will investigate the design of an analog-digital simulation approach for this problem.
C. An example with realistic parameters Similar to the case of the Yukawa theory, near-term experimental implementations of the Schwinger model with realistic analog-digital gate parameters will be viable for small to intermediate systems and can be scaled up straightforwardly similar to what is expected for digital implementations. It is therefore useful to consider small instances of the problem that can be simulated classically to obtain the range of expected parameters in experiment for which interesting phenomena in the model can be observed. Among many interesting features of the Schwinger model are the fermion-antifermion pair creation and the string-breaking dynamics. These quantities can be evaluated from expectation values of fermionic or bosonic observables in a state evolved from a trivial initial state, |ψ(0) , such as the vacuum of the model in the limit g → ∞, i.e., a state with no fermion, no antifermion, and M bosons. Denoting the time-evolved state at time t as |ψ(t) , we simply consider the Loschmidt echo, i.e., | ψ(0)|ψ(t) | 2 . |ψ(t) involves a superposition of states with any number of fermions and antifermions allowed by symmetries, and involving varying boson occupation whose average is shown in the inset of the plots in Fig. 7. The coarsest Trotter digitization is chosen such that quantities can still be described accurately in the time window specified, and the corresponding values are marked with unfilled points. Λ denotes the cutoff on the excitations of the bosons above the M bosons in the initial state. It is notable that the boson cutoff effects are significant even for such a small lattice, motivating further the need for mapping a larger bosonic Hilbert space to the quantum simulator, which is done naturally in the analog-digital proposal of this work using the local phonon modes in the trap.
Since the trap is in the tight-binding regime along the transverse directions, the local-mode frequencies at the location of each ion are equal to a good approximation, simplifying the parameters of the gates applied to each ion. For the phonon-phonon gate, it is ensured that the parameter F η 2 /ω x is small by choosing η and F appro-   Table II. The time t is in units of lattice spacing b, which is set to one.
priately. Furthermore, to generate a specific value of the ratio of the phonon terms a x † a x and (a x † a x ) 2 as required by the HOBM, the red and blue sideband transitions in the spin-spin and spin-phonon gates must be implemented with a frequency shift −δω x , corresponding to the necessary modification to the interaction-picture Hamiltonian. The shift value, along with other relevant trap parameters, are tabulated in Table X of Appendix A 2. As noted in Tables VIII and IX, the spin-phonon gate time is similar to what was obtained for the Yukawa theory. The phonon-phonon gate time needs to be ∼ 100 microseconds to prevent the need for a large value of the standing wave amplitude. Last but not least, the Lamb-Dicke parameters associated with the Raman beams as well as the standing-wave beam must satisfy η √ M 1 and η √ M 1, respectively. This is because the system must be prepared in a state with phonon occupation M according to the HOBM. These conditions are satisfied with the realistic experimental parameters chosen.

V. A QUALITATIVE COST ANALYSIS
In this section, our hybrid analog-digital approach will be compared with the fully-digital realization of the Yukawa theory and the Schwinger model. In the long term, when fault-tolerant digital quantum computation will be a reality, the non-Clifford gate (T-gate) count determines the cost of the simulation. In the near-term and particularly over the next decade, however, quantum computations will be limited by the slow and error-prone application of entangling (CNOT) gates. In the analogdigital protocol of this work, while the hardware-specific spin-phonon gate is an entangling gate in the combined Hilbert space of the phonons and qubits, its operation is governed by dynamics that occur at O(η), where η is the Lamb-Dicke parameter introduced in Sec. II. This is in contrast to the entangling spin-spin operation that is the core of the CNOT implementation in a trapped-ion digital computer, which is governed by dynamics that occur at O(η 2 ). Since η 1 in the Lamb-Dicke regime in which experiments operate, the spin-spin gates are one or two orders of magnitude slower than single-spin and spin-phonon gates. Phonon-phonon self-interaction gates based on standing-wave beams are governed by O( η 4 ) dynamics but their strength can be compensated by laser power, as demonstrated by the examples studied before. Without experimental realizations, the fidelity of this new set of gates cannot be accurately predicted. Nonetheless, the qualitative expectations for the gate times as described here can provide a rough guidance to the relative performance of the simulator in an analog-digital mode compared with a fully-digital mode.

A. The Yukawa theory
Since the implementation of the dynamics associated with the fermionic hopping and mass is the same in both digital and analog-digital circuits, here we focus on implementing dynamics involving the scalar field, namely evolution with the free scalar-field Hamiltonian as well as the fermion scalar-field interacting Hamiltonian.
In the analog-digital protocol, the time evolution of the free scalar-field Hamiltonian H (II) Yukawa in Eq. (24) is implemented at no cost since it amounts to adjusting laser beatnote frequencies in the sideband operations, as explained in Sec. III B. In the fully-digital implementation, the exact circuit design depends on the mapping of the scalar fields to qubit registers. The field basis [13], the harmonic-oscillator basis [29], and the single-particle basis [44] are among the representations considered in literature to perform such a mapping, each with their own benefits and disadvantages. Here, we give a qualitative CNOT-gate count (or in turn MS entangling-gate count) in the harmonic-oscillator basis, as it is the representation closely related to what is considered in Sec. III. In this basis, the occupation number of each harmonic oscillator corresponding to each mo-mentum mode is truncated at some cutoff value Λ and is expressed in a binary representation, which is then mapped to log Λ qubit registers. The evolution operator e −iH (II) 2 )δt can then be realized by a number of single-qubit rotations around the z axis of the qubits' Bloch sphere (see e.g., Ref. [39] for a similar example of implementing the time evolution of the electric-field Hamiltonian in the lattice Schwinger model). This step, therefore, can be considered cost-free in both digital and analog-digital implementations.
Next, the time evolution with the fermion scalar-field Hamiltonian H (III) Yukawa in Eq. (25) is implemented with a single ancilla qubit and with O(N (N + 1)) spin-phonon gates, which are considered to be cost-free compared with the entangling spin-spin gates. In a fully-digital implementation, e −iH (III) Yukawa δt can be implemented by writing H The implementation of these terms follows the procedure explained in Ref. [39] in the case of realizing the dynamics of fermion-gauge interaction term in the Schwinger model. First, Eq. (50) indicates that the operations proportional to A ≡ (d k + d † k ) and B ≡ i(d k − d † k ) operators are only performed if the fermion occupies a given site, necessitating a controlled operation on the qubit register of the fermion. The A and B operators are two near-diagonal matrices whose exponential can be implemented using the shift operators, that are realized using quantum Fourier transform circuits and single-qubit rotations in the Fourier space. As a reminder, for a binary number with log Λ digits, each quantum Fourier transform requires O((log Λ) 2 +log Λ) CNOT operations. Relating A and B operators to the shift operators requires a periodic wrapping of the matrices, i.e., identifying the least and most values of harmonic-oscillator occupation. This unphysical modification can subsequently be removed by application of appropriate log Λcontrolled operations, amounting to O(log Λ) additional CNOT gates [39]. Putting everything together, including the controlled operations required on the fermionic register, and taking into account all the terms in the lattice sums in Eq. (50), the time evolution operator e −iH (III) Yukawa δt can be implemented using O(N 2 (log Λ) 2 ) CNOT operations. Therefore, assuming that spin-phonon gates of the hybrid scheme are free compared with spin-spin entangling gates, the digital approach is inferior to the hybrid approach. Even if the spin-phonon gate performs comparably to the spin-spin (CNOT) gate, the digital scheme require O((log Λ) 2 ) more entangling operations which can be significant when Λ 1.

Yukawa theory
Fermion hopping Fermion mass Free scalar fields Fermion scalar-field interaction

Schwinger model
Fermion-gauge interaction Fermion mass Electric-field term a O(N ) standing-wave phonon-phonon gates which are expected to have comparable fidelity to the spin-spin gates. Such an advantage is at the core of the power of the hybrid approach: phonons are represented naturally and as many phonon excitations as permitted in the dynamics can be generated without the need to cut their spectrum off. Of course, an excessive number of phonons in the system can lead to Kerr cross-coupling [165] and loss of coherence in the simulator, and therefore a balance should be established between accuracy of the simulated theory given a truncated boson spectrum and the experimental error in the simulator. For this reason, experimental benchmarks are necessary in confirming these qualitative theoretical expectations. A summary of the entangling-gate count of both schemes for evolving each term in the Hamiltonian of the Yukawa theory is provided in Table III. B. The Schwinger model Except for the time evolution of the fermion mass term, both the fermion-gauge field interaction and the electricfield term (the boson self-interactions in the HOBM) are implemented differently in the hybrid and fully-digital schemes.
The circuit in Fig. 6 reveals that the time evolution of the interacting fermion-boson field in the HOBM requires O(N ) spin-spin gates and O(N ) spin-phonon gates, with the latter expected to be not too costly. In the fully-digital scheme, e −iH (I) U (1) δt with the Hamiltonian in Eq. (43) can be implemented following the circuit construction described earlier in the case of fermion scalarfield interactions of the Yukawa theory. The only differences are that now the bosons are defined locally, and associated with each site (link) there is one such boson (as opposed to N bosons associate with all momentum modes in the Yukawa theory), and that the fermions correspond to nearest-neighbor sites (as opposed to a local fermion occupation operator in the corresponding term in the Yukawa theory). As a result, the total entangling-gate count of the digital-circuit implementation of this evolution operator is O(N (log Λ) 2 ), where Λ denotes the cutoff on the boson excitations (translating to the electric-field excitations in the original U(1) theory). Therefore, the hybrid implementation will save a factor of O((log Λ) 2 ) in CNOT count compared with the digital implementation assuming again that the spin-phonon gate performs at a much higher fidelity than the spin-spin gate.
Time evolving the electric-field term, that is Eq. (46), requires N standing-wave gates, which can be counted to be comparable to O(N ) entangling gates. The same term in the fully-digital computation requires either free Zrotations (for evolving the d † j d j term similar to what was discussed for the free scalar Hamiltonian in the Yukawa theory) or O(N log Λ(log Λ − 1)) CNOT gates (for evolving the (d † j d j ) 2 term). This latter count can be understood by noting that the application of d † j d j phase depends on the occupation of the d excitation itself. In summary, the analog-digital implementation requires O((log Λ) 2 ) fewer entangling operations (assuming that the standing-wave gate can perform comparably to the spin-spin gate). A summary of the entangling-gate count of both schemes for evolving each term in the Hamiltonian of the HOBM is provided in Table III.

VI. CONCLUSION AND OUTLOOK
Our work combines the benefits of analog and digital trapped-ion simulators to enable both short-and long-term simulations of quantum field theories involving bosonic fields with sizable Hilbert spaces. Among the models considered are the Yukawa theory, i.e., scalar fields coupled to fermions in 1+1 D and the lattice Schwinger model within the highly-occupied bosonic model. By introducing a set of known and new gates, including spin, spin-spin, spin-phonon, and phonon-phonon gates, the Trotterized evolution in each of these models is expressed as a quantum circuit. Small instances of the models with realistic experimental parameters for nearterm proof-of-principle demonstrations of non-trivial dynamics are identified. Experiments with 20 ions using a set of normal or local modes, each occupied with only up to two phonon excitations, already push the limits of what is possible on a classical computer. Most importantly, higher cutoffs on the bosonic excitations can be imposed in the model as the bosons are naturally represented by the phonons in the ion trap, leading to an O (N log Λ) (O (log Λ)) reduction in the number of qubits and at least an O (N 2 log Λ) 2 (O (log Λ) 2 ) reduction in the number of the entangling spin-spin gates compared with a fully-digital simulation in the Yukawa theory (Schwinger model). Here, Λ is the cutoff on the boson-field excitation and N is the number of lattice sites in the simulated theory. These scalings assume that the spin-phonon gate has a higher fidelity compared with the spin-spin gate, while the phonon-phonon gate performs with a fidelity comparable to that of a spin-spin gate. Experimental implementations will test these theoretical expectations in the coming years and establish the analog-digital mode of the simulator as a more powerful paradigm than either of fully-analog and digital schemes in enabling simulations of interacting fermionic-bosonic field theories of phenomenological interest.
A number of promising theoretical and experimental directions may result from extensions of this work to other quantum simulators and more complex stronglyinteracting field theories. These include the following: The bosonic modes are common in a range of quantum simulators. A notable example is a circuit-QED platform in which cavity photons are coupled to superconducting qubits [171,172]. The physics of photon-qubit effective interactions is very similar to that of phonon-ion interactions in an ion trap. For example, cavity photons can be taken advantage of to induce non-local couplings among the qubits [173]. One can envision encoding the dynamics of scalar or gauge fields in the photon modes and that of fermions in the superconducting qubits, hence devising photon-based gates as discussed in this work. Nonetheless, the degree of control and the feasibility of experimental realizations must be examined carefully before the potential of an analog-digital scheme can be evaluated in such platforms.
A highly desirable capability, that was not required for the current work given the nature of the models studied, is the engineering of non-local phononphonon couplings. Such a capability would allow the simulations of a self-interacting scalar field theory, as well as the plaquette (magnetic) interactions in lattice gauge theories in higher dimensions, in an analog-digital fashion. It will be worth exploring the possibility of taking advantage of the virtual spins to mediate the phonon interactions in this context, inspired by the Mølmer-Sørensen scheme, in which spin-spin interactions are mediated by the virtual phonons. Promising experimental demonstrations in circuit-QED have emerged [174] and similar ideas may be applicable to the trapped-ion systems.
To make progress toward the goal of simulating quantum chromodynamics, one must develop feasible simulation proposals for non-Abelian gauge theories. The SU (2) LGT coupled to fermions is a valuable non-Abelian gauge theory that has become the focus of intense studies in recent years [18,19,21,22,52,58,59,68,73,74]. Its prepotential formulation in terms of Schwinger bosons [71] appears the most suitable representation to potentially take advantage of the analog-digital setup presented here. Basically, despite the U(1) theory in 1+1 D which has only an approximate representation in terms of bosonic degrees of freedom, the SU (2) LGT in any dimension has an exact mapping to bosons. The complication is that each link on the lattice is a collection of four sets of oscillators and the physical degrees of freedom are those made up of gauge-invariant combinations of these bosonic operators (along with the fermionic operators). The Hamiltonian matrix elements retain the knowledge of the SU(2) algebra and hence contain non-trivial factors [71]. Ongoing progress in developing fully-digital algorithms for this theory, nonetheless, will be beneficial in developing the analog-digital counterparts. Here, complications associated with the non-local phonon-phonon interactions in the trapped-ion simulator must be dealt with as well.
It will be interesting to consider efficient statepreparation strategies for QFTs when bosonic fields can be encoded into bosonic registers. Besides the reduction in the number of qubits and entangling gates required for the preparation routine in an analog-digital mode, one can envision more efficient state-preparation protocols as well, such that nontrivial initial states involving entangled boson and fermion modes can be achieved more straightforwardly. This direction is, however, less developed and requires additional investigation.
Last but not least, it is important to investigate whether the experimental imperfections and non-ideal gate fidelities are more significant in an analog-digital compared to a fully-digital setting. This is a rather crucial question when it comes to gauge-theory simulations, as local gauge invariance imposed on an initial state needs to be efficiently retained throughout the evolution, or the simulated physics will be fundamentally different than the desired gauge theory. Recent ideas in protecting gauge invariance in the simulation [175][176][177][178] can be investigated in the context of analog-digital simulations as well.

ACKNOWLEDGMENTS
We acknowledge valuable discussions with Mohammad Hafezi and Alireza Seif. ZD is supported in part by the U.S. Department of Energy's Office of Science Early Career Award, under award no. DE-SC0020271, for theoretical developments for map-ping QFTs to quantum simulators, and by the DOE Office of Science, Office of Advanced Scientific Computing Research (ASCR) Quantum Computing Application Teams program, under fieldwork proposal number ERKJ347, for algorithmic developments for scientific applications of near-term quantum hardware. GP and NML acknowledge support by the DOE Office of Science, Office of Nuclear Physics, under Award no. DE-SC0021143, for designing hardware-specific simulation protocols for applications in nuclear physics. GP and NML are further supported by the Army Research Office (W911NF21P0003) and the Office of Naval Research (N00014-20-1-2695).  (35). The values are associated with the model parameters given in Table I Tables VIII and IX.