Hi @kareem_essafty — Thanks for the question!
I believe sf is giving you the correct output, in that a Gibbs state, i.e., a state of the form \exp(-\beta H) is only pure at infinite inverse temperature \beta = 1/(k_B T) where k_B is the Boltzmann constant and T is the temperature.
A quick solution to this problem is to start with a thermal state with a small mean number of photons and then extrapolate to the zero temperature limit.
2 Likes