Simulating Molecular Spectroscopy with Circuit Quantum Electrodynamics

Spectroscopy is a crucial laboratory technique for understanding quantum systems through their interactions with electromagnetic radiation. Particularly, spectroscopy is capable of revealing the physical structure of molecules, leading to the development of the maser - the forerunner of the laser. However, real-world applications of molecular spectroscopy are mostly confined to equilibrium states, due to computational and technological constraints; a potential breakthrough can be achieved by utilizing the emerging technology of quantum simulation. Here we experimentally demonstrate that a superconducting quantum simulator is capable of generating molecular spectra for both equilibrium and non-equilibrium states, reliably producing the vibronic structure of the molecules. Furthermore, our quantum simulator is applicable not only to molecules with a wide range of electronic-vibronic coupling strength characterized by the Huang-Rhys parameter, but also to molecular spectra not readily accessible under normal laboratory conditions. These results point to a new direction for predicting and understanding molecular spectroscopy, exploiting the power of quantum simulation.

Quantum simulation represents a powerful and promising means to overcome the bottleneck for simulating quantum systems with classical computers, as advocated by Feynman [4]. One of the major applications for quantum simulation is to solve molecular problems [5][6][7][8][9]. In recent years, much experimental progress has been achieved in simulating the electronic structures of molecules using quantum devices. Particularly, the potential energy surface of the hydrogen molecule was simulated experimentally [10][11][12]. However, it remains a challenge to scale up this type of experiments for larger molecules, as the phase-estimation method involved requires an enormous amount of computing resources for implementation.
An alternative and potentially more economical approach for quantum molecular simulation has been achieved by using a quantum variational approach [13][14][15][16] that aims to improve the eigenstate approximation through local measurements of the Hamiltonians. So far, most (if not all) of the molecular simulation experiments performed are all confined to the study of static properties of molecules. It is still an experimental challenge to utilize quantum simulators for studying molecular dynamics, in particular, molecular spectroscopy. Furthermore, classical methods in predicting vibrationallyresolved absorption spectra are mostly limited in the gas phase. However, most chemical processes occur in solution, where the molecular vibrational motion depends heavily on the environment; predicting molecular spectroscopy for non-equilibrium states represents a major challenge in quantum chemistry [17].
In this work, we develop and demonstrate a quantum simulation approach for studying molecular dynamics and absorption spectroscopy using a superconducting simulator. Besides simulating molecules in equilibrium, this approach of quantum simulation also allows us to obtain non-equilibrium molecular spectra that are not directly accessible under normal laboratory conditions. In addition, the problem of sampling the absorption spectra of molecules [19] has been found to be related to the problem of Boson Sampling [20], which represents a promising approach to justify that quantum simulators cannot be simulated efficiently with any classical means. Our approach is complimentary with the existing approach [19], where the absorption spectra are obtained by sampling the transition probabilities for each pair of input-output Fock states. The key difference is that we focus on the dynamics of the phonons, instead of the structural shift due to the Duschinsky transformation [21].
More specifically, our approach can be applied to obtaining the temporal correlation function of the electronic transition dipole [3], which yields the information about the absorption spectrum of the molecule, after applying the Fourier transformation. In our superconducting simulator, there are many adjustable control knobs for simulating the spectra for a variety of scenarios. In particular, we are able to simulate molecules in a wide range of values of the Huang-Rhys parameter D, which characterizes the electron-phonon coupling strength.
In this work, we focus on the model approximating the electronic degrees of freedom by a two-level system (Fig. 1a). This model has been applied to study vibronic wavepacket dynamics, chemical reaction rate, Marcus theory for non-adiabatic electron transfer, etc. For molecular spectroscopy, the absorption spectra strongly depend on the initial state of the phonon degree of freedom in the manifold of the electronic ground state. In our experimental demonstration, we have performed simulations by preparing the phonon mode in pure Fock states, as well as simulations for a thermal state and a vacuum  Basic principle of the superconducting simulator. (a) Two identical energy surfaces of a molecule, with one curve displaced from the other along a nuclear coordinate. Near the minimum, the energy surfaces can be approximated by a harmonic potential, with an energy separation of ω0. Here ωeg is the 0-0 energy splitting; in most cases ωeg ω0. (b) The kernel quantum circuit diagram of our method. The circuit consists of an ancilla qubit and a bosonic system. The bosonic system represents nuclear motion mode and is in an initial state |ψ . Similar to DQC1 [18], a composite evolution gate U ≡ e iHg t/ e −iHet/ is applied to the system following a Hadamard gate on the qubit. Finally, measurements along X and Y axis are performed to obtain the correlation function Cµµ = σx + i σy . (c) Device layout and pulse sequence for the superconducting simulator. A "vertical" transmon qubit (dark blue in the enlarged device schematic) on a sapphire chip (light blue) in a waveguide trench couples to two 3D Al cavities. The qubit is first prepared in the ground state |g and the storage cavity (the bosonic system in b) is initialized to different states for various simulations (see the main text for details). As shown in b, the simulation scheme consists of three processes: a qubit rotation R −φ−π/2 (π − γ), a controlled displacement D g (α) of the cavity conditional upon the qubit state |g , and finally a σx or σy measurement of the qubit. Here Rϕ(θ) represents a θ−rotation along ϕ−axis in X-Y plane on Bloch sphere. Note that the rotation angle π − γ in our simulation is not limited to π/2 (see the main text). state with damping. In all cases, we are able to experimentally observe the progression of absorption peaks separated by the vibronic frequency, which is a characteristic feature of molecular spectrum due to vibronic transitions. This flexibility of our superconducting simulator makes it a useful tool for validating theoretical prediction when scaled up (see Methods section).
The architecture of the superconducting simulator is constructed through a three-dimensional (3D) circuit quantum electrodynamics (QED) system [22,23], where a "vertical" transmon qubit is dispersively coupled to two 3D aluminum cavities for storage and readout, as shown in Fig. 1c. The qubit with a transition frequency ω eg /2π = 5.345 GHz is used to model the electronic state {|g , |e } of the molecule. The storage cavity (hereafter referred as the "cavity" for simplicity) is used to model the quantization of the nuclear vibrational motion, i.e., phonons {|0 , |1 , |2 , ...}, with a frequency ω 0 /2π = 8.230 GHz. Note that the energy gap of the qubit is comparable with that of the cavity frequency, i.e., ω eg ∼ ω 0 . However, for a typical molecule, the phonon frequency is much smaller than that of the electronic excitation gap. Therefore, a direct analog molecular simulation with superconducting qubits is not feasible; such a challenge can be overcome by a digital approach of quantum simulation covered in this work.
The working mechanism of our superconducting simulator is summarized as follows (see Figs. 1b and 1c). First, the qubit is initialized to the ground state |g while the phonons (cavity) are prepared in certain given state |ψ for the purpose of simulating the molecular system initially at different nuclear states. In our experiment, we have prepared different phonon states: (i) a vacuum state at zero temperature, (ii) a Fock state |1 at zero temperature, (iii) a thermal equilibrium state, and (iv) a vacuum state with damping. As an example, the pulse sequence for the case of a Fock state |1 is presented in the Supplementary Materials. The qubit is then through a classical microwave pulse turned into a superposition state (|g +|e )/ √ 2, after applying a π/2 rotation (a Hadamard transformation).
Next, a controlled-operation U ctrl is applied to the qubit-phonon system, which drives the evolution of the phonons only if the qubit is in |g , i.e., U ctrl = |g g| ⊗ U + |e e| ⊗ I, where the unitary operator, U (t) ≡ e iHgt e −iHet , first evolves the phonons for a time interval t with Hamiltonian H e , followed by an inverse time evolution with H g for the same time interval. The oper- Here we only present the real part of σ abs . The x axis represents the normalized spectral frequency of electronic transition, which describes the necessary energy for transitions from the electronic ground state to the excited states. We compare experimental results (left) with theory (right) using the same color scale which represents the transition probability. (b) The cross section at D=0, 1, and 4, and the corresponding schematics of the electronic transitions (bottom row). As expected, the absorption peaks of zero-temperature molecular spectrum arises from ωeg, separated by ω0 with a Poisson distribution of intensities. The experimental data are lower by a constant reduction factor f = 0.83 than theory, as expected dominantly due to the decoherence of the qubit during the simulation process.
ation U can be simplified as follows: in the second quantized form, we have the Hamiltonian, H e = ω 0 b † b + ω eg , describing a harmonic oscillator with an equilibrium position shifted by d relative to H g = ω 0 a † a, where b = D(−d)aD(d) = a +d withd = d mω 0 /2 , and D a displacement operator. Consequently, the operator U can be implemented as a displacement operator, U = e −iφ(t) D(d(e iω0t − 1)), apart from a phase factor e −iφ(t) , where φ (t) ≡ ω eg t +d 2 sin ω 0 t (see Supplementary Materials for derivation details).
Note that this phase factor cannot be ignored, as it yields a relative phase instead of global phase with U ctrl . Experimentally, the phase φ is realized in the previous π/2 rotation as an azimuth angle in the X-Y plane on the Bloch sphere (Fig. 1c). The controlled displacement operation D g (α), effective only when the qubit is at |g state as indicated by an extra superscript g, is implemented by a broad selective pulse with a Gaussian envelope truncated to 4σ = 1.34 µs (Fig. 1c). Here the displacement vector α =d(e iω0t − 1). It is worth noting that the decoherence of the system during this long selective pulse lowers the subsequent qubit measurement contrast by a factor of about 0.83 compared to the ideal case.
Finally, as a result the dipole correlation function defined as C µµ (t) = ψ| U (t) |ψ is encoded in the offdiagonal elements of the reduced density matrix of the qubit, i.e., C µµ (t) = σ x + i σ y . σ y and σ x of the qubit can be measured by applying an extra π/2 rotations along X and Y axis (R XorY ) respectively followed by a Z-basis measurement. This general procedure is applicable for any initial state of the phonon, pure or mixed. The absorption spectrum σ abs is finally obtained by a Fourier transform of C µµ (t).
We follow the above procedure to simulate the molecular system initially at a vacuum state and a Fock state |1 at zero temperature. However, in order to simulate molecular spectra with the phonon mode initialized in a thermal state, ρ ≡ e − ω0a † a/kT /Tr(e − ω0a † a/kT ), it is not practical to increase the physical temperature, as the performance of the experimental system would decrease significantly. To overcome this challenge, we can modify the above procedure at Step 1: instead of an equal superposition (after a Hadamard gate), the qubit is ini-  tialized to e −iφ(t) sin γ(t) 2 |g + cos γ(t) 2 |e , where the angle γ(t) is chosen such that sin γ = e 2d 2n (cos ω0t−1) and n = (e ω0/kT −1) −1 (see Supplementary Materials). Similarly, for the case of a vacuum state with damping, we choose sin γ(t) = e −t/τ , where τ is the characteristic time (also see Supplementary Materials). In both cases, following the same remaining procedure as described above, one can obtain the correlation function C thm µµ (t) for an initial thermal state and the damped correlation function C damp Our experimental results are as follows. In our digital simulation, we have set ∆t = 1, t max = 900, ω eg = π/5, ω 0 = π/90. The spectrum lineshape of molecule illustrates the relative probability of electronic transition between different vibrational states in nuclear space. In Fig. 2, we present the progression of absorption peaks for the case where the phonon state is initialized at vacuum and at zero temperature, |ψ = |0 , for various Huang-Rhys parameter D =d 2 . When D = 0, there is only a sharp peak located at the frequency ω = ω eg . This case represents the limit where the electronic transition and the nuclear motion are decoupled. In other words, the molecule is essentially the same as a two-level atom, as far as the spectrum is concerned. When D is increased from zero to, e.g., D = 1, several peaks emerge, and these peaks are equally spaced by the phonon frequency ω 0 . When D is increased further to D = 4, we can observe more equally-spaced peaks. However, the amplitude of the direct transition at ω = ω eg is no longer the largest. In all experimental trials, except for a reduction factor f = 0.83 mainly due to the qubit decoherence, the spectral peaks are in good agreement with the expected Poisson distribution (see Methods).
The absorption spectrum of the other three different initial nuclear states in the molecular system are shown in Fig. 3: (i) Fock state |1 with different D; (ii) thermal equilibrium state at different temperatures characterized by the occupation numbern; and (iii) damped vacuum state with different dissipation rates described by the characteristic time τ . The corresponding electronic transitions for each case have been depicted in the top diagrams of Fig. 3. For clarity, we only show the typical spectra. Except for a reduction of experimental peak values, the experimental results show good agreement with theoretical expectation.
One of the key features of our quantum simulator is that the parameters, such as Huang-Rhys parameter D, can be varied continuously. To better illustrate the progression of the spectrum in Fig. 2 and Fig. 3 as a function of various parameters, we present the peak values at ω = ω eg as an example in Fig. 4. Dots are experimental data by our quantum simulator while the solid curves represent theoretical expectation. After taking into account the reduction of experimental peak values, again mainly due to the system decoherence, by dividing a con-stant reduction factor f (f = 0.83 for Figs. 4a, 4c, and 4d; f = 0.75 for Fig. 4b), the experimental results are in good agreement with theoretical expectations. The smaller f for the case of Fock state |1 is mainly due to the finite Fock state preparation fidelity F = 0.94 (measured Wigner function shown in Supplementary Materials) while all other three cases start from a nearly perfect vacuum state.
To conclude, we demonstrated experimentally a new method to simulate electronic absorption spectra of a molecule, where the nuclear vibrational states may or may not be in thermal equilibrium. Our quantum simulator is based on a superconducting circuit QED architecture with flexible parameter tunability. The simulation results indicate that the resulting molecular spectra are in good agreement with theoretical expectation. Finally, we note that this method can be readily extended to other quantum simulation platform, including photonic [24] or trapped-ion [25] systems. Therefore, our experiment represents the beginning of a new approach of predicting molecular spectroscopy using quantum simulators.
Device parameters.
The transmon qubit has an energy-relaxation time T 1 = 13 µs and a pure dephasing time T φ = 16 µs. The storage cavity has a lifetime τ 0 = 80 µs. The readout cavity has a transition frequency ω m /2π = 7.291 GHz and a lifetime τ r = 42 ns. Together with a Josephson parametric amplifier [26,27] operating in a doublepumped mode [28,29], the fast readout cavity is used for a high fidelity and quantum non-demolition detection of the qubit state (see Supplementary Materials for details). Experimental setup details can also be found in Ref. [30]. The qubit-state-dependent frequency shift of the storage cavity is χ qs /2π = −1.44 MHz, allowing for the qubit-controlled operation on the cavity state as used in our experiment.

Molecular Hamiltonian.
Under the standard Born-Oppenheimer framework, the Hamiltonian H mol of a molecule depends on the nuclear configuration (i.e., position coordinates) q as parameters, H mol (r, q) = K e + U ee (r) + U eN (r, q), where K e is the kinetic-energy term for the electrons, U ee (r) and U eN (r, q) are the electron-electron interaction term and electron-nuclei interaction term respectively. In the low-energy sector, the molecule typically contains an electronic ground state |g and an excited state |e , where the molecular Hamiltonian becomes: with H g = K N + V g (q) and H e = K N + V e (q). Here K N is the nuclear kinetic energy, V g (q) and V e (q) are the potential energies, which are typically approximated as harmonic functions (Fig. 1a), i.e., H g = 1 2m p 2 + mω 2 0 2 q 2 and H e = 1 2m p 2 + Here ω eg is the electronic gap between the minima of both potentials (i.e., 0-0 energy splitting).
The coupling strength between the electronic transition and the nuclear motion is characterized by the Huang-Rhys parameter, D =d 2 , whered = d mω 0 /2 . Similarly, the electronic transition dipole operator is given by µ(q) = µ eg (q) |e g| + µ ge (q) |g e|. However, the dependence of electronic transition moment on nuclear is usually insensitive to the nuclear motion; one can therefore approximate (known as Condon approximation) it with a constant, i.e., µ eg (q) = µ ge (q) = 1 for simplicity.

Absorption lineshape.
The absorption line shape, σ abs (ω) = ∞ −∞ dt e iωt C µµ (t), can be obtained by the Fourier transform of the dipole correlation function C µµ (t) ≡ µ (t) µ (0) , where µ (t) = e iH mol t/ µ(0) e −iH mol t/ . In order to mimic the effects on the molecular spectra due to the influence of the environment [1], one can append a damping factor e −t/τ to the above correlation function, i.e., C damp µµ (t) = e −t/τ C µµ (t), which yields a spectrum with line broadening σ damp abs (ω). Our main task is to apply our superconducting simulator to obtain the correlation functions for the molecules to be simulated. For example, if the initial state is a vacuum state, the correlation function C µµ (t) = e −iωegt e D(e −iω 0 t −1) . The absorption lineshape is σ abs . Thus the spectral peaks are separated by ω 0 with a Poisson distribution of intensities.

Scalability.
Our approach can be scaled up for molecules with multiple vibronic modes.
In this case, the dipole correlation function comes from the contributions of the individual modes, i.e., for n modes, C µµ (t) = |µ eg | 2 e −i(Ee−Eg)t/ F n (t), where F n (t) = Tr e iHgt/ e −iHet/ ρ 1 ⊗ ρ 2 · · · ⊗ ρ n . In other words, the superconducting qubit needs to be coupled with multiple cavity modes. This direction has been realized experimentally [31]. There, a superconducting qubit is coupled to two cavity modes to realize an entangled pair of singlecavity cat states. With a similar geometry, the superconducting qubit can easily be extended to couple to more cavity modes.
We thank R. Vijay and his group for the help with the parametric amplifier.

KEY IDEAS OF THE MOLECULAR SPECTROSCOPY SIMULATION
The main idea of this work is related to the problem of simulating the absorption spectrum of molecules associated with the Huang-Rhys factor D. The molecular spectrum can be obtained by applying the Fourier transformation on the temporal correlation function of an electronic transition dipole operator. Therefore the key part of the experiment involves a simulation of the correlation function through a quantum circuit.
As shown in Fig. 1b in the main text, the core of the quantum circuit consists of three components: 1. Hadamard gate, applied to the ancilla qubit at the beginning of our quantum circuit. Here |g represents the ground state and |e the excited state of the qubit.

A controlled unitary operation
where U ≡ e iHgt/ e −iHet/ is the "Forward and Reversal" time-evolution operator and I is the identity operator.
3. σ y and σ x measurements through π/2 rotations along X and Y axis: Similar to DQC1, We get the real and imaginary part of the correlation function C µµ (t) = ψ| U (t) |ψ = σ x + i σ y , where |ψ is the phonon initial state. Then the absorption spectrum of the molecule can be obtained by a Fourier transformation of C µµ (t).

TEMPORAL CORRELATION FUNCTION
Under the Condon approximation, the dipole operator is given by, where |g and |e represent the ground and excited electronic wave function. Given a molecular Hamiltonian, the time-correlation function C µµ (t) for the dipole operator is given by where σ x (t) = e iHt/ σ x e −iHt/ , and σ x (t = 0) = σ x . Since σ x |g = |e , we can write C µµ (t) = Tr( g| e iHg/ |g g|t σ x e −iHe/ |e e|t |e ) = Tr(e iHgt/ e −iHet/ ).
Here Tr(·) represents the expectation of the operator in the nuclear space. For simplicity, Tr(·) is replaced by · hereafter.

HARMONIC APPROXIMATION
In terms of the creation and annihilation operators, p = i m ω 0 /2(a † − a) and q = /2mω 0 (a + a † ) .
The Hamiltonians of the ground and excited electronic states are given by Here D(d) = ed a † −d a is the displacement operator for a single mode. The Huang-Rhys parameter D is defined as

EVOLUTION PROCESS IN EXPERIMENT
Our quantum simulator simulates the absorption spectra of a molecule in which electronic states (initialized at the ground state |g ) coupled with a single nuclear mode (denoted by |ψ ). In the controlled unitary gate U ctrl = |g g| ⊗ U + |e e| ⊗ I, we have The standard state evolution steps of our quantum simulator are |g |ψ By measuring the expectation of Pauli operators σ x and σ y , we get σ y = Im(e −iωegt−id 2 sin ω0t D(d(e iω0t − 1))) .
In order to simulate molecular spectra with the phonon mode initialized in a thermal state, ρ ≡ e − ω0a † a/kT /Tr(e − ω0a † a/kT ), it is not practical to increase the physical temperature, as the performance of the experimental system would decrease significantly. To overcome this challenge, we can modify the above procedure at the first step: instead of an equal superposition (after a Hadamard gate), the qubit is initialized to e −iφ(t) sin γ(t) 2 |g + cos γ(t) 2 |e , where φ (t) ≡ ω eg t +d 2 sin ω 0 t, the angle γ(t) is chosen such that sin γ = e 2d 2n (cos(ω0t−1) , andn = (e ω0/kT − 1) −1 .
Then we can re-calculate σ x and σ y to get C thm µµ (t) = σ x + i σ y as follows.
This result is consistent with the correlation function for a real thermal equilibrium state (also see [1])

SIMULATING CORRELATION FUNCTION WITH DAMPING
The method to mimic the influence of the environment is similar to the case of a thermal equilibrium state except for sin γ(t) = e −t/τ , where τ is the characteristic time which describes the decay time of the correlation function. This method works for any nuclear state |ψ |g |ψ ⇒ cos γ(t) 2 |e |ψ + e −iωegt−id 2 sin ω0t sin γ(t) 2 |g |ψ (22) It is easy to verify that where C damp µµ (t) (C µµ (t)) is the correlation function with (without) damping.

DEVICE AND READOUT PROPERTIES
The transmon qubit in our experiment is fabricated using the standard Dolan technique [2] with a double-angle evaporation of aluminum after a single electron-beam lithography step on a c-plane sapphire substrate. The experiment is performed in a cryogen-free dilution refrigerator with a base temperature of about 10 mK. A Josephson parametric amplifier (JPA) [3,4] is connected to the output of the readout cavity at the base temperature as the first stage of amplification before the high-electron-mobility-transistor amplifier at 4K. The JPA is operated in a pulsed doublepumped mode [5,6] to minimize pump leakage to the readout cavity. The readout pulse is calibrated to contain only a few photons (see the calibration session below), enough for a high-fidelity single-shot readout of the qubit state. The schematic of the measurement setup can be found in Ref. 7.
The readout property of the qubit is shown in Fig. 1. The readout histogram is clearly bimodal and well separated. A threshold V th = 0 is chosen to digitize the readout signal to +1 and −1 for the ground state |g and the excited state |e , respectively. The qubit has an excited state population of 3.9% in the steady state, presumably due to stray infrared photons or other background noise leaking into the cavity, although the exact source remains unknown. We use measurement-based post-selection to purify the qubit to |g state. This requires the measurement to be high quantum non-demolition (QND). Inset of Fig. 1 shows the qubit readout matrix with the cavity left in vacuum and the corresponding experiment sequence. After the first measurement and post-selection of |g , the subsequent qubit The absence or presence of a π pulse at the beginning together with M1 is used to prepare the qubit's initial state |g or |e by a post-selection. (b) Readout matrix. The numbers outside the parenthesis correspond to a π pulse with σ = 6 ns for α = 1 (n = 1). The numbers inside the parenthesis correspond to a π pulse with σ = 2 ns for α = 4 (n = 16). measurement shows that the probability of the qubit being populated in |g is as high as 0.999 (dashed histogram), demonstrating the high QND property of the qubit readout. If we post-select the qubit at |e state, the subsequent measurement finds that the qubit remains in |e with a fidelity of 0.953. These errors dominantly come from the T 1 process during the waiting time after the initialization measurement (250 ns) and during the readout time (240 ns). Figure 2 shows the properties of non-selective π pulses with a coherent state |α in the cavity. In our simulation, there is always a coherent state presented in the storage cavity (see Fig. 3). The associated measurement pulses are shown in Fig. 2a. The absence or presence of a π pulse at the beginning together with a measurement M1 is used to prepare the qubit's initial state |g or |e by a post-selection. The numbers outside the parenthesis correspond to a π pulse with σ = 6 ns for α = 1 (n = 1) in the storage cavity. The numbers inside the parenthesis correspond to a π pulse with σ = 2 ns for α = 4 (n = 16) in the storage cavity. The fidelity loss mainly comes from two parts: the qubit T 1 process during the measurement time and the waiting time between the two consecutive measurements, and the qubit frequency shifts due to photon number occupations in the cavity.

MEASUREMENT PULSE SEQUENCE
The pulse sequence for the experiment is shown in Fig. 3a. Here we use the case of the nuclear state in a Fock state |1 as an example. The whole experiment can be divided into three main parts: 1) initialization of the system to |g, 0 by post-selecting results of the qubit measurement and the cavity parity measurement; 2) creation of |1 state by first displacing the cavity to a coherent state |α = 1 , then applying a selective π rotation (σ = 360 ns) of the qubit corresponding to N = 1 photon in the cavity, and finally post-selecting the excited state of the qubit; 3) simulation as described in the main text.
The classical microwave pulse to displace the cavity state has either a square envelope with a width of 100 ns or a Gaussian envelope with σ = 335 ns. The amplitudes have been calibrated carefully as shown in the section of "calibration of driving amplitude". All qubit drive pulses have a Gaussian envelope pulses truncated to ±2σ. To eliminate the possible qubit leakage to higher qubit levels, we also apply the so-called "derivative removal by adiabatic gate" technique [8] for pulses with σ =2, 4, and 6 ns. Figure 3b shows the Wigner tomography of the created |1 state and Fig. 3c shows the moduli of the reconstructed density matrix ρ by least-square regression using a maximum likelihood estimation [9,10]. The measured fidelity of this state is F = 1| ρ |1 = 0.94. The Wigner tomography of the cavity state is performed by a cavity's displacement operation D(−β) followed by a parity measurement [7,[11][12][13][14]. The parity measurement is achieved in a Ramsey-type measurement of the qubit with a conditional cavity π phase shift C(π) sandwiched in between two unconditional qubit rotations R Y (π/2) followed by a projective measurement of the qubit.

CALIBRATION OF DRIVING AMPLITUDE
To have a high-QND and high-fidelity single-shot readout is essential for any experiment that requires postselections. There should be enough readout photons for a high signal-to-noise ratio of the readout, but not too many to take a long time for those readout photons to leak out, otherwise slowing down the subsequence operations. In Fig. 4, the readout photon number is calibrated through a measurement-induced dephasing process. Figure 4a shows the measurement pulse sequence which is a Ramsey-type measurement of the qubit inserted by a readout pulse with various amplitudes at a fixed width T m . Figure 4b shows the decaying signal, coming from the dephasing due to readout photons, as a function of the measurement pulse amplitude. An exponential fit gives a calibration of the readout pulse. In our experiment, we typically readout the qubit with aboutn = 5 photons in the cavity and wait for nearly six times the photon decay time to make sure the average remaining photon number in the readout cavity is only about 1%.
The calibration of the controlled cavity displacement D g (α) (Gaussian envelope with σ = 335 ns) is critical for our simulation scheme. This calibration is realized by measuring the probability of the first nine Fock states N = 0, 1, 2, ...8 through a selective π pulse (Gaussian envelope with σ = 360 ns) as a function of the displacement pulse amplitude, as shown in Fig. 5. A nearly perfect global fit to the measurement results with a Poisson distribution not only gives the required calibration of the cavity displacement amplitude, but also indicates good control of the coherent state in the cavity. The calibration of the displacement D(1) with a 100 ns square envelope is performed by measuring the photon number parity instead as a function of displacement pulse amplitude (inset of Fig. 5). Figure 6 shows the typical spectral peak values for the nuclear system at vacuum and at zero temperature (Fig. 2  of the main text). The peaks are all normalized by dividing a constant reduction factor f = 0.83 (mainly due to the system decoherence), and are in good agreement with the expected Poisson distribution. (a) A Ramsey-type measurement sequence for the qubit inserted by a readout pulse with various amplitudes at a fixed width Tm. The second π/2 pulse has a rotating phase relative to the first one in order to have an oscillating interference signal. (b) The probability of the qubit at the ground state as a function of the measurement pulse amplitude. The Ramsey amplitude A ∼ e −TmΓm gives a calibration ofn in the readout pulse, where the measurementinduced dephasing rate Γm = 2nκsin 2 (tan −1 (χqr/κ)), κ is the readout cavity's decay rate, and χqr is the dispersive frequency shift between the readout cavity and the qubit. FIG. 5. Calibration of the controlled cavity displacement D g (α) (Gaussian envelope with σ = 335 ns). The probability of the first nine Fock states N = 0, 1, ..., 8 is measured as a function of the displacement pulse amplitude through a selective π pulse (σ = 360 ns) on the qubit at each resonant transition frequency corresponding to a specific photon number N in the cavity. There is no normalization of the measurement data and the loss of probability dominantly comes from the T1 process in the long duration of the π pulse. Lines are from a global fit with a Poisson distribution P (|α , N ) = A|α| 2N e −|α| 2 /N !, where A is a scaling factor accounting for the probability loss. The excellent agreement indicates good control of the coherent state in the cavity, giving a calibration DAC level=326 corresponding to α = 1. Inset: calibration of the displacement D(1) with a 100 ns square envelope. This calibration is performed by measuring the photon number parity instead as a function of displacement amplitude, giving DAC level=2649 for α = 1. The loss of parity measurement contrast mainly comes from the qubit decoherence during the parity measurement.