It is common to state that a proton is bounded state of three quarks, and that the QCD energy binding (associated to a "cloud" of gluons joining together the three quarks) is responsible for the 99% of the inertial mass of a proton. For example, Durr et al. (Ab-initio Determination of Light Hadron Masses), confidently state that: According to the mass-energy equivalence, $E = m c^2$, we experience this energy as mass. This is an informal statement, and it is not clear to me why any amount of energy in the form of massless quanta of a quantum field contribute to "inertial mass".
In addition, I have checked other articles computing the proton mass, for example, Yang et al. (2018), Proton Mass Decomposition from the QCD Energy Momentum Tensor, presents a detailed and beautiful computation for the component $\langle T_{44}\rangle$ of stress-energy tensor of a proton. Taking the "QCD-mass of a proton" as $M_p = \langle T_{44}\rangle/c^2$, Yang et al. obtain that $M_p = f(m_\pi^2)$ where $m_\pi$ is the mass of a pion, taking the experimental value of $m_\pi c^2 \approx 139$ MeV, they obtain $M_p c^2 \approx 938$ MeV a beautiful numerical agreement.
But: (1) How can we demonstrate that the bounded quark-gluon state which constitutes a proton presents a resistance to be accelerated, or ''inertial mass'', equivalent to $E_p/c^2$? (I understand the informal affirmations, but I am looking for a clear/formal mathematical proof/argument/computation).
My guess is this one: A box contains a quantity of energy, an observer at rest with respect to the center of the box will measure an stress-energyimpulse tensor given by:
$$T^\mu_\nu = \begin{pmatrix} \epsilon & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 \end{pmatrix}$$
(assume for simplicity the space-time is Minkowskian with metric -+++). If the box is accelerated until it moves with uniform velocity (directed along the X-axis) the final stress-energy tensor will be:
$$T^\mu_\nu = \begin{pmatrix} \gamma\epsilon & \gamma\beta\epsilon & 0 & 0\\ \gamma\beta\epsilon & \gamma\beta^2\epsilon & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 \end{pmatrix}$$
with $\gamma = 1/\sqrt{1-v^2/c^2}$ and $\beta = v/c$. During the acceleration process we need a force density $f$, and the equations of motion will be:
$$\frac{1}{c}\frac{\partial T^0_0}{\partial t} + \frac{\partial T^1_0}{\partial x} = \gamma \frac{fv}{c}$$ $$\frac{1}{c}\frac{\partial T^0_1}{\partial t} + \frac{\partial T^1_1}{\partial x} = \gamma f$$
Assuming that the force is continuous and only acts for a finite time, it is easy to see that the force will always be proportional to:
$$f\propto \frac{\text{d}}{\text{d}\tau}\left( \frac{(\epsilon/c^2)v}{\sqrt{1-v^2/c^2}} \right)$$
Thus the term $\epsilon/c^2$, behaves indistinguishably from a mass density at rest $\rho$.
