I’m trying to simulate a sequential circuit process that uses PNR detection, a BS gate, and a displacement gate at each step. After each step I trace out the detected mode (mode 1) and use the state in the non-detected mode (mode 2) as the input state in mode 1 for the next step and repeat. For instance at the start of the sequence:
eng = sf.Engine("fock", backend_options={"cutoff_dim": 15})
prog = sf.Program(2) # 2 mode circuit
with prog.context as q:
Dgate(0.5) | q[1] # displace vacuum
BSgate(np.pi/6, 0) | (q[0], q[1])
MeasureFock() | q[0]
result = eng.run(prog)
I then take the output of mode 2 by doing:
dm = result.state.reduced_dm([1])
and then perform the next step in the sequence as so:
My question is: is using DensityMatrix(dm) | q[0] the correct way to go about doing this or is there a more proper way to accomplish this in strawberry fields?
For completeness: you can always concatenate programs, i.e., continue one program with another. For example, your code can be rewritten as
eng = sf.Engine("fock", backend_options={"cutoff_dim": 15})
prog1 = sf.Program(2) # 2 mode circuit
with prog1.context as q:
Dgate(0.5) | q[1]
BSgate(np.pi/6, 0) | (q[0], q[1])
MeasureFock() | q[0] # measurements reset quantum registers to vacuum!
result = eng.run(prog)
dm1 = result.state.reduced_dm([1])
# continue previous program by initializing
# a new program with the previous one
prog2 = sf.Program(prog1)
with prog2.context as q:
Dgate(1) | q[0]
BSgate(np.pi/3, 0) | (q[0], q[1])
MeasureFock() | q[0]
result = eng.run(prog)
dm2 = result.state.reduced_dm([1])
Remember that measurements reset quantum registers. If there is some structure in your circuit architecture, you can put all this in a for loop.
I forgot to ask this in the original post as well but I wanted to know if strawberry fields also allows one to calculate the probability of detection events when performing a measurement?
For example if one does
with prog.context as q:
MeasureHomodyne(phi=0) | q[0]
I know we can get the resulting x value with
result = eng.run(prog)
x = result.samples[0][0] # homodyne sample
However is there a way to also receive the probability of measuring this particular sample as well? I can see a way of doing this by using something like x_quad_values but I was wondering if there was a simpler way.