First, the Newton-Wigner "operators" are actually not quantum at all, but classically-rooted. What's referred to as "operators" are just quantized version of the classical quantities. Historically, however, the ideas first cropped up in the setting of quantum theory - thus the mistaken attribution of them being "quantum" in essence and the continual mistaken attribution of them being referred to as "operators", rather than as quantized forms of classical quantities.
They may be obtained by solving the equations for the angular momentum (or rotation generator) $$ and boost generator $$ for the internal angular momentum $$ and position $$. The position vector $$ is what you're asking for.
The equations apply to all sub-light systems in Relativity and to all non-zero mass systems in non-relativistic theory:
$$ = × + ,\quad = M - t + α\frac{×}{m+M}.$$
They are connected to the Pauli-Lubanski scalar $W_0 = ·$ and 3-vector $ = M + ×$. In Relativity, they together form the Pauli-Lubanski 4-vector $\left(W_0, \right)$. Explicitly:
$$W_0 = ·,\quad = M + α\frac{×(×)}{m+M}.$$
Together, they satisfy the identities:
$$M^2 - α ||^2 = m^2,\quad ||^2 - α {W_0}^2 = m^2 ||^2.$$
The $\left(m,||^2\right)$ eigenvalues define the representation family, when quantized.
For the non-relativistic version $α = 0$ and $M = m$, with $m$ being the mass. For the relativistic version $α = 1/c^2 > 0$, with $m$ being the rest mass and $M = m\sqrt{1 + α||^2/m^2}$ the "relativistic" or "moving" mass. In the contemporary literature, $M$ is normally written as the "total" energy $E$, via the relation $E = M c^2$. However, $E$ has no non-relativistic limit, i.e. as $c → ∞$, so it's actually more natural to write it in terms of $M$, instead.
Finally, $$ is the momentum 3-vector. In Relativity, you can treat $(,M)$ as the components of the momentum 4-vector, just as well as you could $(E,)$. It's the same 4-vector, apart from the scaling of the components.
Solving for $$ and $$, one gets:
$$
= \frac{}{m} - α\frac{}{m}\frac{W_0}{m+M},\quad = \frac{}{M} + \frac{}{M} t - α \frac{×}{mM(m+M)}.
$$
To quantize this, and to make the representation applicable uniformly, the operator forms should be written as $\hat{} = -i ħ ∇$ and $\hat{M} = m + α i ħ ∂/∂t$. For $E$, you'd thus write $\hat{E} = m c^2 + i ħ ∂/∂t$, instead of $\hat{M} = α i ħ ∂/∂t$ and $\hat{E} = i ħ ∂/∂t$. The replacement of the latter set by the former set is the Foldy-Wouthuysen transform in disguise. With the altered form of $\hat{M}$, you can pass directly over to the non-relativistic limit by taking $α = 0$. With the unaltered forms, you end up getting representations for $\hat{M} = 0 = m$.
In the quantized form, the operator-ordering issue arises with $(,M)$, since $[\hat{},\hat{M}] = i ħ α \hat{}$. So, effectively the ordering ambiguity in $\hat{}/\hat{M}$ gives you an undetermined multiple of $\hat{}$ - alongside the arbitrary $t$ multiple for $\hat{}/\hat{M}$ which positions you at an undetermined point on a worldline with a velocity $\hat{} = \hat{}/\hat{M}$. In Weyl ordering, one writes $\widehat{/M} = ½\left(\hat{}\hat{M}^{-1} + \hat{M}^{-1}\hat{}\right)$.
For the operator form, $(\hat{}, \hat{M}, \hat{})$ are all mutually commuting, and $(m,t)$ are c-numbers, so there is no operator-ordering ambiguity with the remainder of the expression, and you could just freely intermix everything.
The operator $\hat{}$ comes from whatever spin-representation you have. The operators $\hat{W_0}$ and $\hat{}$ come from the expression of the Pauli-Lubanski components $\left(W_0, \right)$ in terms of $(,M,)$ and $(m,α)$, the former set being converted to $(\hat{},\hat{M},\hat{})$.
So, the resulting expression is:
$$\hat{} = \frac12\left(\hat{}\frac1{\hat{M}} +\frac1{\hat{M}}\hat{}\right) + \frac{\hat{}}{\hat{M}} t - α \frac{\hat{}×\hat{}}{m\hat{M}(m+\hat{M})}.
$$
The contrast between the relativistic and non-relativistic versions has been made transparent all throughout by keeping the $α$ deformation parameter in. The relativistic versions are the deformation of the non-relativistic versions by $α = 0 → α = 1/c^2 > 0$.
As far as photons go: you can set up a position operator for them. However, it has a singularity much like the singularity that occurs in spherical coordinates. So, you need different operator representations for different parts of its phase space. Constructing the photon position operator.
The non-relativistic limits of both the photon and (a certain subclass of) the tachyon would be a similar subclass of the mass-zero representations of the Bargmann group (which is the central extension of the Galilei group). They don't have a name, so I've called them "synchrons". Like the photon, the sub-family is "helical": the helicity is invariant. All the massless particles in the Standard Model and all the ones studied for scalar and gauge fields in quantum field theory are helical. That's not generally the case for the "luxon" Wigner family (nor for the "synchron" family).
The synchrons have the property that they convey impulse instantaneously across space at an instant. That's what used to be called a "line of force", except that it acts instantaneously (so it has a "time of occurrence" operator). Some of these features are inherited by tachyons and pass over to them, instead of to luxons. In particular, instead of having a positive mass-squared invariant $m^2 = M^2 - α ||^2 > 0$, a tachyon has a positive impulse-squared invariant $||^2 = ||^2 - M^2 c^2 > 0$ and they convey impulses instantaneously (relative to at least one frame of reference) across space. They're more akin to instantaneous "lines of force" than to particles.
