I am running a circuit with two modes in which I prepare two logical 0 GKP states, using the GKP class with epsilon=0.125 on both modes, and have no gates acting on them.
I noticed that when I plot the Fock probability distribution of either mode using cutoffs larger than 40, I get a set of positive and negative probabilities, some of them with magnitudes larger than 1; which seem to be close to the cutoff independently of the cutoff I use. I attached figures to show this.
This behavior also shows up in circuits where I actually have gates acting on the modes. And I just noticed, that if I were to increase the epsilon value, this would start happening at greater cutoffs (observed with epsilon=0.15 at cutoffs larger than 45)
I am thinking that because of this, I can only trust simulations for GKP states up to a cutoff of 40 (when using an epsilon of 0.125); but nothing beyond that.
Would you think there is a way around this? Or is there some background knowledge to have in mind here?
Thanks for your question! This is a known issue stemming from instabilities on Gaussian operations (e.g. displacement, beam splitter and squeezing) for large cutoff values which start to become noticeable for values around 40 — as you have spotted as well. We are well aware of this issue and are currently working on methods to tame this numerical instabilities.
If you’d like to keep an eye on this, I invite you to open a new issue on The Walrus (where the aforementioned displacement, beam splitter and squeezing operations are implemented) so that you can follow the progress. You can also keep an eye upcoming releases of Strawberryfields and The Walrus.
StrawberryFields has been under maintenance mode for several years and it has now been fully deprecated. The Walrus has received very few updates in the past few years so it should be easy to check if this has been resolved.
You can check the latest commits here so hopefully you’ll find the one you need there!
and The Walrus 
.