Changing hbar convention produces different Fock probabilities

[tguillaume]

Before posting a bug report

• I have searched exisisting GitHub issues to make sure the issue does not already exist.

Expected behavior

Consider a circuit in the Gaussian backend consisting of nModes where the first mode is initialized in a thermal state with average photon number 0.01. I expect that:

(i) Computed Fock probabilities should not depend on the hbar convention
(ii) A thermal state (in Gaussian backend) should always be mixed, regardless of number of modes in circuit and hbar convention

Actual behavior

(i) Fock probabilities depend on hbar convention
(ii) Purity depends on hbar convention and size of circuit

Reproduces how often

All the time

System information

Python version:            3.10.10
Platform info:             Windows-10-10.0.19045-SP0
Installation path:         C:\Python\lib\site-packages\strawberryfields
Strawberry Fields version: 0.23.0
Numpy version:             1.23.5
Scipy version:             1.10.1
SymPy version:             1.11.1
NetworkX version:          3.0
The Walrus version:        0.19.0
Blackbird version:         0.5.0
XCC version:               0.3.0
TensorFlow version:        None

Source code

import strawberryfields as sf
from strawberryfields.ops import *

    
def SF_Thermal(nModes, nbar, hbar_val):

    sf.hbar = hbar_val
    prog = sf.Program(nModes) 
    eng = sf.Engine("gaussian")
    
    with prog.context as q:
        
        Thermal(nbar)              | q[0]

    # Run SF engine
    results = eng.run(prog)
    state = results.state
    
    return state


# Test with hbar = 2

print('Using default hbar=2.....')

out1 = SF_Thermal(nModes=1, nbar=0.01, hbar_val=2)
out2 = SF_Thermal(nModes=15, nbar=0.01, hbar_val=2)

print('1-mode state is pure: ' + str(out1.is_pure))
print('15-mode state is pure: ' + str(out2.is_pure))

equal = out1.fock_prob([1])==out2.fock_prob([1,0,0,0,0,0,0,0,0,0,0,0,0,0,0])

print('1 photon probabilities are equal: ' + str(equal))

print('%%%%%%%%%%%%%%%%%%%%%%%')
        
# Test with hbar = 1
print('Using hbar=1')

out1 = SF_Thermal(nModes=1, nbar=0.01, hbar_val=1)
out2 = SF_Thermal(nModes=15, nbar=0.01, hbar_val=1)

print('1-mode state is pure: ' + str(out1.is_pure))
print('15-mode state is pure: ' + str(out2.is_pure))

equal = out1.fock_prob([1])==out2.fock_prob([1,0,0,0,0,0,0,0,0,0,0,0,0,0,0])

print('1 photon probabilities are equal: ' + str(equal))

Tracebacks

No response

Additional information

I dug a little bit, and this issue may stem from the fact that the state is wrongly flagged as pure (some precision limit?) Because the SF object is considered pure, the method fock_prob uses the function thewalrus.twq.pure_state_amplitude rather than thewalrus.density_matrix_element. But it's unclear to me why the state is wrongly flagged as pure when hbar=1 (but correctly flagged as mixed when hbar=2). Also, this issue might potentially be similar to #488. However, unlike that issue, I only see the incorrect pure flag when I set nModes >= 15.

[CatalinaAlbornoz]

Hi @tguillaume, thank you for reporting this bug!

My colleague @sylviemonet has made two PRs to fix this: #746 and #747.

Unfortunately, there seem to be all sorts of problems with the automated checks that need to pass for these PRs to make their way into master. My recommendation in that case is that you make a Fork of the Strawberry Fields repository and add the changes in your own Fork. Hopefully this can help you bypass the failing checks and work with the fixes!

[tguillaume]

Thank you for looking into this so quickly! Setting hbar=2 seems to solve the issues for my purposes. But it's useful to see these PRs.

[CatalinaAlbornoz]

That's great to hear @tguillaume! Thanks for letting us know.