What is the Interpretation of This Field Operator: $\vec{X} = \int \pi \vec{x} \phi \operatorname{d}^3 x$? In scalar quantum field theory the field momentum operator is constructed from the canonical field operators, $\phi$ and $\pi$, in the equation:
$$P_j = -\int \pi \partial_j \phi \operatorname{d}^3x.$$
As long as the field operators obey the equal time commutation relations, $\left[ \phi\left(\vec{x}\right),\pi\left(\vec{y}\right)\right] = i \delta\left(\vec{x}-\vec{y}\right)$, it is possible to show that the following operator:
$$X_j = \int \pi x_j \phi \operatorname{d}^3x ,$$
obeys the commutation relation $[X_j,P_k] = i\delta_{jk}$.
In detail:
$$\begin{align}
X_j P_k & = -\int \pi(y) y_j \phi(y) \operatorname{d}^3y \int \pi(x) \frac{\partial}{\partial x_k} \phi(x) \operatorname{d}^3x \\
& = -\int \operatorname{d}^3 x \operatorname{d}^3 y \left(\pi(y) y_j \left[i \delta\left(\vec{x}-\vec{y}\right) +\pi(x) \phi(y)\right] \frac{\partial}{\partial x_k} \phi(x)\right) \\
& = -i \int \operatorname{d}^3 x\, \pi(x) x_j \frac{\partial}{\partial x_k} \phi(x) \\
&\hphantom{=} - \int \operatorname{d}^3 x \operatorname{d}^3 y\, \left(\pi(x) \frac{\partial}{\partial x_k}\left[-i\delta\left(\vec{x}-\vec{y}\right) + \phi(x) \pi(y) \right] y_j \phi(y)\right)\\
& = -i \int \operatorname{d}^3 x\, \pi(x) x_j \frac{\partial}{\partial x_k} \phi(x) + i \int \operatorname{d}^3 x \operatorname{d}^3 y\, \pi(x) \frac{\partial}{\partial x_k} \delta\left(\vec{x}-\vec{y}\right) y_j \phi\left(\vec{y}\right) + P_k X_j\\
& = -i \int \operatorname{d}^3 x\, \pi(x) x_j \frac{\partial}{\partial x_k} \phi(x) - i \int \operatorname{d}^3 x \operatorname{d}^3 y\, \pi(x) \frac{\partial}{\partial y_k} \delta\left(\vec{x}-\vec{y}\right) y_j \phi\left(\vec{y}\right) + P_k X_j\\
& = -i \int \operatorname{d}^3 x\, \pi(x) x_j \frac{\partial}{\partial x_k} \phi(x) + i \int \operatorname{d}^3 x \operatorname{d}^3 y\, \pi(x) \delta\left(\vec{x}-\vec{y}\right) \frac{\partial}{\partial y_k} \left[ y_j \phi\left(\vec{y}\right)\right] + P_k X_j\\
& = i \delta_{jk} + P_k X_j .
\end{align}$$
What is the physical interpreation of $X_j$?
 A: *

*As Cosmas Zachos comments, the fact that your $X_j,P_i$ fulfill the same commutation relations as the usual non-relativistic position and momentum operators is not surprising. The map $\mathcal{o}\mapsto \mathcal{O} := \int \pi \mathcal{o} \phi\mathrm{d}x$ for an operator $\mathcal{o}$ in the position representation is a infinite-dimensional variant of the Jordan map, $\mathcal{o}$ being considered as an "infinite-dimensional matrix".

*$X_j$ is a "position operator" with respect to the $P_i$ at the chosen timeslice only in so far as it fulfills the correct commutation relation in that timeslice. However, it does not behave correctly under Poincaré transformations since you are missing a zeroth component to make it a proper 4-vector. For the $P_i$ the Hamiltonian $H = P^0$ naturally fulfills this role, using $p^0 = \sqrt{\vec p^2 + m^2}$ to write $P^\mu = \int p^\mu a^\dagger(\vec p) a(\vec p)\mathrm{d}^3 \vec p$. No such thing is possible for the $X_j$ - there is no time operator. Therefore, your definition of the $X_i$ is frame-dependent and not properly covariant; therefore they have no "physical interpretation" in a fully relativistic field theory. Everything you construct in the naive Hamiltonian formalism of relativistic QFT must be checked for its covariance properties; it is crucial and a priori non-trivial that the usual derivations done in textbooks yield covariant results while using non-covariant steps.

*Attempting to make the construction of position operators $[X_i,P_j] = \mathrm{i}\delta_{ij}$ Lorentz covariant is generally doomed to failure, see e.g. this answer by  Valter Moretti and links therein.
A: Interpreting $X_j$ as the position operator seems to contradict the interpretation of $\phi^2(x)$ as a probability density — note that his is a real scalar field, so complex conjugation is unnecessary, but a parallel development is possible for a complex scalar field. The answer to that is: it does, and that's because the correct probability density isn't given by $\phi^2(x)$ for a free scalar field. 
If we define the spatial Fourier transform of a function as $\tilde{f}(\mathbf{k}) = \frac{1}{(2\pi)^{3/2}} \int \operatorname{d}^3 x\, f(\mathbf{x}) \operatorname{e}^{-i \mathbf{k}\cdot\mathbf{x}}$ (i.e. symmetric convention), then the free Klein-Gordon field the momentum space annihilation operators are given by:
$$\begin{align}
a(k) &\equiv \sqrt{\frac{\omega}{2 }} \left(\tilde{\phi}(\mathbf{k}) + \frac{i}{\omega} \tilde{\pi}(\mathbf{k})\right)\\
\omega &\equiv \sqrt{\mathbf{k}^2 + m^2},
\end{align}$$
giving a Hamiltonian:
$$H = \int \operatorname{d}^3 k\, \omega \left(a^\dagger(\mathbf{k})\, a(\mathbf{k}) + \frac{1}{2} \left[a(\mathbf{k}),\, a^\dagger(\mathbf{k})\right]\right).$$
The Hamiltonian leads directly to the definition of the particle number operator as:
$$N = \int \operatorname{d}^3k\, a^\dagger(\mathbf{k})\, a(\mathbf{k}).$$
Backing out to write this in terms of the real space field operators:
$$\begin{align}
N & = \int \operatorname{d}^3k\, \frac{\omega}{2}\left(\tilde{\phi}^2(\mathbf{k}) + \frac{1}{\omega^2} \tilde{\pi}^2(\mathbf{k}) + \frac{i}{\omega}\left[\tilde{\phi}(\mathbf{k}),\, \tilde{\pi}(\mathbf{k})\right]\right).
\end{align}$$
This expression for $N$ is identical, classically, to $\int \tilde{\pi}(\mathbf{k}) \tilde{\phi}(\mathbf{k}) \operatorname{d}^3k = \int \pi(\mathbf{x}) \phi(\mathbf{x}) \operatorname{d}^3x$. The relashionship follows from the classical dropping of commutators and applying the equations of motion:
$$\begin{align}
\frac{\partial^2 \phi}{\partial t^2} &= \nabla^2 \phi - m^2 \phi \Rightarrow \\
i \frac{\partial \phi}{\partial t} & = \sqrt{m^2 - \nabla^2} \phi,\ \mathrm{and} \\
\pi & = \frac{\partial \phi}{\partial t}.
\end{align}$$
Thus, because the probability interpretation of vanilla "first quantization" quantum mechanics is classical with respect to quantum field theory, it is not a stretch to say that the probability density for a Klein-Gordon particle at $\mathbf{x}$ is given by $\pi(\mathbf{x}) \phi(\mathbf{x}) = \dot{\phi}(\mathbf{x})\phi(\mathbf{x})$, and not $\phi^2(\mathbf{x})$.
