Actually, notwithstanding the no-go result, there is a position vector for photons; but it is singular in much the same sense that spherical coordinates are singular.
The issue can be best addressed by looking at the Wigner classification - but within the framework of symplectic geometry, rather than Hilbert spaces.
The real meaning and import of the no-go theorem is that the Wigner class which photons belongs to (which I term, below, the helical subfamily of the luxons, or the "helions") has no spin-orbit decomposition, so that the usual expressions for spin and position cannot be developed for helions. The symplectic geometry for the helion subclass shares many features in common with the symplectic geometry for magnetic monopoles (the latter which is discussed in LNP 107), except that the roles of the (q,p) coordinates are reversed.
Like all symplectic geometries, the coordinates for a symplectic leaf pair off into (q,p) pairs, and the helions have 3 Darboux pairs, which can be arranged (with a little manipulation and adjustment) into the usual form (𝐫,𝐏) for position and momentum. But unlike the Newton-Wigner position vector, 𝐫 is singular, when expressed as a functions of (𝐉,𝐊,𝐏,E) = (angular momentum, moving moment, momentum, energy). It has a coordinate singularity of the above-mentioned type.
The Wigner classes for the Poincaré group consist of the following:
(0) Homogeneous classes (unnamed by Wigner) (𝐏 ≡ 𝟎, E ≡ 0),
(1) Tardions (P² < αE²), where I will use α = 1/c² here and in the following,
(2) Luxons (P² = αE²), with 𝐏 ≢ 𝟎,
(3) Tachyons (P² > αE²).
where ≡ refers to conditions that hold on the symplectic leaf that characterizes the given representation,
(Most of what I describe here and below, by the way, also applies also to non-relativistic theory, by taking α = 0; except that the Luxons and Tachyons merge into a single unnamed family: the mass 0 representations for the Bargmann group - a class I named the "Synchrons". I also coined the term "Vacuon" for class (0).)
Over all classes, there are two invariants:
m² = M² − αP² = constant: mass shell constraint,
W² − αW₀² = constant: "spin/helicity shell" constraint
(the latter name being for lack of a better term),
where, for convenience, I will also use M = αE for "moving mass" here and below; where
(W₀,𝐖) = (𝐏·𝐉, M𝐉 + 𝐏×𝐊)
is the Pauli-Lubanski vector. For tardions, the second invariant reduces to
W² − αW₀² = m² S² (tardions only)
where S is the spin; and there are decompositions for:
Angular Momentum (Spin-Orbit): 𝐉 = 𝐫×𝐏 + 𝐒
Moving Mass Moment: 𝐊 = M𝐫 − 𝐏t + α𝐏×𝐒/(m + M)
where t may be arbitrarily selected, and 𝐫 adjusted accordingly. This can be inverted to express (𝐫,𝐒) in terms of (𝐉,𝐊), the result yielding that is known as the "Newton-Wigner" position vector for tardions.
For all families (1), (2), (3), there is a sub-family given by (W₀,𝐖) = (0,𝟎) Pauli-Lubanski vector - called "spin 0". For this class, too, there is a similar decomposition:
Angular Momentum: 𝐉 = 𝐫×𝐏
Moving Mass Moment: 𝐊 = M𝐫 − 𝐏t
and one can write
𝐫 = 𝐊/M + 𝐯t, 𝐏 = M𝐯
The indeterminacy in t - the same as what occurs generally for tardions - characterizes the trajectory for a worldline:
{ (𝐫,t) ∈ ℝ³×ℝ: 𝐫 = 𝐊/M + 𝐯t }.
For this subclass, 𝐖 ≡ 𝟎, and W₀ ≡ 0, which results as a secondary constraint.
For the quantized form of the symplectic decomposition, 𝐊 and M are represented by operators that do not commute with one another (their brackets are [𝐊,M] = iħα𝐏), so the quotient is only determined up to "factor ordering ambiguity" - which here means: up to an undetermined multiple of 𝐏, i.e. t. So the − 𝐏t term in the expression for 𝐊 already comes out automatically, in the quantized form of the classification.
For spin non-zero tardions, the expression for 𝐫 is 𝐫 = 𝐫₀ + 𝐯t, where 𝐫₀ is:
The Newton-Wigner Position Vector: 𝐫₀ = 𝐊/M − α 𝐏×𝐒/(m(m + M)).
The expression for 𝐒 is
Spin Vector: 𝐒 = 𝐖/m − αW₀𝐏/(m(m + M))
The most important features of the classes and subclasses are that:
(a) they are each characterized by the invariants and by what conditions apply to them,
(b) subsidiary invariants may also occur for subfamilies,
(c) the number of free parameters left over after removing the constraints from the set (𝐉,𝐊,𝐏,M) (or (𝐉,𝐊,𝐏,E)) is even,
(d) the remaining free parameters pair off into (q,p) variables - which is the essential statement of the Darboux Theorem,
(e) upon quantization, these pairs yield Heisenberg pairs - and this is where the Heisenberg relations come from.
For classes (1)-(3), the spin-0 systems have 4 constraints (0 Pauli-Lubanski vector) and, thus, 6 free variables, which combine to give you the 3 Heisenberg pairs (𝐫,𝐏). The extra parameter t can be normalized to 0 ... which is how it is normally done with the Newton-Wigner vector ... and so is inessential. (In the quantized version of the symplectic classification, one normalizes 𝐊/M − 𝐏t to the symmetric product ½(𝐊M⁻¹ + M⁻¹𝐊).)
For class (0) there are subsidiary invariants K² − αJ² and 𝐉·𝐊 that emerge, so that only 4 parameters at most are left free. The subclasses may have 2 pairs of Darboux coordinates (a "vacuum with spin and moment") or 0 (the "vacuum"); in the latter case the additional constraints are just K² = αJ² and 𝐊 ≡ 𝟎.
For class (1), the spin non-zero subclasses (i.e. where S² > 0) have 4 Darboux pairs. The fourth pair corresponds to the azimuthal component of angular momentum and the longitude and is normally quantized by the "m" number for spin states.
I won't describe class (3) in any detail, since it is a mess. The spin non-zero subfamilies all have 4 Darboux pairs.
Class (2), the Luxons, has 3 subclasses,
(a) spin 0: (𝐖, W₀) ≡ (𝟎, 0),
(b) helical: 𝐖 ∥ 𝐏, i.e. 𝐖×𝐏 ≡ 𝟎 (or equivalently, W² ≡ αW₀²), with 𝐖 ≢ 𝟎,
(c) general (or "continuous spin"), W² − αW₀² > 0
Note that the identity 𝐖·𝐏 = MW₀ follows from the definition of the Pauli-Lubanski vector, so from the constraint M² = αP², must follow that W² − αW₀² ≥ 0. Equality can only occur if 𝐖 ∥ 𝐏, which is why the constraints 𝐖×𝐏 ≡ 𝟎 and W² ≡ αW₀² are equivalent for Luxons.
The most important properties of these subclasses are that:
(a) the spin 0 subclass has only 3 Darboux pairs, which can be represented as (𝐫,𝐏),
(b₀) helicity (i.e. the component of 𝐉 parallel to 𝐏) is a subsidiary invariant for the helical subclass,
(b₁) the helical subclass, therefore, also has only 3 Darboux pairs(!),
(c) the continuous spin class has 4 Darboux pairs, and they are not represented by any spin orbit decomposition(!!).
Photons fall into the helical subfamily. The same is true for all fundamental particles ... in their true massless states before they are endowed with the appearance of mass by interaction with the Higgs. The reason for this is that weak nuclear charge is a multiple of left helicity for matter and right helicity for anti-matter and - by virtue of being a charge - it must first and foremost be an invariant property of the particle, which means the particles can only be helions or spin 0. That's why a Higgs mechanism is required for electroweak theory.
There is no spin-orbit decomposition, per se, for the helical subfamily, simply because there are only 3 Darboux pairs, rather than 4. Photon helicity is not spin! Classically, this corresponds to the fact (as Hehl has frequently pointed out) that the free electromagnetic field has no spin current and presents a symmetric stress tensor. For the interacting electromagnetic field (i.e. the field in a medium) the spin current would be proportional to 𝐃×𝐄 + 𝐁×𝐇, which is only non-zero if the constitutive laws for (𝐃,𝐁) versus (𝐄,𝐇) ... or (𝐄,𝐁) versus (𝐃,𝐇) ... are non-isotropic.
For electromagnetic fields inside a medium (like water) light goes slower than light speed in vacuuo, so the corresponding dressed quanta would fall into the tardion class and would have spin-orbit decompositions. In the quantized version of this, one would probably represent such "fields inside media" by effective Lagrangians, integrating out the external modes comprising the medium, and the dressed photons would acquire - in addition to the two values m = ±1 that come out of helicity - an extra mode for m = 0 and the dressed photons would "acquire mass". This is directly related to the very phenomenon in solid state physics that inspired the idea of the Higgs mechanism itself.
The question you're asking is: what about the helical subfamily? Since there are 3 Darboux pairs, then they do admit a quantization that has 3 Heisenberg pairs, notwithstanding the so-called no-go theorem. What it is really saying is that there is no spin-orbit decomposition and no analogue of the Newton-Wigner position operator that can be derived in that way.
However, there is a position operator, simply by virtue of the fact that the symplectic representation has 3 pairs of Darboux coordinates! The situation, like that of mapping coordinates for the sphere, is that at some point, the coordinates will go singular.
The sphere does not admit a globally non-zero linearly independent pair of vector fields on it. A similar situation occurs with the symplectic geometry that characterizes the helions. The similarity of its symplectic geometry to that of the magnetic monopole has been noted in the literature. The situation is analogous, except for the (q,p) reversal.
To write down a position operator, you can start by simply writing down a decomposition analogous to the "spin-helicity" decomposition for tardions:
𝐉 = 𝐫×𝐏 + η𝐏/M, 𝐊 = M𝐫 − 𝐏t ⇒ W₀ = ηP²/M, 𝐖 = η𝐏
the helicity being ηP/M = ηc.
It does, indeed, work - except that the 𝐫-𝐫 Poisson bracket relations acquire a deficit that is proportional to η. It's possible to adjust the definition of 𝐫 to eliminate this deficit, resulting in a bona fide Heisenberg pair set for (𝐫,𝐏), but the expression for 𝐫 will be singular in the components of 𝐉 and 𝐊. It's a coordinate indeterminacy, like that which the spherical coordinates (r,θ,φ) have at the poles when expressed as functions of Cartesian coordinates (x,y,z).
Would you like to see what it is? (Chomping at the bit, after all this long discussion, hmm?) Should I tell you? (Tease, tease!) No I think I'll end the reply here and leave it hanging...
Well, on second consideration...
They're in my notes somewhere and I'll have to look and check (and review it closely).
Here it is. There is no one solution. Instead, you need to pick a unit vector 𝐧. Then you can write down the decomposition:
𝐉 = 𝐫×𝐏 + ηP²/M 𝐧×𝐏×𝐧/|𝐧×𝐏|², 𝐊 = M𝐫 − 𝐏t + η 𝐧·𝐏 𝐧×𝐏/|𝐧×𝐏|².
This is obtained by taking the unadjusted 𝐫 and making an adjustment (𝐉,𝐊) → (𝐉 + δ𝐫 × 𝐏, 𝐊 + M δ𝐫) for a suitable δ𝐫 that fixes the deficit in the 𝐫-𝐫 brackets, while preserving (W₀,𝐖).
The representation goes singular in directions 𝐏 ∥ 𝐧, so you need a second 𝐧-vector to cover this region of the symplectic geometry. Two coordinate maps and regions, at minimum, are required to cover the symplectic geometry.
It's the same situation that occurs with magnetic monopoles, and η plays a role analogous to the electric-magnetic charge product.
To find 𝐫, you'll have to solve the above relations for 𝐫, which I'll leave to you and to the interested reader.
If you examine the little group for this subclass, using (𝛚,υ,𝛆,τ) to denote infinitesimal (rotations, boosts, spatial translations, time translations), you will find that it includes
(1) rotations 𝛚 ∥ 𝐏,
i.e. rotations along the axis collinear with 𝐏 or "helical" rotations,
(2) spatial translations 𝛆 ∥ 𝐏
combined with time translations τ such that ε = cτ,
(3) transverse boosts/rotations, 𝛚,υ ⊥ 𝐏,
combined with a compensating translations 𝛆,
such that 𝛚 = (𝐏/P)×υ/c and 𝛆P² + η𝛚 = 𝟎.
Properties (1) and (2) single 𝐫 out as a center-of-mass worldline, while property (3), which is just a "null boost" (combined with a translation perpendicular to both the boost and 𝐏), shows that there is a compensating relocation of the worldline, under a transverse boost.