{'fraction_non_alphanumeric': 0.04128850824491883, 'fraction_numerical': 0.010865396906557587, 'mean_word_length': 4.3687521445728015, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Digital quantum simulation of scalar Yukawa coupling: Dynamics following an interaction quench on IBM Q

Motivated by the dearth of studies pertaining to the digital quantum simulation of coupled fermion-boson systems and the revitalized interest in simulating models from medium- and high-energy physics, we investigate the nonequilibrium dynamics following a Yukawa-interaction quench on IBM Q. After adopting -- due to current quantum-hardware limitations -- a single-site (zero-dimensional) version of the scalar Yukawa-coupling model as our point of departure, we design low-depth quantum circuits that emulate its dynamics with up to three bosons. In particular, using advanced circuit-optimization techniques, in the one-boson case we demonstrate circuit compression, i.e. design a shallow (constant-depth) circuit that contains only two CNOT gates, regardless of the total simulation time. In the three-boson case -- where such a compression is not possible -- we design a circuit in which one Trotter step entails 8 CNOTs, this number being far below the maximal CNOT-cost of a generic three-qubit gate. Using an analogy with the travelling salesman problem, we also provide a CNOT-cost estimate for quantum circuits emulating the system dynamics for higher boson-number truncations. Finally, based on the proposed circuits for one- and three-boson cases, we quantify the system dynamics for several different initial states by evaluating the expected fermion- and boson numbers at an arbitrary time after the quench. We validate our results by finding their good agreement with the exact ones obtained through classical benchmarking.


I. INTRODUCTION
The reinvigorated research interest in digital quantum simulation (DQS) [1][2][3] has chiefly been motivated by the tantalizing recent progress in the development and deployment of quantum hardware based on superconducting- [4], trapped-ion- [5], or neutral-atom systems [6]. This recent flurry of research activity in DQS was preceded by a large body of investigation pertaining to analog quantum simulators [3] of many-body systems that had been accrued over the period of nearly two decades. Both of these research strands have in large part been inspired by the pioneering assertion of Feynman that a generic quantum system could efficiently be simulated using a device whose operation is also governed by quantum-mechanical laws [7]. Admittedly, the field of DQS has heretofore been dominated by the development and application of quantum algorithms for simulating purely fermionic systems [8][9][10][11][12][13][14][15][16][17][18][19]; this is mostly related to the fact that -owing to the existence of the Pauli exclusion principle -simulations of such systems typically require relatively modest quantum-hardware resources.
Motivated by the dearth of research activity pertaining to DQS of coupled fermion-boson models, as well as the emerging trend of simulating models from highand medium-energy physics [34][35][36][37][38][39][40], in this paper we present a DQS of the nonequilibrium dynamics following a quench [41] of scalar Yukawa interaction [42]. Such coupling -first proposed in the context of interaction between nucleons mediated by pions [43] -involves fermion-, antifermion-, and boson degrees of freedom. Due to technological limitations of current noisy intermediatescale quantum (NISQ) hardware [44], here we investigate the low-energy limit of scalar Yukawa coupling on a single lattice site, with low truncation numbers for the real scalar (boson) field. We show that, despite its inherent simplicity, the resulting single-site (zero-dimensional) model shows nontrivial quantum dynamics.
By first making use of the local charge conservation in this system, which allows us to encode the fermionantifermion sector of the problem using a single qubit, we subsequently design low-depth quantum circuits that emulate the system dynamics with up to three bosons. Furthermore, we implement those circuits as a sequence of quantum gates on state-of-the-art IBM Quantum (IBM Q) hardware, based on quantum processors that rely on the use of controlled-NOT (CNOT) as a representative two-qubit gate. Based on the designed circuits, we characterize the system dynamics for several different initial states by computing the expected fermion-and boson numbers at an arbitrary time after the quench on IBM Q. We also provide a benchmark for the results obtained on IBM Q by carrying out a numerically-exact evaluation of the aforementioned quantities on a classical computer. We validate our DQS by finding a reasonably good agreement between the two sets of results.
Using sophisticated circuit-optimization techniques, in the one-boson case -which in our realization corresponds to a two-qubit system -we demonstrate an example of circuit compression. More precisely, we design a constant-depth circuit that contains only two CNOT gates, regardless of the total simulation time. In other words, in this case we find a much more efficient realization than the conventional one in which the circuit depth scales linearly with the number of Trotter steps. Such circuit compressions are only possible for certain special types of the system (coupled-qubit) Hamiltonian and have quite recently been discussed in the context of the transverse-field Ising-and XY models by means of Liealgebraic methods [45], as well as using reflection symmetry and the quantum Yang-Baxter equation [46]. The shallow circuits resulting from such compressions have depths that are completely independent of the Trotterstep size and depend only on the number of qubits.
On the other hand, in the more nontrivial three-boson case -which in our adopted state-encoding scheme corresponds to a three-qubit system -a compression to a constant-depth circuit is no longer possible. We show that in this case it proves beneficial to make use of the second-order Trotter-Suzuki product formula [47,48], which entails a symmetrized Trotter step. Using this symmetrized Trotter decomposition, we design a circuit in which one Trotter step entails 8 CNOTs. This last number of CNOTs is far below the maximal CNOT-cost (14) of a generic three-qubit gate [49], which speaks in favor of the efficiency of the proposed circuit [49,50].
We also address cases with larger boson-number truncations by providing the CNOT-cost estimate for quantum circuits that emulate the system dynamics in those cases. We do so based on an analogy of our circuit optimization with the travelling salesman problem (TSP) [51]. To this end, we make use of the exact solution of the TSP based on the Bellman-Held-Karp dynamic-programming algorithm [52,53] -which due to its unfavorable time-complexity is limited to rather low truncations -as well as the Christofides algorithm, a polynomial-time heuristic that approximately solves the TSP on a metric graph [54]. Being polynomial in character, the latter can be utilized for much larger bosonnumber truncations.
The remainder of this paper is organized as follows. In Sec. II we introduce the fermion-boson system under consideration and its underlying Hamiltonian. Section III is devoted to a short description of our scheme for describing the dynamics following an interaction quench via DQS, followed by a brief recapitulation of our adopted approaches for encoding the relevant Fock states of fermions and bosons on dedicated registers of qubits. Section IV is devoted to the design of low-depth quantum circuits that emulate the dynamics of the system un-der consideration, both in the two-qubit-and three-qubit cases. We complement this by providing a CNOT-cost estimate for higher boson-number truncations in Sec. V. The principal results for the relevant quantities describing the system dynamics, along with their comparison with the exact ones obtained through classical benchmarking, are presented in Sec. VI. We summarize the paper and underscore its main conclusions in Sec. VII. Some cumbersome derivations, as well as certain relevant mathematical and computational details, are relegated to Appendices A, B, C, and D.

II. SYSTEM AND HAMILTONIAN
To set the stage for further considerations, we first introduce the system at hand and its underlying fermionboson Hamiltonian. We start with some general remarks on Yukawa-type interactions (Sec. II A) then introduce the scalar Yukawa-coupling model (Sec. II B), and, finally, specialize to the single-site (zero-dimensional) version of the latter model (Sec. II C).

A. Yukawa-type interactions: Generalities
Being one of the basic interaction mechanisms involving both fermionic and bosonic particles, Yukawa-type interactions represent a paradigm for particle production and absorption. The simplest form of Yukawa-type interaction is described by a coupling term of the type where ψ(x, t) is a fermion field, φ(x, t) a boson field, while g is the Yukawa (fermion-boson) coupling strength. This interaction term is a scalar in the fermion space and involves a scalar boson field. For the sake of completeness, it is worthwhile mentioning that there exist much more complicated Yukawa-type interaction mechanisms [55]. For example, ψ † σψ·∇φ is an interaction term that couples a pseudovector in the fermion spinor space to a pseudoscalar boson, this term being also of scalar character. Likewise, a scalar in the configuration space that involves a boson of isovector character through internal degrees of freedom is given by a term of the type ψ † τ ψ · φ, in which a fermion isospinor combination is an isovector and φ ≡ {φ 1 , φ 2 , φ 3 } is a three-component field in the isospin space; this last interaction is an isoscalar that guarantees charge independence of the system. The Yukawa-interaction mechanism was originally introduced to describe the nuclear force between nucleons mediated by pions [43]. It later attracted interest in the context of the Standard Model of partice physics, where it describes the coupling between the Higgs field and massless quark and lepton fields (i.e. the fundamental fermion fields); the latter fermion fields acquire masses, via the Higgs mechanism, after electroweak symmetry breaking [42]. Finally, it is worthwhile mentioning that interaction mechanism somewhat analogous to Yukawa coupling are of relevance in the context of strongly-correlated systems in condensed-matter physics [56].

B. Scalar Yukawa coupling
In the following, we will be interested in the Yukawatype interaction of a scalar field φ and a Dirac field ψ, described by the Hamiltonian where g is the dimensionless coupling strength. While pions, which mediate the nucleon-nucleon interaction, are pseudoscalar mesons, we will assume a scalar interaction for the sake of simplicity. In our subsequent discussion, we adopt natural units, where = c = 1.
To perform a DQS of the real-time dynamics of scalar Yukawa coupling, we have to discretize the theory and select a convenient basis. While approaches that are directly based on the algebraic properties of the field operators do exist, in what follows we make use of an expansion in terms of creation and annihilation operators. We further select the momentum eigenstates |p as our preferred basis since it diagonalizes the free Hamiltonian. Using the definition of the Fourier transform φ(x) = (2π) −3 d 3 p e ip·x φ(p), one may rewrite the interaction Hamiltonian as with the field operators in momentum space being defined as with the relativistic dispersion relations ω 2 p = m 2 + p 2 (for bosons) and Ω 2 p = M 2 + p 2 (for fermions) [42]. By discretizing the integrals in Eq. (3) and using Eq. (4), one obtains the momentum-space version of H int in terms of creation/annihilation operators where a l is the lattice spacing.

C. Single-site model
Due to technological limitations of current quantum hardware, it appears prudent to start investigating scalar Yukawa interaction on a single lattice site with low truncation numbers for the real scalar. While for large lattices the gate cost of the time evolution scales more favourably in real space, the momentum-space representation is still suitable for small lattices -especially in the extreme case of a single-site model. The single grid point in momentum space corresponds to a momentum below a certain small cutoff value for all particles, implying that all particles can be considered to be approximately at rest. While zero-dimensional, single-site models of this type are widely used in the context of the DQS of high-/medium energy models [57,58], it is obvious that they cannot be expected to yield accurate results for any relevant physical observable. However, such "toy models" can still display nontrivial quantum dynamics and can thus be be seen as a useful first step towards more complex simulations that should be achievable in the not-toodistant future.
As derived in Appendix A, the effective interaction Hamiltonian of our single-site model is given by with η ≡ 4mgβ 3/2 being the effective coupling strength. At the same time, the single-site version of the free Hamiltonian H 0 = H Dirac + H Klein-Gordon is given by The sum of the Hamiltonians in Eqs. (6) and (7) -i.e. the Hamiltonian H tot = H 0 + H int -will be our point of departure for simulating the dynamics following a quench of the scalar Yukawa coupling in what follows. It is worthwhile pointing out that the single-site (zerodimensional) coupling Hamiltonian in Eq. (6) features a fermion-boson interaction of the density-displacement type, i.e. it describes coupling of the fermion (antifermion) density c † c (a † a) with the boson displacements b+b † . One widely known model that features this type of interaction is the time-honoured Holstein model [59] from condensed-matter physics, which describes a local coupling of a single electron (hole) to dispersionless phonons on each site of a discrete lattice, these phonons being described as zero-dimensional bosons. While this last model also involves zero-dimensional bosons, its principal difference from our effective model is that a fermion in the Holstein model -which can be considered spinless -is itinerant in nature (i.e. has a nonzero hopping amplitude between different lattice sites). Consequently, this model is discussed in scenarios that involve at least two lattice sites. Therefore, the analogy to the Holstein model might be more relevant for a multi-site generalization of our effective model.

III. DQS OF HAMILTONIAN DYNAMICS AND ENCODING PARTICLE STATES ON QUBITS
To faciliatate the discussion of the DQS of the Yukawacoupled system in the remainder of this paper, we begin with some general considerations of the DQS of the dynamics following an interaction quench and quantumcircuit synthesis (Sec. III A). We then recapitulate the basic aspects of encoding the relevant particle states on qubits, both for fermions (Sec. III B) and bosons (Sec. III C).

A. Simulating the dynamics following an interaction quench
In what follows, the dynamics of the fermion-boson system at hand following a scalar-Yukawa interaction quench at t = 0 will be considered within the DQS framework. We will evaluate the expected particle (fermion, boson) numbers at an arbitrary time t after the quench starting from several different initial states, i.e. for different choices of the initial fermion-and boson populations in a single-site realization of the fermion-boson system at hand. To this end, we will first evaluate the timeevolution operator of the system e −iHtott by first representing it in the form of a quantum circuit, i.e. as a sequence of single-and two-qubit gates.
It is useful to recall that the unitary time-evolution operator U (t) corresponding to the Hamiltonian H(t) satisfies the time-dependent Schrödinger equation in the operator form ∂U/∂t = −iH(t)U (t), with the initial condition U (t = 0) = ½. We will hereafter assume that we are dealing with a time-independent Hamiltonian H(t) = H sys . By discretizing time into n time steps of duration ∆t, such that t ≡ n∆t is the total evolution time, the timeevolution operator U (t) at time t can be expressed as The inherent tradeoff pertaining to such a decomposition is that a smaller time step corresponds to a longer circuit. The last expression for U (t) is typically approximated using first-order Trotter-Suzuki-type decomposition [47,60]. For a time-independent Hamiltonian H sys = l H l the latter approximates exp(−iH sys t) by ( l e −iH l t/n ) n , with the corresponding error being upper bounded by where N is the number of qubits [1]. The time-evolution operator of the system at hand is represented through quantum circuits that emulate the system dynamics. Those circuit are decomposed into single-qubit and two-qubit gates. Typical single-qubit gates include the Pauli gates [61] the single-qubit rotation gates as well as the Hadamard-(H), phase-(S) and π/8 (T ) gates The essential two-qubit gate is controlled-NOT (CNOT), whose matrix representation is Generally speaking, simulating the dynamics of a quantum many-body system on a classical computer is typically a highly nontrivial problem due to the large Hilbert-space size that necessitates the allotment of a large amount of memory. The task of designing a circuit that emulates these same dynamics on a quantum computer [61] -a process often referred to as unitary synthesis -is exponentially hard. This is intimately related to exponential growth of a generic Hilbert space with the system size; for instance, for interacting fermion systems it is known that generating such a circuit is an NP-hard problem. What further complicates the already challenging task of unitary synthesis is the fact that near-term quantum hardware is both very limited in size and inevitably noisy [44]. Therefore, optimizing the quantum circuit that is supposed to mimic the dynamics of a relevant quantum system -i.e. reducing its depth in order to obtain as shallow a circuit as possible -is a task of paramount importance in the realm of DQS.

B. Fermion-state encoding
In the familiar Jordan-Wigner (JW) encoding [62] of fermion states, due to the Pauli exclusion principle, the fermionic occupation numbers are restricted to the set {0, 1}. This allows for a one-to-one mapping from the occupation-number basis to the computational basis. In other words, a qubit in the state |0 (|1 ) corresponds to an empty (occupied) fermionic orbital where f j is the occupation number of orbital j (j = 0, . . . , n − 1) and |q j the corresponding qubit state It should be emphasized that in writing Eq. (13) we assumed that the rightmost single-qubit state corresponds to the qubit 0, a convention that we will hereafter also utilize for operators. The JW-type encoding requires O(n) qubit operations to represent a fermionic operator, where n is the number of fermionic orbitals (or, in the problem at hand, momentum states of spin-up and spin-down particles/antiparticles). While the alternative, Bravyi-Kitaev encoding [10,63] reduces this cost to O(log n), in what follows we will rely on the JW encoding as we will be concerned with Hamiltonians that are diagonal on the fermionic subspace in the occupation-number basis.
The JW qubit mappings of the fermion creation-and annihilation operators are, respectively, given by Here the operators (X j ± iY j )/2 change the occupation numbers of target spin orbital, while the action of the string of Z operators amounts to computing the parity of the state.

C. Boson-state encoding
Bosonic creation/annihilation operators satisfy the commutation relation In contrast to fermions, bosonic operators on different sites commute, thus they can be constructed such that they act locally on a bosonic mode. The main challenge in encoding bosonic states pertains to their unbounded occupation number, which forces us to impose an occupation-number cutoff Λ. The encoding of a bosonic Fock state with cutoff Λ requires O(n log Λ) qubits, where n is the number of bosonic modes. This scaling is achieved by representing the occupation numbers of the basis states of the truncated Fock space as binary strings and mapping the digits to single qubits [21]. For a Fock state |k , where k is an integer whose binary decomposition reads the mapping is given by The mapping of the truncated creation/annihilation operators to their pseudospin-1/2 (qubit) counterparts is then achieved by finding their corresponding Pauli-basis representations for arbitary truncations. For a truncation at Λ = 2 N − 1, a bosonic creation operator can be represented as where F j,k is an operator defined as (18) Another often used operator is the particle-number operator b † b, which in the Pauli basis is given by

IV. CIRCUIT DESIGN
In the following we describe the design of quantum circuits emulating the dynamics of the system at hand. We start with the constant-depth two-qubit circuit that corresponds to the one-boson exchange case (Sec. IV A), followed by its three-qubit counterpart in the case of three-boson exchange (Sec. IV B).
Before embarking on the construction of specific circuits, we point out how a specific property of the system at hand -namely, the local charge conservation -allows us to encode the fermion-antifermion sector of the problem using only one qubit. Starting from a general expression originating from JW mapping, namely we restrict the system to the subspace with a total charge of Q = 0. This allows us to rewrite a † a + c † c − 1 using only a single-qubit operator: This explains why the fermion-antifermion sector of the problem can be encoded using a single qubit. A. Exchange of up to one boson (two-qubit circuit) The next step is to apply both the fermion and boson mappings from Secs. III B and III C, respectively, to obtain the qubit Hamiltonian. Here we truncate the bosonic occupation number at Λ = 1 such that b + b † = X. By making use of Eq. (21), the two contributions to the total qubit Hamiltonian H tot = H 0 + H int are given by The spectrum of the free Hamiltonian H 0 is symmetrized by dropping the constant term (2M + m)/2. This changes the time evolution only by an irrelevant global phase.
The general approach to perform the time evolution of such a Hamiltonian with non-commuting terms is based on the Trotter-Suzuki decomposition, also known as Trotterization. To be more specific, for the Hamiltonian in Eq. (22) our starting point is the standard second-order Trotter-Suzuki product formula [47,48] where the free Hamiltonian H 0 is used for symmetrization as its cost in terms of CNOT gates is zero. The circuit for a single Trotter step, based on the "star" configuration of CNOT gates [see Fig. 8(a)], is depicted in Fig. 1.
In the conventional scenario of using the product formula in Eq. (23), one obtains a quantum circuit whose depth scales linearly with the number of Trotter steps. In other words, the circuit depth grows with the total simulation time, which in many systems limits feasible simulations to relatively short times. We show that in the system at hand -at least in the two-qubit case (i.e. in the case of up to one boson exchange) -we can defy this conventional scenario and design a constant-depth circuit. To this end, it is instructive to start by analyzing the transition zone between subsequent Trotter steps, which is represented by the following circuit: The simplification of the circuit in Fig. 2 is made possible by simple properties of the CNOT gate -namely, its commutation with a z-rotation gate on the control qubit and with an x-rotation gate on the target qubit, as well as the fact that the CNOT gate is self-inverse. It is also worthwhile pointing out that the occurrence of x-rotation gates on the target-qubit wire in the last circuit originates from the identity R X (θ) = HR Z (θ)H, which follows from the fact that X = HZH.
From this last circuit, one can easily go back to the exact time evolution by taking the limit ∆t → 0, n → ∞ while keeping n∆t ≡ t constant. This approach yields the exact constant-depth quantum circuit that corresponds to infinitely many small time steps, i.e. to the finite evolution time t; this circuit is shown in Fig. 3, where the top wire corresponds to the fermion, the bottom one to the boson. It should be emphasized that it is the specific form of the total Hamiltonian H tot = H 0 + H int in Eq. (22) that allows one to perform circuit compression between different Trotter steps such that only two CNOT gates are required to perform an arbitrary number of steps. At this point, the problem of decomposing the unitary two-qubit time evolution operator has been reduced to the decomposition of a single-qubit gate. We choose to break down exp [i(mX + ηZ)t/2] by performing a ZXZ Euler and ω 2 = m 2 + η 2 . The final circuit, which emulates the system dynamics over a finite-evolution time t using only two CNOT gates, then assumes the form depicted in Fig. 4. In connection with the form of this circuit, it is worthwhile pointing out that circuit compressions of the type utilized here have quite recently been discussed in the context of interacting qubit arrays described by the transverse-field Ising-and XY models [45,46], paradigmatic models in condensed-matter physics. These models can be mapped to free fermionic models and are also known to be classically simulatable with polynomial resources [64]. It is thus interesting that the fermion-boson system under consideration offers another, much less common example -namely, that of an XZ-coupled pair of qubits, with mutually unequal external fields in the z-direction acting on either qubit [cf. Eq. (22)] -where such a circuit compression is also possible. In order to describe an exchange of up to three bosons, two qubits are required to encode the bosonic Fock space. Using the bosonic qubit mapping, we modify the Hamiltonian as To perform the time evolution of this Hamiltonian, we start with Trotterization. As this Hamiltonian does not allow for the same efficient circuit compression as in the one-boson case, we have to be more careful how to trotterize. To reduce the error of a Trotter step to O(∆t 3 ), we once again make use of the second-order Trotter-Suzuki product formula [cf. Eq. (23)], also using the free Hamiltonian H 0 for symmetrization. While a first-order Trotterization of exp(−iH int ∆t) requires 8 CNOT gates, we will instead construct an exact decomposition with the same gate count. The key idea is that H int can be written as a tensor product of operators acting locally on the fermionic or bosonic registers. To calculate the matrix exponential, we may diagonalize these operators separately. The fermionic operator Z is already diagonal, leaving us with the task of diagonalizing b + b † . We find b + b † = SDS † , where in the matrix forms D and S are given by with λ ± ≡ 3 ± √ 6 andλ ± = λ ± / √ 3. The next step is to perform the Kraus-Cirac decomposition of S, a transformation intimately related to the Cartan decomposition of the Lie algebra su(4) [65,66]. In this manner, we find that where the matrices K 1 , . . . , K 4 ∈ SU(2) are given by with ϕ = arctan √ 2/(1 + √ 3) . The minimal circuit representation of S is depicted in Fig. 5. The next step is to express D [cf. Eq. (27)] in the Pauli basis: Then, the Hamiltonian H int can be recast in the form where η ± = η (3 ± √ 3)/2. The exact circuit for exp(−iH int, diag ∆t) is then given by Fig. 6. Finally, by combining the circuits shown in Figs. 5 and 6 with those corresponding to the free Hamiltonian H 0 , we obtain the circuit that emulates a single Trotter step of the system under consideration. This circuit is depicted in Fig. 7.

V. CNOT-COST ESTIMATION FOR HIGHER BOSON-NUMBER TRUNCATIONS
For simulations with higher truncation numbers (or more grid points), a more systematic approach towards circuit optimization is required due to the fast-growing number of Pauli strings. Even when considering only a single grid point with a boson-number truncation at Λ = 7 (three bosonic qubits), circuit optimization by hand in the problem under consideration becomes rather tedious. Namely, for different possible orderings of Pauli strings, we obtain different CNOT counts.
In what follows, we will focus on first-order Trotterization, which means that each Pauli string appears exactly once per Trotter step. All considerations in this section are based on the star + ancilla layout [61,67], which is depicted in Fig. 8(b). Using this technique, we can exactly simulate any Hamiltonian of the form where σ k c(k) ∈ {I, X, Y, Z} is a Pauli operator acting on the k-th qubit and H is given by a single Pauli string. The design of a quantum circuit emulating the time evolution governed by such a Hamiltonian is explained in Appendix C. For Hamiltonians with multiple non-commuting Pauli strings, we append the same circuit structure for each string and justify the circuit decomposition through the Trotter product formula.
Qubits with identity operations can be disregarded and qubits with X or Y are transformed to Z using X = HZH and Y = SHZHS † . Whenever we deal with multiple non-commuting Pauli strings, we use Trotterization to recover the circuit structure from Fig. 8(b) for each string.
For the sake of simplicity, throughout this section, we ignore the fermionic qubit in the interaction Hamiltonian H int , and focus entirely on the bosonic term Our goal is therefore to estimate the CNOT-cost of the circuit representation of exp (−iH b ∆t) for an arbitray truncation Λ. We restrict ourselves the case where the truncation can be expressed as Λ = 2 N − 1, where N is in integer corresponding to the amount of qubits forming the bosonic register. Let S N = {P 1 , P 2 , . . . , P k } denote the set of Pauli strings required to encode b + b † . In Appendix B, we show that the number of Pauli strings k for a truncation to N qubits is given by FIG. 9: Time evolution circuit for the Hamiltonian H = XX + ZZ using the Star+Ancilla layout. We define the zones on the left (right) of the first (last) RZ-gate as the start (end) layer. A zone between two subsequent RZ-gates is referred to as a transition zone.
To find the optimal order of these strings, we use an approach based on an analogy with the TSP; a similar technique was proposed in Ref. [67]. The formulation of this optimization task as a TSP works as follows. We first define a complete graph G on k nodes and set the distance between each pair of nodes i and j to |P i − P j | CNOT , where |P i − P j | CNOT denotes the CNOT-cost of the transition zone between the circuits representing exp (−iP i ∆t) and exp (−iP j ∆t). A graphical illustration of the definition of the transition zone is provided in Fig. 9. In Appendix D we prove that this CNOT-cost is precisely given by the Hamming distance Compared to the distance metric proposed in Ref. [67], this reduces the CNOT-cost by up to 50%. TSP aims to minimize the cost of a closed path visiting all nodes exactly once. For a circuit made of k Pauli strings, there are k − 1 transition zones. A closed path on the graph described above would however consist of k transition zones. This issue can be fixed by inserting an ancilla node which is fully connected with a cost of zero [67]. We build on this approach by taking the additional CNOTcost from the start-and end layers of the circuit (Fig. 9) into consideration. We find the cost of such a layer to be given by the Hamming weight The closed path on our modified graph now contains k + 1 edges, of which two correspond to the start-and end layers, meaning that we have k − 1 transition zones, as desired. Last, to build our graph, we need to know the Pauli strings explicitly. We generate the sets of Pauli strings recursively using with S 1 = {X}, where the tensor product is performed element-wise (e.g., S N ⊗ Z = {P 1 ⊗ Z, . . . , P k ⊗ Z}) and S(O) denotes the set of Pauli strings building the operator O. A naive expansion of the last set S(σ ⊗N + ⊗ σ − + σ ⊗N − ⊗ σ + ) gives us 2 N +1 different strings where each character can be either X or Y . But as the operator is Hermitian, all strings with imaginary coefficients cancel out. These strings are precisely the ones with an odd number of Y s. This leaves us with a subset containing 2 N strings with an even number of Y s. An example of the completed graph for a truncation at Λ = 3 is shown in Fig. 10.
To solve the TSP exactly, we use the Bellman-Held-Karp dynamic-programming algorithm [52,53], which scales in time with O(k 2 2 k ). Due to the exponential growth of the number of Pauli strings, we are limited to truncations with N ≤ 3 (Λ ≤ 7). As we are dealing with a metric graph, TSP can be 1.5-approximated  in polynomial time O(k 3 ) using the Christofides algorithm [54]. We make use of this heuristic to approximate the CNOT-cost for truncations with N ≤ 8 (Λ ≤ 255). Last, we compare to an upper bound on the cost derived in Appendix D. Figure 11 shows the results obtained using all three methods. We find, as expected, that the exact TSP produces the best results. The heuristic outperforms the upper bound we propose, suggesting that a tighter bound can be found.

VI. RESULTS AND DISCUSSION
In what follows, we present and discuss the obtained results for the cases of one-boson-(Sec. VI A) and threeboson exchange (Sec. VI B). We benchmark the accuracy of our DQS by comparing the results obtained on IBM Q processors with those resulting from a numerically-exact treatment of the quantum dynamics of the model under consideration on a classical computer.
To show that our single-site model of scalar Yukawa coupling exhibits nontrivial quantum dynamics, we perform DQS of this model for a variety of initial states of the fermion-boson system. After preparing the desired initial state, we construct the states |ψ(t) of the system at different times t by applying the designed two-qubitand three-qubit circuits that emulate the system dynamics following a Yukawa-interaction quench at t = 0. Finally, for each of those constructed states we perform measurements that give us access to the desired quantities, bearing in mind the relevant expressions in the Pauli basis [cf. Eqs. (14) and (19)].
Given that the DQS of the single-site fermion-boson model under consideration only requires two-and threequbit systems [cf. Sec. IV], for benchmarking purposes we made use of two freely available IBM Q processors. More precisely, we performed DQS in the case of up to one boson exchange on the five-qubit ibm oslo processor, while for the case of up to three bosons we used the seven-qubit ibmq quito. The qubit-connectivity graphs (coupling maps) of these two quantum processors are de-picted in Figs. 12(a) and 12(b), respectively.
Generally speaking, the errors in numerical experiments on IBM Q devices strongly depend on the number of CNOT gates in the corresponding quantum circuit. The principal reason for this is that the CNOT-gate error is an order of magnitude larger than that of singlequbit gates. Another reason is that the CNOT-gate time is much longer than that of its single-qubit counterparts, which leads to the accumulation of errors due to energy relaxation and dephasing, this two processes being quantified by the respective decoherence times T 1 and T 2 . Therefore, in order to minimize the error, in our DQS we make use of the connected subsets of qubits with the smallest average CNOT error.

A. Exchange of up to one boson
For the DQS of the quench dynamics in the case with up to one-boson exchange (two-qubit case) we utilized the IBM Q processor ibm oslo (average CNOT error on this device is 1.321 × 10 −2 , average readout assignment error 2.596×10 −2 , each data point averaged over 20, 000 shots), whose corresponding qubit-connectivity graph is shown in Fig. 12(a). In line with the strategy indicated above, for the qubit pair that we selected (physical qubits 3 and 5) the CNOT error is smaller -only 5.606 × 10 −3 . At the same time, the readout assignment errors are 1.850×10 −2 (qubit 3, used for encoding bosons) and 9.9 × 10 −3 (qubit 5, used to encode fermions).
The obtained results in the two-qubit case are illustrated in Figs. 13 and 14, which show the expectation values of the particle-number operators -both for fermions and bosons -at an arbitrary time t after the quench for various initial states. In particular, initial states with a particle number of 1/2 were prepared such that |ψ 0 = (|0 + |1 )/ √ 2. What can be inferred from both figures is that the simulation of the compressed circuit with only two CNOT gates [cf. Sec. IV A] shows an excellent agreement with the exact (continuous-time) quantum dynamics. Importantly, the small discrepancy between the data points and the exact results can mostly be ascribed to the readout assignment errors.
Another important observation is that for some of the chosen initial states, examples of which are presented in Fig. 13, the fermion-boson dynamics is trivial. In other words, these states do not cause any interaction between fermions and bosons, which is manifested by the fact that the ensuing boson dynamics are completely independent of the fermion number. We can straightforwardly see why that is the case by finding an explicit expression for the expected boson number at time t based on Eq. (19): By considering separable initial states |ψ 0 = |ψ 0 b ⊗ |ψ 0 f , we can infer that the fermion-boson interaction vanishes for states with ψ 0 |Z|ψ 0 f = 0 or ψ 0 |X|ψ 0 b = ψ 0 |Y |ψ 0 b = 0. By contrast to Fig. 13, nontrivial quantum dynamics of fermions and bosons are illustrated in Fig. 14, where the boson dynamics clearly depends on the fermion number. In particular, Fig. 14(b) shows a special case where both fermion-and boson numbers are conserved.
The states from the previous simulation in Figs DQS results show a reasonably good agreement with the exact results for the first few Trotter steps, but subsequently start to deviate visibly. These deviations are mainly caused by the CNOT-gate error. For instance, given that we employed up to n = 10 Trotter steps and that our designed three-qubit circuits entail 8 CNOTs per one step, obtaining the last data points requires one to carry out 80 CNOTs. While the Trotter error also has to be taken into account, by performing idealized simulations of the Trotter circuit we have verified that this error is low compared to the accumulated CNOT error.
Initial states with a fermion number of 1/2 were once again prepared such that |ψ 0 f = (|0 + |1 )/ √ 2. The initial states for bosons in Figs. 15 and 16 were prepared the same way as for the two-qubit circuit, with the natural extension to the higher-dimensional Hilbert space by setting the additional qubit to state |0 . While the initial state in Fig. 17 also corresponds to an occupation number of 1/2 as in Fig. 16, it lies in the extended Hilbert space and is encoded as |ψ 0 b = ( √ 3 |00 + |10 )/2. Finally, the boson state in Fig. 18 is given by |ψ 0 b = (|00 + |01 + |10 + |11 )/2. Figures 15 and 16 show states with a non-interacting behaviour, while the states in Figs. 17 and 18 exhibit fermion-boson interactions. Note that the states from Fig. 16 did interact in the twoqubit circuit.
By inspecting different initial states, it can be inferred that the agreement between the simulation-and exact results is much better in certain special cases, such as, for example, that of Fig. 18(b). This can likely be attributed to the fact that the readout assigment error is small for states like |0 + |1 , while at the same time being large for states such as |0 and |1 .

VII. SUMMARY AND CONCLUSIONS
To summarize, in this paper we investigated the nonequilibrium dynamics of a coupled fermion-boson system following a quench of Yukawa-type interaction within the framework of digital quantum simulation. We showed that a single-site abstraction of scalar Yukawa coupling can successfully be simulated on the existing IBMQ hardware in the limit of low boson-number exchange. Using advanced circuit-optimization techniques -exemplified by the Kraus-Cirac decomposition -we designed efficient, low-depth quantum circuits for simulating the exchange of up to three bosons. Based on these circuits we computed the expected boson-and fermion numbers at an arbitrary time after the quench and demonstrated a good agreement of the obtained results with their classical counterpart.
In the one-boson case we designed a constant-depth circuit with only two CNOT gates, regardless of the total simulation time. This constitutes an example of circuit compression, which is only possible for certain special types of the underlying system Hamiltonian and allows one to defy the usual linear scaling of the circuit depth with the number of Trotter steps. In the three-boson case -where such a compression is no longer possiblewe designed a circuit in which one Trotter step requires 8 CNOT gates, which is far below the maximal CNOTcost of an arbitrary three-qubit gate (14 CNOTs) [49]. Finally, by making use of an analogy with the travelling salesman problem, we derived a CNOT-cost estimate for quantum circuits emulating the system dynamics for higher boson-number truncations.
The findings of the present study motivate future work in several different directions. For example, it is of interest to extend our single-site digital quantum simulation to multi-site realizations, including much larger systems for which benchmarking may require access to quantum devices with a smaller two-qubit gate error and all-toall connectivity [68]. For such systems, tensor networks may provide an efficient approximation for certain classes of quantum states [69]. Finally, accurate simulations of larger systems might be enabled through the use of circuit-optimization techniques based on machine learning [70,71] and quantum-error mitigation [72][73][74][75][76][77][78].
where the momentum argument index/argument is dropped. For the spinor products in the non-relativistic limit we haveū By taking into account these last identities, the single-site Hamiltonian in Eq. (A5) reduces to Finally, we ignore the spin degrees of freedom and restrict ourselves to the states with the vanishing total charge, which yields At this point we can express the lattice spacing a l in terms of the boson mass via a l /(2m) = βγ. Because β ought to be small for our single-site abstraction (the nonrelativistic limit), we use γ ≈ 1. In this manner, we finally obtain the zero-dimensional interaction Hamiltonian with η ≡ 4mgβ 3/2 being the effective coupling strength.

Appendix B: Generation of Pauli strings
In this appendix, we derive an algorithm to efficiently calculate the Pauli strings for a truncation using N qubits. Our approach is based on finding a recurrence relation between b + b † truncated with N and N + 1 qubits [79]. Writing out b + b † in the Fock state basis, we find Notice that the top left block is nothing but (b + b † ) N . In order to succinctly express the recurrence relation, we define I ± = (I ± Z)/2 and σ ± = (X ± iY )/2. The expression we find is given by where (b + b † ) 1 = X and M BR,N is the bottom right structure in Eq. B1. The reason we do not write it out explicitly is that M BR,N consists of the same Pauli strings as (b + b † ) N . Bearing this in mind, we write down the recurrence relation for the set of Pauli strings S N [cf. Eq. (37)], with S 1 = {X}.
Here, the Pauli operators I and Z come from the expansion of I ± . As discussed in Sec. V, the set S(σ ⊗N + ⊗ σ − + σ ⊗N − ⊗ σ + ) consists of 2 N Pauli strings, which implies that the number of Pauli strings |S N | satisfies the difference equation with |S 1 | = 1. The explicit solution to this equation is given by Appendix C: The distance metric In this section, we base our observations on the star+ancilla layout, depicted in Fig. 8(b). While the pure star layout from Fig. 8 (a) makes use of a smaller number of CNOT gates, the ancilla qubit allows for easier gate cancellation in the transition zones, as no basis transformations are applied on the target qubit.
Let us consider the transition zone of two consecutive Pauli strings P 1 and P 2 . We will briefly describe the approach of Ref. [67] before introducing our improved cancellation technique. On indices where P 1 [i] = P 2 [i], the CNOT gates cancel out. But for indices where P 1 [i] = P 2 [i], one has to differentiate between two cases: (a) neither P 1 [i] nor P 2 [i] is I or (b) one of them is I. Case (a) requires two CNOTs, while (b) only needs one. Combining these results, the CNOT distance metric is defined as We will now modify this approach by optimizing case (a). We start by considering transitions of the type X → Z or Z → X, represented by the circuit in Fig. 19, which can be reduced to a single CNOT gate. Next, we investigate the Z → Y and Y → Z transitions. Of course, Y → Z is just the inverse of Z → Y and can be obtained as the Hermitian conjugate. The circuit for Z → Y is shown in Fig. 20. In the first step, we used that S commutes with a control, and in the second one, we reuse the circuit from Fig. 19 for the X → Z transition. It remains to simplify the circuit for X → Y , depicted in Fig. 21. In the first step, we use that HSH = S † HS † up to a global phase and the commutation of S † with the control. In the second step, we once again recycle the circuit for X → Z and undo the global phase. All these results together prove that any transition with P 1 [i] = P 2 [i] can be implemented with only one CNOT gate, which implies that we have where |P 1 − P 2 | Ham is the Hamming distance counting the number of disagreeing indices. Compared to the previous cost function in Eq. (C1) this cuts the cost by up to 50%.
Appendix D: Upper bound on the CNOT-cost Using the recurrence formula for the Pauli strings (37) and the functions for the CNOT-cost of the start-and end layer as well as the transition zones, we will now derive an upper bound on the CNOT-cost of the circuit representing exp (−iH b ∆t) in first-order Trotterization. We first recall the recurrence formula for the Pauli strings [cf. Eq. (37)] where S 1 = {X}. In the following, C(O) denotes the CNOT-cost of implementing the Pauli strings of the operator O in first-order Trotterization. Using that C(S N ⊗ I) = C(S N ), C(S N ⊗ Z) = C(S N ) + 2 and the fact, that we can always cancel 2N CNOT-gates when transitioning from S N ⊗ I to S N ⊗ Z, we find For the second part, namely S(σ ⊗N + ⊗ σ − + σ ⊗N − ⊗ σ + ), we recall that this set contains all (N + 1)-digit Pauli strings made from X and Y with an even number of Y s. The start-and end layers of this set give us 2(N + 1) CNOTs. To maximize the gate cancellation, we arrange the set using every other element of an (N + 1)-digit gray code. This way, only two characters change per transition. For 2 N − 1 transition zones, we then find a total cost of 2 N +1 − 2. Combining these results, we have with C(S 1 ) = 2. The explicit solution to this recurrence relation is given by