Discrepancies between SF vibronic spectra calculation and the classical harmonic oscillator overlap approach

Photonic Software Strawberry Fields
1

Hello,

In order to strengthen my understanding of the Strawberry Fields vibronic spectra calculation using Gaussian Boson Sampling, I’ve implemented a simple vibronic spectra solver that computes the Franck-Condon matrix directly from overlaps of harmonic oscillator eigenstates. After weeks of debugging, I still cannot get the two methods to agree, even for very simple toy examples.

For the simplest case, a single mode is considered. Classically, this corresponds to a 1D harmonic oscillator. For a temperature of 0 K, the procedure is as follows:

  1. Compute ground and excited vibronic states, as 1D harmonic oscillator eigenstates), taken directly from Wikipedia.
  2. Build the Franck-Condon matrix F(m, n) by computing the overlap (i.e., inner product) between excited state m and ground state n; where m and n are integers ranging from 0 to a cutoff.
  3. Compute the spectrum. Being at 0 K, only the state n = 0 is occupied for the ground state, so the spectrum corresponds to a sum of Lorentzians (or Gaussians) centered at the transition frequencies (i.e., difference between excited state energy and ground state n = 0 energy).

Note that this procedure has been verified by other means (i.e., from comparing results with proprietary algorithms).

Now, for SF implementation, it should be relatively straightforward:

  1. Retrieve the Duschinsky matrix U_D and displacement vector delta. In this case, for 1D, U_D is a scalar with value = 1.
  2. Retrieve the interferometer parameters required for boson sampling.
  3. Acquire the spectrum.

Although it only takes a few lines of code here, the results are not the same. I’ve had issues with the fact that I was using different units than SF, but now I think that’s solved.

Item 1 of the above list is

U_D, delta = sfapps.qchem.duschinsky(
        np.array([[1]]),
        np.array([[1]]),
        np.array([ri]),
        np.array([rf]),
        wf,
        np.array([m])
    )

where ri and rf are scalars, representing the center of the harmonic oscillator eigenstates, wf is the excited vibronic state frequency and is given in cm^-1, and m the mass of the mode and is a scalar, given in atomic mass units. Note that if ri = rf, the two methods agree, whatever the frequencies and mass, but if the displacements differ, then the spectra differ too.

Finally, let me show an example of two different spectra. The parameters are the following:

  • m: 1 a.m.u
  • ri: 0.5
  • rf: 0.0
  • wi: 500 cm^-1
  • wf: 1000 cm^-1
  • T: 0 K
  • Hilbert space cutoff: 5
  • SF number of samples: 1000

Results:

image
image936×326 32.5 KB

The quantum (i.e., SF implementation) results on the left clearly disagree with the analytical (i.e., classical implementation) results shown on the right. If you have any suggestions, please let me know. I apologize for the long post but I wanted to give a lot of details. Thank you!

Andre

2

Hi @andre,

Welcome to the forum!

It’s great to see users like you who are interested in deeply understanding quantum algorithms. After an initial read of your post, I’m not exactly sure what the issue may be, but I have an idea.

In the case of computing the vibronic spectrum of a molecule, which is the context assumed by StrawberryFields, the parameters r_i and r_f denote the coordinates of the atoms in a molecule. These change between ground and excited electronic states. The case of a single mode corresponds to a diatomic molecule (a single atom doesn’t vibrate), in which case r_i and r_f contain the 3D coordinates of both atoms. When flattened, this is an array with six parameters.

So I believe that the problem may be due to the fact that the simple 1D harmonic oscillator case you are studying may not correspond to the context that the StrawberryFields functions assume, namely the vibrational structure of molecules. In particular, I’m not sure that it makes sense to consider r_i\neq r_f and to treat them as a one-dimensional parameter. In my mind this explains why you get correct results by setting r_i=r_f. In other words, I’m afraid you potentially may be misinterpreting the meaning of the parameters that enter the GBS algorithm.

In a vibronic transition, the initial and final states correspond to different states of the electrons in a molecule. Because the electronic density changes, the nuclei respond to the new electric potential and settle in a new equilibrium configuration. So the coordinates of the nuclei change (displacement). Additionally, because the electronic structure changes, so do the normal modes of the molecule (squeezing and rotation). The Franck-Condon factors are overlaps between states of different vibrational quanta of the ground and excited electronic states. But since the modes change, this effect has to be accounted for in terms of the Doktorov operator. The point is that all these quantities are defined in the context of molecules undergoing changes in their electronic states.

When we implemented these functionalities in StrawberryFields, we checked against calculations for the vibronic spectrum of Formic acid, as described in this tutorial. This is how we ultimately became convinced the calculations were correct. Data for this molecule can be found in the data module.

I know it’s more complicated than the 1D harmonic oscillator example, but I would recommend that you compare your own classical calculations to the results obtained with StrawberryFields for the molecules in the data module: (formic acid, pyrrole, and water). If the problem persists, then the issue may be something else.

Let me know if that helps! I’m happy to keep discussing if needed.

Juan Miguel

1 Like
3

Hi @jmarrazola ,

Thank you very much for the extensive answer. I’ll mark your post as being the solution to my problem, although I still haven’t solved the issue on my side.

The solution, I believe, would be to transform the molecular geometry displacements to normal modes coordinates, the system used in the classical implementation based on Franck-Condon factors.

Thanks!

Andre

1 Like