quantum mechanics - Pegg-Barnett phase implementation does not seem to work

I attempt to monitor the phase of a wavevector $|\psi\rangle$.

As I found (e.g. here ), a matrix representation for the Pegg-Barnett phase operator in Fock base can be obtained as

$$\Phi=\sum_{m,n,n'=0}^s\frac{2\pi m}{(s+1)^2}\,\exp\left[\frac{2\pi i m(n-n')}{s+1}\right]|n\rangle\langle n'|$$

for a particle number cutoff $s$. Making the indices start at one instead of zero (for implementation in matlab), this becomes

$$\Phi=\sum_{\tilde{n},\tilde{n'}=1}^\tilde{s}\sum_{m=0}^{\tilde{s}-1}\frac{2\pi m}{\tilde{s}^2}\,\exp\left[\frac{2\pi i m(\tilde{n}-\tilde{n'})}{\tilde{s}}\right]|\tilde{n}\rangle\langle \tilde{n'}|$$ such that $|\widetilde{n+1}\rangle=|n\rangle$.

When I construct the wavefunction though of e.g. a coherent state with $\alpha=\sqrt{30}$; $\langle\alpha|\Phi|\alpha\rangle$ with $\tilde{s}=200$ gives me a phase ~2.7, whereas this coherent state is located in phase-space on the positive X-axis (so that phase zero would be expected.)

What went wrong?

NOTE ADDED: by trial and error, I found that we can obtain a phase prediction of zero if the summation over m also goes up to $\tilde{s}$ and in addition the reference phase $\theta_0$ is chosen $-\pi$ instead of zero. I don't understand why though, and if this is to be trusted as a measure of phase. At first sight, it does not seem to fullfil the small angle approximation $\phi\approx \frac{P}{X}$. Any insight is still welcome.

CODE USED with the original parameters

%% Constructing the initial state
alpha=sqrt(30);Nmax=200;
a=sparse(Nmax,Nmax);for it=1:Nmax-1,  a(it,it+1)=sqrt(it); end;adag=a';  
%annihilation and creation operators for use in the displacement operator.

vacuumpsi=zeros(Nmax,1); vacuumpsi(1,:)=1;
largepsiinit=exp(-0.5*abs(alpha)^2)*(expm(alpha*adag)*vacuumpsi);

%% The PB operator (tildes omitted)
s=Nmax;
PBphaseop=zeros(s);
thetazero=0;
for n=1:s
 for nprime=1:s
  for m=0:s-1

  PBphaseop(n,nprime)=PBphaseop(n,nprime)+(-
  thetazero+2*pi*m/(s))/Nmax*exp(1i*(n-nprime)*2*pi*m/(s));
   end
 end
end

%% Calculate expectation value of phase
phaseexp=largepsiinit'*PBphaseop*largepsiinit

Answer

What works quite nicely is defining the operator from eq.33 in https://arxiv.org/pdf/hep-th/9304036.pdf .

Also, it was to be expected that the PB-phase of a coherent state on the positive real axis is not zero for reference phase $\theta_0=0$, namely the eigenvalue spectrum is bounded by 0 and $2\pi$ and as a coherent state has a finite phase uncertainty, it crosses the cut in phase.

Aside from these issues, this PB-operator used has proven useful in monitoring the phase of a numerical wavevector and gives quite intuitive results. For further improvement, one can also iteratively update the used $\theta_0$ to avoid crossing of the phase-cut.