I was looking for papers on the energy-time, momentum-length, and mass-tau trio of uncertainty relations when I ran across your question.
The mass-tau uncertainty relation $\{mc, \tau\}$ certainly exists.
In answer to the second part of your question, it is the source of the subpopulation of virtual particle pairs whose combined net momentum during their brief existence is zero as measured from the frame of the observer. Virtual particle pairs with non-zero combined momentum during their existence must instead be described using the energy-time uncertainty relation $\{\frac{E}{c}, t\}$, which is itself a composition of the two linearly independent uncertainty relations $\{mc, \tau\}$ and $\{\frac{E}{c}, t\}$.
I think -- I am not sure, that is why I was looking -- that little is heard about mass-tau uncertainty because it tends to end up being confused with energy-time uncertainty. If so, that is sloppy, since mass-tau is linearly independent of momentum-length, while energy-time is not.
An experiment that differentiates mass-tau and energy-time uncertainty would have to distinguish between virtual particle pairs with zero and non-zero net momentum. While both cases are implied by QED and Feynman diagrams, I have no idea if anyone has ever done experimentation that specifically differentiates zero-momentum virtual pairs (the mass-tau uncertainty pairs) from non-zero-momentum virtual pairs (the energy-time uncertainty pairs). It is unlikely such experiments would find anything new, since the Feynman QED formulation has proven extraordinarily precise and predictive in its description of the entire population of virtual particle pairs.
If you do are wondering why I keep describing energy-time uncertainty as a composition of the linearly independent mass-tau and momentum-length relations, please see below.
Uncertainty relationships all have units of action, so to see how they relate to each other it is convenient to express the relevant equations in momentum-like units and length-like units.
Start with the energy-momentum relation $E^2 = p^2c^2 + m^2c^4$. Divide by $c^{2n}$ to get:
$\begin{array}{cccccll}
(Ec^{0})^2 & = & ({p_x}c^{1})^2 & + & (mc^{2})^2 & & \leftarrow\mbox{energy units}
\\
(Ec^{-1})^2 & = & ({p_x}c^{0})^2 & + & (mc^{1})^2 & & \leftarrow\mbox{momentum units}
\\
(Ec^{-2})^2 & = & ({p_x}c^{-1})^2 & + & (mc^{0})^2 & & \leftarrow\mbox{mass units}
\end{array}$
To emphasize the momentum units of the second form above, define $p_E=\frac{E}{c}$ and $p_m=mc$:
$p_E^2 = p_x^2 + p_m^2$
All three terms have directional interpretations, so the above quadratic equation can be re-interpreted as a linear vector sum:
$\overrightarrow{p_E} = \overrightarrow{p_x} + \overrightarrow{p_m}$
The right side terms two orthogonal axes of a plane, so the vector sum can also be expressed in terms of two linearly independent unit vectors $\hat{p_x}$ and $\hat{p_m}$:
$\overrightarrow{p_E} = p_x\hat{p_x} + p_m\hat{p_m}$
The resulting Pythagorean vector sum and planar space look like this:
Completing the set with the unit vectors for the plane $p_{y}p_{z}{\bot}p_{x}$ gives the Euclidean 4D set of momentum unit vectors $\{\hat{p_{m}},\hat{p_{x}},\hat{p_{y}},\hat{p_{z}},\}$, which serves the same linearization purpose as the Dirac $\gamma$ matrices. The Dirac matrices are more complicated because Dirac chose to factor a mixed-signature $[-+++]$ Minkowski space, which requires matrices. I do not know if Dirac was aware of (or cared about) the similarity of his approach to vector addition in an all-plus Euclidean space. Mass remained very firmly a scalar in Dirac's approach, and that decision guided many of his subsequent choices. Mass as a momentum-like vector replaces negative energy with negative mass, which is a considerably simpler concept in that it avoids the asymmetry of positive-mass electrons moving with negative energy.
The other half of the uncertainty trio is the length-unit relations. From SR, the corresponding relation for time, length, and proper time (again arranged to give an all-positive Euclidean signature) is:
$t^2 = x^2 + {\tau}^2$
As with the momentum equation, relabeling with $x_t=t$ and $x_{\tau}=\tau$ emphasizes their length-like units:
$x_t^2 = x^2 + x_{\tau}^2$
The observer-to-observed relationship, with the vector origins in the observer frame, provides a straightforward and experimentally meaningful interpretation of these lengths as vectors, which again allows the quadratic equation to be simplified into a linear vector sum:
$\overrightarrow{x_t} = \overrightarrow{x} + \overrightarrow{x_{\tau}}$
The two right side terms again form two orthogonal axes of a plane, with the linearly independent unit vectors this time being $\hat{x}$ and $\hat{x_m}$:
$\overrightarrow{x_t} = x\hat{x} + x_{\tau}\hat{ x_{\tau }}$
Diagrammatically, the result is this:
Adding the unit vectors for $yz{\bot}z$ gives a set of four length unit vectors $\{\hat{x_{\tau}}, \hat{x}, \hat{y}, \hat{z}\}$, again with strong parallels to the Dirac $\gamma$ matrices, and an exact parallel to the earlier set of momentum unit vectors. This set uses $\tau$ instead of $t$, however, so it does not describe regular spacetime. For example, this tau spacetime (I'm not sure if it has a common name) has a velocity angle range of $0^{\circ}{\le}\theta_{t}{\le}90^{\circ}$, versus the velocity angle range of $0^{\circ}{\le}{\alpha}{\le}45^{\circ}$ for regular or t-spacetime.
Placing the two vector sum equations in parallel produces the three uncertainty relations:
$\begin{array}{ccccc}
\{p_E,x_t\} & & \{p_x, x\} & & \{p_m,x_{\tau}\}
\\
\overbrace{\overrightarrow{p_E}} & = & \overbrace{p_x\hat{p_x}} & + & \overbrace{p_m\hat{p_m}}
\\
\overrightarrow{x_t} & = & x\hat{x} & + & x_{\tau}\hat{ x_{\tau }}
\end{array}$
Particles reside simultaneously as chiral waveforms in both of the 4D spaces (momentum and length), with symmetric Fourier transforms linking the corresponding pairs of momentum and length axes. It is the Fourier relationships between these four cross-space pairs $\{x_{\tau}, p_{m}\}$, $\{x, p_{x}\}$, $\{y, p_{y}\}$, and $\{z, p_{z}\}$ that guarantees reciprocal uncertainty within each pair.
In terms of these pairs, the prevalent use of energy-time uncertainty over mass-tau is most likely just a bit of observer bias, since for the observer frame (only) energy is mass ($E=mc^2$) and time is tau ($t=\tau$), making the two uncertainty pairs fully interchangeable. Since the clock-based $\tau$ time is always real, measurable, and causal, just like spatial length $x$, what we call "time" or $t$ ends up being measured in terms of $\tau$ anyway.
A side note: The main significance of a span of $t$ when interpreted as a Euclidean vector sum is that its total length remains invariant for all observed frames. E.g., if the observer frame sees 60 seconds of $t$ time as measured via the $\tau$ of a clock at rest relative to the observer, then the vector sum of $x$ (distance traversed) and $\tau$ (change in an associated moving clock) for any causally observed frame must add up to 60 seconds, regardless of the velocity of that frame. At $v=0$ the composition of $t$ will be 100% $\tau$ seconds (no motion), while at $v=1$ (in $c$ units) $t$ becomes 100% spatial length, with change halted entirely (photons).
This was (is) one of my biggest bugbears whilst learning QFT. The reasoning behind using a Fock space is actually really simple and intuitive for a scalar field (provided you are comfortable with standard QM) but it's always masked by the horrible concept of 'canonical quantisation'.
Take the Fourier transform of the Klein-Gordon equation:
$$
\partial_t^2\hat{\phi} = -\omega_{k}^2\hat{\phi}
$$
this is the classical equation of motion for a harmonic oscillator. There's one of these for every possible momentum (this corresponds to $\phi$ being a field over space transforming into $\hat{\phi}$ which is a field over momentum).
It's harder to transform the Lagrangian (since it includes terms like $\left(\frac{\partial\phi}{\partial t}\right)^2$), but one can assume it similarly describes $\hat{\phi}$ as an infinite spectrum of harmonic oscillators. Given this, it's reasonable to assume that the quantum Lagrangian/Hamiltonian for $\hat{\phi}$ similarly corresponds to an infinite spectrum of quantum harmonic oscillators. This we do know the form of:
$$
H = \int \frac{d^{3}p}{2E_{p}}a_{p}^{\dagger}a_{p},
$$
where I've dropped the zero point energy because you can (it corresponds to reordering the fields, which is an ambiguity that exists in going to the non-commuting quantum case from the commuting classical case) and it avoids the usual infinite-energy problems.
Now we just claim that $\phi$ has the same quantum Lagrangian as in the classical case and that the above Hamiltonian is the Fourier transform (of the Legendre transform) of the Lagrangian. If you work through you find that you get out the canonical form of the field. David Tong does this on page 24 of his notes, though he does it by essentially proposing the canonical form as an ansatz.
Then you just use your infinite set of annihilation and creation operators that arose naturally to generate the infinite set of QHO number states (one for each momentum). This is identical to the Fock space generated by the momentum operator, so you just treat it as a Fock space.
Best Answer
I) Well, the Legendre transformation can be e.g. seen as the leading classical tree-level formula of a formal semiclassical Fourier transformation.
This fact is e.g. used in QFT when relating the quantum action $S[\varphi]$, the partition function $Z[J]$, generating functional $W_c[J]$ for connected diagrams, and the effective action $\Gamma[\Phi]$.
II) To see the correspondence in detail, let $x$ and $p$ be the two dual/conjugate variables (in both senses!). Let
$$\tag{1} f(x;\hbar)~\equiv~\sum_{n=0}^{\infty}\hbar^n f_n(x)\quad\text{and}\quad g(p;\hbar)~\equiv~\sum_{n=0}^{\infty}\hbar^n g_n(p)$$
be two formal power series in $\hbar$ with function coefficients. Consider their semiclassical exponentials
$$\tag{2} F(x;\hbar)~:=~e^{if(x;\hbar)/\hbar}\quad\text{and}\quad G(p;\hbar)~:=~e^{ig(p;\hbar)/\hbar}.$$
Now assume that
$$\tag{3} G(p;\hbar)~=~\int \! dx~ e^{-ipx/\hbar} F(x;\hbar)$$
is the Fourier transform of $F(x;\hbar)$. We can use the WKB/stationary phase approximation to deduce that the classical parts $f_0(x)$ and $-g_0(p)$ are then Legendre duals of each other, i.e.
$$\tag{4} g_0(p) ~=~ -px+f_0(x)\quad\text{where}\quad p~=~f_0^{\prime}(x),$$
for a sufficiently nice function $f_0(x)$.