Quantum simulation of weak-field light-matter interactions

Simulation of the interaction of light with matter, including at the few-photon level, is important for understanding the optical and optoelectronic properties of materials, and for modeling next-generation non-linear spectroscopies that use entangled light. At the few-photon level the quantum properties of the electromagnetic field must be accounted for with a quantized treatment of the field, and then such simulations quickly become intractable, especially if the matter subsystem must be modeled with a large number of degrees of freedom, as can be required to accurately capture many-body effects and quantum noise sources. Motivated by this we develop a quantum simulation framework for simulating such light-matter interactions on platforms with controllable bosonic degrees of freedom, such as vibrational modes in the trapped ion platform. The key innovation in our work is a scheme for simulating interactions with a continuum field using only a few discrete bosonic modes, which is enabled by a Green's function (response function) formalism. We develop the simulation approach, sketch how the simulation can be performed using trapped ions, and then illustrate the method with numerical examples. Our work expands the reach of quantum simulation to important light-matter interaction models and illustrates the advantages of extracting dynamical quantities such as response functions from quantum simulations.


I. INTRODUCTION
The fundamental physics of light-matter interactions was conceived more than half a century ago with the formulation of quantum electrodynamics [1,2], followed by the distillation of a theory of quantum optics [3,4].Despite this, new surprising phenomena in the realm of how light interacts with matter are still being uncovered today.This is especially true in regimes of light-matter interaction where the electromagnetic field and the material system it interacts with must both be treated quantum mechanically -i.e., where semiclassical approximations break down.Examples of recent results in this area are the revelation that a single photon can be jointly absorbed by two atoms given the right conditions [5], and the establishment of the fundamental limits and trade-offs present in building detectors for single or few photons [6][7][8][9].
There are also many applications that benefit from accurate modeling and simulation of the interaction of weak fields with complex material systems, including: (i) understanding and mimicking light absorption by photosynthetic organisms, where it is a challenge to understand the mechanisms by which the initial stages of photosynthesis can be efficient, even in the weak illumination conditions that many organisms live in [10,11], (ii) interpretation of recently developed nonlinear spectroscopies that use entangled states of few photons, which have the potential to provide unprecedented resolution of electronic, molecular, and condensed phase dynamics [12,13], and (iii) designing next-generation nanoscale engineered photodetectors that can tailor the interaction dynamics between light and matter [6][7][8][9].In such settings, where exotic light fields like single photon wavepackets interact with nanoscale structured materials and molecules, one typically has a quantum many-body model description of the physics that cannot be solved exactly and is also intractable to solve numerically on a computer.Tractable semiclassical approximations that are suitable in other settings often produce inaccurate predictions in the settings described above.This model complexity is an obstacle to many important phenomena arising from weak-field light-matter interactions.Additionally, experiments themselves can be difficult to perform due to the need to prepare exotic quantum states of light.
These obstacles provide the motivation for the quantum simulation technique we present in this paper.Quantum simulators are nascent hardware platforms that have the potential to transform the landscape of what is tractable for simulation of physics and chemistry models.This technique can be applied to any problem involving the interaction of matter with weak field states of light and can be used to study the coherent dynamics of such systems.Additionally, this technique can be generalized to model other types of non-locally and/or time-dependently coupled baths such as phononic wavepackets.We focus on analog quantum simulation platforms, where the underlying physics of the system to be simulated is encoded into the Hamiltonian of a tunable system, who's dynamics then naturally carries out the simulation task.Trapped ions are one of the leading platforms for analog quantum simulation due to the long coherence times, high degree of controllability, and scalability afforded by the platform [14].In fact, recent demonstrations using trapped ions have included simulation of models of quantum magnetism [15], of many-body localization phenomena [16], of discrete time crystals [17], and of energy transfer phenomena [18].
In this work we study how light-matter interactions between exotic weak light fields and fundamental matter degrees of freedom can be simulated using ions trapped in linear radio-frequency (RF) traps, see Fig. 1.A key aspect of the quantum simulation model we propose is that the electromagnetic field is modeled by the quantized vibrational degrees of freedom in a trapped ion system.Due to the high degree of controllability of the trapped ion system, almost arbitrary states of the vibrational modes can be engineered, thus allowing simulation of interactions with exotic electromagnetic field modes that would be difficult to prepare, especially at optical frequencies.
While electromagnetic field modes and vibrational modes are both bosonic degrees of freedom, an immediate obstacle to a simulation of the former with the latter is that while a general electromagnetic field is described by a continuum of harmonic modes, any simulation platform only contains a small, discrete number of vibrational modes.To overcome this difficulty, we formulate a novel form of simulation that proceeds via reconstruction of response functions from ensembles of quantum simulation experiments, and analyze the simulation cost in terms of the model being simulated.We also provide a detailed description and analysis of the implementation of our simulation protocol on a trapped-ion quantum simulator.
In the more general context, the approach we sketch is an example of using quantum simulation for calculation of physically relevant response functions [19][20][21][22].Although we present a formulation of the approach specific to analog trapped-ion quantum simulators, it is suitable for other hardware platforms with bosonic degrees of freedom (e.g., circuit-QED [23] or cavity-QED [24]), and there is a natural digitization of the approach, e.g., through Trotterization, which we will comment on in Sec.VI.Quantum simulation is widely believed to have potential to enable simulations that are intractable on conventional computers.Our methods introduce the simulation of weak-field light matter interactions into the quantum simulation toolbox, and as quantum simulation capabilities scale, our methods could be applied to model the interaction of weak electromagnetic fields with materials and molecules that are too complex for conventional simulation and modeling tools.
We note that there have been previous proposals to simulate vibrational degrees of freedom < l a t e x i t s h a 1 _ b a s e 6 4 = " 1 a B u S v 0 h u 9

()
FIG. 1.We develop a quantum simulation approach that enable platforms with controllable bosonic degrees of freedom, such as the trapped ion system shown on the right, to simulate the dynamics of matter, such as the chromophoric complex on the left, in response to illumination by extremely weak light fields.Such dynamics is relevant for understanding the behavior of photodetectors, photovoltaics, biochemical systems such as photosynthetic light harvesting complexes, and for modeling new weak-field spectroscopy experiments.
in molecules with bosonic degrees of freedom in trapped-ion systems [25] and to utilize quantum optical networks to emulate and sample vibronic spectra of molecules [26,27].To our knowledge, our proposal is the first to show that simulation of quantized states of electromagnetic fields is also possible with trapped-ion bosonic modes.
The remainder of the paper is structured as follows.In Sec.II we present a brief overview of trapped-ion physics and the type of Hamiltonian models that can be engineered on this platform.
Sec. III presents the light-matter interaction and dynamics that we wish to simulate, and then in Sec.IV we detail our response function approach for simulating these dynamics on a platform like trapped ions with a fixed number of stationary bosonic modes.Sec.V presents examples illustrating the scheme, with numerical simulations.Then in Sec.VI we discuss important considerations when scaling the proposed simulation scheme to large systems, and in Sec.VII we discuss application of our response function approach to quantum simulation to models beyond light-matter interactions.Finally, Sec.VIII concludes with a summary of contributions made in this work.

II. TRAPPED ION HAMILTONIANS
The trapped-ion platform allows for realization of a rich set of quantum models with localized degrees of freedom, encoded in internal states of the ions, and distributed bosonic degrees of freedom, encoded in quantized motional modes of ion motion.In the following, we briefly review the building blocks of quantum simulation using trapped ions, with a particular focus on the physics and achievable Hamiltonian models that are relevant to our setting.We focus on linear ions chains of Ca + ions for concreteness, although neither is an important restriction.We emphasize that this is a narrow review and refer the reader to Refs.[14,28,29] for more comprehensive treatments.
We encode the ground state and optically connected excited state of point-like absorbers (e.g., atoms, molecules) in two Zeeman sublevels of a Ca + ion.In particular, a commonly used sublevel 2 ) and |e → D 5/2 (m J = 1 2 ) , that uses states in the stable S 1/2 orbital and the metastable D 5/2 orbital [18].These states will form the eigenbasis for localized degrees of freedom and in the following, when we write Pauli operators acting on localized degrees of freedom, they will be with respect to this basis.Note that it is possible to consider point-like absorbers with more than one excited state by encoding into other states in Zeeman sublevels of the ion, but we will restrict attention to single excited states here.
An ion chain with n ions has 3n quantized motional modes (corresponding to collective motions of the n ions under the influence of the overall trap potential and their mutual Coulomb repulsion) that can serve as modes of the electromagnetic field.This is less than the continuum required to accurately model electromagnetic (EM) fields but in the following we will construct a scheme for overcoming this obstacle.All of the motional modes can in principle be coupled to the internal degrees of freedom of the ions through laser-induced interactions.The basic interaction that enables most of trapped ion quantum simulation is given by the following Hamiltonian that describes a single ion interacting with a laser field nearly resonant with the energy difference of the |e and |g states of the ion ( = 1 here and in the rest of the paper) [30]: where ∆ is the detuning of the laser from the atomic transition, ϕ is the phase of the laser with respect to the atomic polarization, Ω is the amplitude of the laser field (expressed in terms of the Rabi frequency), ω T is the trap frequency.σ+ is the raising operator for the internal degrees of freedom, and â is the annihilation operator for the quantized motion of the ion around its equilibrium position along, for example, the z axis.η ≡ k z 1 2mω T is the Lamb-Dicke parameter, with k z being the projection of the laser's wavevector along the z direction.In the Lamb-Dicke limit, which describes the regime of small light-induced changes in momenta, or η (â + â † ) 2  1, this Hamiltonian can be approximated as [30] H SI (t) ≈ Ω σ+ e −i(∆t−ϕ) + σ− e i(∆t−ϕ) Then, through choices of ∆ and ϕ one can engineer a wide variety of ion-mode interactions.The above description of the dynamics of a single ion generalizes easily to multiple ions, with the motion now being interpreted as modes of the collective motion of the ion chain.
In addition to the ion-mode interactions, we can also mediate interactions between the internal states of ions through their common interaction with motional modes.There are a few schemes for engineering such interactions, and the commonly used Mølmer-Sørenson scheme -which illuminates two ions with a bichromatic laser with frequencies ω 0 ± (ν + δ), where ω 0 is the energy difference between |g and |e , ν is the frequency of a vibrational mode (usually an axial mode), and δ is a detuning -generates an effective interaction between two ions of the form σx Putting these ingredients together we write the family of Hamiltonians that can be engineered with trapped ions -which we denote with the superscript TI -as Here, ĤTI a is the Hamiltonian for the internal states of the n ions, representing localized degrees of freedom.ω j 0 is the local transition energy and can be tuned by AC Stark shifting the magnetic sublevels of each ion with laser beams that locally address the ions.The second term represents coupling between some subset of ions and are implemented via the Mølmer-Sørenson interaction as mentioned above.By tuning the Mølmer-Sørenson interaction this coupling between ions can be tuned to the rotating wave version where only the terms σ+ The magnitude of the couplings, J ij , can be tuned via the intensities and detunings of the laser beams implementing these interactions.We note that engineering fully tunable J ij interaction terms through the Mølmer-Sørenson (MS) interaction requires coupling to one motional mode per interaction term, which must be accounted for when accounting for quantum simulation resources.ĤTI m is the Hamiltonian representing the m ≤ 3n motional modes used in the simulation.The harmonic frequencies ν j are set by the trapping potential and number of ions in the trap but the effective mode frequencies that the ions interact with can often be tuned via how the interactions are engineered [18].ĤTI I represents the interactions between the ions and motional modes, and accounts for two types of interactions.The first sum represents a coherent exchange of energy and is implemented via the interaction discussed in Eq. (1) (by setting ∆ = −ω T , or a red sideband drive).The second represents a shift in the energy of the ions that depends on the position of the motional mode and can be implemented via a bichromatic local addressing of the ion involved [18].The index α in these sums represent the number of interactions, and j α , k α index the ions and modes involved in an interaction.The interactions parameters, κ α and λ α , can be tuned by adjusting the intensity of the local addressing laser beams.Importantly, one ion cannot participate in both types of interactions.We note that other types of Hamiltonian terms are also possible but the family of Hamiltonians described above will be sufficient for our purposes.
Having described the types of Hamiltonians that can be engineered on the trapped-ion platform we will now summarize the other two requirements for a quantum simulation, state preparation and measurement.There are a variety of mechanisms in trapped ion platforms for extracting entropy and preparing desired states of the internal and motional degrees of freedom with high fidelity.The internal states of the ions can be prepared in the |g state through optical pumping [31].The states of the motional modes can be prepared in thermal states, including very low temperature thermal states, through a variety of cooling mechanisms.The most common such cooling mechanism, resolved sideband cooling, proceeds by coupling a mode to an ion using the same interaction as shown in Eq. ( 1) and cooling the mode by dumping excess energy into the ion, which emits into an optical mode usually controlled by coupling the long-lived excited state of the ion to a fast decaying excited state [31].Each mode that is used in the simulation can be cooled sequentially.If many modes need to be initialized cooling via electromagnetically induced transparency might be an attractive alternative.We note that other cooling mechanisms (e.g., Doppler cooling) must first be applied to the ion chain before the ion motion is sufficiently cold for resolved sideband cooling to be effective.Measurement of the ion internal degrees of freedom is accomplished via the electron shelving technique, which selectively excites one of the encoding states (usually |e ) to a short-lived higher lying state whose emission is monitored.
In the following sections we will develop a scheme for simulating light-matter interactions using the quantum simulation building blocks discussed above.

III. INTERACTION OF MATTER WITH WEAK FIELDS
Matter interacting with light must be considered as interacting with all modes allowed by the confining geometry, so that the Hamiltonian of an arbitrary system in the rotating-wave approximation can be written as where ĤM is an arbitrary Hamiltonian for the matter subsystem, âj is the annihilation operator for allowed mode j, and L describes the action of the light-matter interaction on the matter system (e.g., promoting a two-level system from the ground state to an excited state and vice versa, in which case L = σ− ).In the event that the geometry of the space containing the matter is unconfined; i.e., the matter exists in free-space in at least one direction, then the allowed modes form a continuum.For the 1D case (for simplicity and clarity) modes in this continuum can be represented by their frequency ω and we can write Solution of this Hamiltonian using a Markov approximation results in the emergence of the phenomenon of spontaneous emission, as first identified by Wigner and Weisskopf, as a direct consequence of the availability of a continuum of emission channels.
Consequently, this Hamiltonian cannot be straightforwardly described by one, or even a few, bosonic modes, as in the case of ion trap system of Sec.II.Nonetheless, as we will now show, it is possible to use the ion trap system to simulate light-matter interactions of the former kind.

IV. SIMULATION THROUGH RESPONSE FUNCTIONS
First, we will rewrite Eq. ( 6) according to the quantum noise formalism.This can be done by constructing new field operators.With the field only in the interaction picture, where ω 0 represents a frequency that will be associated with the incoming pulse.The new operators behave like noise operators, and obey the same statistics, such that [32] [ dB t , dB t † ] = δ t,t dt A monochromatic single photon pulse with frequency ω 0 traversing the matter system with a broad temporal profile ε(t) is constructed as where |0 is the field vacuum.
The form of the transformed Hamiltonian, ĤM−F (t), makes it seem like the system is only coupled to a single field mode.However, the behavior of the dB t operators will result in different dynamics.In particular, due to Eq. ( 7) the interaction Hamiltonian provides an additional, nonunitary term at linear order that must be included when propagating the dynamics; to see this explicitly, note that the time evolution operator for a short interval is Due to the nonunitary contribution we will write the evolution in terms of the density matrix We note that incoherent processes (e.g., due to couplings to environments other than the EM field) can be included as additional terms, e.g., of the Lindblad form X † X, ρ , however we will omit these presently for clarity.We will also "vectorize" the density matrix, performing the transformations ρij = This allows us to express superoperators acting on ρ, such as commutation and the Linbladian D, as linear operators acting on ρ.We can now integrate Eq. ( 9) directly.Let Then Here, the interaction between the system and field (represented by F dB (t)) is treated perturbatively and each term in this expansion corresponds to an interaction of a certain order.For a system with the matter subsystem initially in the ground state and the field in a single photon state, i.e., this series truncates and we get with L+ (t, t ) = Ḡ0 (t − t )ε(t )e −iω 0 t L † 1 for compactness in order to emphasize the structure (note that by construction L † 1 † = 1 L † ).The reason for the early truncation of this series is that once an absorbed photon is re-emitted into the continuum field it cannot interact with the matter system again; i.e., d Bt † |0 ε = |1 ε1 .This is the essence of spontaneous emission, given by the third term of the sum.
We consider now the case where the system interacts with a single stationary bosonic-mode as in the trapped-ion context.The Hamiltonian, in an interaction picture with respect to the mode's free Hamiltonian, is: where the superscript "sm" on this quantity and subsequent quantities indicates that a single mode is modeling the EM field.This is simply related to the general trapped ion Hamiltonian described in Eq. ( 2); ĤM is the matter Hamiltonian engineered through ĤTI a in Eq. ( 2), we have restricted to one mode, so k α is the same for all α, and L = α κ α σjα − is the engineered interaction between the matter subsystem (modeled by the ions) and the mode.Finally, we have set λ α = 0.In addition, for reasons that will become clear shortly, we will include a Markovian decoherence term to the dynamics of the form D[ Ŷ ], where Ŷ is an operator on the internal states of the matter subsystem.
Given this setup the unitary contribution to the infinitesimal time evolution is governed by Then, with the dissipative contribution included, and so that the dynamics can be expanded in the light-matter interaction for ρ(t with L expressions defined similarly to the continuum mode case (e.g., Lsm − (t, t ) = Ḡsm 0 (t − t )γ(t )e iω 0 t L 1 ).We can now compare Eq.( 11) and Eq.( 15), and find two essential differences.
Firstly, the series in the continuum mode case, Eq. ( 11), cuts off after (spontaneous) emission, but continues with an infinite number of coherent absorption/emission cycles in the single mode case, Eq. (15).Unlike in the continuum case, an excitation emitted into the stationary mode can be reabsorbed by the matter system.A second, related, difference is that in the continuum case, the light-matter interaction leads to a decoherent effect captured by the real terms in the argument of Ḡ0 and the second term in the definition of F dB (t), whereas in the single mode case decoherent effects are generated by the added term Ŷ .Thus, in order to reproduce the continuum mode case with the single mode system, we must suppress the additional terms of Eq. ( 15) and set Ŷ = L. Unfortunately, neither of these tasks is straightforward to carry out.For the first, the absorption/reemission dynamics of the single mode case are inherent in the system.For the second, since Ŷ is engineered via an off-resonant drive between the relevant internal states of the ions and higher lying state(s) with short lifetimes, or through other dissipative mechanisms [33], it must be a local operator, unlike L, which couples multiple ions to the field.We shall deal with both of these issues in turn, and show how simulation of relevant continuum dynamics can still be achieved with the single mode system.

A. Moving to response functions
A way to overcome the difficulty of not having a continuum of modes and the resulting phenomenon of spontaneous emission is to take advantage of the time dependence of the coupling, γ(t), in Eq. ( 14) and move to a response function framework.
First we note that, the population of an excited matter state i M at time t, which is often the quantity of interest, is given by the second term of Eq. ( 11) after tracing over the field state: where so that the dynamical properties are captured by this two-time Green's function.We write such Green's functions as G i M (t 1 , t 2 ), where the subscript i denotes the state whose population is being measured.If this quantity is known, then the dynamics in response to any wavepacket profile shape ε(t) can be determined.The problem can then be reduced to finding the response to delta-like inputs.
To simulate the delta function response with just a single stationary mode, we can consider modulating the field coupling to be nonzero only briefly and at given times.Consider a square "pulse" of width t γ and area n γ , so that γ(t) We note that these pulses do not need to be strictly square in shape; we assume as such for illustrative purposes.
When Ḡ0 (t) can be approximated as constant over the interval t γ , we obtain for a pair of pulses at t and t ).
We can see that when || L|| 2 n 2 γ 1 the terms to linear order dominate.Setting Ŷ = L for now, and comparing to Eq. ( 17), letting t m ≡ t − t and t int ≡ t − t we find that Therefore, the Green's function relevant to the multimode case, Eq. ( 17), can be extracted from the single mode system.We note that G i M (t m , 0) can be found by using a single pulse.Thus by choosing a proper dissipative process ( Ŷ ) and appropriately short pulses, we can reconstruct G i M and simulate the matter response from the continuum mode case.

B. Sampling error
In practice P M sm must be estimated from multiple trials resulting in a measurement of state i, which will either be populated or not.Thus the statistics governing the measured P M sm will be that of a binomial distribution and the sampling error in the estimate of P M sm is given by where N is the number of trials.

Then the sampling error in
Note that the impact of larger n γ is to reduce the sampling error estimates σ G by an overall factor of approximately 1/n γ ; this must be weighed against the error present in the underlying P M sm (and any additional experimental sources of error), with the optimal value being case dependent.We will return to this point later on when considering example systems.
Further, can develop an error bound on the population esimates by discretizing the integral in Eq. (16).Consider the response to a particular wavepacket with temporal profile ε(t), the overall error in the simulated population is where t i are times spaced ∆t apart, and the first step above assumes that ∆t is small enough that the integral in Eq. ( 16) is well approximated by the discrete sum.It is worth examining the case where the error in G i M (t m , t int ) can be taken as roughly the same for all t m , t int ; i.e., σ G (t m , t int ) ≈ σ G .In that case, the normalization of the wavepacket means that Since ∆t ∝ 1/n int and σ G ∝ 1/ √ N , the total number of experiments needed to obtain a given error in the wavepacket response can be actually be minimized in principle by taking a larger number of intervals and rather than a larger number of samples per interval.

C. Simulating dissipation
In the previous subsection, we showed that the response of the matter system to single photon pulses could be recovered with a trapped ion simulator with a modulated ion-vibration coupling, and where a vibrational mode simulates a temporal mode of the photon field.However, to do this we had to introduce an additional Markovian dissipative process D[ Ŷ ] acting on the ion internal states, and set Ŷ = L.We now examine challenges posed by this term.
As described above, additional dissipative processes can be introduced into a trapped ion platform (e.g., through optical pumping), but these act locally and independently on each ion.However, the L operator is non-local operator acting on all excited states in the model and thus, in a correlated manner on several ions -i.e., if the optically coupled levels in ĤM are encoded in the excited states of N ions, L = N j l j σ− j .This mismatch makes setting Ŷ = L non-trivial.The most straightforward solution is to transform the basis so that L is local, however, operators associated with any additional decoherence processes will also be transformed, and may become nonlocal.In the general case there will not be a basis in which all operators are local; additionally, the state to be measured may become a superposition of several ionic states, compounding the difficulty.To overcome these difficulties, we develop a method for simulating non-local decay-type processes of the same form as L using tailored couplings to an auxiliary ion that undergoes a local, fast decay.
We introduce an auxiliary ion and isolate two internal states within it governed by a two-level system with Hamiltonian Ĥaux = ωaux 2 ( 1 − σz aux ).Then we assume that the excited state is subject to a fast decay process: D[ X], with X = χσ − aux .This could be engineered through optical pumping to a higher lying state, for example [33].In Appendix A we show that by coherently coupling the system ions to this fast decaying auxiliary ion with a Hamiltonian of the form and adiabatically eliminating the auxiliary ion, results in a Markovian dissipative process on the system ions that takes the form D[ 2 Ĵ χ ].Therefore, we can engineer the non-local decay terms necessary to model the dissipative effect of coupling to a continuum of modes by tuning the coherent couplings in Ĵ such that 2 Ĵ χ = L. Thus, using a single bosonic mode and an auxiliary ion, we can simulate the response to a single-photon traversing the system in a free-field setting.
A graphical representation of the actual and simulated process is depicted in Fig. 2.
In practice, this engineered non-local dissipation, which is the "desired" relaxation dynamics, will be in competition with local, single qubit relaxation processes.For trapped ions intrinsic relaxation occurs on very long timescales and therefore || L|| is likely to be a much faster rate than local relaxation rates.However, in the unlikely event that single qubit relaxation rates dominate the non-local desired relaxation prescribed by L, then all simulation parameters should be scaled up such that the slowest rate is faster than single qubit relaxation and dephasing rates.

D. Experimental protocol
Putting together the ingredients from the previous subsections, we can specify an experimental protocol to follow to measure the response functions to weak-fields and estimate the response to an arbitrary wavepacket using a quantum simulator (Fig. 2(a)): 1. Initialize "matter" system in its ground state, the auxiliary system in the ground state, and the bosonic mode in the first excited state with the coupling set to zero.18).For an arbitrary wavepacket profile ε(t), the response at a time t is simulated by multiplying Gi M (tm, tint) by ε(t − tm − tint)ε(t − tm) and summing over relevant tm, tint (Eq.( 16)).(b) A depiction of a system interacting with a field characterized by a continuum with a single-photon wavepacket (red oscillatory lineshape), traversing the system.At early times the matter system (black circle), represented here as a two-level system with a decay mode (red lines), interacts only with the vacuum and remains in the ground state.At some point during the traversal of the wavepacket the matter system may absorb the photon and enter the excited state, leaving the field in the vacuum state.Later, the matter system may emit a photon back into the field, decaying to the ground state; afterwards the photon is not present at the matter system and may not be reabsorbed.(c) A system comprising a single mode for the field (red circle) whose coupling can be rapidly adjusted with time and a highly incoherent bath (dashed line circle) coherently coupled (gray double arrow) to the matter system (black circle).Initially the boson mode is in the first excited state and uncoupled to the matter system.At some time the interaction (red double arrow) is rapidly turned on and off, allowing for possible transfer of the excitation.At some point later this population will be transferred to the bath and almost immediately decay, mimicking a spontaneous emission event.Shown is the response for a single coupling pulse; as discussed in the text, in general multiple pulses are necessary to capture the impact of interference between excited matter states created by interaction with an external field at different times.
2. Modulate coupling with the profile γ(t), providing a "pulse" at time −t int .
3. At time 0 modulate the coupling again to provide a second "pulse".
4. At time t m measure the relevant population ρ ii of the matter system's state to obtain a sample of P M sm (t m , t int ).
5. Perform 1-4 for different t int , t m , repeat N times (as needed for the desired statistical error on the estimates, see Eq. ( 19)), to generate the Green's function G i (t m , t int ).
This protocol is resolving the two-dimensional function G i (t m , t int ), and the total number of experiments required is a function of the resolution required in t m and t int .The maximum required resolution is set by the fastest timescale in the system evolution (R ≡ ), which sets the required sampling resolution and discretization of t m and t int , while the optical decay rate (|γ|), sets the maximum time interval over which there is interesting dynamics.These are conservative estimates, especially the temporal resolution of 1/R, since the timescale of system evolution is typically much slower than R. We will see an example of this in Sec.V.

E. Multiple photons
For the case of n photons, the above scheme can be easily (if tediously) generalized by considering an initial bosonic mode in its nth excited state with coupling pulses at 2n different times, leading to a 2n-time Green's function.To isolate the desired Green's function, a similar expression to Eq. ( 18) must be used so that "duplicate" same-time terms can be subtracted out.For n = 2 where the . . .indicate the additional same-time terms (e.g., G i M (t m , 0, t int , 0)) omitted for brevity; in general there will be (2n)!− 1 of these.
The remaining unknowns are the settings of simulation parameters -the shape of pulse profiles, parameters of the auxiliary system, and sampling of time intervals -required to reach the necessary limits for desired accuracy.We will now consider a small example system to analyze the dependence on these parameters, as well as a larger example system akin to the kind this method is intended for.
V. ILLUSTRATIONS

A. Example 1: coupled chromophores
To illustrate the ideas presented in the previous sections and gain intuition, we consider the simplest example of a composite matter system interacting with a weak-field wavepacket.Consider two chromophores, each with optically active excited states interacting with a single-photon Gaussian wavepacket, see Fig. 3(a).The Hamiltonian and L operator for this example is given by: , where the states are ordered as the common ground (HOMO) state of both molecules with zero energy, the excited (LUMO) states of molecule 1 with energy ω 1 , and the excited state of molecule 2 with energy ω 2 .Note that we have assumed there is no direct Coulomb coupling between the molecules for simplicity.We also restrict the above operators to the single excitation manifold since that is all that is needed to model the interaction with a single photon.In the multiphoton case we would need to model the multi-excitation energy states also.
Suppose we want to determine the population of state 1 in response to single photon pulses at the resonance frequency ω 1 .The two-time Green's function in Eq.( 17) for this system can be written We construct a trapped ion simulation of this two chromophore system with two ions, each with the following Hamiltonian and L operators: and bosonic mode and auxiliary qubit for the field are described by The couplings are There are three parameters that must be adjusted to ensure that the simulation operates in the appropriate limit for simulating G 1 : n γ ,χ, and t γ .Additionally, the times t m and t int for which G 1 are simulated must be chosen from a sufficiently dense grid to ensure the pulse profiles of interest can be accurately reconstructed.In Fig. 4 we show the dynamics for exact and simulation systems, given ω 1 = 1.0, l 2 1 = 0.0036, ω 2 = 0.8, and l 2 2 = 0.0064, as well as the impact of n γ , and χ on the population response to two pulses used to generate G 1 (t m , 250).Fig. 4(a) shows the exact result for a single t int .The main features are the step changes at the pulse times, the long decay due to spontaneous emission, and oscillations at frequency ω 1 − ω 2 due to the different energies of states 1 and 2. We emphasize that, as given by Eq. ( 18), in order to obtain G 1 (t, 250) we must subtract out the t int = 0 from the overall response to the two pulses.Shown in Fig. 4(b), in this case the necessary parameter values required to reach the asymptotic limit.χ 2 must essentially be about two to three orders of magnitude faster than the rates l 2 i (which are relatively small).To understand the impact of n γ , recall that this parameter corresponds to the area of the pulses of coupling between the ions and the mode.Therefore, larger n γ means greater population exchange between the mode and internal states of the ion, and thus increased reabsorption.This explains the deviation of the computed response and true response with increasing n γ .n γ ∼ 1 (resulting in the expected population shown in Fig. 4(a)) is sufficient for most purposes.
The impact of different choices of t γ is more nuanced.While a larger t γ , which results in wider pulses, might seem to result in a poor approximation of a delta impulse (which is what the pulse mimics), the fact that the Green's function is convolved with a smooth wavepacket profile results in considerable tolerance to the value of t γ .To demonstrate this, consider Fig. 5(a), which shows the G 1 computed for different pulse widths.For narrow pulses, oscillations due to the difference ω 1 −ω 2 are quite pronounced, while for wider pulses these are substantially reduced.However, since according to Eq. ( 17) these computed Green's function are further convolved with the wavepacket profile, if the wavepacket is slowly varying in time the impact on the final result of larger t γ is minimal until it begins to exceed the wavepacket's temporal width.This is demonstrated in Fig. 5(b) for a Gaussian wavepacket with σ t = 100.The inset shows the exact response P M (t), while the contour plot shows the deviation from the exact value for P M (200) for different values of t γ and sampling interval ∆t of t m and t int .The computed response is robust to fairly large t γ ; in fact, larger values are somewhat preferable if one wishes to minimize the impact of choice of the ∆t sampling interval, which exhibits pathologies at particular values.These are at multiples of 2π/(ω 1 − ω 2 ) = 31.412and are due aliasing from sampling intervals that align with the period of oscillation of G 1 .
Note that the timescale of system evolution in this example is 2π/(ω 1 −ω 2 ), which is significantly larger than 1  ||H M || since the coherent dynamics is only between the excited states in the system, and the optical decay rate is much slower than the coherent timescales.This means that the resolution required of the two-dimensional function G 1 (t m , t int ) is naïvely ∆t < 2π/(ω 1 − ω 2 ), and indeed, the error in Fig. 5(b) (for constant t γ ) is smallest for ∆t in this range.However, note that the error is also small for larger ∆t values as long as the aforementioned pathological values are avoided.This small error at larger ∆t values results from a cancellation of errors due to the integral in Eq. ( 16); i.e., the fast oscillations of G 1 and slow variation of (t) imply that the integral evaluates to almost zero, even if G 1 is coarsely sampled.This is specific to this system, and not generically expected.

B. Example 2: energy transfer with non-Markovian environment
The system considered in the last section is of course trivial to simulate conventionally.Even in more complicated cases, single photon response can often be efficiently simulated by restricting the calculation to an energetically accessible subspace.However, this is not always especially helpful, as in many systems -especially those featuring many body interactions -the single (or few) photon accessible manifold can nonetheless be quite large and involve subsystems where multiple energy states must necessarily be included.For example, an important consideration and phenomenon of interest is the dynamics and behavior of such systems under the influence of non-Markovian (e.g., 1/f ) noise due to impurities such as charge traps at interfaces.These can be modeled by adding auxiliary -frequently two-level -systems driven by Markovian noise coupled to the main system of interest.Scaling due to these additional degrees of freedom added by considering explicitly structured noise sources cannot be reduced, quickly rendering conventional simulation intractable.The present scheme, however, is not impacted by such considerations; as long as the necessary elements can be incorporated into the experimental setup, the number of samples required experiences constant scaling.
To illustrate this, we now consider a model that is representative of the kind of complex mate- rial whose response one might want to calculate with our approach.We model an optically active molecule or quantum dot that is Coulomb coupled to several optically inactive molecules/dots.
The system is also coupled to environmental fluctuators that cause non-Markovian decoherence of the system.We are interested in the energy transfer dynamics induced by single photon illumination in such a system and thus monitor the population at a optically inactive site.Due to the non-Markovian dynamics, simple reduced treatments of the optically activated dynamics are not possible.This minimal model is representative of the structures found in biochemical molecular complexes or quantum dot nanostructures.In Fig. 6(a) we depict such a system: an optically active subsystem on the left, coherently coupled to two dark systems with an incoherent decay to a stable state at the end, mimicking an electron transport pathway.The middle state of the coherent chain is coupled to a pair two-level fluctuators acting as a non-Markovian noise source.
The matter Hamiltonian for this example is given by The field coupling is only to the bright dot and therefore, Additionally, the system is driven by two incoherent processes, generated by Lindblad terms D[ Ẑ], with the following generators The first of these represents the incoherent transfer from the third (dark) dot to a final stable state f on the last dot in the chain (which is monitored by the electron transfer pathway/measurement), and the second represents the incoherent flipping of the fluctuator state.The reduced dynamics of the quantum dots is non-Markovian and complex due to the coupling to the fluctuators and the EM field.
This system is simulated with the trapped-ion platform using the scheme shown in Fig. 6(b).
Individual ions represent the active and dark quantum dots and the fluctuators.In addition, a vibrational mode represents the EM field and we introduce an auxiliary ion to model the field damping as described above.The terms in the matter Hamiltonian ĤM are simulated using either local Stark shifts to ion energies or through Mølmer-Sørenson interactions in the rotating wave regime (for the coherent interactions in Ĥij (Eq.( 25)).The coupling between the quantum dot 2 and the fluctuators, Ĥ2j in Eq. ( 26), is not of any of the standard forms discussed in Sec.II.
It describes a shift in energy of ion 2, that depends on the state of the fluctuators.This type of interaction, although not as popular in trapped-ion quantum computing as the Mølmer-Sørenson interaction, can be engineered and has been demonstrated [34,35].Therefore, this example also requires such σ z ⊗ σ z interactions between ion 2 and the ions modeling the fluctuators.
The incoherent dynamics generated by Eq. ( 29) are generated by driving each of the fluctuator ions with microwave fields superimposed with broadband noise [36].Finally, the incoherent coupling to the stable state f , generated by Eq. ( 28), is achieved by encoding this state as a metastable hyperfine state of ion 3.Then, by optically pumping the excited state of ion 3 to a high-lying state that decays into the state f, we achieve incoherent coupling desired.The strength of the coupling can be tuned by the detuning of the optical pump.
As with Example 1, the spontaneous emission (decay) from the optically active quantum dot is modeled through a coupling to an auxiliary ion that is damped.The coupling between the ion 1 and this auxiliary ion, and the vibrational mode that models the optical field, and the internal Hamiltonians for both of these systems take the form: where L, as in Eq. ( 27), is a local operator acting on ion 1.As described in previous sections, the auxiliary ion is optically pumped in order to simulate dissipation.
In Fig. 6(c,d) we show the response characteristics for a measurement of the stable state f following the passage of a Gaussian photon pulse of the same kind considered above.In Fig. 6(d) the function G f (∞, t int ) is plotted.The decay processes result in a sharper response with respect to the interval, reducing the overall window over which G f is significant and must be computed.In Fig. 6(c) we show the expected population after the simulated pulse has passed and the population of f has stabilized for varying t γ and n γ .As before, the final result is robust for rather large choices of t γ and ∆t int .Interestingly, the longer pulse widths act to mitigate this error; by spreading out the excitation over a longer period, the incoherent process are given more time to damp out the excitation, limiting accumulation in the 1 state that can be coherently transferred back to the bosonic mode, resulting in less error as n γ is increased.
Finally, we consider the number of trials N necessary to estimate G.For both cases the values of P M sm are around 0.01.This suggests that, in order to obtain σ G < 0.001 one should take N ∼ 10 4 samples.If we consider a sampling interval for t m , t int of ∆t = 10 for a total time domain of 1000, then according to Eq. ( 20) an error of 0.01 in the population corresponds to around one million experiments.As mentioned previously, depending on the experiment objectives and the system in question, there is some freedom to "reallocate" trials from the sampling of G to the sampling of t m and t int , reducing the overall number of experimental runs needed.

VI. SCALING CONSIDERATIONS
Now we address the scalability of the proposed approach to quantum simulation of weak-field light-matter interactions.Ultimately, the goal is to simulate interactions with material systems with hundreds or thousands of localized degrees of freedom (e.g., atoms).This corresponds to hundred or thousands of qubits since each localized degree of freedom requires at least one qubit to model.Furthermore, one requires a constant number of auxiliary qubits (or other degrees-offreedom) to capture the effects of spontaneous emission.While it is an engineering and technical challenge to scale the trapped-ion (or any other quantum computing) platform to such sizes, it is within the roadmap of the technology.The aspect of the simulation protocol that requires more thought is the scalability of engineering all the interactions necessary for the general model, given by the Hamiltonian (in the rotating frame with respect to mode â's frequency ω 0 ): where we have assumed the matter Hamiltonian H M to be of the form that can be engineered in trapped ions (i.e., H a in Eq. ( 3)).In addition to this Hamiltonian, fast dissipation at rate χ 2 must be engineered on the auxiliary qubits.
In the following, we discuss the primary concerns when implementing this model at scale on the trapped ion platform.
1. We begin with the need to engineer a fast decay process on one or several auxiliary ions.As mentioned above, this can be accomplished by optical pumping of an ion, i.e. driving the ions coherently with laser light followed by a spontaneous decay.The spontaneous decay process, however, generates a recoil, i.e. there is a probability of η 2 α (n α + 1) that a phonon in mode α is generated.While this probability is small for small motional mode occupations n α and typcial Lamb-Dicke parameters η α of a few percent.Still if a large number of spontaneous emission process needs to be simulated, this heating process needs to be countered.Probably the most promising avenue to accomplish this is to employ another ion species or isotope which allows cooling using light not resonant with any of the transitions of the primary ion species.Thus, the temperature of the ion crystal can be controlled while preserving the coherence of the ions.
2. As mentioned in Sec.II, one also must be careful about accounting for the number of motional modes required.Not only do we need a mode for the electromagnetic field (with annihilation operator â), but each qubit-qubit interaction must be mediated by a separate motional mode.While the most general model can contain n(n − 1)/2 non-zero J ij coupling values, in physically relevant models, the coupling structure is more sparse.For example, if these couplings arise from screened Coulomb interactions, as is common in biochemical models, the coupling values decay rapidly with separation distance and consequently, in a complex of many constituents (e.g., molecules) only the couplings between the closest ones have to be taken into account to get an accurate picture of the dynamics.Suppose there are m < n(n − 1)/2 such non-zero coupling terms.Finally, implementing all the qubitqubit interactions in H M −aux requires p modes, where p is the number of optically active subsystems in the model (the number of subsystems/qubits that enter the definition of the optical coupling L).Therefore, in sum, we need m + p + 1 motional modes to implement this model.A chain of n ions has 3n total vibrational modes, and the allocation of these modes to mediate the interactions in a model will be a complex optimization that depends on hardware and trap constraints and the model.However, the minimum necessary condition that we need to satisfy is m + p + 1 3n.
3. If the number of motional modes required by the model is feasible, the next concern is the complexity of having all interactions turned on at once to implement the model.Even assuming the technical challenges of having so many addressing laser beams can be met, one could be worried about interference or crosstalk effects between terms.Not only do many ion-mode interactions have to be on at once, but the same ion needs to interact with multiple modes at once.The main concern is that light coupling a particular ion to a particular mode with strength κ will also couple this very same ion off-resonantly to other motional modes α detuned by ∆ α with strength n α (κ/∆ α ) 2 .Anticipating a typical coupling strength of order kHz and aiming at a coupling of less than 1 % at nearby modes, we require ∆ α n α ×10 kHz.
In order to mitigate crosstalk from strongly coupled modes to weakly coupled ones, it will be advisable to group them such that weakly coupled modes have large frequency differences to the strongly coupled ones.Typical motional frequencies of the transverse modes of an ion strings are of order 5 MHz spanning about 1 MHz.Shaping the axial potential such that those mode frequencies are uniformly spread, there is room to control the coupling to close to 100 modes.If less crosstalk should be desired or the ion crystal is not near the motional ground state, all interaction strengths could be reduced to slow down the simulation at the expense becoming more sensitive to decoherence such as motional heating and qubit decoherence, both of which can be larger than 100 ms.
The ultimate path forward to scaling this approach might be to go beyond analog simulation and instead use a digitized model that performs Trotterization of the above model.This would enable one to implement each interaction between subsystems/qubits separately, thereby allowing one to recycle the modes and thus reducing the required number.Although digital simulation offers many advantages, especially for scaling up the approach, it is important to realize that digitized evolution can take longer, and thus require a platform with better operation fidelities and coherence the electromagnetic field did not have be bosonic since we only considered one photon states; such simulations could be performed solely by encoding in internal state of ions.However, for multiphoton interactions, the bosonic/harmonic nature of the vibrational modes becomes essential.In this work we have focused on the trapped-ion platform for concreteness.However, any quantum simulation platform with controllable bosonic modes with tunable coupling to other localized (e.g., qubit) degrees of freedom could implement our scheme.
T y o 6 p K 1 T L t m N + Q e L w P u 8 n Z i 3 n 8 c n 7 0 7 m h 2 f D L s y j 5 5 Q p 6 S 5 y Q m r 8 g x e U t O y Y I w U p I v 5 C v 5 F n w O v g c / g p 9 b a T A Z P I / J l Q h + / Q U o 5 / + X < / l a t e x i t >

FIG. 2 .
FIG. 2. (a)Example pulse/measurement schedules (top) and a sketch of their inclusion in the response reconstruction (bottom).The pulse profiles γ(t) (red) are chosen with spacing tint and measurements are performed at tm after the second pulse, allowing for the generation of samples of GiM according to Eq.(18).For an arbitrary wavepacket profile

FIG. 3 .
FIG. 3. Schematic for the example illustrated in Sec.V. (a) The physical scenario being simulated.Two chromophores represented by black circles are excited by a weak-field wavepacket.Each chromophore has a ground state, corresponding to a lowest unoccupied molecular orbital (LUMO), and an excited state, corresponding to the highest occupied molecular orbital (HOMO); these are drawn in red and cyan to indicate optical coupling and are shown with wavy arrows to indicate the optically mediated decay process.(b) Schematic of trapped ion simulator capable of simulating the physical setup in (a).The internal states of each chromophore are encoded into states of distinct ions (indicated by the dashed circles).The excited states no longer have internal decay processes, but are instead coherently coupled to an auxiliary ion that is damped.In addition, each of ions encoding chromophore states is coupled to a vibrational mode representing the field (red circle).The monitored state of the red chromophore, encoded in the top ion, is the shown attached to an output channel that represents measurement.

FIG. 4 .
FIG. 4. (a) The population of the excited state of 1 in response to two δ-pulses spaced t int = 250 apart (gray) along with the independent contributions from each pulse (dashed line) and G 1 (t, 250) (blue), the latter of which captures the interference between coherent excitations due to the pulses.(b) shows the impact of n γ and χ 2 on the simulated value of G 1 (50, 250) in the limit of δ-pulses (t γ = 1).The dashed line marks the exact value.χ 2 = 5 for the n γ sweep, and n γ = 0.25 for the χ 2 sweep.

FIG. 5 .
FIG. 5. (a) A plot of G 1 (200, t) for different t γ .(b) The population of excited state 1 in response to a single photon with a Gaussian-shaped temporal profile.The top panel shows the exact response P M (t) with the simulated time of measurement marked by a vertical dashed line.The gray shaded area shows the profile shape of the intensity of the photon pulse.The contour plot in the bottom panel shows the relationship between the error (computed population less the exact value indicated in the top panel) in the simulated population at t = 200 for different values of t γ and sampling intervals ∆t of t m and t int , for n γ = 1 and χ 2 = 5.

FIG. 6 .
FIG. 6.(a) Schematic of the Example 2. This model has one optically active element, e.g., a quantum dot, that couples to a chain of optically inactive (dark) elements, which could be other quantum dots.We are interested in the transfer of a photoexcitation down the chain to a state (f ) in the final dark element.Note that the two-headed arrows indicate coherent coupling between elements and the single-headed arrow indicate incoherent transfer.Making this transfer difficult to model is that one of the intermediate dark elements (2) is coupled to classical fluctuators (shaded circles), which induce non-Markovian stochastic dynamics in the remainder of the system.(b) The trapped ion simulation is setup similarly to the small system.Again, we model the internal states of each of the quantum dots with internal states of ions (shown in dashed circles).The coherent coupling between quantum dots and fluctuators are implemented through Hamiltonian engineering (see main text).Ion 3 encodes both states in quantum dot 3 and also the f state.The incoherent coupling to state f is implemented through optical pumping via an intermediate state (similar to ion readout).(c) The measured population of the end state for different n γ as a function of t γ = 2∆t int .(d)The G f (∞, t) that describes the response of the system to input pulses of varying width.In both cases, ω 0 = ω 1 = 1.0, ω 2 = 0.8, ω 3 = 0.9,ω f = 0.5,ω α = ω β = 0.1,J 12 = 0.096,J 23 = 0.1,J 2α = J 2β = 0.02, l = 0.01, Γ 2 f = 0.4, Γ 2 α = 0.25,Γ 2 β = 0.5, n γ = 1 and χ 2 = 5.

where 1 , 2 , 3 ,
f index the quantum dot subsystems and α and β are the two fluctuators.Ĥi represent local energies of the quantum dot excited states and fluctuator energies, Ĥij represent the coherent Coulomb coupling between the quantum dots with i, j indexing the dots that are coupled (see double headed arrows in Fig. 6(a)), and Ĥ2j representing the coupling between the quantum dot excited states and fluctuator states.