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:
- Compute ground and excited vibronic states, as 1D harmonic oscillator eigenstates), taken directly from Wikipedia.
- 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.
- 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:
- Retrieve the Duschinsky matrix U_D and displacement vector delta. In this case, for 1D, U_D is a scalar with value = 1.
- Retrieve the interferometer parameters required for boson sampling.
- 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([]), np.array([]), np.array([ri]), np.array([rf]), wf, np.array([m]) )
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
wi: 500 cm^-1
wf: 1000 cm^-1
T: 0 K
- Hilbert space cutoff: 5
- SF number of samples: 1000
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!