The question can be addressed and resolved in a way that spans the divide between classical and quantum physics, between relativity and non-relativistic theory, and applies across the board uniformly to them all; and this I will do here, now.
In the quantum case, the analysis will lead to a set of position-momentum Heisenberg triples - one for each non-interacting body in the system - and the usual method for resolving probability distributions under constraints on the momentum may then be applied. To get sensible results, you can't have a sharp constraint on the momentum, but rather something smeared, like with a Gaussian. Then, in that case, you are using coherent states with the position and momentum operators.
Here's the general answer.
A physical system has a state space that is acted on by the kinematic symmetry group associated with the underlying paradigm. In the case of non-relativistic physics, this is the central extension of the Galilei group, known as the Bargmann group. For Relativity, it is the Poincaré group. However, to make it compatible with the correspondence limit, and to have the Bargmann group as its non-relativistic limit, it should actually be raised to a (trivially) central extended form of the Poincaré group. We will see where that comes into play in a moment.
The kinematic groups define symmetry transforms for spatial translations, time translations, rotations and "boosts". The latter are what a transform from a stationary to a moving frame is called. By the Noether Theorem, for an isolated non-interacting system, associated with each type of symmetry transform is a conserved quantity - its Noether charge. The charges, themselves, are acted on by the symmetry transforms as the "co-adjoint" representation of the symmetry group. The Noether charges are conventionally identified as:
- linear momentum $𝐏$, for spatial translations,
- angular momentum $𝐉$, for rotations,
- (negative) energy $-H$, for time translations,
- adjusted mass moment $𝐊$, for boosts,
- mass $μ$, for translations with respect to a fictitious fifth coordinate.
The last transform is needed, in the non-relativistic case, to provide a matching coordinate representation for the Bargmann group. This is exactly the extension that takes place when going from the Galilei group to the Bargmann group, and in the corresponding central extension, $μ$ is the central charge. I call it $μ$, instead of $m$, for reasons that will become shortly apparent.
The Galilei group is what results from setting $μ = 0$. Its representations are "zero mass" systems, which can either be instantaneous: like an instantaneous action-at-a-distance impulse, at a specific time, acting over an interval in space in a specific direction, or homogeneous: uniform over space, with no center of mass. A special case of this is the homogeneous system that is isotropic and boost-invariant: the "vacuum".
Under the co-adjoint action, all the Noether charges are invariant with respect to translation on the extra coordinate, so that it acts as a gauge degree of freedom. By the Noether theorem, this makes the associated Noether charge, the mass $μ$ a conserved quantity.
For a non-relativistic system with these Noether charges, and with non-zero mass, it is possible to use the co-adjoint action to transform the momentum by a boost to $𝐏 → 𝟬$. This lands in the "rest frame" of the system. Under this boost, the angular momentum transforms as $𝐉 → 𝐒$ to the internal angular momentum of the system, and the energy $H → U$ to the internal energy of the system. The mass and moment remain invariant: $μ → μ$, $𝐊 → 𝐊$ as they are both boost-invariant. In Relativity, matters are a little more complicated: the moment in that case is not quite boost invariant.
In the rest frame, $𝐊$ is a bona fide mass moment. Under a spatial translation, it may be transformed $𝐊 → 𝟬$. All other Noether charges remain unchanged. This lands you at the center of mass of the system.
Applying these processes in reverse with the following transforms leads to the following results for the step-by-step transform of $(𝐉,𝐊,𝐏,H,μ)$:
- Spatial translation by $𝐫$:
$(𝐒,𝟬,𝟬,U,μ) → (𝐒,μ𝐫,𝟬,U,μ)$,
- Boost by $-𝐯$:
$(𝐒,μ𝐫,𝟬,U,μ) → \left(μ𝐫×𝐯 + 𝐒, μ𝐫, μ𝐯, ½μ|𝐯|^2 + U, μ\right)$,
- Time translation by $t$:
$\left(μ𝐫×𝐯 + 𝐒, μ𝐫, μ𝐯, ½μ|𝐯|^2 + U, μ\right) → \left(μ𝐫×𝐯 + 𝐒, μ(𝐫 - 𝐯t), μ𝐯, ½μ|𝐯|^2 + U, μ\right)$,
- Final result:
$\left(μ𝐫×𝐯 + 𝐒, μ(𝐫 - 𝐯t), μ𝐯, ½μ|𝐯|^2 + U, μ\right) = (𝐉,𝐊,𝐏,H,μ),$
which leads to the following decomposition:
$$𝐉 = m𝐫×𝐯 + 𝐒, \quad 𝐊 = m(𝐫 - 𝐯t), \quad 𝐏 = m𝐯, \quad H = \frac{m|𝐯|^2}{2} + U, \quad μ = m.$$
This describes a system of mass $m = μ$, a center of mass position $𝐫$ and associated velocity $𝐯$. You will also notice that the mass moment is projected back to time $t = 0$. It has an explicit dependency on $t$ and - by this account - a degree of ambiguity. The locus of all positions, as you vary time $t$ is the trajectory of the system's center of mass.
The appearance of internal values for angular momentum and energy are justified, after the fact, by the requirement that a composite system, having such elementary systems as components, should also have a set of attributes of its own that are obtained by requiring all the Noether charges to be additive. Here, we can see that for two bodies, with the following rule of composition:
$$
μ = μ_0 + μ_1, \quad 𝐫 = \frac{μ_0𝐫_0 + μ_1𝐫_1}{μ}, \quad 𝐯 = \frac{μ_0𝐯_0 + μ_1𝐯_1}{μ},\\
𝐒 = 𝐒_0 + 𝐒_1 + \bar{μ}\bar{𝐫}×\bar{𝐯}, \quad U = U_0 + U_1 + \frac{\bar{μ}|{\bar{𝐯}}|^2}{2},\\
\bar{μ} = \frac{μ_0μ_1}{μ}, \quad \bar{𝐫} = 𝐫_0 - 𝐫_1, \quad \bar{𝐯} = 𝐯_0 - 𝐯_1.
$$
The internal dynamics of the two component bodies making up the larger body is captured by a fictitious body with effective mass $\bar{μ}$, position $\bar{𝐫}$ and velocity $\bar{𝐯}$. The non-linearity of the expressions for $𝐉$ and $H$ is compensated for by extra additions made to $𝐒$ and $U$, respectively.
There is nothing in principle that requires either of these internal values to vanish for non-composite systems. Therefore, it is meaningful to speak of "spin" and not just "internal energy" - even for classical physics, and even for non-relativistic theory. Also, the decomposition helps bolster the description of $𝐫$ and $𝐯$, respectively, as center of mass position and velocity.
For any system, be it elementary or composite, the relations may be inverted, and position may be defined in terms of the Noether charges as
$$𝐫 = \frac{𝐊 + 𝐏t}{μ}.$$
All of what was described here applies to both classical and quantum theory. However, when passing over to quantum theory, the Noether charges become "q-numbers" and may have non-zero commutators, which can lead to operator ordering ambiguity. Fortunately, for the decomposition just presented, there is no non-commutativity between the factors at all. In Relativity, however, there will be a minor exception in the case of $𝐊$.
How are the commutators determined? In the classical theory, Poisson brackets may be associated with the Noether charges by replicating the Lie brackets on the Noether charges, themselves. For instance:
$$\left\{K_i, P_j\right\} = δ_{ij} μ, \quad \left\{J_x, J_y\right\} = J_z,$$
with cyclic permutations of $(x,y,z)$ in the latter case. In quantized form, this becomes, respectively:
$$\left[K_i, P_j\right] = iħ δ_{ij} μ, \quad \left[J_x, J_y\right] = iħJ_z.$$
The parameter $t$ is treated as a classical coordinate and does not directly enter any of the relations.
From this follows the usual Heisenberg relations, e.g.
$$\left[𝐫, P_x\right] = \left[\frac{𝐊 + 𝐏t}{μ}, P_x\right] = \frac{\left[𝐊, P_x\right]}{μ} = \frac{iħμ(1,0,0)}{μ} = iħ(1,0,0),$$
where $𝐫 = (x,y,z)$, as well as those for internal angular momentum, e.g.:
$$\left[S_x, S_y\right] = iħS_z.$$
The two sets of operators for $𝐫$ and $𝐒$ arise from quantities that have zero Poisson brackets with each other, so that they are mutually commuting.
Your question restricted its focus to free systems - many-body systems whose component bodies are mutually non-interacting. In such a system, each body possesses bona fide Noether charges of their own as conserved quantities, so that all of what was described applies to each body individually, and not just to the whole system.
This is true in the classical case and in the quantum case of a free field described as a many-body system with a Fock space.
If the bodies are interacting, then you lose the ability to make a clear description of the individual bodies, since their "Noether charges" are no longer conserved quantities and are not actually Noether charges at all anymore. In the non-relativistic case, you can still use them as time-varying quantities and describe the system dynamics in terms of them. That would be a lot more difficult to do in the Relativistic case, since you have a road-block known as the no-interaction theorem (or Leutwyler Theorem) in the classical case, and the Haag Theorem in the quantum case.
Now, for the relativistic case, if we use the central charge, then the brackets for the moment and momentum become the following:
$$\left\{K_i, P_j\right\} = δ_{ij} \left(μ + \frac{H}{c^2}\right) \quad (i,j = x, y, z),$$
while for the moment, the brackets become
$$\left\{K_x, K_y\right\} = -\frac{J_z}{c^2},$$
with similar brackets obtained by cyclically permuting $(x,y,z)$. These are actually the only differences between the extended Poincaré group and Bargmann group. So, we can treat it as a deformation of the Bargmann group with an extra "deformation" parameter that goes from $0$ to $1/c^2$.
What makes this an extension of Poincaré is that it has 11 Noether charges, instead of 10, and the respective Lie algebra has 11 dimensions, instead of 10. The reduction to 10 is obtained by combining the energy $H$ - which is just kinetic energy plus internal energy - and the mass $μ$ into a single parameter $E$ for "total energy". Equivalently we can also express this as the "relativistic mass" $M = μ + H/c^2$. The relation between the two, of course, is $E = Mc^2$.
Both $H$ and $μ$ can then be replaced by $E$ or by $M$. Conventionally, the choice is with $E$, though $M$ actually gives rise to more intuitive descriptions and - more importantly - has a correspondence limit to non-relativistic theory. There is no non-relativistic limit for $E = Mc^2$, since it blows up to infinity as $1/c² → 0$.
In non-relativistic theory, there are two invariants that can be formed from the central charge, energy and momentum:
$$μ, \quad |𝐏|^2 - 2μH.$$
Under deformation to the extended Poincaré group, these respectively become:
$$μ = M - \frac{H}{c^2}, \quad |𝐏|^2 - 2μH - \frac{H^2}{c^2} = |𝐏|^2 - 2MH + \frac{H^2}{c^2}.$$
The only combination (up to functional dependence) that can be formed from them that removes the explicit dependency on $(μ, H)$ is:
$$μ^2 - \frac{1}{c^2}\left(|𝐏|^2 - 2μH - \frac{H^2}{c^2}\right) = M^2 - \frac{|𝐏|^2}{c^2} = m^2,$$
which - for systems satisfying the condition $M^2 c^2 ≥ |𝐏|$ - may be identified as the square of the rest mass $m$. The term "rest" is justified - at least in the case $|m| > 0$ by the ability to transform to a rest frame, where $M → m$.
If we keep the central charge intact and stay within the extended Poincaré group, then it is possible to carry out a set of transforms similar to the non-relativistic case, but slightly more complicated, that lands at:
$$(𝐉,𝐊,𝐏,H,μ) → (𝐒,𝟬,𝟬,U,μ),$$
just like before. Finding the right combination of transforms to achieve this, however, turns into an exercise much like solving the Rubik's Cube.
Going in the opposite direction, upon transform back to a general state, what we get is the following decomposition that is a deformation of the non-relativistic case:
$$
𝐉 = M𝐫×𝐯 + 𝐒, \quad 𝐊 = M(𝐫 - 𝐯t) + \frac{1}{c^2}\frac{M𝐯×𝐒}{m+M},\\𝐏 = M𝐯, \quad H = \frac{M|𝐯|^2}{m + M} + U, \quad μ = m - \frac{U}{c^2},
$$
where
$$M = \frac{m}{\sqrt{1 - |𝐯|^2/c^2}}.$$
The reduction from 11 charges to 10 that arises from replacing $H$ and $μ$ with $E$ or $M$ corresponds to the elimination of $U$, resulting in the following relations in their place
$$
𝐉 = M𝐫×𝐯 + 𝐒, \quad 𝐊 = M(𝐫 - 𝐯t) + \frac{1}{c^2}\frac{M𝐯×𝐒}{m+M},\\
𝐏 = M𝐯, \quad H = \frac{M|𝐯|^2}{m + M} = E - mc^2, \quad μ = m\\
M = \frac{m}{\sqrt{1 - |𝐯|^2/c^2}} = \frac{E}{c^2}.
$$
This makes both $μ$ and $H$ redundant and unnecessary, so that they drop out.
As a result, the internal energy $U$ drops out of the picture. I don't know if retaining it would provide a means to overcome the no-interaction theorem, but without its presence, there is nothing to absorb the non-additivity for the energy, and this creates part of the problem that leads to the no-interaction theorem. But the theorem may be robust enough to continue to block non-trivial interacting many-body systems, even if internal energy is included for each body.
You will notice that the moment has acquired a dependence on the spin or internal angular momentum $𝐒$. This arises from the $𝐊-𝐊$ Poisson brackets, now being a non-zero multiple of $𝐉$. This is directly connected to a feature unique to Relativity: that two successive boosts in non-collinear directions results also in a small rotation.
Another feature that arises is that the components of the Noether charge don't all commute (or have non-zero Poisson brackets) anymore. Though they continue to do so for the Noether charges $𝐉$, $𝐏$, $H$, $μ$ (and for $M$ and $E$), they don't for $𝐊$. The Poisson bracket for the mass and moment deforms from the non-relativistic case to the relativistic case as:
$$0 = \{μ, 𝐊\} ⇒ \{M, 𝐊\} = \frac{1}{c^2}\{E, 𝐊\} = \frac{M}{c^2}\{E, 𝐫\} = \frac{M}{c^2}𝐯 = \frac{𝐏}{c^2}.$$
Thus,
$$\{M, 𝐫\} = \frac{\{M, 𝐊\}}{M} = \frac{𝐯}{c^2}.$$
Correspondingly, for the commutators, we have:
$$[M, 𝐫] = iħ\frac{𝐯}{c^2}.$$
The operator ordering ambiguity may be absorbed into a shift on $t$ in the expression $M(𝐫 - 𝐯t)$, which is already a free-running parameter (though the shift is by an imaginary value). This results in a time-uncertainty in the position of the body on its center of mass trajectory.
The prevailing convention is to adopt Weyl operator ordering, i.e.
$$\widehat{M𝐫} = \frac{1}{2}(\hat{M}\hat{𝐫} + \hat{𝐫}\hat{M}).$$
Similarly, when inverting the relation between $𝐊$ and $𝐫$ to obtain an expression for $𝐫$, Weyl operator ordering may be used in the quantum case. The resulting expression is known as the Newton-Wigner position operator.
$$𝐫 = \frac{𝐊 + 𝐏t}{M} + \frac{1}{Mc^2}\frac{𝐒×𝐏}{m+M}
$$
As before, the Poisson brackets on the Noether charges reduce to the usual Poisson and Heisenberg brackets on $𝐫$ and $𝐒$, as a result. Upon quantization, with Weyl ordering, the part of $𝐊$ with the non-commuting factors is set as:
$$\widehat{\frac{𝐊}{M}} = \frac{1}{2}\left(\frac{1}{\hat{M}}\hat{𝐊} + \hat{𝐊}\frac{1}{\hat{M}}\right).$$
There is a dependence on the internal angular momentum or (for elementary systems: spin), $𝐒$. But the extra term commutes with $𝐏$, so it has momentum eigenstates as its own eigenstates.
For a non-interacting many-body system, this may be applied to each body's set of Noether charges, as well as to the system as a whole. The way seems to be blocked, however, to making this work - in both the classical and quantum cases - for interacting bodies, where we allow the "Noether charges" to be time-varying. The Haag Theorem blocks non-trivial many-body dynamics on Fock space in the quantum case, while the Leutwyler Theorem does the same for the classical case.
