Order Operation Vibronic Spectra

[Hendrik]

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)

[isaacdevlugt]

Hey @Hendrik! Welcome to the forum

In the source code of sf the order is …

I think you’re referring to sf.apps.qchem.vibronic.VibronicTransition. Just wanted to check that your circuit matches the output of just using this operator (your result is on the bottom, and just using VibronicTransition is on the top):

There are some differences in peak intensity, but the locations of the peaks are the same. Here is my code:

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:
    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]

    Interferometer(U1) | q

    for i in range(2):
        Sgate(r[i]) | q[i]

    Interferometer(U2) | q

    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

a = plot.spectrum(e, xmin=-1000, xmax=10000)
a.show()

############################################
# now use SF's VibronicTransition operator #
############################################

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:
    qchem.vibronic.VibronicTransition(U1, r, U2, alpha) | q
    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)

a = plot.spectrum(e, xmin=-1000, xmax=10000)
a.show()

I try to code like in equation 19, but it is not the same spectrum anymore

One thing that I’m not certain on in your code is where you got the hard-coded values for Ud, w, wp, and delta. Maybe that’s why things don’t match what you’re seeing in the paper and in the demo?

[isaacdevlugt]

As a follow up, I tried switching alpha_neu with alpha in your code and the relative peak intensities aren’t the same — they match VibronicTransition more closely. Your code is on top and using VibronicTransition is on the bottom.

[Hendrik]

Thank you for answering my question
Yes I´m referring to “VibronicTransition” and try to code the same function but with changed order of operations. The hard-coded values are for a different molecule (sulfur-dioxide) and I use them because 2 modes are way easier to control than 7 (from the original) paper.

The differences in peak intensities are what matters actually. The positions are the same, because the resulting fock state (after the Doktorov operation) are weighted with the mode frequencies. The whole code is about evaluating the intensities. If one looks into the paper of Huh they change the order of operations so that a coherent state is produced which passes through the squeezing and rotation. In “VibronicTransition” the order is changed so that the squeezed vacuum is displaced at the end. This is fine theoretically.

What bothers me is that the intensities change when I try the rearange the order according to Huhs paper even though I changed the argument accordingly. It seems like strawberry fields doesn´t like displacing and then doing stuff. Could it work in a different backend?

[Hendrik]

As a follow up, I decided to dig deeper into the gaussian backend and squeezing. I tried to compare the photon number statistics of squeezed vacuum with the analytical results (see Wikipedia section about single mode squeezing). What I did is to prepare a squeezed vacuum state and sample. The probability to find the resulting state in a number state is then approximately given by the number of how often is the resulting state measured in that number state divided by the number of samples.

The anayltical result for the same probability is given by the formula (from the link above). Since the squeezed state there is given in number state representation, you have to “bra” it with the number state of interest and the sum collapes. The absolut square of this is then the anayltical probability.

import numpy as np
import strawberryfields as sf
from strawberryfields.ops import *
from collections import Counter
import math
from gmpy2 import isqrt


prog = sf.Program(1)
eng = sf.LocalEngine(backend="gaussian")
z = [2]
nr_samples = 1000

with prog.context as q:
    Sgate(np.abs(z[0]), np.angle(z[0])) | q[0]
    
    MeasureFock() | q
    
samples = eng.run(prog, shots=nr_samples).samples

c = Counter(list(samples[:,0]))
print("n=0 strawberry fields: " + str(c[0]/nr_samples))


def get_analytical(n, z_ana):
    r, theta = np.abs(z_ana), np.angle(z_ana)
    if(n%2 == 1): return 0
    n = int(n/2)
    cosh = 1/(np.sqrt(np.cosh(r)))
    tanh = (-1*np.exp(complex(0, theta))*np.tanh(r))
    bruch = (isqrt(math.factorial(2*n)))  /  (2**n*math.factorial(n))
    sqrt_probability = cosh * tanh**n * bruch
    return float(abs(sqrt_probability)**2)

print("n=0 anyltical: " + str(get_analytical(0, z[0])))

The resulting probabilities differ drastically. My guess why this might be the case is, that summing over all n converges super slowly and that this slow convergence rate leads also to problems within strawberry fields (especially for higher squeezing parameters).

[isaacdevlugt]

Hey @Hendrik! Just having to consult someone here and will get back to you shortly!

[isaacdevlugt]

Alrighty — I think this is a numerical (approximation, finite-precision, cutoffs, etc.) error . I can’t seem to find a reason why your code is incorrect. It’s probably something similar to what you mentioned in your last comment in that squeezed states cover many photon number states, which is hard to approximate well.

[Hendrik]

Alright, thank you @isaacdevlugt for you intense research . I mark this as the solution of the problem. Do you think there might be an option to quantify the intensitiy of this approximation-error?

Cheers

[isaacdevlugt]

My pleasure!

Do you think there might be an option to quantify the intensitiy of this approximation-error?

That’s a tough one. Since there are many compounding effects that occur when numerical issues arise it’s tough to quantify them on some order of magnitude. I suppose you can compare results that don’t have numerical issues to those that do !