The trick is given in equation 4.4 of the attached article:
First couple the theory to gravity, (by introducing a metric tensor in the integration measure and for each index raising) obtaining the action:
$S = \int_M d^4x \sqrt{-g} \mathcal{L}$
Then vary the action with respect to the metric tensor:
$T_{\alpha\beta} = \frac{1}{\sqrt{-g}} \frac{\delta S}{\delta g^{\alpha\beta}}$
Then replace the metric tensor by a flat one to return to the original theory. The resulting stress energy tensor is called Belinfante -
Rosenfeld stress energy tensor.
The Belinfante - Rosenfeld stress-energy tensor is not equal to the canonical stress-energy tensor corresponding the translations (derived from Noether's theorem) , but differs by the divergence of an anti-symmetric 3-tensor, thus will be conserved whenever the canonical one is conserved.
The 3-tensor is the canonical conserved current corresponding the Lorentz symmetry. Thus for spinless fields both tensors will be equal.
The Belinfante - Rosenfeld stress-energy tensor is usually considered to be better, since it is always symmetric and gauge invariant.
The reasoning why both methods give the same stress- energy tensors (“up to a total divergence) is as follows:
When one covariantizes the theory, it becomes invariant under
diffeomorphisms. Thus the variation of the action with respect to the
diffeomorphisms vanishes. Now, this variation is composed from the
variation due to the (matter) fields and the variation due to the metric. In flat space the variation due to the fields is the canonical Noether current due to translations, thus the variation with respect the metric should cancel this contribution. Therefore if we vary with respect to the metric, keeping the fields fixed, we obtain the same variation with an opposite sign.
The subtlety of not obtaining the same stress-energy tensor is explained in the following
article by: Gotay and Marsden. This article contains a few examples of
the derivation of the symmetric (Belinfante - Rosenfeld ) stress-energy
tensor. Please, see also the following article by Forger and Romer for further explanations and examples.
Symmetry of the canonical energy-momentum tensor can be related to the spin of the object(s) that contribute to it (in other words, the representation of the Lorentz group under the fields transform). Note that the canonical EM tensor is obtained by using the Noether's procedure for translational symmetry
$$
T_{\mu\nu} = \sum\limits_r \frac{\delta {\cal L}}{\delta \left( \partial^\mu \phi_r \right)} \partial_\nu \phi_r - g_{\mu\nu} {\cal L}
$$
This expression is clearly not symmetric. However, we can comment about its (non-)symmetry by looking at the conserved quantity corresponding to Lorentz transformation. By the Noether's procedure, we can show that this is
$$
M_{\mu\nu\alpha} = T_{\mu\alpha} x_\nu - T_{\mu\nu} x_\alpha - \frac{\delta {\cal L}}{ \delta \left( \partial^\mu \phi_r \right) } \left( J_{\nu\alpha} \right)^{rs} \phi_s
$$
Here, $\left( J_{\nu\alpha} \right)^{rs}$ is the representation of the Lorentz algebra under which the set of fields $\phi_r$ transforms. Conservation of this quantity implies
$$
0 = \partial^\mu M_{\mu\nu\alpha} = T_{\nu\alpha} - T_{\alpha\nu} - \partial^\mu \left( \frac{\delta {\cal L}}{ \delta \left( \partial^\mu \phi_r \right) } \left( J_{\nu\alpha} \right)^{rs} \phi_s \right)
$$
This implies
$$
T_{\nu\alpha} - T_{\alpha\nu} = \partial^\mu \left( \frac{\delta {\cal L}}{ \delta \left( \partial^\mu \phi_r \right) } \left( J_{\nu\alpha} \right)^{rs} \phi_s \right)
$$
Non-symmetry of the stress-energy tensor is an indication that the fields that are contributing to it transform non-trivially under the Lorentz group. In particular, the canonical EM tensor is symmetric only if the theory contains only scalars.
The way I like to think about the process of "symmetrizing the EM tensor" is the following. The canonical EM tensor does not contain any "spin information" and one needs the angular momentum tensor for that information. However, a symmetrized EM tensor is essentially defined to "absorb" in the spin information of the field content so that the angular momentum tensor is no longer needed (it contains more information?). The reason I think of this like this is that in terms of the symmetric EM tensor we can define another conserved quantity
$$
{\tilde M}_{\mu\nu\alpha} = {\tilde T}_{\mu\alpha} x_\nu - {\tilde T}_{\mu\nu} x_\alpha
$$
Since ${\tilde T}_{\mu\nu}$ is symmetric the tensor above is trivially conserved and does not contain any new information. However, this modified angular momentum tensor still generates all the conserved quantities as $M_{\mu\nu\alpha}$. It then seems to me that ${\tilde T}_{\mu\nu}$ already contains information about the conserved angular momentum.
This also reconciles with the fact that the symmetrized EM tensor is often the same as one obtains by varying the metric, usually defined as
$$
T_{\mu\nu} = \frac{2}{\sqrt{-g}} \frac{\delta S}{\delta g^{\mu\nu} }
$$
Since the metric (gravitation) couples to all particles in a universal fashion, the above definition of the EM tensor should involve spin information as well, and therefore, must be equal to (or atleast closely related to) ${\tilde T}_{\mu\nu}$ described above.
Best Answer
You should think at the way the Noether current is obtained. When an infinitesimal symmetry transformation is made spacetime dependent, that is the parameters $\omega^a$ that control the symmetry are taken as functions of the spacetime point $\omega^a=\omega^a(x)$, the action is not longer left invariant $$ \delta S=-\int d^D x\, J^{a\,\mu}(x)\partial_\mu \omega^a(x) $$ but rather it provides the definition of the current $J^{a\,\mu}$ that is conserved on-shell.
Now, let's look at the case of the energy momentum tensor: in this case, the translations $x^\nu\rightarrow x^\nu+\omega^\nu$ are made local $x^\nu\rightarrow x^\nu+\omega^\nu(x)$ so that $$ \delta S=-\int d^D x\, T_\nu^\mu(x)\,\partial_\mu \omega^\nu(x)\,. $$ Actually, one looks for a symmetric tensor $T^{\mu\nu}=T^{\nu\mu}$ so that one can rewrite the expression above in the following form $$ \delta S=-\frac{1}{2}\int d^D x\, T^{\nu\mu}(x)\,(\partial_\mu \omega_\nu(x)+\partial_\nu\omega_\mu)\,. $$ Now, here is the catch: if we were to transform the spacetime metric $g_{\mu\nu}$ (equal to $\eta_{\mu\nu}$ in the case at hand) as if $x^\nu\rightarrow x^\nu+\omega^\nu(x)$ was just an infinitesimal change of coordinates, that is $$g_{\mu\nu}\rightarrow g_{\mu\nu}-(\partial_\mu \omega_\nu(x)+\partial_\nu\omega_\mu)\,,$$ then the action (rendered coordinates independent by the inclusion of the metric in the usual way such as $d^D x\rightarrow d^D x \sqrt{|g|}\,,\ldots$) would be left invariant $$ \delta S=-\frac{1}{2}\int d^D x\, (\partial_\mu \omega_\nu(x)+\partial_\nu\omega_\mu)\left. \left(\sqrt{|g|}\,T^{\mu\nu}(x)+2\frac{\delta S(x)}{\delta g_{\mu\nu}}\right)\right|_{g_{\mu\nu}=\eta_{\mu\nu}}=0\,. $$ From this equation it follows that the current associated with spacetime translations can be written as $$ T^{\mu\nu}=\left. -\frac{2}{\sqrt{|g|}}\frac{\delta S}{\delta g_{\mu\nu}}\right|_{g_{\mu\nu}=\eta_{\mu\nu}} \qquad \mbox{evaluated on the bkg }g_{\mu\nu}=\eta_{\mu\nu}\,. $$ It should be apparent that this definition gives a symmetric energy-momentum tensor that matches the one appearing the Einstein's equations. From the derivation above it should also be clear that the alternate versions of $T_{\mu\nu}$ arise because the definition of $T^{\mu\nu}$ via the variation of the action when the translation is made spacetime dependent does not uniquely fix it. For example, given a valid $T_{\mu\nu}$, one can always define another one $T^{\mu\nu}\rightarrow T^{\mu\nu}_B=T^{\mu\nu}+\partial_\rho B^{\rho\mu\nu} $ with an arbitrary $B^{\rho\mu\nu}=-B^{\mu\rho\sigma}$ which also gives $$ \delta S=-\int d^D x\, T^{\mu\nu}(x)\partial_\mu\omega_\nu(x)=-\int d^D x\, T_B^{\mu\nu}(x)\partial_\mu\omega_\nu(x) $$ up to an integration by parts. Einstein's equations break this degeneracy and nicely identify ``the'' energy-momentum tensor.