Continuous matrix product states for coupled fields: Application to Luttinger Liquids and quantum simulators

A way of constructing continuous matrix product states (cMPS) for coupled fields is presented here. The cMPS is a variational \emph{ansatz} for the ground state of quantum field theories in one dimension. Our proposed scheme is based in the physical interpretation in which the cMPS class can be produced by means of a dissipative dynamic of a system interacting with a bath. We study the case of coupled bosonic fields. We test the method with previous DMRG results in coupled Lieb Liniger models. Besides, we discuss a novel application for characterizing the Luttinger liquid theory emerging in the low energy regime of these theories. Finally, we propose a circuit QED architecture as a quantum simulator for coupled fields.

A way of constructing continuous matrix product states (cMPS) for coupled fields is presented here. The cMPS is a variational ansatz for the ground state of quantum field theories in one dimension. Our proposed scheme is based in the physical interpretation in which the cMPS class can be produced by means of a dissipative dynamic of a system interacting with a bath. We study the case of coupled bosonic fields. We test the method with previous DMRG results in coupled Lieb Liniger models. Besides, we discuss a novel application for characterizing the Luttinger liquid theory emerging in the low energy regime of these theories. Finally, we propose a circuit QED architecture as a quantum simulator for coupled fields. Quantum Information and Quantum Technologies are providing both a new language and a new experimental landscape for the study of large quantum many-body systems. The study of entanglement in extended lattice models has made it possible to tackle the successful numerical renormalization group (NRG) 1 and the density matrix renormalization group (DMRG) 2,3 and provide them with a solid theoretical background based on the distribution of bipartite entanglement in 1D systems. This understanding made it possible to introduce new methods based on the matrix product states (MPS) formalism that allow studying both static 4 and time-dependent phenomena [5][6][7][8][9] , together with generalizations for critical 10 and two-dimensional systems 11 . As examples of the success of these methods we can remark the extremely good accuracy of DMRG studies in studying quantum phase transitions of lattice models 3 , as well as the success in the quantitative modelisation of novel experiments with cold atoms [12][13][14] , ions 15 and photonic systems 16,17 .
The above examples rely on lattice models. Sometimes, however, physics is best described via continuum 1D field theories. This includes 1D Bose-Einstein condensates under strong confinement and interaction, long Josephson junctions or nonlinear materials [18][19][20][21][22] . In the seminal work of Verstraete and Cirac 23 the MPS formalism was extended to treat continuum 1D quantum mechanical systems. The continuous matrix product state (cMPS) was formulated as a variational ansatz for obtaining ground states of continuum one dimensional and non relativistic fields 24 . More recently, the cMPS formalism has been used for tackling excited (1 particle) states 25 and 1 + 1 relativistic theories 26,27 .
In this work we introduce a natural extension of cMPS states to study coupled fields. We show that, thanks to the cMPS formalism, the interaction between fields does not have to be treated perturbatively, developing the appropriate algorithms to compute ground state properties. This is the main result of the paper and, by itself, it has potential applications to describe systems present experiments with interacting 1D Bose gases 28 , as well as the fermion-fermion or fermion-boson interactions, occurring in the so-called ladders 18,19 .
A well known peculiarity for one dimensional models is their low-Energy description as Luttinger liquids 29,30 . These liquids, no matter of the original model, are effective theories of bosonic character and are described by Sine-Gordon-like models. The parameters in the effective theory must be extracted from the original (microscopic) model. In the case of field theories, the parameters are given in terms of the ground state 31,32 . The second important result in this work is that cMPS can be used to derive those Luttinger parameters, both for the single and the coupled field case.
Finally, we also relate the coupled cMPS ansatz to the simulation of coupled quantum fields. We provide a recipe for building cMPS in the lab: engineering discrete quantum systems coupled to transmission lines, as in circuit QED setups. Those lab-layouts are nothing but prototypes for quantum simulators of field theories within cavity QED 33 .
The rest of the paper is organized as follows. The following section II is an (almost) self contained summary of the cMPS theory. Next, Sect. III is our first application. We use the cMPS, still single field, for obtaining the parameters in the Luttinger liquid theory, explained in III A and applied to the Lieb Linniger model, Subsect. III B. section IV explains our extension for coupled fields and V reports in our numerical results for coupled bosonic species. We finish, in section VI, commenting on the application of cMPS for constructing quantum simulators and summarize our results. In a similar way, we propose that coupled fields can also be constructed by coupling independent ancillas.

II. OVERVIEW OF CMPS
We review here the basics of continuous matrix product states (cMPS) for a single field. The cMPS are trial states for a variational estimation of ground states in onedimensional quantum field theories. We start by considering a quantum system described in second quantization by means of the field operatorsψ(z). According to the spin-statistics theorem, these operators must satisfy (anti)commutation relationsψ(z)ψ † (z ) ±ψ † (z )ψ(z) = δ(z − z ) according to whether they are fermions or bosons. Our system will be defined in a length L with periodic boundary conditions. The explicit form of the state can be written, as introduced in the seminal work of Verstraete and Cirac 23 ( = 1 is used through the text) where P denotes path-ordering (we follow the prescription in which, for the argument of the exponential, the value of z increases as we move to the right). Q(z) and R(z) are complex D × D matrices acting on an auxiliary Hilbert space H aux . The partial trace tr aux is taken over H aux . The suffix in 1 F emphasizes that it is the identity for the field. Finally, the state |Ω is the vacuum of a free theory,ψ (z)|Ω = 0 From now on, we will restrict ourselves to translational invariant setups in which the matrices R and Q become independent of z. The cMPS are complete 34 , i.e., any one dimensional quantum field can be casted in the form (1). This class of states can be obtained as the continuum limit of MPS, with bond dimension D. The bond dimension can be understood as a measure of the block entanglement. In one dimension the block entanglement saturates, thus, D is expected to be sufficiently small. If so, we are able to reach any quantum state with a relatively small number of variational parameters (2D 2 ). This combined with the variational method results in a very powerful technique for finding ground-states of onedimensional theories.
In a relativistic scenario, the block-entanglement has a UV logarithmic divergence. This can be understood since the ground state of a relativistic theory, also in 1 + 1, will contain zero-point fluctuations from all energy scales, which are the ones contributing more to the entropy 35 . A related argument due to Feynman is quoted in Ref. 26. The ground state will be dominated by the high energy contributions. As a consequence, in the variational procedure, the accuracy for describing the low-E sector is lost. Therefore, a cutoff must be introduced. Though challenging, the description of relativistic field theories has been succesfully described via cMPS introducing a regularization scheme 26,27 . Here we will face the most favourable case of one dimensional non-relativistic theories.
Inherited from their discrete countenparts, the cMPS is not unique but the gauge transformation Q → gQg −1 and R → gRg −1 leaves the state |χ invariant 36 . It turns out that the gauge is quite convenient. In this gauge, the cMPS state (1) can be rewritten as, and, K = K † Hermitean: which implies the gauge condition (3). The unitary operator U , in Eq. (5), is formally equivalent to a evolution in z-time for the fieldψ(z) and a D-level (auxiliary) system with Hamiltonian K. Field and ancilla are coupled via iR ⊗ψ † (z) − iR † ⊗ψ(z). The ground state is described in terms of the matrices K and R, i.e., in terms of an auxiliary zero-dimensional system. This suggests an holographic interpretation for the cMPS 37 . See Fig.  1(a) for a pictorial interpretation. It remains to provide operational rules for computing within the cMPS formalism. To be precise, we must be able to write any field observable χ|O(ψ,ψ † )|χ in terms of the matrices R and Q. As detailed in Ref. 38, the following relations are found: here, The Kronecker products in the ancilla space occurs because some identities, e.g., tr{A}tr{B} = tr{A ⊗ B}, have been used.
To avoid those products in the auxiliary space, the isomorphism |a |b → |a b * | is introduced. This allows us to map vectors in H aux ⊗ H aux into operators acting on H aux . This can be understood from the fundamental property The former also implies that operators acting on H aux are mapped into superoperators. Therefore, the action of T on a ket |ρ will be mapped into T [ρ], where T is a superoperator acting on the state (matrix) ρ. Under the isomorphism, it is straightforward to show that This is nothing but the dissipator governing a Linbladlike evolution d z = T for the irreversible dynamics of a system coupled to a reservoir. In this case, the role of the system is being played by the ancilla and that of the bath by the field (see Fig. 1(a) and the discussion above on the holographic interpretation). The Linbladian is a positive-semidefinite operator, T ≤ 0, having at least one zero eigenvalue 39 . With this at hand, Eq. (7) can be rewritten as χ|χ = l|e T L r = Tr e T L |r l| = Tr l · e T L r .
Here, l| and |r are the left and right eigenvectors of T (respectively) associated with its zero eigenvalue. We have assumed implicitly the limit L → ∞ where this eigenvalue yields the principal contribution to the exponential. In the third equality, the above introduced isomorphism has been used. Note that the zero eigenvectors of T , under the isomorphism, are mapped into the stationary solutions of the Linblad equation (left and right equations). Accordingly, the action of T into the bra l| can also be mapped into the action of a superoperator on a matrix: l|T ⇔ Q † l + l Q+R † lR. It is easy to see that, under the gauge (3), l * = 1 is a solution of the stationary Linblad-like dynamics (d z l * = 0). Combining all of this, we end up with χ|χ = tr(r * ) = 1 (14) where d z r * = 0. In a similar way, we can re-express the expectation value of any operator in terms of the steady-state solution * of the right Linblad equation With this we conclude our overview of the cMPS formalism. In the limit L → ∞, we will be concerned with the ground state energy density e 0 = χ|Ĥ(ψ, ψ † )|χ (wherê H is the Hamiltonian density operator). The latter can be computed by minimizing with the matrices R and Q as input and using the latter relations (and similar ones).
Once the minimization procedure has finished, observables can be computed with the same relations using the optimized matrices.

A. Bosonization
At low temperatures, a large class of one dimensional theories exhibit excitations of bosonic nature and their correlation functions are characterized by power laws. An interesting feature of 1D is that this class makes almost no distinction between bosons and fermions. Haldane 29,30 termed this class of theories Luttinger liquids. The bosonic nature of the low-energy excitations in 1D is due to the enhanced role quantum fluctuations acquire in low dimensional systems.
For a given microscopic model, the so-called bosonization prescription, consists in expressing the original degrees of freedom in terms of new fields which capture the collective behaviour characterizing the low-energy regime. For the case of a bosonic field, we will introduce the density-phase representation 18 whereρ(x) := ψ † (x)ψ(x) is the particle density field and θ(x) the phase field. Close enough to the ground state we can safely approximate the density operator bŷ where ρ 0 is the ground state density and the operator φ(x) characterizes the fluctuations over the ground state.
The commutation relations for bosonic fields will translate into a canonical commutation relation for theθ and φ fields B. Calculation for the Lieb-Liniger model We are going to apply the previous ideas to the Lieb-Liniger model 40 . The former describes a 1D nonrelativistic bosonic gas interacting via a repulsive zerorange potential The Lieb-Liniger model is exactly solvable by means of a Bethe ansatz. In fact, the solution shows that at low-energies, this model displays a Luttinger liquid behaviour 41 . An excellent agreement between the exact ground state energy density and the cMPS solution has already been provided 23 . Finally, note that this model conserves the particle number density. This quantity will represent a minimization constraint when finding the ground state numerically.
Following the bosonization scheme, the effective Hamiltonian describing the low-energy behaviour of the Lieb-Liniger model iŝ Hence, the low-energy regime can be completely characterized by means of two parameters (Luttinger parameters): the velocity v and the dimensionless parameter K. These, in turn, can be related to the ground state energy density e 0 (ρ) of the microscopic Hamiltonian (21). The corresponding relations are 31 It is possible to obtain asymptotic (analytic) expressions for the former parameters in terms of the dimensionless coupling constant γ = M c/ρ 0 . In Fig. 2 we compare those asymptotic limits (small and large repulsion, see Ref. 31) for v and K with (23) and (24)  energy density was calculated for up to twelve different densities. By interpolation, we constructed the continuous function e 0 (ρ) and the derivatives were calculated from it. Results show that for a moderately small bond dimension (D = 6), it is possible to match the predicted asymptotic behaviour up to high values of γ. Results for D = 2 are not shown for the sake of clarity (such a small bond dimension does not capture correctly the ground state of the Lieb-Liniger model).

IV. EXTENSION/GENERALIZATION FOR COUPLED FIELDS
The cMPS formalism can be naturally extended to treat a multi-species system. Let us consider a system of length L in which coexist q bosonic and/or fermionic particle species which are annihilated by the operatorsψ α , α = 1, ..., q. These operators satisfy (anti)commutation relationŝ where η αβ = +1 if at least one of the fields α or β is bosonic and η αβ = −1 if both fields are of fermionic nature.
The q-species cMPS state is defined as 24 here, matricesR α have been introduced for each one of the fieldsψ α and a single HamiltonianK for the auxiliary system. We will employ the tilde notation for the variational parameters of the multi-species cMPS state to differentiate them from their single field counterparts, cf. Eq. (1). The matrixQ is now defined as At difference with the single field case, a regularity condition must be imposed on theR α matrices in order that the expectation value of the non-relativistic kinetic energy, as computed with (27), will not become divergent. This condition reads In other words, the matricesR α inherit the (anti)commutation relation of their corresponding fields. With these ideas in mind we can extend the operational rules for computing expectation values with cMPS. For example, where the transfer operator (11) has been generalized tõ and translational invariance has been assumed for simplicity. Special care must be taken into account for systems where two or more fermionic species coexist. Let us discuss correlators like ψ † α (x)ψ β (y) . Expanding the pathordered exponential in (27), which acts on the vacuum |Ω of the field theory, and taking the annihilation operators to the right (normal ordering prescription) we obtain 24 (x > y) where the generalized transfer operatorT α deals with the exchange statistics For the case of bosonic systems,T α =T . The transfer operatorT governs the evolution of states in the ancillary space. Similarly to the single field case, this evolution can be mapped to a dissipative dynamics corresponding to the following Linblad quantum master equation Thus, we have again the picture of the ancilla coupled to a bath (the fields) by means of the operatorsR α .
Consider now the case of two bosonic fieldsψ 1 andψ 2 . We are interested in studying how the matrices (R α and K), which define the cMPS state in this two-species system, can be constructed from the matrices which characterize a single field. The simplest scenario considers two uncoupled fields. We have seen how the problem of computing expectation values in the ground state can be reduced to a dissipative dynamics going on in the auxiliary space -where the state of the total auxiliary system is described in terms of the density matrixρ. In the absence of a coupling between the fields, we should be able to recover our single field solutions. This is nothing but to demand the density matrix to be separable, that is,˜ = 1 ⊗ 2 . Both fields do not need to be identical, therefore, each of them will have associated a different set of matrices R α and K α which act on the corresponding auxiliary space A α . For simplicity, we assume that both A 1 and A 2 have the same bond dimension D 1 = D 2 = D. The total auxiliary space for the two fields will be the tensor product of the individual spacesÃ = A 1 ⊗ A 2 . Due to the tensor product structure, the bond dimension of the total auxiliary space is nowD = D 2 . The ancillas evolve independently according to the total Hamiltoniañ Similarly, each auxiliary system will couple to its quantum field by means of the matrices R α . The extension of these to the product space is Notice that the matricesR α satisfy the bosonic commutation relation [R 1 ,R 2 ] = 0 as it is demanded for a multi-species system (29). As desired, our construction let us recover the results for single fields. For instance, ψ † 1ψ 1 = tr(ρ * R † 1R 1 ) = tr(ρ * 1 R † 1 R 1 )trρ * 2 = tr(ρ * 1 R † 1 R 1 ) (as trρ * 2 = 1 for a density matrix). How is this picture modified in the presence of a coupling betweenψ 1 andψ 2 ? An arbitrary operator C, mappingÃ into itself, can be represented asC = i c i A i ⊗ B i , where A i acts on A 1 and B i acts on A 2 . Therefore, this defines the most general structure for the matricesR α andK. Those general matrices must satisfy the regularity conditions (29), which complicates their construction.
A possible solution is the following. We use the intuitive interpretation for the cMPS in terms of a systembath, see Fig. 1(b) and the discussion on the holographic interpretation below Eq. (13). Starting from the decoupled solution (35), (36) and (37), we switch on the coupling adiabatically and expect that our solutions will start to modify. This is depicted schematically in Fig.  1(b). Here, as one introduces the coupling between the physical fields, the individual auxiliary spaces will also start to interact. Inspired by this procedure, we propose the following construction in the presence of a coupling. First of all, the matricesR 1 andR 2 will continue to be described by (36) and (37) respectively. In this way, we guarantee that they commute, satisfying (29) trivially. In order to render the state non-separable, the matrixK is written in a general way but containing the uncoupled solution as a limit (35). This is done as follows In order to keepK Hermitean, we will demand that the matrices Z (p) i are Hermitean too. The number P of pairs of Z matrices is arbitrary. In principle, we will expect it to grow with the strength of the coupling.
We have seen for the single field case, that the cMPS ansatz is able to map the properties of a continuous onedimensional field theory by means of 2D 2 variational parameters (with of course, a relatively small bond dimension D). Doubling the number of fields, as well as introducing P pairs of the already defined Z matrices, increases the total number of variational parameters to (4 + 2P )D 2 .

V. TWO-SPECIES BOSONIC SYSTEM
The system we have in mind to test the cMPS method for coupled fields is a two-component bosonic system. Binary systems of this kind (as well as bosonic + fermionic mixtures) are Luttinger liquids with a rich phase diagram 42,43 . We will consider two Lieb-Liniger gases with a density density coupling. This is described by the following Hamiltonian: In order to obtain the low-energy behaviour of this model we will use the bosonization technique introduced in Sect. III A. As already explained, this consists in rewriting the bosonic fields in terms of the collective fieldsθ α andφ α which characterize the bosonic lowenergy excitations. Hamiltonian (39) conserves the individual particle densities, [Ĥ,ρ α ] = [Ĥ,ρ] = 0 (α = 1, 2). Therefore, we can fix these two densities as minimization constraints in the cMPS procedure. In Eq. (19), we have considered the lowest order term in a harmonic expansion of the density operator. A more careful treatment 18 shows that, the correct expansion for the density operator in terms of the fieldφ α is of the formρ α (x) = [ρ 0α − ∂ xφα (x)/π] p e i2p(πρ0αx−φα(x)) . Our former simplification is justified due to the fact that, at long distances (low-energies), the phase terms oscillate very fast and will average to zero upon integration. In performing the bosonization, we must retain the most dominant terms at low-energies. For the case of our inter-species coupling, this supposes to consider also the first harmonic p = 1. This leads, at low temperatures, to a coupling contribution of the form: 1/2π dx [2g x ∂ xφ1 ∂φ 2 + g c cos(2(φ 1 −φ 2 ) + πδx)] (with δ = ρ 01 − ρ 02 ). Of particular interest for us will be the case of equal filling ρ 01 = ρ 02 (δ = 0). Species 1 and 2 in the low-energy effective Hamiltonian can be decoupled by introducing the normal modesφ . In terms of these we have that the low-energy excitations of (39) can be described by the effective Hamiltonian Similarly to the single field case, the Luttinger parameters v ± and K ± can be related to the ground state energy density e 0 (ρ + , ρ − ) (as a function of the normal densities) of Hamiltonian (39).
Coupled species have been thoroughly studied 42,43 . In this work we study coupled bosonic species described by (39). The range of parameters considered coincides with the one in Refs. 44 and 45 where a DMRG study is reported, hence a direct comparison is possible.
In Fig. 3 (a) we plot the ground state density-density correlations |∆ρ 2 | = | ρ 1 (x)ρ 2 (x) − ρ 1 (x) ρ 2 (x) | as a function of the interspecies coupling g for different values of P . Both, the repulsion strength c and the bond dimension D are kept fixed (c = 1.5 and D = 6). As expected, P = 0 renders the state separable and no correlations are observed. Making P = 0 the correlations between the two fields build up. They grow with the coupling strength. In this range of parameters, P = 2 seems to be sufficient for account with the physics.
The ground state energy density as a function of g is shown in Fig. 3 (b). Only the last term of (39) depends on g, therefore, (for a fixed value of c) the first two terms yield a constant contribution. The case P = 0 yields a mean field treatment where the interaction is replaced by g ρ 1 ρ 2 . It was already mentioned that [Ĥ,ρ 1 ] = [Ĥ,ρ 2 ] = 0. Thus, with P = 0, the energy as a function of g has a linear dependence with slope ρ 01 ρ 02 . Including Both as a function of the interspecies coupling strength g for equal densities ρ01 = ρ02 = 0.63 and c = 1.5. The ancilla space for each field has bond dimension D equal to 6. We couple the ancillas by means of P pairs of Z matrices. Results are shown for P = 0 (circles), P = 1 (triangles), P = 2 (squares) and P = 3 (diamonds). Inset: ground state energy density for a single Lieb-Liniger chain as a function of c (D = 6). The vertical line denotes the value of c at which we are coupling two fields.
quantum correlations (P = 0) the energy is no longer a linear function of the coupling, as seen in 3 (b).
As it was already discussed, the low-energy description of our model is characterized by the Luttinger parameters (41) and (42). In particular, the difference in normal velocities v ± yields the charge-spin separation -a typical experimental characteristic in mixtures. In Fig. 4 we compare for both the v ± and K ± our cMPS results with the DMRG values extracted from Refs. 44 and 45. Let us remark the excellent agreement in K + and v + and the minor discrepancy in K − or v − . The very small differences could be attributed to a number of issues: small bond dimension in the cMPS (D = 6) to be compared with the DMRG (several hundreds), or the fact that the DMRG theory is discretized an the cMPS is fully continuous.

VI. QUANTUM SIMULATION OF COUPLED CMPS
There exist two approaches towards the quantum simulation of continuous or discrete field theories 46 . The conventional one consists on taking a flexible quantum system, such as a Bose-Einstein condensate, ultracold atoms in an optical lattice or a superconductor, and working with it to implement the full field theory, or an approximate version of it, in the experiment. This "analogue" quantum simulator therefore evolves and equilibrates as the original model dictates and all observables may be directly studied on the experiment itself.
A second possibility for quantum simulation arises from the physical interpretation of cMPS. The idea is that there exists a mapping between a continuous Matrix Product state and a physical process operating on a small quantum mechanical object. This mapping between states and channels was already evidenced for discrete MPS 47,48 and has been recently generalized for cMPS 33 , by means of their physical interpretation in The ancillas consist on two cavity-qubit setups and the fields are the input and output of the EM field. The Hamiltonian of the cavity-qubit setups simulatesK and the cavity operators couple to the external field, that is, they correspond to the matricesRα. terms of a system (the ancilla) coupled to a bath (the field). The beauty of this mapping is that it is quite general and applies to a variety of quantum optical systems. The prototypical system is an atom-cavity setup (the system or ancilla in the language of this paper) that interacts with external input and output fields through the bath (the field in cMPS) in this case the electromagnetic field. However, any other quantum discrete system coupled to an outer field, where different order correlations of the latter can be measured, such as circuit QED 49,50 would do the job.
Let us now summarize the proposal in Ref. 33. The atom-cavity system is described through the well known Jaynes-Cummings (JC) model, H ancilla =Ĥ JC = Ωâ †â + σ +σ− + g(â †σ− + h.c) (43) hereâ (â † ) are bosonic annihilation (creation) operators describing the main stationary mode of a cavity. The atom (with two relevant states splitted by ) is coupled to this fundamental mode of Ω-frequency with a strength g. Theσ − (σ + ) are lowering (raising) operators for the two-level system. The atom-cavity is coupled to an EM-environment, that in second quantization is given by the free Hamiltonian,Ĥ EM = dω ωb † (ω)b(ω). Taking an interaction picture with respect to the EM field, the system-bath (ancilla-field) coupling can be written aŝ Limiting the integration region to frequencies ω near Ω we can safely assume the RWA. Also, assuming a pointlike interaction in space, the coupling function is flat in momentum. It is customary to introduce the time dependent operatorsÊ + (t) = i/ √ 2π dω e −iωtb (ω) and Hermitian conjugate. They correspond to the electric field components of the EM field. In this way, we can finally write the total Hamiltonian aŝ The electric field operators can, in turn, be decomposed into in-out components 51 . The in component corresponds to the field that impinges on the system while the out component consists of a reflected part plus a radiated one due to the interaction of the EM field with the system. If we take the in state of the EM field to be the vacuum, it can be shown 52 that the evolution governed by (45) can be reduced to that of the non-Hermitian Hamiltonian This is the same kind of evolution which generates the cMPS ansatz (see Eqs. (1) and (6)) once we trace over the degrees of freedom of the ancilla. We thus make the following identification: While we do not have control over R, we can modify the variational parameter K by properly tuning the couplings (Ω, ,g) of the cavity-atom system. The continuous fieldψ(x) will map into the output field operators of the electromagnetic field:ψ(x) =Ê + (t). Being the EM field in a cMPS state, computing expectation values of operators will translate into measuring correlations of the EM field itself, i.e., measuring the normalized correlation functions g (1) (t, t ), g (2) (t, t ) g (1) (t, t ) = Ê − (t)Ê + (t ) Ê− (t)Ê + (t) Ê− (t )Ê + (t ) (48) g (2) (t, t ) = Ê − (t)Ê − (t )Ê + (t)Ê + (t ) Ê− (t)Ê + (t) Ê− (t )Ê + (t ) (49) and higher orders depending on the model we wish to simulate. Following our previous identification, the correlators Ê − (t)Ê + (t ) and Ê − (t)Ê − (t )Ê + (t)Ê + (t ) map to ψ † (x)ψ(x ) and ψ † (x)ψ † (x )ψ(x)ψ(x ) respectively. It was shown numerically that the atom-cavity setup could simulate the Lieb-Liniger model giving correlations acceptably well 33 .
With this work at hand, our proposal has also a natural realization. In our case, we envision two superconducting cavities interacting each one with one [53][54][55] or several superconducting qubits (Fig. 5): That are coupled to different baths (different fields) througĥ H coupling (t) = α √ κ α dωâ † αbα (ω)e −iωt + h.c. (51) Therefore in our case (α = 1, 2), the identifications are the following,R α = √ κ αâα (52) and with and Z (1) and Z Note that in Sect.IV we demanded that the matrices Z should be Hermitian. Eqs. (55) and (56) can always be brought into a sum of tensor products of Hermitian operators. The former equations are of the form C = A ⊗ B † + A † ⊗ B. We can split any operator in terms of its Hermitian components. In the case of A, the decomposition reads: A = A r + iA i (and similarly for B). Here, A r = 1/2(A † + A) and A i = i/2(A † − A). It is straightforward to show that C can be rewritten as: Finally, as for the single field, EM field correlations need to be computed. In addition, cross-correlations, for instance, Ê − i (t)Ê + j (t ) will be necessary. In circuit QED this is possible as reported in the literature [56][57][58][59][60] .

VII. SUMMARY AND CONCLUSIONS
In this work we have proposed an extension of continuous Matrix Product States (cMPS) to study the ground state properties of 1D coupled fields. Our treatment has been confronted to previous DMRG numerical results, showing good convergence properties even for moderately large coupling strengths. Finally, we have discussed how it could be possible to realize computations for coupled fields using a quantum simulator to implement the cMPS ansatz and optimizing over the ansatz parameters 33 . We believe that extensions of this ansatz, together with new ideas on time evolution and the study of quasiparticle excitations 36,61 can provide a valuable insight on existing experiments with 1D atomic Bose-Einstein condensates 21,62 .