Experimental quantum simulation of fermion-antifermion scattering via boson exchange in a trapped ion

Quantum field theories describe a variety of fundamental phenomena in physics. However, their study often involves cumbersome numerical simulations. Quantum simulators, on the other hand, may outperform classical computational capacities due to their potential scalability. Here we report an experimental realization of a quantum simulation of fermion–antifermion scattering mediated by bosonic modes, using a multilevel trapped ion, which is a simplified model of fermion scattering in both perturbative and non-perturbative quantum electrodynamics. The simulated model exhibits prototypical features in quantum field theory including particle pair creation and annihilation, as well as self-energy interactions. These are experimentally observed by manipulating four internal levels of a 171Yb+ trapped ion, where we encode the fermionic modes, and two motional degrees of freedom that simulate the bosonic modes. Our experiment establishes an avenue towards the efficient implementation of field modes, which may prove useful in studies of quantum field theories including non-perturbative regimes.

Q uantum simulators are devices designed to predict the properties of physical models associated with target quantum systems 1,2 . Their intrinsic physical behaviors are fully governed by the laws of quantum mechanics, making it possible to efficiently study complex quantum systems that cannot be solved by classical computers 3,4 . Trapped ions and superconducting circuits have proved to be promising for experimentally simulating a variety of paradigmatic quantum models such as various spin models [5][6][7][8][9] , relativistic Dirac equations [10][11][12][13] , embedding quantum simulators [14][15][16][17][18] , and fermionic models 19,20 . More recently, a digital quantum simulation of a fermionic lattice gauge theory has been performed in trapped ions 21 . However, it would be desirable to realize a quantum simulator that involves interacting fermionic and bosonic fields as described by quantum field theories (QFTs) 22 . In this sense, fermionic modes could be mapped in the ion internal levels, whereas bosonic modes could be naturally encoded in the motional degrees of freedom.
Here we report an experimental quantum simulation of interacting fermionic and bosonic quantum field modes, where fermions are encoded in four internal levels of an 171 Yb + ion and the bosonic modes in the motional degrees of freedom, following the proposal by Casanova et al. 23 . Therefore, this analog quantum simulation constitutes a step forward towards a digital-analog quantum simulator 19,[24][25][26][27] of perturbative and non-perturbative QFTs. In this sense, a remarkable feature of our experiment is that it contains all orders in perturbation theory, which is equivalent to all Feynman diagrams for a finite number of fermionic and bosonic modes. Moreover, our approach could in principle be scaled up by progressively adding more ions allowing the codification of additional fermionic and bosonic field modes, which may lead to full quantum simulations of QFTs such as quantum electrodynamics (QEDs) 22 .

Results
Hamiltonian for quantum simulation of QFT. The common way to analyze QFTs is via a Dyson series expansion in perturbation theory and Feynman diagrams 22 . If we consider larger coupling parameters, standard perturbative methods become cumbersome for a finite-mode Dyson expansion, mainly because only a reduced number of perturbative Feynman diagrams can be calculated. On the other hand, a trapped-ion quantum simulator with its high degree of quantum control 28 could overcome these limitations and simulate QFTs more efficiently than classical computers 29 . Based on the proposal of ref. 23 , our experimental quantum simulation of finite-number interacting quantized field modes includes all terms of the Dyson expansion. We experimentally implement a fundamental QFT model in a single trapped-ion considering (i) one fermionic and one antifermionic field mode, (ii) one or two bosonic field modes, which already reveals interesting QFT features such as self-interactions, particle creation and annihilation, and perturbative and non-perturbative regimes. The general Hamiltonian involving the continuum of fermionic and bosonic fields reads where b p and d p are fermionic and antifermionic annihilation operators, respectively, whereas a k are the bosonic annihilation operators. Here, ω (ω k ) is the fermion and antifermion free energy (boson free energy), whereas ψ(x) denotes the fermionic and A(x) the bosonic fields 23 . As a stepped experimental demonstration, we first consider the simplest situation with only one bosonic mode, which can be implemented by a single vibrational mode of the ion. The fermion and antifermion modes are considered as two comoving modes describing incoming Gaussian wave packets, which are centered in the average momentum and have distant average initial positions 23 . These modes describe self-interacting dressed states by emission and absorption of virtual bosons. They also represent the lowest-order in perturbation theory of the scattering of the incoming wave packets that will collide in a certain region of spacetime. The pair creation and annihilation is local and takes place only when the two wave packets of the fermion and antifermion overlap, namely, when the particles scatter. A diagram of these interactions, in the spirit of a Feynman diagram, is shown in Fig. 1(a). It is noteworthy that the loop of this diagram includes all terms in a finite-mode Dyson expansion, which is different from the standard perturbative approach with only a reduced number of Feynman diagrams. By considering slow massive bosons, as described in ref. 23 Fig. 1 Fermion-antifermion scattering process and its mapping to an 171 Yb + ion system. a Diagram of the interactions between a fermion, an antifermion, and bosons. The fermion emits and absorbs virtual bosons through the self-interaction process. In the fermion-antifermion scattering process, the middle dashed loop represents the summation of all terms in a finite-mode Dyson series expansion. b Diagram of the encoding and operations to implement the interaction Hamiltonian H I with an 171 Yb + trapped ion. The vacuum state and the fermionic states are mapped onto four internal states through the Jordan-Wigner mapping. The bosonic mode is directly implemented with the vibrational mode along the X radial direction. The self-interaction is implemented by a displacement operation, which shifts the center of the harmonic oscillator without changing the internal states. The fermion and antifermion scattering is simulated by the combination of the red-and the blue-sideband transitions, which change the internal states together with the vibrational mode Hamiltonian we would like to realize turns into where the associated time-dependent coupling strength is and δ ¼ ω f þ ω f À ω 0 . Here, ω f , ω f , and ω 0 represent the energy of the fermionic field mode b, the antifermionic field mode d, and the bosonic field mode a 0 , respectively. The ratio g 2 /g 1 gives the relative strength between pair creation and self-interaction. T is the total time of the pair-creation process, whereas σ t is the temporal interval of the interaction region. Our formalism, explained in detail in ref. 23 , involves considering incoming comoving fermionic and antifermionic modes at the lowest order in perturbation theory. The time dependence of the interaction of the incoming particles, as they collide, maps onto a time dependence of the interaction Hamiltonian coupling. Applying a Jordan-Wigner mapping 23 from the fermionic modes to a 2 × 2 Hilbert space, b y ¼Î σ þ ;b ¼Î σ À ; ð4Þ the vacuum state and fermionic states are represented by where 1 f 1 f denotes the state containing one fermion and one antifermion. Thus, the interaction Hamiltonian reads Trapped ion implementation. We point out that, due to the asymmetric role of fermionic annihilation and antifermionic creation operators in the fermionic field, the one-antifermion state is a dark state of the Hamiltonian in equation (2) and therefore the antifermion does not have self-energy at first order (it has, when considering more modes and higher orders). We implement the Hamiltonian on a single 171 Yb + ion trapped in a three-dimensional harmonic potential as shown in Fig. 1(b). The radial harmonic potential is generated by an oscillating electric field V RF in the radial direction with the two trapping frequencies along X and Y directions being (ω X , ω Y ) = (2π) (2.4, 1.9) MHz. The bosonic modes are mapped onto these radial vibrational modes and we choose the X mode for experiments involving a single bosonic mode. The vacuum state 0 f 0 f is mapped to the hyperfine state |F = 0, m = 0〉. Fermionic states are mapped onto Zeeman states as With the mapping of the bosonic mode and the fermionic states onto the 171 Yb + ion system, the Hamiltonian (7) is naturally divided into three parts: displacement, red-sideband, and blue-sideband operations.
The operations of the self-interaction and scattering processes of the fermion and the anti-fermion are realized by σ + -polarized Raman laser beams [30][31][32] counter-propagating along the direction of the magnetic fieldB. The strength of the magnetic field at the position of the ion is around 7 G, which produces ω ZM = (2π)10 MHz Zeeman splitting. The magnetic field is generated by a pair of Helmholtz coils and is aligned along the angle bisector  ( 2 ) The two AOMs are driven with different frequencies ω R1 and ω R2 , where ω R1 is fixed at (2π)231 MHz and ω R2 is adjusted in the range of (2π)233~253 MHz. Quarter-wave plates are used for polarization adjustment of the laser beams. b, c Frequency combs of the pico-second pulsed lasers and choice of the effective Raman beat-note frequency. The frequency interval of the "comb" is the repetition rate of the laser pulse, which is stabilized at ω rep = (2π)76.51 MHz. The frequency difference between the two AOMs is tuned near to the trap frequency ω X (b) for the displacement operation or (c) to produce the hyperfine frequency with the addition of 165 intervals, ω HF . d, e The basic level structure and transitions of 171 Yb + system coupled by σ + polarized Raman laser beams. The beat-note frequency of the Raman beams (d) for the displacement operation is Thick lines represent the two times stronger displacement operation on state 1 f 0 f . e The frequency difference between the carrier transition and the red-sideband (blue sideband) operation is NATURE COMMUNICATIONS | DOI: 10.1038/s41467-017-02507-y ARTICLE NATURE COMMUNICATIONS | (2018) 9:195 | DOI: 10.1038/s41467-017-02507-y | www.nature.com/naturecommunications direction of the X and Y axes, which allows the laser beams to couple both of the vibrational modes, as shown in Fig. 2(a). The laser beams are modulated with acousto-optic modulators (AOMs), which are driven with different frequencies ω R1 and ω R2 . For the Raman transitions, the mode-locked picosecond laser is used with a wavelength of 375 nm, which is Δ = (2π)12 THz red detuned from the optical transition 2 S 1/2 ↔ 2 P 1/2 . The train of laser pulses in the time domain can be considered as an equally spaced frequency "comb" 33 , which in our experiment had a repetition rate of ω rep = (2π)76.51 MHz. As shown in Fig. 2(b, c), we use the frequency "comb" to select a Raman beat-note frequency according to the relation ω R = Δω + n × ω rep , where Δω = ω R2 − ω R1 and n = 0, ± 1, ± 2, …. For transitions between different motional levels of the same electronic state, we simply use n = 0 and make Δω close to ω X or ω Y . For transitions between two different electronic states, we use n = 165. Figure 2(d, e) shows the Raman schemes needed to implement Hamiltonian (7), which are naturally divided into three parts, namely, ÀgðtÞ ÀgðtÞ Here, the first part is ω 0 -detuned displacement operation, the second part is δ-detuned red-sideband operation between 0 f 0 f $ 1 f 1 f , and the last part is (2ω 0 + δ)-detuned bluesideband operation between 0 f 0 f $ 1 f 1 f .
The first part corresponds to a displacement operation with 1 : 2 : 1 : 0 relative ratios among the strength coefficients of states 0 f 0 f , 1 f 0 f , 1 f 1 f , and 0 f 1 f . Figure 2(d) shows how to implement the displacement operation through the counterpropagating Raman laser beams shown in Fig. 2(a). The σ +polarized Raman beams produce the exact ratios in the strength of displacement operations, as state 1 f 0 f is coupled to two levels in the 2 P 1/2 manifold, states 0 f 0 f , and 1 f 1 f to one level, and state 0 f 1 f to no level. The strength coefficient of a Raman path is given by Ω R = g 1 g 2 /2Δ R , where g 1 and g 2 are Rabi frequencies of the two Raman beams coupled to the transition between 2 S 1/2 and 2 P 1/2 and the detuning Δ R ≈ Δ = (2π)12 THz. The coefficients of all possible Raman paths are added, as all optical transitions between 2 S 1/2 ↔ 2 P 1/2 states have the same coefficients in absolute values. We note that the frequency difference ω HF = (2π)12.6 GHz between states 0 f 0 f and 1 f 1 f is much smaller than the detuning Δ of the Raman laser beams acting on the manifold 2 P 1/2 , which produces around a 0.1% difference in the strength of the displacement operations. Finally, we measure the strength of the displacement operations and observe the ratios (see Methods). In principle, we can also implement other ratios of displacement operations by applying additional σand π-polarized Raman beams (see Methods).
The second and third parts are realized by the red-and the blue-sideband transitions as shown in Fig. 2(e). The timedependent strength-coefficient g(t) in equation (3) is implemented by the change of laser intensity, which is proportional to the RF power on the AOMs of Fig. 2(a). We generate the timedependent RF signal from an arbitrary waveform generator (AWG) and apply it to the AOM R2. By using the AWG, we can generate all the necessary RF frequencies and powers, which realizes the full Hamiltonian (7) containing the displacement operation, red-, and blue-sideband transitions.
Experimental procedure of the quantum simulation. In the experiment, we initialize the motional and internal state of the ion to the state 0 f 0 f ; n ¼ 0 by standard Doppler cooling, resolved sideband cooling, and optical pumping 34,35 . The residual average phonon number and the heating rate are measured to be 〈n〉 = 0.016 ± 0.025 and 3.8 ± 1.2 quanta s −1 , respectively. The heating effect can be neglected in the typical duration of a single simulation, which is of < 2 ms. Then we implement the target Hamiltonian (7) and let the system evolve for a time t. Finally, we measure the average boson number 〈n〉 and the populations of various fermionic states, as well as the correlation between the bosonic mode and the fermionic state. A detailed discussion of the measurement procedures can be found in the Methods section. We compare the experimental results with the ideal theoretical calculations. In our simple situation of single bosonic, fermion, and anti-fermion modes, we are able to numerically calculate the exact evolution with the full Hamiltonian and find a perturbation method that works for a short time dynamics. The whole evolution is then computed by accumulation of the latter. We note that such numerical methods would not be allowed as the system size grows. Typically, one considers the size corresponding to 50 qubits to be intractable. For example, a realistic situation with 16 ions, 16 modes, and 8 considered levels per mode would be beyond the capabilities of classical computers. Self-interaction and particle creation and annihilation. We first study the fermion self-interaction processes by setting g 2 = 0, starting from the initial state 1 f 0 f ; 0 b . Then the self-interacting dynamics occurs via the couplings 1 f ; 0 f ; n b ! 1 f ; 0 f ; n b ± 1 . Figure 3(a) shows experimental data for the time-dependent bosonic vacuum populations and the average boson numbers for different self-interaction strengths g 1 /ω 0 = 0.1 and 0.15, which quantitatively coincide with the theoretical calculations within experimental errors. We clearly observe the expected emission and reabsorption processes of virtual bosons and the growth of the average number of virtual bosons with the self-interaction strength g 1 .
Subsequently, we realize the annihilation of a fermion-antifermion pair and the creation of bosons with parameters g 1 = 0.01ω 0 , g 2 = 0.21ω 0 , and σ t = 3/ω 0 . We choose the initial state to be the state with a fermion-antifermion pair and no bosons: 1 f 1 f ; n ¼ 0 . Figure 3(b) shows the dynamics of the fermion-antifermion scattering process via the population of the fermionic-pair state and the average boson number. It can be clearly seen that the initial fermion-antifermion pair disappears, creating a single boson.
Next, we realize the process of scattering with parameters g 1 = 0.1ω 0 , g 2 = ω 0 , σ t = 4/ω 0 , where g 2 ≥ ω 0 . In this regime, the interaction Hamiltonian ((2)) cannot be regarded as a perturbation. In such a strong coupling situation, we cannot easily discriminate the contributions from the self-interaction and pair production processes. When the initial fermion-antifermion pair disappears, more than a single boson is created in the process as shown in Fig. 3(c), which is qualitatively different from the dynamics shown in Fig. 3(b). As the size of the Hilbert space is not too large, we numerically calculate the dynamics of the Hamiltonian by direct numerical integration (see Methods), which is in agreement with the experimental results shown in Fig. 3(c). However, as the number of fermion-antifermion pairs and bosons increases, the exact numerical calculation will be intractable by classical means. We have also developed a perturbation method based on the observation that, for a reasonably small time g 2 t ( 1, the effect of the coupling term g 2 does not produce a divergence in the dynamics. We divide the total time of the process by 100 and apply the perturbation method (see Methods) to the unitary evolution operator in each time slice. We find that after including terms up to the 7th order in the perturbation parameter, the deviation of the perturbative dynamics from the complete one is below 10 −4 . However, even this approach, based on a perturbative expansion within time slices, would be difficult to use for large Hilbert space dimensions with more fermions and bosons.
Finally, as a demonstration of scalability, we realize fermion self-interaction processes extended to two bosonic modes by using both X and Y phonon modes of a single trapped ion. We set g 1 = 0.15ω 0 , the first boson mode frequency ω 1 = ω 0 , and the second boson mode frequency ω 2 = 0.9ω 0 . We note that the g 1 (g 2 ) coupling to the mode Y (X) is negligible, as the detuning to the mode Y (X) is larger by a factor of 50, which effectively suppresses the strength by the same amount. We choose the initial state to contain one fermion and no bosons, 1 f 0 f ; n 1 ¼ 0; n 2 ¼ 0 . Then the self-interacting dynamics is given by the transition 1 f 0 f ; n 1 ; n 2 $ 1 f 0 f ; n 1 ± 1; n 2 ± 1 . As the bosonic modes have different frequencies, we observe that the considered fermion emits and reabsorbs bosons differently from the single-boson case. Instead of a sine curve, we see a clear beatnote behavior of the fermionic population as shown in Fig. 3(d). We also clearly observe the dynamics of both bosonic modes in a good agreement with the theoretical expectation. By increasing the number of bosonic modes, we would simulate the continuous regime of bosonic modes, which would be related to scattering experiments. In such large number of bosonic modes, the nonperturbative behavior of fermionic or antifermonic mode could be intractable. On the way of increasing bosonic modes, a technology of correlation measurement of multiple phonon modes could be applied 36,37 .

Discussion
In conclusion, this work considers an experimental quantum simulation of interacting fermionic and bosonic quantum field modes. Our approach could be in principle scaled up by progressively incorporating more fermionic and bosonic field modes, which may lead to a full-fledged digital-analog quantum simulation of QFTs such as QED [22][23][24] or the Holstein model 38 , where correlations between multiple fermions and phonons have critical relevance. In our current experimental system, an extension to multi-fermion and multi-phonon (bosonic) modes could in principle be implemented by loading a number of ions in a single trap, where the spins of ions map the fermionic modes through Jordan-Wigner transformation 24 and the vibrational modes of ions directly map the bosonic modes. The many-body operators or spin-spin interactions appearing after mapping of the fermionic modes onto spins can be efficiently implemented via a combination of two Mølmer-Sørensen gates and a local gate as shown in ref. 24 . Other than the spin-spin interactions in the Holstein model, e.g., the couplings between fermionic modes and  bosonic modes can be implemented by the same Raman laser beams that are individually addressing single ions and tuned to specific mode frequencies. In this respect, it has been shown that the number of gates grows polynomially as the number of fermions and bosons 38 . Ref. 38 also discussed the estimated infidelities from the gate errors in realistic experimental decoherence condition up to four sites, which clearly showed the degree of control is more demanding when the coupling strengths between the modes increase. As demonstrated in our experiment, we do not observe any clear degradation of the simulation when using two modes, although here we do not have the technical problem of individual addressing. We may implement the model in a fully analogue way together with proper spin-spin interactions [39][40][41] , which would allow us to study the pairing or polaron physics occurring in many unconventional superconducting systems 42,43 with the controls of various parameters. In particular, we remark that already with 16 two-level ions and 8 phononic levels per ion, one could perform quantum simulations of interacting quantum field modes that are beyond the reach of classical computations, i.e., a Hilbert space dimension of 16 16~264 , which would otherwise require a lengthy quantum algorithm with 64 qubits 44,45 . This experiment opens an avenue that aims at out-performing the limitations of classical computers, with in principle scalable quantum simulations.
We also point out that there are no known efficient classical algorithms for simulating interacting fermionic models in arbitrary spatial dimensions, whereas with our approach, with a trapped-ion quantum simulator, fermionic models in arbitrary dimensions could be analyzed with polynomial resources 24 . The verification of our proposed scalable experiment requires polynomial resources, as for the detection of the number of bosonic excitations produced, or the population of the fermionic or antifermionic states, only a polynomial number of measurements is required.

Methods
Uniform red sideband. The "uniform red sideband" 46-48 is implemented as an adiabatic transition where the transfer speed between 0 f 0 f ; n and 1 f 1 f ; n À 1 is the same for all n = 1, 2, …. It is realized by adding a time-dependent amplitude A (t) = sin(πt/d) and a time-dependent phase φ(t) = −1/β sin(πt/d) to the normal redsideband operation, and some additional terms to compensate for the AC Stark shift. Here, d = c π red is the duration of the transition, and β = ((l + 1) (h + 1)) −1/4 /c is an empiric parameter depending on the lower bound l and the upper bound h of the phonon number n. We typically choose c to be 10, such that the transition duration d is c/2 = 5 times the red-sideband operation period. Therefore, we achieve more than 99% of theoretical fidelity for all phonon numbers between l and h.
Displacement strength adjustment. We experimentally measure several strength coefficients to check the strength ratios depending on the electronic states. We first prepare the initial state m; n ¼ 0 j i , where m j i is either 0 f 0 f or 1 f 0 f . Then we apply the displacement operation for a small period τ. After this, we should obtain a coherent state m; α j i, α = Ωτ, where Ω is the desired strength coefficient. Subsequently, with several different τ, we fit the parameter Ω by measuring each time the remaining population on state m; n ¼ 0 j iwith the "uniform red sideband" method, which should be e Àα 2 ¼ e ÀΩ 2 τ 2 . After careful beam alignment and quarter wave plate adjustment, the measured strength coefficients of 0 f 0 f and 1 f 0 f are (2π)7.2 and (2π)14.4 kHz, respectively, which are consistent with the theory ratio 1:2.
As the magnetic quantum number is conserved during the displacement operation, two virtual optical transitions in a Raman path should have the same polarization. If both polarizations are purely σ − , then the relative strengthcoefficient ratio between states 0 f 0 f , 1 f 0 f , 1 f 1 f , and 0 f 1 f is 1 : 0 : 1 : 2. If both polarizations are purely π, then the relative strength-coefficient ratio between states 0 f 0 f , 1 f 0 f , 1 f 1 f , and 0 f 1 f is 1 : 1 : 1 : 1. In general, if the ratio between σ + , σ − , and π polarization is a : b : c, then the relative strength-coefficient ratio between states 0 f 0 f , 1 f 0 f , 1 f 1 f , and 0 f 1 f is a + b + c : 2a + c : a + b + c : 2b + c.
Fermionic state measurement. To measure P 1 f 1 f À Á , we simply apply a π rotation between states 0 f 0 f and 1 f 1 f to swap their populations, with a microwave horn. Then the measured population of state 0 f 0 f is equal to the original P 1 f 1 f À Á .
To measure P 1 f 0 f ; n ¼ 0 À Á , however, we need a phonon projective measurement [46][47][48] . Instead of using fluorescence detection together with a postselection scheme, which may introduce significant heating errors because of photon scattering, here we use an auxiliary state as a swap buffer. It is noteworthy that the interaction Hamiltonian (7) does not have terms related to state 0 f 1 f . Therefore, we employ state 0 f 1 f as the auxiliary state and always initialize it to zero. We first apply three consecutive π swap gates between 0 f 0 f and 0 f 1 f , and between 1 f 1 f and 1 f 0 f . After these operations, 1 f 1 f is swapped with 1 f 0 f , and 0 f 1 f is swapped with 0 f 0 f . Then we apply a "uniform red sideband" π rotation to swap the population in 0 f 0 f ; n>0 with that of 1 f 1 f ; n À 1 . Then, we measure the remaining vacuum-state population, P 0 f 0 f ; n ¼ 0 À Á , which is equal to the original population, P 1 f 0 f ; n ¼ 0 À Á . The uncertainty of the measurement mainly comes from the quantum projection noise of binary result of single measurements 49 .
For the experiment involving two boson modes, we first measure P 1 f 0 f ; n 1 ¼ 0 À Á using the same method as that of the single-boson case. Next we consecutively apply a "uniform red sideband" to the first mode and another "uniform red sideband" to the second mode. Then, we measure the population of the upper state, which should be P 1 f 0 f ; n 1 ¼ 0; n 2 >0 À Á . Therefore, we obtain the desired population from the relation: It is noteworthy that this scheme is clearly scalable in the number of bosonic modes.
Average boson number measurement. For the average boson number measurement, we first use optical pumping to trace out electronic states 48 and then apply a blue sideband time sweep from t = 0 to t = 12 π blue 50 . We get the phonon number distribution by fitting the result signals through the maximum likelihood method with parameters of the Fock state populations 46,47 . The main uncertainty in the average phonon measurement comes from fitting and we include one standard deviation as an uncertainty throughout the manuscript.
Ideal theoretical calculations. The exact dynamics of the Hamiltonian (7) can be obtained by solving the time-dependent Schrödinger equation i h ∂ ∂t ψðtÞ j i¼ H I ðtÞ ψðtÞ j i. We numerically solve the equation with the built-in function of Mathematica, NDSolve, which finds a numerical solution to the ordinary differential equation mainly based on Runge-Kutta method. In the numerical calculation, we include 4 internal levels and up to 10 phonons per mode for the Hamiltonian (7), which changes the Schrödinger equation to the ordinary differential equation with 40 and 400 components for single mode and two modes of the state ψ j i, respectively. With the option of infinite Maxsteps in Mathematica, the numerical calculations converge and do not show any error messages. We point out that all the parameters in the simulation are experimentally determined, not obtained via fitting. The main limitation of the numerical calculation would be the size of the Hilbert space when we scale up the system with multiple fermions and bosons.
Feynman diagram calculation. In the interaction picture, the evolution operator U I (t, t 0 ) satisfies the following differential equation which can be exactly solved as the so-called Dyson series, By introducing the time-ordering operator T , the above solution can be written in a formally succinct way Truncating at some finite N in equation (12) provides a straightforward perturbation treatment of the evolution operator U I (t, t 0 ). However, when the evolution time increases, the unitarity of the perturbation expansion becomes difficult to guarantee, because the truncation error is proportional to (t − t 0 ) N+1 . In order to deal with the long-time dynamics, we make use of the composition property of the evolution operator and interleave M − 1 equally spaced points between t 0 and t. Then, the evolution operator U I (t, t 0 ) is identically written as the product of M evolution operators, each of which governs the dynamical evolution over a short period of time, with t M ≡ t. For any dynamics with finite duration, saying t − t 0 is finite, we can always assign a sufficiently large M so that Δt ≡ t m − t m−1 is a small but finite quantity. Consequently, U I (t m , t m−1 ) is readily to be treated perturbatively.
Denote the n-th order perturbation expansion of U I (t m , t m−1 ) as U Then the whole dynamics can be treated perturbatively as follows, The deviations are related to the number of sliced sections M in time and the order of perturbations N. In our numerical calculations of the Dyson series, we divide the total time by M = 100 and apply the perturbations up to N = 9th order. Figure 4 shows the deviations of the norm from 1 1 À U I t; t 0 ð Þψ t 0 ð Þ j i j j 2 and the infidelity 1 À ψ exact ðtÞ U I ðt; t 0 Þ j j ψðt 0 Þ h i j j j j depending on the order of the perturbations for the case of Fig. 3(c). Here ψ exact (t) is the result by the ideal numerical calculation. From the 7th order, the deviations in the perturbation calculation are below 10 −4 from the ideal norm of 1.
Data availability. The data that support the findings of this study are available from the corresponding author on request.
Received: 12 January 2017; Accepted: 6 December 2017  Fig. 4 Convergence of Dyson series. In order to evaluate the validity of the perturbation calculations, we use (a) the norm of the state from 1, 1 À ψðtÞjψðtÞ h i j j j j and (b) the infidelity of the state, 1 À ψ exact ðtÞjψðtÞ h i j j j j , where ψ exact (t) is the result of the ideal numerical calculation, for the case shown in Fig. 3(c) with M = 100 divisions of time. After including up to 7th order perturbation, the deviation of the norm from 1 and the infidelity reduce to below 10 −4 . We note that for the case of Fig. 3(b), even the first and the second order perturbations provide the deviation of the norm and infidelity below 5 × 10 −2 and 2 × 10 −4 , respectively NATURE COMMUNICATIONS | DOI: 10.1038/s41467-017-02507-y ARTICLE NATURE COMMUNICATIONS | (2018) 9:195 | DOI: 10.1038/s41467-017-02507-y | www.nature.com/naturecommunications