Exploring Interacting Quantum Many-Body Systems by Experimentally Creating Continuous Matrix Product States in Superconducting Circuits

Improving the understanding of strongly correlated quantum many body systems such as gases of interacting atoms or electrons is one of the most important challenges in modern condensed matter physics, materials research and chemistry. Enormous progress has been made in the past decades in developing both classical and quantum approaches to calculate, simulate and experimentally probe the properties of such systems. In this work we use a combination of classical and quantum methods to experimentally explore the properties of an interacting quantum gas by creating experimental realizations of continuous matrix product states - a class of states which has proven extremely powerful as a variational ansatz for numerical simulations. By systematically preparing and probing these states using a circuit quantum electrodynamics (cQED) system we experimentally determine a good approximation to the ground-state wave function of the Lieb-Liniger Hamiltonian, which describes an interacting Bose gas in one dimension. Since the simulated Hamiltonian is encoded in the measurement observable rather than the controlled quantum system, this approach has the potential to apply to exotic models involving multicomponent interacting fields. Our findings also hint at the possibility of experimentally exploring general properties of matrix product states and entanglement theory. The scheme presented here is applicable to a broad range of systems exploiting strong and tunable light-matter interactions.

Improving the understanding of strongly correlated quantum many body systems such as gases of interacting atoms or electrons is one of the most important challenges in modern condensed matter physics, materials research and chemistry. Enormous progress has been made in the past decades in developing both classical and quantum approaches to calculate, simulate and experimentally probe the properties of such systems. In this work we use a combination of classical and quantum methods to experimentally explore the properties of an interacting quantum gas by creating experimental realizations of continuous matrix product states -a class of states which has proven extremely powerful as a variational ansatz for numerical simulations. By systematically preparing and probing these states using a circuit quantum electrodynamics (cQED) system we experimentally determine a good approximation to the ground-state wave function of the Lieb-Liniger Hamiltonian, which describes an interacting Bose gas in one dimension. Since the simulated Hamiltonian is encoded in the measurement observable rather than the controlled quantum system, this approach has the potential to apply to exotic models involving multicomponent interacting fields. Our findings also hint at the possibility of experimentally exploring general properties of matrix product states and entanglement theory. The scheme presented here is applicable to a broad range of systems exploiting strong and tunable light-matter interactions.
Progress in revealing relations between solid state physics and quantum information theory has constantly extended the range of quantum many body problems which are tractable with classical computers. One such successful approach is the density matrix renormalization group (DMRG), which was introduced by White in 1992 [1] and since then developed into a leading method for numerical studies of strongly interacting one dimensional lattice systems [2]. Later it was realized that the DMRG can be interpreted as a variational optimization over matrix product states (MPS) [3,4]. The class of matrix product states [5,6] naturally incorporates an area law for the entanglement entropy [7] and is thus ideally suited to parameterize many-body states with finite correlation length [8]. An interesting connection between the MPS formalism and open quantum systems has recently been discovered [9], which led to the suggestion of using the high-level of experimental control achievable over open cavity QED systems to create continuous matrix product states [10] as itinerant radiation fields for the purpose of quantum simulations [11,12]. In this letter we provide experimental evidence that this concept is indeed capable of determining the properties of strongly correlated quantum systems and offers promising perspectives, complementary to existing digital and analog quantum simulation approaches [13] explored with trapped atoms and ions [14][15][16][17][18].
In clear distinction to previous experiments, here, we simulate a continuous quantum field theory rather than a lattice model. In particular, we study the ground-state of the Lieb-Liniger modelĤ = dxĤ = dx(N +T +Ŵ ) [19], which describes a gas of interacting bosons confined in a one-dimensional continuum [20], as schematically depicted in Fig. 1a. Here,N = −µψ † xψx is the potential energy,T = ∂ xψ † x ∂ xψx is the kinetic energy of particles, and * Current address: Department of Physics, Princeton University, Princeton NJ 08544, USA. Correspondence to ceichler@princeton.edu  Repulsion between particles mediates an interaction energyŴ . The particle density ρ of the gas is controlled by the chemical potential µ. (b), We experimentally simulate the ground-state ofĤ by employing a variational minimization procedure. A tunable cavity QED system is used to generate radiation fields emulating continuous matrix product states |φ(λ) in a 1D transmission line. The average energy E λ = φ(λ)|Ĥ|φ(λ) of the simulated Hamiltonian is experimentally determined from measured correlation functions. External control fields λ are used as variational parameters. (c),(d), Examples of first-order G (1) (τ ) and second-order g (2) (τ ) correlation functions measured (dots) in superconducting circuits and simulated (solid lines) using a master equation approach for the indicated drive rates Ω, effective anharmonicity α/2π = 5.2 MHz, and cavity decay rate κ/2π = 2.2 MHz.  proportional to the potential energy N , (b) dωω 2G(1) (ω) proportional to the kinetic energy T , and (c) second-order correlator g (2) (0) proportional to the interaction energy Ŵ . The smallest drive rate is Ω0/2π = 0.37 MHz.
x is the interaction energy expressed in second quantization by the field operatorψ x . The Lieb-Liniger model has only one intensive parameterṽ = v/ρ, where ρ = ψ † xψx is the average particle density and v the interaction strength. The ground-state energy of this model can be calculated analytically using the Bethe ansatz [19]. The calculation of two point correlation functions requires the use of numerical methods such as quantum Monte Carlo or DMRG [10]. The fact that this model is well understood makes it an ideal test case to benchmark the as yet unexplored quantum variational algorithm recently proposed by Barrett et al. [12].
In the experiments presented here, we prepare continuous matrix product states |φ(λ) [10,12] as microwave radiation fields propagating along a one-dimensional transmission line, see Fig. 1b. The radiation fields are generated by an ancillary quantum system -in our case a tunable circuit QED system [21] -which is coupled with rate κ to the transmission line. Notably, any radiation field generated in this way is described by a continuous matrix product state with a bond dimension depending on the number of participating ancillary energy levels [9]. We vary the quantum state |φ(λ) by tuning a set of two external variational parameters λ = (α, Ω). Here, Ω is the drive rate of a coherent field applied resonantly to the upper eigenmode of the coupled system and α is the effective anharmonicity of the driven mode. Tunability of the effective anharmonicity α is experimentally achieved by employing a qubit of which both the frequency and its coupling to the cavity are adjustable in-situ during the experiment [22], see Appendix A 5 for details.
The variational ground-state of the Lieb-Liniger Hamilto-nianĤ is found by measuring the expectation value E λ = φ(λ)|Ĥ|φ(λ) = N + T + Ŵ and minimizing E λ with respect to different variational states |φ(λ) created in our experiment. The simulated Hamiltonian is thus solely determined by the measurement observable. Any model of which the corresponding expectation values can be measured, is therefore accessible with this approach. Given the photonic realization of |φ(λ) , the measurement of E λ translates into the measurement of photon correlation functions. Spatial correlations in the fieldψ x are mapped onto time correlations in the cavity output fieldâ out (t) by identifyinĝ ψ x =â out (t = x/s)/ √ s, where the scale parameter s = x/t acts as an additional variational parameter [12]. Entanglement in the matrix product states thus corresponds to entanglement between photons emitted from the cavity at different times. According to this correspondence, E λ for the Lieb-Liniger Hamiltonian is given by the first-and second-order correla- [23]. More specifically, the average kinetic energy T = s −3 dωω 2G(1) (ω) is calculated from the Fourier transform of the first-order correlation functionG (1) (ω), the interaction energy is Ŵ = s −2 vG (2) (0) and the potential energy is given by the average photon flux N = −s −1 µG (1) (0).
The presented variational approach thus crucially relies on the ability to generate and probe a wide range of different (quantum) radiation fields with high efficiency. For fast and reliable correlation measurements [24] we have developed a quantum-limited amplifier which allows for phase-preserving amplification at large bandwidth and high dynamic range [25]. The examples of measured correlation functions shown in Fig. 1c-d illustrate their dependence on the drive rate Ω at constant α. While G (1) (τ = 0) equals the total average photon flux a † out a out , the limit G (1) (τ → ∞) is proportional to the square of the coherence of the field | a out | 2 . Therefore, G (1) increases with drive rate due to the enhanced photon production rate. The normalized second-order correlation functions g (2) (τ ) ≡ G (2) (τ )/(G (1) (0)) 2 show anti-bunched behavior (g (2) (0) < 1) for weak drive and Rabi type oscillations when the drive rate Ω becomes larger than the decay rate [24]. Both the measured first-order and second-order correlation functions are in agreement with the results obtained from master equation simulations (black solid lines).
We have measured ∼ 10 3 such correlation functions over a wide range of variational parameters Ω and α by acquiring data for one week. Without the employed parametric amplifier the time for measuring this set of data would have been on the order of years because of the exponential scaling be-tween statistical error and correlation order [26]. Based on this collection of time-resolved correlation functions we have evaluated the three relevant terms G (1) (0), dωω 2G(1) (ω) and g (2) (0) that enter the calculation of E λ , see Fig. 2. As expected, the average photon flux G (1) (0) (Fig. 2a) increases with drive rate Ω and is suppressed for increasing anharmonicity α. The kinetic energy term dωω 2G(1) (ω) in panel b is determined by the power spectral densityG (1) (ω). Only the spectral weight of photons generated at finite detuning from the drive frequency (|ω| > 0) contributes to the integral. In the Bose gas picture such photons correspond to particles in finite momentum states and therefore carrying kinetic energy. The rate of scattering events from drive photons into photons with finite ω increases with drive strength and has a non-trivial dependence on α. Finally, the second-order correlator g (2) (0) in Fig. 2c clearly reveals the crossover from antibunched radiation (g (2) (0) → 0) for large α and small Ω to coherent radiation (g (2) (0) → 1) when lowering the anharmonicity.
Based on the three measured quantities shown in Fig. 2 and a chosen interaction parameter v we evaluate the energy E λ for all prepared states |φ(λ) . We identify a local minimum in variational space (blue region), which corresponds to the variational ground-state of the Lieb-Liniger model, see Fig. 3. When changing the interaction parameter v we find the energy E λ to be minimized by a different set of parameters (α, Ω) in  After identifying the variational ground-state for each interaction strength v as the respective minimum in the energy landscape, we further investigate its properties. We compare the experimental results with the numerically exact results obtained from a variational matrix product state algorithm executed on a classical computer with bond dimension D = 14 [10]. We followed the usual convention and rescaled all quantities so that they correspond to a particle density of ρ = 1 [10]. The Lieb-Liniger ground-state energy density E LL , which is the sum of kinetic energy and interaction energy, increases with interaction strength and ideally converges towards the Tonks-Giradeau limit lim v→∞ E LL = π 2 /3 indicated as a dashed line in Fig. 8a. Given the small number of variational parameters the experimental data (blue dots) reproduces the characteristic dependence of E LL onṽ of the exact solution (red solid line) quite well.
Importantly, having physical access to the ground-state wave functions |φ(λ) we can also probe quantities beyond the ground-state energy, such as two-point correlation functions. First-order correlation functions ψ † xψ0 are obtained from G (1) (τ ) by converting time into spatial coordinates (Fig. 8b). As expected, we observe a decrease in correlation length with increasing interaction strengthṽ. Due to the absence of spontaneous symmetry breaking in one dimension [27], the exact ground state of the Lieb-Liniger model does not exhibit Bose-Einstein condensation. The observed finite limit lim x→∞ ψ † xψ0 is a characteristic feature of matrix product states which do not support the U (1) symmetry of the model for finite bond dimensions.
The nontrivial nature of the ground-state in the presence of interactions also becomes manifest in the particle-particle correlator ψ † 0ψ † xψxψ0 shown in panel c. With increasing v particles are more likely to repel each other leading to antibunching. Our experiments clearly resolve this crossover from a weakly into a strongly interacting Bose gas by accessing variational wave functions for interaction parametersṽ over two orders of magnitude. While this general behavior is qualitatively well reproduced, an accurate quantitative agreement with the numerical results would require a larger number of independent variational parameters in the experimental realization. This becomes apparent when comparing the experimental results with numerical calculations based on continuous matrix product states of different bond dimensions, D = 2 and D = 14, where the number of variational parameters is 2D 2 . Correlation functions simulated with low bond dimension D = 2 deviate from the exact results (D = 14) similarly as the measured ones ( Fig. 8d-g).
In summary, we have experimentally revealed connections between open quantum systems, the matrix product state variational class, and quantum field theories that can be used for practical quantum simulations. The presented quantum variational algorithm is general in the sense that it can be applied 0.8 0.8 to any one-dimensional quantum field theory. Exploring interacting vector field models seems particularly appealing, since they are difficult to simulate on classical computers. Experimentally, this could be achieved by coupling tunable quantum systems to multiple transmission lines each representing one of the vector field components [28]. Higher accuracy in the simulation will require more variational parameters and quantum systems with more degrees of freedom, which is achievable with tunable superconducting circuits. In this context, the collective dissipation of multiple emitters coupled to the same transmission line at finite distance [29] may turn out very useful. Extensions of the presented approach may be envisaged to also explore dynamical phenomena and discrete lattice models using experimentally created matrix product states.
We thank Christoph Bruder, Ignacio Cirac, Tilman Esslinger and Atac Imamoglu for discussions and comments. This work was supported by the European Research Council (ERC) through a Starting Grant, by the NCCR QSIT and by ETHZ. CE acknowledges support by Princeton University through a Dicke fellowship. TJO was supported by the ERC grants QFTCMPS and SIQS, and by the cluster of excellence EXC201 Quantum Engineering and Space-Time Research. The experiments presented in the main text are performed with a device consisting of a superconducting circuit, see Fig. 5a for details about the experimental setup. The sample ( Fig. 5b) consists of a λ/2 transmission line cavity with a resonance frequency ω res /2π ≈ 7.3425 GHz of the fundamental mode. The resonator has one output port which dominates the total decay rate κ/2π ≈ 2.2 MHz and one weakly coupled input port κ in /κ ≈ 0.01 which is used for coherent driving of the cavity field. The resonator is fabricated using photolithography and reactive ion etching of a Niobium thin film sputtered on a sapphire wafer. We have fabricated a superconducting qubit (Fig. 5c) close to one end of the resonator. Both its transition frequency ω ge and coupling strength g to the resonator are tunable by varying the magnetic fluxes threading the two SQUID loops [22,30]. Flux control is achieved by a combination of a superconducting coil mounted on the backside of the sample holder and a local flux line which couples predominantly to one of the two SQUIDs. The currents feeding the coil and the flux line are generated by voltage biased resistors at room temperature. The qubit is fabricated using double-angle evaporation of aluminum on a mask defined by electron beam lithography. The qubit decay and dephasing times are measured to be T 1 ≈ 1.8 µs and T 2 ≈ 1.2 µs. The anharmonicity of the qubit is (ω ef − ω ge )/2π ≈ −80 MHz, where ω ef is the transition frequency from the first excited state |e to the second excited state |f of the qubit.
The sample is mounted on the base plate of a dilution refrigerator cooled down to a temperature of about 20 mK. The qubit and the resonator are coherently driven through attenuated charge control lines. The microwave radiation emitted from the cavity is guided through two circulators to a Josephson parametric dimer (JPD), which provides quantumlimited amplification at large bandwidth and dynamic range [25,31,32]. A directional coupler is used to apply and interferometrically cancel the reflected pump field. The ampli- fied signal reflects back from the JPD, passes a bandpass (BP) filter, is further amplified by a high electron mobility transistor (HEMT) amplifier, and is down-converted to an intermediate frequency (IF) of 25 MHz. After low-pass (LP) filtering and IF amplification the down-converted signal is digitized using analog-to-digital conversion (ADC) and further processed with field programmable gate array (FPGA) electronics.

Controlling the qubit frequency and the coupling strength
We characterize the coupled cavity-qubit system by probing the transmission coefficient of the cavity and fitting the data to the absolute square of the expression which we obtain from input-output theory for the Jaynes-Cummings model [33]. Here, ω is the probe frequency, γ is the qubit decoherence rate and A is a scaling factor. The probe power is chosen such that the average number of excitations of the coupled resonator qubit system is much smaller than one. In this case the qubit may be approximated by an harmonic oscillator. We determine the qubit detuning ∆ = ω res − ω ge and its coupling strength g to the cavity by fitting spectroscopically obtained transmission data to the above model, see Fig. 6. The magnetic fluxes through the qubit SQUID loops and with that the qubit parameters g and ∆ are controlled by a pair of voltages V 1 and V 2 applied to coil bias resistors. For small g and ∆ we approximate the relation between (g, ∆) and (V 1 , V 2 ) by linear equations of the form g ∆ = m 11 m 12 m 21 m 22 We determine the coupling matrix elements m ij and the offset voltages V 1,0 and V 2,0 by recording transmission data for pairs (V 1 , V 2 ) k . For each set of data we extract the corresponding parameters (g, ∆) k and perform a least-square fit to determine the model parameters m ij , V 1,0 and V 2,0 . By inverting Eq. (A2) we calculate the voltages V 1 and V 2 for a given set of desired qubit parameters (g, ∆). In order to further fine-tune the parameters, we have developed an automated calibration algorithm which minimizes the deviations from desired target values by iteratively measuring, fitting and readjusting the control parameters V 1 and V 2 . With this control procedure we are able to independently set the qubit frequency and the interaction strength. We also use this procedure to monitor and correct for slow qubit frequency drifts occurring during long runs of the experiment. We demonstrate individual control of qubit parameters by either tuning the qubit frequency for approximately constant coupling strength g (Fig. 6a) or by varying g for fixed qubit frequency (Fig. 6b). For all sets of data we have turned on the JPD amplifier and have divided out its frequency dependent gain. For the measurements in Fig. 6b we have kept the qubit resonant with the cavity (∆ = 0) and have varied g. These measurements demonstrate the ability to tune the system from the fast cavity (κ g γ) into the strong coupling regime (g γ, κ) [22].

Josephson parametric dimer amplifier
To measure higher order correlation functions efficiently while retaining a high level of linearity we employ a Josephson parametric dimer (JPD) amplifier. For details about the operational principles of the JPD we refer the reader to [25]. The gain G measured vs. signal frequency f is approximately described by a Lorentzian function (Fig. 7). The amplifier bandwidth at full width half maximum is approximately 35 MHz. In order to increase dynamic range, we have chosen a moderate maximum gain of 17 dB. A measurement of the gain as a function of signal power P in results in a 1 dB compression point at about -107 dBm which corresponds to a photon flux of 3000 µs −1 , see Fig. 7b. The largest photon flux which is generated in the described experiments is less than 50 µs −1 . The JPD amplifier is thus far away from its compression point for all measured correlation functions. The improvement in detection efficiency becomes apparent when measuring the noise power spectral density S δ in units of photons per Hz per second when the JPD amplifier is turned ON and when it is turned OFF. The effective noise level, which is referenced back to the input of the JPD amplifier is decreased by more than an order of magnitude when it is turned on. The scaling of S δ is based on a comparison between the frequency dependent gain G δ and the JPD amplifier noise [34]. The deviation from the quantum limit is due to the additional noise from the following HEMT amplifier which is comparable to the noise at the output of the JPD. The equivalent detection efficiency of the amplification chain is η amp = 1/S δ ≈ 50%.

Calibration of drive rate and output power
We have calibrated the total gain of the detection chain including all cable losses in order to reference the measured photon flux back to the output of the cavity. Calibrating the total gain of the detection chain is equivalent to calibrating its detection efficiency η tot . The detection efficiency is typically limited by losses between the cavity and the first amplifier (η loss ) and by noise added during the amplification process (η amp ), see Fig. 8. We compare the nonlinear response of the  8: (a) The detection efficiency of the cavity output field is typically limited by radiation loss, schematically represented as a beamsplitter with finite transmittivity η loss , and by noise S δ = η −1 amp added in the amplification chain. The total detection efficiency is given by the product ηtot = η loss ηamp. (b) Measurement (points) and fit (solid line) of the cavity photon number and absolute square of the coherent amplitude for the coupled cavity-qubit system driven through the cavity input port at rate Ω. The smallest drive rate is Ω0/2π = 0.37 MHz. coupled cavity-qubit system with master equation simulations to perform this calibration. We bias the qubit with parameters (∆, g)/2π = (3, 5.7) MHz and apply a drive field to the input port of the cavity at rate Ω and resonant with the frequency ω + = ω ge + ∆/2 + g 2 + ∆ 2 /4 = 2π × 7.35 GHz of the upper Jaynes-Cummings doublet state. For these settings we measure the coherent photon flux κ| a | 2 and the total photon flux κ| a † a | 2 emitted from the cavity for different drive rates Ω. We fit these data sets to the results obtained from master equation simulations leaving the the total gain factor of the measurement chain and the absolute drive rate incident to the sample as free parameters (Fig. 8b). The equivalent detection efficiency resulting from this fit is equal to the inverse of the scaled noise level and found to be η tot = η loss η amp = 0.27. Together with the estimate for the efficiency of the amplification chain η amp stated in the previous section, we extract a radiation loss of 1 − η loss = 0.46 between the cavity and the JPD, which is reasonable given the components and cables connecting the two stages.

Drive scheme and variational parameters
In the original proposal for simulating the Lieb-Liniger model with cavity QED it has been suggested to keep the qubit resonant with the cavity and use the qubit drive power Ω q and the coupling strength g as two variational parameters. We have experimentally realized this scheme and found that in the limit of small coupling strengths g the total emission rate becomes extremely small which in turn limits the signal to noise ratio. This is because the effective emission bandwidth scales like g 2 /κ when g < κ and thus decreases quadratically with g. We have therefore developed an alternative scheme for which the photon emission rate remains proportional to κ even in the limit of small g. We have therefore made use of the ability to tune both the qubit frequency and the coupling strength. Rather than keeping the qubit at fixed frequency we adjust for each value of g the detuning ∆ such that ω + remains at constant frequency resonant with the drive frequency ω + = ω d = 2π × 7.35 GHz. The spectroscopy data for the qubit bias points used in the quantum simulation experiment are shown in Fig. 9a and demonstrate constant ω + over the entire range of coupling strengths. In order to keep ω + constant, we compensate the larger splitting when increasing g by tuning the qubit further away from the cavity, see Fig. 9b. For this specific tuning scheme we find that the effective anharmonicity α of the upper Polariton ladder decreases with increasing coupling strength g, as shown in Fig. 9c. To illustrate this effect we show a schematic drawing of the energy levels |n± of the Jaynes Cummings model for the resonant case (∆ = 0) and for the case of finite qubit detuning ∆ < 0. The inverse proportionality between anharmonicity α and coupling strength g illustrates the appearance of antibunched behavior for the small values of g and the observed coherent radiation for large g values for this bias scheme, see Fig. 2 of the main text. The drive rate Ω of a coherent field applied to the cavity input port acts as a second variational parameter for the quantum simulation.

Measurement of correlation function and master equation simulation
We employ fast real-time signal processing performed with an FPGA at a clock rate of 100 MHz for the measurement of time-resolved correlation functions. The cavity output field is processed, as described in section 1A. After digitization we multiply the sampled voltages with digital sine and cosine waves of frequency ω IF /2π = 25 MHz to obtain the quadrature components I(t) and Q(t), respectively. We then apply an FIR filter with an effective bandwidth of Γ/2π = 10 MHz to the quadratures in time-domain to obtain the filtered quadraturesĨ(t) andQ(t), which in the following we write as the single complex valued amplitude S(t) = I(t) + iQ(t). To measure the first-order correlation function in S(t) we take the discrete Fourier transform (F) of M time traces S i (t) of length 10.24 µs, multiply with their complex conjugate, and average Similarly we extract the second-order correlation function by calculating the absolute square of S(t) before Fourier transforming We record each of these quantities with the drive field turned on, giving Γ ON (τ ), and with the drive turned off, giving Γ OFF (τ ). To avoid effects due to slow drifts we alternate between all four measurements every 12.5 µs. As explained in detail in reference [26] and as demonstrated in many experiments since then [24,35,36], we can use these four measurements to extract the correlation functions G (1) (τ ) = κ a † (τ )a(0) and g (2) (τ ) = a † (0)a † (τ )a(τ )a(0) / a † a 2 of the output field of the cavity. In these expressions, a (a † ) is the annihilation (creation) operator of the intra-cavity field. For the cases in which the average photon number is small ( a † a ) 1 we find g (2) (0) values which are systematically smaller than the corresponding master equation simulation. We attribute this to a weak thermal background radiation during the off measurements which we correct for [37]. Good agreement between the measured and simulated correlation functions is found when correcting for a thermal photon flux of n th ≈ 0.03/µs in the detection band.
We compare these measurements with correlation functions obtained from master equation simulations. For these simulations we describe the system by the Hamiltonian expressed in a frame rotating at the drive frequency ω d /2π = 7.35 GHz. Here, b and b † are annihilation and creation operators for an excitation of the transmon. In addition to that, we account for qubit decay γ, qubit dephasing γ φ and resonator emission κ with standard Lindblad terms. Simulations are run in a Hilbert space including 6 resonator and 3 transmon levels. In order to account for the finite detection bandwidth when simulating the second-order correlation function we employ the techniques described in [38]. In this approach we introduce an ancillary mode c of frequency ω d , which is weakly coupled with rate /2π = 20 kHz to the cavity and decays with a rate equal to the detection bandwidth Γ. The second-order correlation function in c is then simulated based on the total Liouvillian and taken as an estimate for the filtered correlation function of mode a.
Appendix B: Theoretical aspects and data analysis 1. Calculation of the Lieb-Liniger energy from correlation functions The expectation value to be minimized Ĥ = T + Ŵ + N is composed of the kinetic energy of the bosons T , its interaction energy Ŵ and the potential energy N . According to the correspondence between the field operatorψ x and the time-dependent radiation fieldâ(t) = √ sψ x=st , each of these expectation values is proportional to a specific measured correlation function is the Fourier transform of the first-order correlation function normalized such that dωG (1) (ω) = G (1) (0). The energy terms in Eq. (B1) thus explicitly depend on the scaling parameter s, which is to be treated as an additional variational parameter. We explicitly minimize H with respect to s by solving ∂ ∂s H = 0, which results in Correspondingly we obtain an explicit expression for E λ = Ĥ , which only depends on the measured correlation functions and the model parameters µ, v. We find the variational ground state for a given set of model parameters (v, µ) by minimizing E λ with respect to α and Ω.

Scaling transformation
We follow the usual convention and study the Lieb-Liniger ground state subject to the constraint that its particle density ρ is equal to one. Using the procedure described in Sec. B 1 we find a ground state which generally does not obey this property. We therefore apply a scale transformation by adjusting the chemical potential µ such that ρ → 1. Under the following transformation (µ, v) → (μ,ṽ) = (µy 2 , vy), the ground state remains invariant up to a change in the parameter s →s = s/y, which immediately follows from Eq. (B2). Any variational ground state can therefore be transformed into another ground state satisfying ρ = 1, by choosing y appropriately. We apply the following procedure to perform this scale transformation: • Chose parameter v, set µ = 1, and find the set of variational parameters (s min , Ω min , α min ) minimizing E λ .
• The variational parameters (s min , Ω min , α min ) specify the variational ground state for the model with interaction strengthṽ and unit particle density.
Ground states for different interaction parameters are obtained by starting the above procedure with a different value for v. The Lieb-Liniger energy E LL at interaction strengthṽ (Fig. 4a of the main text) is given by Correlation functions for the Lieb-Liniger model are directly obtained from the measured correlation functions by identifying and analogously for the second-order correlation function ψ † x ψ † 0 ψ 0 ψ x = g (2) (τ = x/s min ). (B5)

Numerically exact solution
The exact solution [19,39], via the Bethe ansatz, of the Lieb-Liniger model at unit density is only possible for one specific value of the interaction parameter, namely v = 2. In order to calculate properties of the model at unit density for other values of v it is necessary to take recourse to numerical methods. We exploited a variational method over translation invariant continuous matrix product states (cMPS) [10,11,40] |Ψ[Q, R] ≡ tr P exp ∞ −∞ Q ⊗ I + R ⊗ψ † x dx |Ω , (B6) where P exp denotes the path ordering of the argument from left to right for increasing values of x. The operators Q and R act on an auxiliary space C D , and |Ω is the Fock vacuum. The variational parameters specifying the cMPS are precisely the two D × D matrices Q and R. The Lieb-Liniger hamilto-nianĤ (in the presence of a chemical potential) is comprised of three termsĤ =T +Ŵ +N , and the expectation values of these three terms can be readily computed [40] in terms of the variational parameters according to where ρ ss is the solution of the matrix equation and K = iQ + i 2 R † R. In order to find the variational minimum of E[v, µ; Q, R] ≡ Ψ[Q, R]|Ĥ(v, µ)|Ψ[Q, R] = tr [Q, R] † [Q, R] + vR † 2 R 2 − µR † R ρ ss (B11) with respect to Q and R a tangent-plane method using the time-dependent variational principle (TDVP) in imaginary time was exploited [41]. This method proceeds as follows. Firstly, v and µ are selected. Then a value D as large as possible is chosen. An initial guess |Ψ[Q 0 , R 0 ] for the ground state results from random choice of Q 0 and R 0 . Also a tolerance η and a step size δ is selected. Set j = 0 and perform the following sequence of operations until the desired convergence is reached. 5. Set j = j + 1, unpack z j+1 into the two matrices Q j+1 and R j+1 , and repeat step (1) until convergence of the energy expectation values reaches the prespecified tolerance η.
After the above algorithm has terminated the cMPS corresponding to unit density (and at the rescaled interaction parameter) is obtained via the rescaling procedure described in the previous section. This method was used to calculate a cMPS representation for the Lieb-Liniger ground state across a range of interaction parameters from v = 0 to v = 1000. The value D = 14 was used throughout as the results so obtained are indistinguishable from the known exact solutions at v = 0, 2, ∞. For the other values of the interaction parameter the accumulated errors were estimated and found to be negligible.