Modeling Non-Covalent Interatomic Interactions on a Photonic Quantum Computer

Non-covalent interactions are a key ingredient to determine the structure, stability, and dynamics of materials, molecules, and biological complexes. However, accurately capturing these interactions is a complex quantum many-body problem, with no efficient solution available on classical computers. A widely used model to accurately and efficiently model non-covalent interactions is the Coulomb-coupled quantum Drude oscillator (cQDO) many-body Hamiltonian, for which no exact solution is known. We show that the cQDO model lends itself naturally to simulation on a photonic quantum computer, and we calculate the binding energy curve of diatomic systems by leveraging Xanadu's Strawberry Fields photonics library. Our study substantially extends the applicability of quantum computing to atomistic modeling, by showing a proof-of-concept application to non-covalent interactions, beyond the standard electronic-structure problem of small molecules. Remarkably, we find that two coupled bosonic QDOs exhibit a stable bond. In addition, our study suggests efficient functional forms for cQDO wavefunctions that can be optimized on classical computers, and capture the bonded-to-noncovalent transition for increasing interatomic distances. Remarkably, we find that two coupled bosonic QDOs exhibit a stable bond. In addition, our study suggests efficient functional forms for cQDO wavefunctions that can be optimized on classical computers, and capture the bonded-to-noncovalent transition for increasing interatomic distances.

Introduction Materials and chemical modeling are considered to be among the most important applications of quantum computing.So far, most quantum computing algorithms have been applied to study short-range chemical bonds, primarily in small molecules, cf.[1] for a review on the topic.However, long-range non-covalent interactions are key to understand many properties of molecules and materials [2,3], motivating the developement of efficient and accurate models capturing these effects.Long-range forces originate from the electromagnetic interaction between electrically neutral atoms or molecules [4][5][6][7] and arise from the coupling of matter to the background quantum electrodynamic gauge field [8][9][10][11][12][13][14][15][16][17], and are at the core of the description of the properties of materials and macromolecules such as their structure [18], stability [19,20], dynamics [21][22][23] and coupling to a background gauge field [24][25][26].The inclusion of dispersion interactions can be done in an economical way by means of many-body methods [27][28][29][30][31][32][33][34][35][36][37][38][39][40], in particular through the so-called many-body dispersion (MBD) framework, whose accuracy was proven in the literature [39,41].In the MBD framework, drawing inspiration from the Drude atomic model, the response of valence electrons in atoms is assumed to be linear and this can be implemented through the introduction of quantum Drude oscillators (QDO), for which a harmonic potential is assumed between an effective electron (called drudon in this context) and an oppositely-charged nucleus, as illustrated in Fig. 1.QDOs define a compact coarse-grained quantum-mechanical model for dispersion forces, in which molecules are then defined as a collection of QDOs interacting through dipole-dipole interaction [42].Though simple, this system was shown to capture long-range phenomena even in large biomolec-FIG.1. Illustration of a system composed of a pair of QDOs.The nuclei are considered to be infinitely massive.The drudons interact with their nucleus through a harmonic potential, but with the other QDO through Coulomb interaction.xi denote the relative position of the drudon with respect to its nucleus in QDO i .r12 denotes the position of QDO 1 with respect to QDO 2 .
ular systems [43].By construction, the MBD framework relies on the dipole-dipole approximation to the Coulomb interaction, and therefore comes with the usual limitations of multipolar-type expansions, leaving aside any contribution coming from higher-order couplings.These limitations can be addressed via multipolar generalizations of the pairwise second-order perturbative approaches [35,36,40,44].We propose here, in the spirit of the Full Configuration Interaction (FCI) approach of [45], to study the MBD model with full Coulomb interaction between its constituents.Going beyong the quadratic dipolar QDO Hamiltonian of course comes at the cost of loosing integrability of the model, and this is where numerical methods, and in particular quantum computing approches become relevant.
Here, we show that quantum computing can be applied to efficiently and accurately model non-covalent interactions.Our study substantially extends the applications of quantum computing to atomistic modeling, by prob-ing the use of Noisy Intermediate Scale Quantum (NISQ) algorithms to the study of concrete quantum chemistry models, beyond the standard electronic-structure problem of small molecules [46][47][48].
We remark that the authors of Ref. [49] implemented a Variational Quantum Eigensolver for the simulation of an effective one-dimensional QDO model.There, the quantum harmonic oscillator Hilbert space, built as a Fock space, was truncated at some fixed level Λ, and mapped to the Hilbert space of a system of ⌈log 2 Λ⌉ qubits.The Hamiltonian was then mapped to a linear combination of Pauli strings.They studied there in great detail the case of a linear coupling between the QDOs, and explored the case of a quartic potential; However it is known on the one hand that the quadratic model has a known analytical solution, and we argued on the other hand that it is beneficial to explore QDO models beyond the multipolar expansion.Moreover the drudons were allowed to move either longitudinally or transversally with respect to the axis connecting the QDOs.However, as we will show, both of these models/configurations do not capture the most critical properties of the full-fledged threedimensional cQDO model, namely existence of a bound state, and smoothness of the binding curve [45,50].The present work adresses both limitations.On the other hand, our work also differs from a choice of hardware point of view.Indeed, in our case we choose the photonicbased continuous-variable quantum computing paradigm [51][52][53][54][55][56][57][58], guided by the intrinsic bosonic nature of QDOs, therefore effectively mapping a bosonic problem onto a bosonic hardware.The Fock space of a single harmonic oscillator is directly identified with the Fock space of one mode (one photon channel, in quantum optics language) of the quantum electrodynamic gauge field.In particular, the position and momentum of the drudon particle are directly identified with the position and momentum quadratures of the electromagnetic field, and a collection of N QDOs is then represented by a collection of 3N photon channels in an optical circuit.Similar ideas were successfully used to simulate bosonic Euclidean quantum fields [59,60] on a lattice.Concerning the quantum simulation of interacting harmonic oscillators in a first quantized framework, let us mention [61], where the authors encode simple supersymmetric quantum mechanical systems into a system of qubits, using again a hard cutoff for the bosonic degrees of freedom.In our case, and following [46,[62][63][64][65], an optical circuit composed of parameterized linear gates like beam-splitters, phase shifters, but also Gaussian gates implementing displacement and squeezing, and quartic gates (Kerr gates) is optimized through a Variational Quantum Eigensolver (VQE) type algorithm in order to accurately prepare the ground state wavefunction of the cQDO model.For the implementation, we leverage on Xanadu's strawberry fields quantum machine learning framework for continuous-variable hybrid quantum-classical optimization [62][63][64]66].Following [64] for state preparation, we probe our approach using strawberry fields simulator API.We success-fully derive the binding energy curve for a pair of QDOs, and study various properties of the corresponding ground state wavefunction.We note in particular the existence of a bound state, as already showed in [45,50], and we observe an interesting behavior of the quantum mutual information as a function of the interatomic distance, in correlation with the formation of the bound state as the two QDOs approach each other.
Definition of the model The Hamiltonian describing a system of N Quantum Drude Oscillators in three dimensions is given by: with the Coulomb interaction containing contributions from interacting drudon-drudon and drudon-nucleus pairs: and where x i denotes the relative position of drudon i with respect to its nucleus, and r ij = r i − r j denotes the position of nucleus i with respect to nucleus j .The usual approach consists in solving the theory in the multipolar expansion framework, in which the potential can be expressed as a power series in the inverse distance separating the two nuclei: with the following scaling behavior V n (x i , x j ) ∝ r −n−3 ij .The potential V 0 corresponds then to the dipole-dipole interaction, and is at the core of the Many Body Dispersion (MBD) model.V 1 corresponds to the dipole-quadrupole interaction, and V 2 to the quadrupole-quadrupole and dipole-octupole interaction.
One obvious limitation of the multipolar expansion, which assumes a large distance between the nuclei with respect to the typical separation between the drudons and their nucleus, is the lower bound it imposes on the interatomic distance.One can easily see that within the MBD (dipole-dipole) model by direct diagonalization of the quadratic Hamiltonian in terms of normal modes.In that case, one of the normal modes (the center of mass mode) develops a purely imaginary frequency at short range.Higher order physical effects are also neglected in the MBD model, motivating the study of the QDO model with full Coulomb interaction potential between its constituents.
Let us define the following dimensionless position and momentum operators associated to QDO i: in terms of which the Hamiltonian reads One can define the 3N creation and annihilation operators in terms of which the Hamiltonian reads Let us from now on restrict the problem to a pair of QDOs (N = 2) for concreteness, though the following developements carry to systems of N coupled QDOs.We denote by d the interatomic distance.
In order to reduce the complexity of the problem, let us define one-dimensional instances of the QDO model as follows: we restrict the movement of the two drudons to be along a common axis (directed by a unit vector êθ ) which form an angle θ ∈ [0, π/2] with respect to the vector r 12 connecting the two nuclei.We therefore have a family of one-dimensional models which can be obtained from the full-fledged 3d model simply by setting to zero the contribution from the oscillator modes belonging to the plane perpendicular to êθ .Let us denote by (X i , P i ) i∈{1,2} the remaining position and momentum degree of freedoms.As limiting cases, we obtain models in which the drudons are constrained to move either in the direction parallel to the axis separating the two nuclei (θ = 0), or perpendicular to the latter (θ = π/2).Those two models were studied in [49] in the dipole-dipole approximation.As mentioned in the introduction, in that paper the authors encode the states in the truncated Fock space of the system into the state of a set of qubits, and run VQE-type algorithms on IBQ quantum processors.However as we will see, the angle θ captures the competition between existence of binding (for small θ) and smoothness (for large θ), and interesting one-dimensional models actually sit at values of θ in the open segment (0, π/2).
In the case of a generic angle θ and interatomic distance d, the one-dimensional Coulomb potential reads: We therefore have a family of Hamiltonians parameterized by (θ, d) Alternatively, in terms of creation and annihilation operators, one has For the numerical simulations, we set ℏ = 4πϵ 0 = 1 as well as m i = q i = ω i = 1 for both QDOs.
Photonic circuit and Variational Algorithm From a molecular physics and long-range intermolecular perspective, we are mainly interested in knowing the ground state of the system (9).Among the various approaches, the hybrid quantum-classical variational algorithms were shown to be particularly efficient at capturing the properties of potentially complicated quantum states.Following [63,64], we apply a continuous-variable version of the VQE algorithm to extract the ground state of the QDO system.The ground state is obtained by optimizing the parameters of a parameterized optical circuit composed of linear multi-modes gates (two-mode beamsplitters and rotation gates), Gaussian gates (single-mode squeezing and displacement gates), as well as non-Gaussian gates (Kerr gates) implementing non-linearity.A generic multimode Gaussian state can be decomposed into a sequence of two-mode beam-splitters and single mode rotations, squeezing and displacement.Inserting a non-Gaussian operation ('non-linearities' in the classical machine learning language) in between each Gaussian transformation, allows to evade the strictly Gaussian realm, and provide enough flexibility to capture arbitrarily complicated wavefunctions by statcking multiple layers, increasing the expressivity of the ansatz space.As discussed later, one can justify a posteriori this specific choice of non-linearity by observing that our ground state is well-captured by cat-like states, whose preparation is known to require the use of Kerr interactions.Within a layer, the number of parameters in the model scales quadratically in the number of photon channels (due to the beam-splitters).We refer the reader to fig. 2 for an illustration of a single layer composing the architecture of the quantum neural network.In the procedure, the relative coordinate of the drudons with respect to their nucleus is directly identified with the position quadrature of the quantized electromagnetic field along the two channels of the optical quantum circuit.Overall, denoting by ω the set of all parameters, the circuit implements a unitary transformation U (ω) acting on an input reference state that we simply take to be the Fock vacuum state |0⟩ ⊗ |0⟩.The state prepared by the circuit is therefore given by |ψ(ω Once the ansatz state |ψ(ω)⟩ has been produced, one extracts the value of the energy in that state, given by ⟨ψ(ω)|H|ψ(ω)⟩ .Let us denote expectations in the state |ψ(ω)⟩ by angular brackets ⟨•⟩.To be specific, let us take our model (9), and let us denote by angular brackets the expectation in state |ψ(ω)⟩.One has by linearity of the expectation.On the second line one has to compute an expression of the form ⟨f (X 1 , X 2 )⟩.One therefore needs to extract the joint statistics of the position quadratures by preparation and measurements of the state |ψ(ω)⟩ in the quadrature basis, as summarized in alg.(1).Once the joint probability density ρ of (X 1 , X 2 ) in the state |ψ(ω)⟩ is known, one can compute where the integral above should be understood as a finite sum over a sufficiently refined grid G X × G X in the position quadratures plane.For the numerical simulation we take a linear grid G X composed of 500 points in the interval [ −6, 6].
Let us note that alternatively, and following [64], the simulator of strawberry fields actually provides direct access to the output statevector expressed in the Fock basis: The Fock space is actually truncated at a some fixed energy level, we choose that cutoff level to be 5.This choice of cutoff is not a real limitation and a larger value could be chosen instead, but 5 is sufficient, as was observed in particular using a FCI approach [45].The amplitude of a specific pair of the quadratures (X 1 , X 2 ) is then given by: in terms of the Hermite polynomials.The joint law of the quadratures in the state |ψ(ω)⟩ is then given by ρ After extracting as well the mean photon numbers ⟨a † 1 a 1 ⟩ and ⟨a † 2 a 2 ⟩, one obtains ⟨H⟩.We then define the cost function and update the parameters ω of the optical circuit in order to minimize that cost.This is summarized in alg.
(2).As a side remark, let us note that given the form of the Hamiltonian in terms of creation and annihilation operators (10), an alternative to the measurement of the position quadratures and photon number operators could be to perform a measurement in the coherent basis through heterodyne measurements.
Binding energy curve We gather here the results of the simulations.We focus on the case of 2 QDOs.In particular we study the profile of the binding energy as a function of the distance between the two nuclei, and make a few observation about the behavior of the entanglement entropy of the system.
We fix a grid G θ × G d ⊂ [0, π/2] × (0, 3.5], with card(G θ ) = 20 and card(G d ) = 200.For each pair (θ, d) in the grid we perform the continuous-variables VQE algorithm to extract properties of the ground state of the Hamiltonian (9).Let us denote by |ψ θ,d ⟩ the corresponding converged ground state.The binding energy is defined as the difference between the ground state energy of the system of interacting QDOs and the ground state energy of a system of uninteracting QDOs H 0 , i.e. with electric charge turned off: In fig. 3 we report the binding energy curve, namely the value of the binding energy as a function of the interatomic distance d, for a fixed value of the angle θ.We choose θ = 0.58 to illustrate most of the results.On fig. 3 (top) we observe a perfect agreement with a fit of Morse type [67][68][69][70][71]: The location of the bound state is given by d b ≃ 0.54, with energy −E b ≃ 0.46.The length scale s is given by s ≃ 2.75.The quality of the Morse fit actully increases monotonically with the angle, cf.fig. 3 (bottom), the quality being defined by the ℓ 2 -distance between the numerical result and the fit.By varying the value of the angle θ, we observe two different regimes.For small values of the angle, the curve deviates more and more from a Morse curve, and in the extreme longitudinal case θ = 0, becomes highly nonsmooth at short interatomic distances.This small angle regime is also characterized by the existence of a negative global minimum of the binding energy, hence by the existence of a bound state.On the other hand, as the angle increases, the ℓ 2 -quality of the fit becomes better.However, passed some value of the angle, the global minimum disappears together with the corresponding bound state.This competition between smoothness and existence of a bound state is illustrated in fig. 3 (middle), and of course signals the presence of a phase transition in the model where θ plays the role of control parameter.Below a critical value θ ⋆ of the order parameter, the system is characterized by the existence of a bound state, which disappear above the critical value of θ ⋆ .
Both regimes can be understood physically as follows: for very small angle θ ≪ θ ⋆ , and in particular for the longitudinal model θ = 0, as the two QDOs are getting closer and closer to each other, unstable configurations in which the two drudons are getting arbitrarily close to each other start appearing.Indeed, space being 1-dimensional and the two drudons moving along a common axis, as the two drudons get close to one another, the associated Coulomb repulsion component of the ground state energy diverges.
On the other hand, for large angle θ ≫ θ ⋆ , and a fortiori for the transverse model θ = π/2, two main configurations of the drudons may occur (thinking classically) depending on the relative position of the drudons with respect to the axis connecting the two nuclei.In the first configuration, the drudons are sitting on opposite sides.In that case, the dominant contribution to the energy between the two QDOs is the Coulomb repulsion between the nuclei.In the second configuration, the drudons are on the same side.In that case, in addition to the repulsive force between the two nuclei, one can also add up the repulsion force between the two drudons, leading to an even more repulsive scenario.Summing up, no binding can occur at θ = π/2, and by smoothness of the binding energy as a function of θ, this should also be the case in an open neighborhood of θ = π/2.The above observation suggests the following recipe.The longitudinal model predicts the existence of negative minima of the binding energy curve.It is however unstable due to the configurations of superposed drudons, as explained above.One can then regularize this 1d model by allowing for a non-zero angle θ, and slightly increase it until reaching a certain level of smoothness, that we have chosen here to be quantified by the quality of a Morse fit.The angle should however not be too large, smaller than the transition point beyond which the bound state disappears.This procedure defines a small range of models characterized by an angle θ ∈ [θ min , θ max ], as illustrated on fig. 3 (bottom).Finally, the QDO parameters usually represent the properties of real atoms, such as their polarizability and dispersion coefficients.Following for instance the paper of Jones et al. [72], the parameters could be fixed by fitting the dispersion coefficients predicted by the cQDO model to the reference (computed ab initio or measured experimentally) dispersion coefficients of the underlying atomic species.In addition, recently an optimized QDO parameterization has been proposed using only well-known dipolar properties of gasphase atoms [73].
Phase space representation In order to have a better intuition about the nature of the ground state of the system, let us turn to its representation in the quadratures phase space.In fig.4, we provide two representations of the system at different values of the interatomic distance (for the model at angle θ = 0.58).On the left side, from top to bottom, we represent the marginal Wigner function of each of the two QDOs for closer and closer distances.At large distance d = 3.16, the two QDOs are far apart and do not feel each other.This can be seen by the fact that the marginal Wigner functions are characteristic of purely Gaussian states.At distance d = 1.36, the two QDOs start feeling each other's presence.This is illustrated by the two marginal Wigner functions starting to be elongated, mostly along the position quadrature axis.Distance d = 0.82 corresponds to the distance of maximal entanglement entropy, as discussed a bit later in the paper.We observe in that case that the marginal Wigner functions are characterized by a signifi-cant spread in both position and momentum quadrature directions.Finally, at the bound state distance d = 0.54, the marginal Wigner functions appear to be mostly elongated along the momentum quadrature axis.This is clear since even though the configuration of the system has reached stability, characterized by the fact that the Wigner functions suggest a fairly neat localization of the drudons in space, the drudons being close to each other lead to an increase in the spread of the momentum density.Observe also the appearance of regions of negativity of the partial Wigner functions, signature of the nonclassicality of the state of the system [74], and witness of the entanglement in bipartite states [75].The right column of fig. 4 instead illustrates the joint probability distribution of the position quadratures.We observe again a Gaussian behavior at large distances, as it should for two independent harmonic oscillators, a maximal elongation at the maximal entanglement entropy point, and a more localized appearance at the bound state distance.Let us mention that a refined analysis around the distance of maximal position quadrature density spread shows that the unimodal joint position quadrature probability distribution develops a bimodal profile, with maximum of the original mode located around the origin (both drudons being in expectation centered at the locus of their respective nucleus), and maximum of the newly appeared mode shifted away from the origin.This bimodal feature of the joint position quadrature density, signature of tunneling, is very easily observable for a smaller angle theta, as can be seen on fig.5, which is consistent with the large curvature of the binding curve for such small angles, at the corresponding range of interatomic distances.After tunneling, the mode initially centered at the origin disappears, and remains only the shifted mode describing a stable bound state configuration in which each of the QDOs has developed a non-zero expected dipole moment, oriented in opposite directions.

Ground state wavefunction and entanglement entropy
In order to gain more intuition about the information possessed by one part of our bipartite system concerning the other part, let us introduce the partial density matrix associated to QDO 1 .The total state of the system, we recall, is expressed in the Fock basis as the pure state: This state is directly accessible in strawberry fields when using the simulator, and could be obtained on a genuine hardware by state tomography techniques [76].Note that the full knowledge of the state is not necessary for the optimization of the VQE parameters.We use this knowledge here solely to analyze the physics of the obtained ground state.Given the pure state describing the full system, the density matrix of the system is simply given by: The partial trace associated to QDO 1 is therefore given by: 2 )/2, diagonal section of the joint distribution and of the potential.Tunneling between the quadratic and Coulomb wells is maximal at the intermediate interatomic distance as can be seen from the bimodal joint distribution, and the neat appearance of two local vacua in the classical potential.This figure corresponds to the small angle θ = 0.17.
Since the state of the total system is pure, QDO 2 can be interpreted as purifying the system composed solely of QDO 1 .The two QDOs therefore have identical von Neumann entropy S(ρ 1 ), the entanglement entropy.The quantum mutual information of the system is therefore given by with the von Neumann entropy being defined as The profile of the entanglement entropy for different values of the angle θ as a function of the interatomic distance is provided in fig.6 (top).
In this figure we can see how for small values of θ the entanglement entropy presents a sharp peak which grad- ually turns into a plateau-like behaviour as θ increases (θ = 0.58 is the highest angle for which a peak is found).
As a comment, we observe fluctuations in the entanglement entropy, mostly for θ = 0.33.One can check that the distance between the density matrices in the sense of the Frobenius norm between subsequent values of the interatomic distance shows fluctuations in the same windows as those observed in the entanglement entropy.The fluctuations are still relatively small, but suggest that there exists a small neighbourhood of acceptable ground states in the ansatz space.This in turn suggests that these fluctuations could be mitigated by imposing a stronger convergence criterium in the VQE training loop.
A kernel fit of the entanglement entropy for angle θ = 0.58 is provided in fig.6

(middle).
This finding, together with the behaviour of the joint probability reported in fig.4, gives us an intuition of what happens to the ground state of the system along the binding process, allowing to open the VQE blackbox.
For d → ∞ the two QDOs will not interact and hence the natural state for them will be the vacuum state |0⟩|0⟩ (ground state of independent harmonic oscillators).On the other hand when d ∼ d bonding we see that the system is shifted towards an antisymmetric configuration where ⟨X 1 ⟩ = −⟨X 2 ⟩ and ⟨P 1 ⟩ = ⟨P 2 ⟩ = 0, which can be approximated by the bipartite coherent state |α⟩| − α⟩ with α = ⟨X 1 ⟩.In the transition region instead the system will pass through a strongly entangled state (hence the peak in von Neumann entropy), which is reflected in the position joint probability by the transition from a single mode to a bimodal distribution, as we discussed previously.This is naturally represented as a superposition of the form where N is a normalization factor, which can be viewed as a displaced Schrödinger cat state.For each value of the interatomic distance, one obtains the closest cat state to the ground state by optimizing the displaced cat state parameters through maximization of the state fidelity (overlap) between the two states.Although there is a slight deterioration in the fit for decreasing values of theta, this minimal ansatz can indeed approximate the ground state of the system up to a good degree for all the considered thetas: it manages to capture the bonding mechanism (states for which there is an entropy peak) with F ∼ 0.96.
A comparison at the level of Wigner functions between the original and fitted states for those bonding points is reported in fig. 7. Let us mention that the production of cat states typically involves the use of Kerr and squeezing gates [77,78], both of which being present in our quantum neural network architecture, Kerr gates being the source of non-Gaussianity and entanglement of the ansatz.A more detailed study concerning cat states in the context of long-range interactions will appear in future work.
Given the joint position quadrature density, another interesting quantity to consider is the quantum correlation between the position quadrature of the two QDOs, which appears to be highly correlated to the behavior of the entanglement entropy, cf.fig.6 (middle and bottom).Denoting again by angular brackets the expectation of an observable in the ground state, we define the correlation coefficient C(X 1 , X 2 ) by: Even though we do not have a neat theoretical understanding of it, we note that the maximum of the entanglement entropy, as well as the minimum of the correlation coefficient, converges towards the inflexion point of the binding energy curve as the value of θ increasing.A zeroth order line of argument could be the following: the binding energy is composed of a Coulomb repulsion and an attractive correlation energy contribution: The inflexion point d ⋆ , which in terms of the Morse fit ( 18) is located at d ⋆ = d b + log(2)s, can be interpreted as the transition between a large distance regime in which it is beneficial for the system to increase correlation in order to lower its ground state energy, and a short range regime in which the Coulomb repulsion becomes dominant.Further understanding of this phenomenon in the language of quantum phase transitions would be very interesting, but is left for a more theoretical investigation.Let us mention that the use of cat states as discussed above, for which the expression of the quantum mutual information turns out to be particularly simple, could also provide an approximate understanding for the presence of the peak.This will be discussed further in future work.

Conclusion
In this work, we showed that continuous variables photonics quantum computing is particularly adapted to the study of the full-Coulomb cQDO model, the fundamental degrees of freedom of the latter being bosonic in nature.Beyond standing as a proofof-concept that NISQ algorithms can be successfully applied to quantum chemistry problems beyond the usual approach using the second quantized formulation of the electronic structure problem for small molecules, we observed that the cQDO model for Many Body Dispersion interaction can be further simplified by reducing it to a single effective spatial dimension, while still capturing the existence of binding.The obtained groundstate wavefunction suggests that Schrödinger cat states can provide an economical ansatz space for systems exhibiting bonding-like behavior.This observation happens to also hold for the fully-fledged 3D cQDO model.Finally, note that Morse type binding curves are usually characteristic of covalent bonding in molecules and materials, suggesting that the cQDO model can potentially be used (up to a reparameterization of its defining parameters) to capture effects beyong those motivating its introduction, namely dispersion interactions.In particular the exponential repulsion behavior is reminiscent of exchange effects, which are a priori not built in the initial definition of the model.The extension of our study to many-body cQDO Hamiltonians are also relatively straightforward.
The quantum neural network architecture used in our study is quite general, and one could wonder whether a more specific architecture tailored to the cQDO model could be engineered.However, the many-body nature of the model makes it difficult to guess any a priori correlations between the various photon channels.A possible approach to get a priori intuition about a possible entanglement graph among the various modes could be, drawing inspiration from the MBD model [39], to first work in the dipole-dipole approximation.In this case the Hamiltonian is quadratic, hence quasi-free, and we can diagonalize it first to get the MBD normal modes.These normal modes do not constitute of course the eigenmodes of the full Coulomb model, but could provide a first coarse intuition about a possible circuit architecture tailored to the model.It would constitute an interesting direction, and we postpone that to possible later exploration.
Finally, concerning the possibility of scaling up to very large systems, one is obviously confronted to the fact that the bigger the system, the higher dimensional the position quadrature grid.
Drawing inspiration from the classical reinforcement learning literature, in which a similar Monte Carlo average of a random variable has to be obtained during the training loop, it could be interesting to see whether a very coarse estimate made out of very few samples could be enough for the VQE training loop to converge to a sensible result.

FIG. 2 .
FIG.2.One layer in the optical quantum circuit.The various gates, namely the beamsplitters, rotation, squeezing, displacement and Kerr gates are collectively parameterized by ω.

FIG. 3 .
FIG.3.On the top.Binding energy curve for the model at angle θ = 0.58, together with its Morse fit.On the middle.Binding energy curve for three different values of the angle θ illustrating the tension between models close to θ = 0 and models close to θ = π/2.For small angle (blue curve), the strong curvature prevents a good Morse fit, while for large angle (green curve), the transverse configuration of the drudons prevents formation of a bound state.The orange curve corresponds an intermediate angle exhibiting both binding and an excellent Morse fit.On the bottom.Blue curve: Quality of the Morse fit as a function of the angle θ, defined as the ℓ 2norm of the difference between the simulation and the Morse fit.Yellow line: Minimal angle above which the models are considered smooth.Red line: Maximal angle above which the models do not exhibit binding anymore.

45 FIG. 4 .
FIG. 4. From top to bottom: d = 3.16: QDOs are far apart, d = 1.36:QDOs start feeling each other, d = 0.82: entanglement entropy is maximal, d = 0.54: deep in the bound state.From left to right: Wigner distribution of the left QDO, Wigner distribution of the right QDO, joint position quadrature distribution of the two drudons.

FIG. 5 .
FIG. 5. From top to bottom: d = 3.16: QDOs are far apart, d = 1.75:Maximal interaction between the QDOs, d = 0.51: deep in the bound state.From left to right: joint position quadrature distribution of the two drudons, classical potential energyV (X1, X2) = V Coul (X1, X2) + (X 2 1 + X22 )/2, diagonal section of the joint distribution and of the potential.Tunneling between the quadratic and Coulomb wells is maximal at the intermediate interatomic distance as can be seen from the bimodal joint distribution, and the neat appearance of two local vacua in the classical potential.This figure corresponds to the small angle θ = 0.17.

FIG. 6 .
FIG.6.On the top.Entanglement entropy as a function of the distance d between QDOs for multiple values of θ.On the middle.Smoothing of the entanglement entropy vs. interatomic distance at angle θ = 0.58.On the bottom.Position quadratures correlation coefficient vs. interatomic distance at angle θ = 0.58.

FIG. 7 .
FIG. 7. Comparison between the states of the system on the entropy peaks (top) and the fitted ansatz in the same point (bottom).We plot the Wigner function sliced along the plane (X = X1 = −X2, P = P1 = −P2).