# What is the current consensus on a well-behaving 'Quantum Phase Operator' (in Fock-state basis)?

I'm interested in a quantum operator that returns the phase of some quantized light in the Fock state basis.

That is, if I perform (in a density matrix picture) the expected value of this 'phase operator': $$\langle \phi \rangle = \mathrm{Tr}(\rho \phi)= \sum_{nm} \rho_{nm} \phi_{mn}$$, I obtain the 'phase' of the light.

In particular, if I give this operator a coherent state $$\alpha = |\alpha|e^{i \theta}$$ that $$\langle \phi \rangle_\alpha = \mathrm{Tr}(\rho_\alpha \phi) = \theta$$.

It seems as though there's been a lot of internal debate in the 90's (Although, I have vaguely heard, nowadays, that this is a fairly resolved question.) I have found two operators for phase in the Fock state basis: the Susskind-Glogower operator and the Pegg-Barnett operator. I've attempted to calculate the phase of a coherent state using these operators.

(Note, I'm aware that for a perfect coherent state it's much easier to calculate the phase; you can simply grab the imaginary log of the coherence term between $$|0\rangle$$ and $$|1\rangle$$. I am simply using the phase of a coherent state as a test for the phase operator, since the phase is very well defined for a coherent state.)

My attempt at calculating the phase using the Susskind-Glogower operator:

The Susskind-Glogower operator is not unitary, and therefore $$\langle \phi \rangle$$ is undefined. To remedy this, other operators related to this operator can be calculated, such as $$\cos(\phi)$$:

\begin{align} (N+I)^{-\frac{1}{2}} a + a^{\dagger}(N+I)^{-\frac{1}{2}} =(a^\dagger a+I)^{-\frac{1}{2}} a + a^{\dagger}(a^\dagger a+I)^{-\frac{1}{2}} \end{align}

But when I try to calculate value for a coherent state, I get results I do not expect. (Which I posted the code for at the end, if anyone is interested.) I've posted a plot for illustrating this (in which I graph the obtained output phase $$\langle \phi \rangle$$ as I very the input phase $$\theta$$ for the coherent state $$\alpha = |\alpha|e^{i \theta}.$$

The output corresponds to the operator for $$\langle \cos(\phi) \rangle$$, and it varies from $$-\pi/2$$ to $$\pi/2$$. Shouldn't this vary from $$0$$ to $$1$$ since it is the cosine of the argument? And either way, what use is it to find the cosine of the phase?

My attempt at calculating the phase using the Pegg-Barnett operator:

The Pegg-Barnett operator has been said to resolve the nonunitary issues of the Susskind-Glogower operator. The paper I'm using as a reference doesn't appear to explicitly write the operator in the Fock state basis, although they write what a Fock state is in this phase-operator basis:

$$|n\rangle = U|\phi(n)\rangle = \sum_m \exp(in(\theta + 2m\pi)/(s+1))$$

And thus I can write the phase operator as: $$|\phi(n)\rangle \langle \phi(n)| = U^\dagger |n \rangle \langle n| U$$ . But when I calculate the expected value of the phase $$\langle \phi \rangle$$ for coherent states, I get incorrect answers (that also vary as I change the value of $$|\alpha|$$. Here's a visual of how $$\langle \phi \rangle$$ varies as I change the input phase of $$\alpha = |\alpha|e^{i \theta}$$:

Here is the Matlab code I used if anyone is interested in going through it (Note: I just recently added a 3rd part that corrects for an error I made in the Pegg-Barnett method):

    % the following script calculates (the cosine of) the phase of the Susskind operator of an input that is a perfect coherent state
% this output phase is plotted as the input phase is varied

N = 30;
a = diag(sqrt(1:N),1);
at = diag(sqrt(1:N),-1);
iden = diag(ones(1, N+1));
perf_p = CoherentStateDensityMatrix(0+pi, 3, 30);
phaseOP = sqrtm(inv(at*a+iden))*a+at*sqrtm(inv(at*a+iden));
close all
phaseout = zeros(5, 1);
phasein = zeros(5, 1);
totaln = 30;
for ii=1:totaln
phasein(ii) = 2*ii*pi/totaln+pi/2;
perf_p = CoherentStateDensityMatrix(phasein(ii), 3, 30);
phaseOP = 1/(2)*((at*a+iden)^(-1/2))*a+at*(at*a+iden)^(-1/2);
phaseout(ii) = trace(perf_p*phaseOP);
end
scatter(wrapToPi(phasein), phaseout)

% the following script calculates the phase of the Pegg-barnett operator of an input that is a perfect coherent state
% this output phase is plotted as the input phase is varied

s = 30;
theta0 = 0;
Umat = zeros(s+1,s+1);
for nnn = 0:s
Umat(nnn+1,:) =  1/(s+1)*exp(1i*nnn*(theta0+2*(0:s)*pi/s+1));
end
phaseOP2= ctranspose(Umat)*diag(0:s)*Umat;
perf_p = CoherentStateDensityMatrix(0+pi, 3, 30);

%I'm assuming this matrix is unitary, consider switching with
%normal inverse if bugs.

for ii=1:totaln
phasein(ii) = 2*ii*pi/totaln+pi/2;
perf_p = CoherentStateDensityMatrix(phasein(ii), 1, 30);
phaseout(ii) = trace(perf_p*phaseOP2);
end
figure
scatter(wrapToPi(phasein), abs(phaseout))

%% Second Attempt at Pegg-barnet
%using: https://arxiv.org/pdf/hep-th/9304036.pdf
% the following script calculates the phase of the Pegg-barnett operator of an input that is a perfect coherent state
% this output phase is plotted as the input phase is varied

s = 30;
theta0 = 0;
rho_n = zeros(s+1, s+1, s+1);
rho_total = zeros(s+1, s+1);
for nnn = 0:s
eig_theta =  theta0+2*(nnn)*pi/s+1;
theta_n =  1/(s+1)*exp(1i*nnn*(theta0+2*(0:s)*pi/s+1));
rho_n(:, :, nnn+1) = ctranspose(theta_n)*theta_n*eig_theta;
rho_total(:,:) = rho_total(:,:)+ rho_n(:, :, nnn+1);
end

perf_p = CoherentStateDensityMatrix(0+pi, 3, 30);

%I'm assuming this matrix is unitary, consider switching with
%normal inverse if bugs.

for ii=1:totaln
phasein(ii) = 2*ii*pi/totaln+pi/2;
perf_p = CoherentStateDensityMatrix(phasein(ii), 1, 30);
phaseout(ii) = trace(perf_p*rho_total);
end
figure
scatter(wrapToPi(phasein), abs(phaseout))

function densitymatrix = CoherentStateDensityMatrix(theta, alpha, maxfock)

rows = 0:maxfock;
a = alpha;
single_row = exp(-a^2/2)*a.^rows./sqrt(factorial(rows)).*exp(1j*theta*rows);
single_column = ctranspose(single_row);
densitymatrix = single_column*single_row;

end

• A good summary of the debate you refer to is found in The Quantum Phase Operator: A Review by Stephen M. Barnett, John A. Vaccaro (CRC Press, 2007). – Emilio Pisanty Oct 8 '18 at 22:54
• and for a cutting counter argument to the Pegg-Barnett approach see Bergou, János, and Berthold-Georg Englert. "Operators of the phase. Fundamentals." Annals of Physics 209.2 (1991): 479-505. The core difficulties (v.g. periodicity) remain largely unresolved despite considerable recent (~10 yrs) literature on the topic. – ZeroTheHero Oct 8 '18 at 23:18

As I understand the situation, the current consensus is that the best way to measure phase is using a POVM, $$\int \frac{\mathrm d\phi}{2\pi}|\phi⟩⟨\phi| = 1 \quad\text{where}\quad |\phi⟩ = \sum_{n=0}^\infty e^{in\phi}|n⟩,$$ and that it's not particularly worth it to go looking for an explicit operator.

Using this convention, if you want the phase distribution for a coherent state $$|\alpha⟩$$, you just plot $$f(\phi) = |⟨\phi|\alpha⟩|^2$$ which looks like this for $$\alpha = r e^{i\theta}$$ with $$\theta=0$$ and $$r=0.5$$, $$1$$, $$2$$ and $$5$$ (in blue, yellow, green and red, resp.):

As you'd expect, the higher the mean photon number, the smaller the phase uncertainty.

I'm unable to locate good references for this formalism at the moment. I'll add them later if I do find them.

• I'm looking specifically for an operator that will tell me the expected value of the phase: $\langle \phi \rangle = \langle \psi | \phi |\psi \rangle$ (Where $| \psi \rangle$ is in the Fock state basis, ie $a|0\rangle + b|1\rangle +$...)You write $|\phi⟩ = \sum_{n=0}^\infty e^{in\phi}|n⟩$ but how is this used to calculate the phase? In the Pegg-Barnett operator this operator can be constructed $| \phi \rangle \langle\phi |$ over some truncated set of $|\phi_n \rangle$ and since $\phi_n = \phi_0 + 2 \pi/(s+1)$ is known, we can find an operator that can be written in the Fock state basis. – Steven Sagona Oct 9 '18 at 17:29
• Maybe to give more context: I am interested in finding $\langle \phi \rangle$ for states that closely resemble coherent states. So first I'll calculate $\langle \phi \rangle$ for coherent states (and make sure that $\langle \phi \rangle$ always matches $|\alpha|e^{i \theta}$ as I vary $\theta$ (In the matlab examples, I how that my calculated output phases aren't anywhere close to matching my input phases, even functionally), and then I'll give a slightly distorted coherent state and observe $\langle \phi \rangle$. – Steven Sagona Oct 9 '18 at 17:38

The Pegg-Barnett operator is in principle a good phase operator if the Fock-space is truncated (which is usually not a problem, to put a wavevector into the computer, you'll need a truncation anyway.)

So, what is the reason that you get this 'wavy'-pattern in your figure instead of a straight line? A first thing to keep in mind is the periodicity of the phase: so the best thing you can hope for is a sharp sawtooth function. In practice, a PB-operator (and basically any other attempt to calculate a phase) must assume an interval e.g. $$[0,2\pi[$$, implying a sharp cut at the edge. I see you did retrieve the periodicity. The reason your edges are smoothening is the following: A coherent state is not a phase eigenstate, so it has a finite variance of the phase. This becomes an issue when your coherent state is located close to the phase-cut: then part of the state (a tail of the distribution) actually crosses this cut! If you have for example a state $$|\alpha\rangle=|0\rangle$$, half of the state has a phase $$\approx0$$ and the other half a phase $$\approx2\pi$$, giving a total result of $$\pi$$!

As a practical approach to find a way out, you (as I briefly explained in an an answer to an old question of my own https://physics.stackexchange.com/a/380324/25292) can exploit that the PB-phase has a free parameter (the reference phase, corresponding to where you lay the phase-cut) and update this one dynamically until you get a proper result. Some MATLAB code I used (currently assumes wavefunction instead of density matrix, but the extention is straightforward as I have indicated in the code):

function [phase,Nit,phaseop]=refinedPBphase(psi)
%REFINEDPBPHASE calculate the phase of a wavevector psi (column). Starting point is
%the Pegg-Barnett formalism (see arXiv:hep-th/9304036). Through the
%iterations, it is possible to avoid nasty effects from the phase-cut by
%dynamically adapting the reference phase (as long as the state considered
%is at least somehow localized in phase space). Outputs the phase, numer of
%iterations used and the final phase operator.

Wouter V, 2018

Nmax=length(psi); indvec=1:Nmax;
thetazero=-pi; %Initiate reference phase

Tol=0.001; %Criterion when to be satisfied with the result.
updphase=0;
diff=Tol+1;

Nit=0; %Keep track of number of iterations
maxit=100; %maximal amount of iterations
while (diff>Tol)

%construct the phase operator
phaseop=2*pi/Nmax*exp(1i*(indvec.'-indvec)*thetazero)./(exp(2i*pi*(indvec.'-
indvec)/Nmax)-1);
phaseop(isinf(phaseop))=thetazero+(Nmax-1)*pi/Nmax;

oldphase=updphase;
updphase=wrapToPi(real(psi'*phaseop*psi)); %updated phase; change to wrapToPi(real(trace(phaseop*rho))) for mixed state

diff=abs(angdiff(updphase,oldphase));
thetazero=updphase-pi; %new reference phase, further away from our state of interest

Nit=Nit+1;

if(Nit>maxit)
warning('refPB:noconv',['refinedPBphase unable to reach tolerance after maxit ',num2str(maxit),' iterations']);
Nit=Inf; %Retrieve failure afterwards
break;
end

end

phase=updphase;

end


I had exactly the same requirement as you: the linear (sawtooth) relation between input and 'measured' phase for a coherent state, and it works perfectly. In practice, I have been applying this mainly to states with a 'banana' shaped W-function, so it should definetely work for your states. Convergence is usually reached without a problem as long as the density is not very close to zero. Of course, the state should not be spread out a full $$2\pi$$ for it to work (then there wouldn't be any phase anymore). You obtain the suitable operator as a result, which can be further used to calculate higher moments of the phase.

UPDATE: the method I described above can also be found at appendix B of https://arxiv.org/pdf/1902.09306.pdf

• Wow! Yes, this I think explains my issue. Your solution is similar to a similar issue I recently solved due to an issue with averaging phase-cuts. I get the idea that the variance of the phase is affected by the phase-cut - but I'm not sure exactly what's going on in the example Matlab code. Some questions: "what is angdiff?" Is this a function you've written that is not included? Maybe it's an array that's defined outside the code? (Try to answer these questions quick so I can give you the bounty!) – Steven Sagona Oct 15 '18 at 19:46
• Also, does this operator work in the general case where you don't know what the phase of the state "should be"? It seems like you're forcing it to agree with the correct answer, and I'm not sure how we can use that in situations where we don't know what the right phase is a priori. – Steven Sagona Oct 15 '18 at 22:13
• @StevenSagona angdiff is a matlab function, apparently from the Robotics toolbox mathworks.com/help/robotics/ref/angdiff.html . To answer your second question: for the function I've shown you don't need to know the phase in advance. As an initial value of the reference phase, I chose -pi, but this is arbitrary, it gets updated later on (I cannot exclude that for a very unlucky choice for the input phase like 'inputphase'=thetazero there will be an infinite loop, but this should not occur if the phase of your 'inputstate' is sufficiently random). – Wouter Oct 16 '18 at 8:41
• The function returns an operator in the end, which will be different for every phase you give as input. This operator will be the best behaving for all phases near the one of your input. But if you are only interested in the phase itself, you don't need the operators anymore. – Wouter Oct 16 '18 at 8:49

The eigenvalues of self-adjoint linear operator are real numbers. But a phase is not a number but either an angle (as you requested it), which is an infinite set of real numbers differing by an integral multiple of $$2\pi$$) or (the more standard definition) a complex number of modulus 1. Therefore no self-adjoint linear operator can have the phase as an eigenvalue, and your request is unanswerable. However, sine and cosine of the angle $$\theta$$ are commuting, self-adjoint operators, with spectrum in $$[-1,1]$$ and joint eigenvectors, and $$P:=\cos\theta+i\sin\theta$$ is a normal operator whose spectrum consists of the phases $$e^{i\theta}$$. This is the closest sensible answer to your question.

For the impossibility of quantizing an angle directly, and for the definition of the sine and cosine operators, see the discussion and references in my answer at Why doesn't the phase operator exist?

This is just a proposal. I have no knowledge whether it has been tried in the literature. (Having said that, I am aware of the phase operators that you are refering to. This proposal is not related to any of that.)

Since what you are interested in is the analogue of the phase of a coherent state, one may use the phase space analogy for coherent states represented in terms of Wigner functions. So, if one can compute the Wigner function $$W(q,p)$$ for your state, one can determine its centre of mass via the first moments in phase space $$q_0 = \int W(q,p) q\ dq$$ $$p_0 = \int W(q,p) p\ dp$$ Then one can construct $$\alpha = \frac{1}{\sqrt{2}} (q_0+ i p_0)$$ and from that compute the phase.

In my view that would give you the closest representation of the phase of a state that would corresponds to its equivalent for a coherent state. Unfortunately, I don't see a simply way to define a quantum operator that can perform the complete process.