Hi all,

I’m trying to understand the example of calculating vibronic spectra with strawberry fields. (Vibronic Spectra)

I know the code is based on the paper of Huh et al. from 2014. In the paper the final order of the operators of the Doktorov transformation differs from what is coded.

The final order (equation 19) is 1. Displacement 2. Rotation 3. Squeezing 4. Rotation

In the source code of sf the order is 1. Rotation 2. Squeezing 3. Rotation 4. Displacement (like equation 14)

I know these two orderings are equivalent if you change your displacement argument.

I try to code like in equation 19, but it is not the same spectrum anymore. I modiefied the argument by multiplying the inverse of J.

Could it be that the change of the order of the operators causes numerical artifacts within sf?

This is the related code. Either use alpha_neu or alpha at 1. postition or 2. position respectively.

```
import numpy as np
import strawberryfields as sf
from strawberryfields.ops import *
from strawberryfields.apps import data, qchem
from numpy.linalg import inv
prog = sf.Program(2)
eng = sf.LocalEngine(backend="gaussian")
Ud = np.matrix([[0.9979, +0.0646],[-0.0646, 0.9979]])
w = np.array([989.5, 451.4]).T
wp = np.array([1178.4, 518.9]).T
delta = np.array([-1.8830, 0.4551]).T
T=0
t, U1, r, U2, alpha = qchem.vibronic.gbs_params(w, wp, Ud, delta, T)
J = np.diag(wp**0.5) @ Ud @ np.diag(w**-0.5)
alpha_neu = np.array(inv(U2 @ np.diag(np.exp(r)) @ U1) @ alpha)[0])
with prog.context as q:
if(1):
Dgate(np.abs(alpha_neu[0]), np.angle(alpha_neu[0])) | q[0]
Dgate(np.abs(alpha_neu[1]), np.angle(alpha_neu[1])) | q[1]
sf.ops.Interferometer(U1) | q
for i in range(2):
sf.ops.Sgate(r[i]) | q[i]
sf.ops.Interferometer(U2) | q
if(0):
Dgate(np.abs(alpha[0]), np.angle(alpha[0])) | q[0]
Dgate(np.abs(alpha[1]), np.angle(alpha[1])) | q[1]
MeasureFock() | q
s = eng.run(prog, shots=10000).samples
if np.any(t == 0):
s = np.pad(s, ((0, 0), (0, 2))).tolist()
e = qchem.vibronic.energies(s, w, wp)
from strawberryfields.apps import plot
plot.spectrum(e, xmin=-1000, xmax=10000)
```