In a theory with spinor fields, what general condition on the action ensures that the stress-energy tensor can be made symmetric? In general relativity, the stress-energy tensor is normally defined by
$$
T^{ab}\equiv \frac{2}{\sqrt{|g|}}\,\frac{\delta S_m}{\delta g_{ab}},
\tag{1}
$$
where $S_m$ is the "matter" action (the part of the action without the Einstein-Hilbert term), $g$ is the metric, and $|g|$ is the magnitude of its determinant. The definition (1) is natural because this is the quantity that appears automatically when the equation of motion for the metric field is written as $\delta S/\delta g_{ab}=0$, where $S$ is the total action $S$ ($=S_m+$ the Einstein-Hilbert term). The tensor (1) is manifestly symmetric, because $g_{ab}$ is symmetric.
But when the "matter" consists of a Dirac spinor field $\psi$, we can't write the equation of motion that way. Instead, we need to express $g_{ab}$ in terms of a frame field (vielbein) $e^A_a$ and write the equation of motion as $\delta S/\delta e^A_a=0$ instead. Each lowercase index $a,b,c$ is a spacetime index, and an uppercase index $A$ labels the internal components of the frame field. After contracting with the frame field so that all remaining free indices are spacetime indices, we again get the usual form of the equation of motion for the metric field, but with the stress-energy tensor given by
$$
\tilde T^{ab}\equiv \frac{1}{|e|}\,
 \frac{\delta S_m}{\delta e^A_a}\,g^{bc}e^A_c.
\tag{2}
$$
In contrast to (1), the tensor (2) is not manifestly symmetric. In fact, it generally really isn't symmetric (reference 1), but it can often be made symmetric by adding terms that vanish on-shell — that is, when $\psi$ and other "matter" fields all satisfy their own equation of motion. In particular, the free Dirac spinor field has this property (reference 2). But checking this case-by-case is tedious and unenlightening.
Question: What simple, physically-motivated condition on $S_m$ ensures that (2) can be made symmetric by adding terms that vanish on-shell? One obvious sufficient condition is when $S_m$ depends on $e^A_a$ only via the metric $g_{ab}$. With that condition, we can easily show that (2) is equivalent to (1), but that condition excludes spinor fields. I'm asking for a sufficient condition that works with spinor fields, too, which means that the $e^A_a$-dependence of $S_m$ cannot be expressed using $g_{ab}$ alone.
By the way, even though the motive for considering spinor fields is to study quantum effects, here I'm only considering a classical spinor field to avoid complications related to operator ordering and renormalization.


*

*Gotay and Mardsen (2001), Stress-Energy-Momentum Tensors and the Belinfante-Rosenfeld Formula (https://authors.library.caltech.edu/19366/1/GoMa1992.pdf)


*The end of section 3.8 in Birrell and Davies (1982), Quantum Fields in Curved Space (Cambridge University Press)
 A: Thanks to a couple of tips provided in comments, I can answer my own question.
As in the question, I'll use lowercase for spacetime indices and uppercase for veilbein indices. Let $\psi$ be a spinor field, and let $\omega$ be the spin connection. The spin connection can be expressed in terms of the vielbein $e_a^A$. The combination $(\partial+\omega)\psi$ transforms covariantly under internal/local Lorentz transformations, which act on the vielbein indices. The action is typically invariant under such transformation by construction, such as the action for a free spinor field in curved spacetime.
For any variation whatsoever, the effect on the matter action $S_m$ is
$$
\newcommand{\opsi}{\overline\psi}
\delta S_m=\int \frac{\delta S_m}{\delta e^A_a}\delta e^A_a
 +\int \frac{\delta S_m}{\delta \psi}\delta \psi
 +\int \delta \opsi\frac{\delta S_m}{\delta \opsi}
\tag{3}
$$
where the integrals are over spacetime. Now consider a special variation, namely one representing an infinitesimal local/internal Lorentz transformation. In this case, we have $\delta e^A_a=h^A{}_B e^B_a$, where $h_{AB}=-h_{BA}$. Note that $h$ may depend on the spacetime coordinates. We normally consider models for which $S_m$ is invariant under such transformations. This includes the example mehtioned in the question (free Dirac spinor field in curved spacetime). For such an action, we have $\delta S_m=0$ for this transformation, so (3) becomes
$$
\newcommand{\opsi}{\overline\psi}
0=\int \frac{\delta S_m}{\delta e^A_a}h^A{}_B e^B_a
 +\int \frac{\delta S_m}{\delta \psi}\delta \psi
 +\int \delta \opsi\frac{\delta S_m}{\delta \opsi}.
\tag{4}
$$
Using equation (2) in the question, we can rewrite this as
$$
\newcommand{\opsi}{\overline\psi}
0=\int |e|\,\tilde T^{ab} h_{AB} e^A_b e^B_a
 +\int \frac{\delta S_m}{\delta \psi}\delta \psi
 +\int \delta \opsi\frac{\delta S_m}{\delta \opsi}.
\tag{5}
$$
The equations of motion for the spinor field are
$$
 \frac{\delta S_m}{\delta \psi} = 0
\hspace{2cm}
\frac{\delta S_m}{\delta \opsi} = 0.
\tag{6}
$$
Since $h_{AB}(x)$ is arbitrary except for the condition $h_{AB}(x)=-h_{BA}(x)$, equations (5) and (6) imply that the antisymmetric part of $\tilde T^{ab}$ must be zero when the spinor field satisfies its equation of motion.
Altogether, if the action $S_m$ is invariant under local/internal Lorentz transformations (acting on the vielbein indices), then the stress-energy tensor defined by (2) is symmetric whenever the spinor field satisfies its equation of motion. That answers my question.
A: I'll try to answer clearly, but the detailed calculations are pretty involved and laborious.
Firstly, the Dirac field lagrangian density is assumed to have the following form (this expression gives the usual Dirac equation in curved space, when we do an arbitrary variation of $\Psi$):
\begin{equation}\tag{1}
\mathscr{L} = i \frac{1}{2} \Bigl( \bar{\Psi} \, \Gamma^{\mu} D_{\mu} \Psi - (D_{\mu} \bar{\Psi}) \, \Gamma^{\mu} \Psi \Bigr),
\end{equation}
where $\Gamma^{\mu} = \gamma^a \, e_a^{\mu}$ are the point dependant Dirac matrices (often called "Fock-Ivanenko coefficients", and also "curvy" Dirac matrices, while $\gamma^a$ are the "flat" Dirac matrices).  $D_{\mu}$ is the spinor covariant derivative (under local Lorentz transformations):
\begin{equation}\tag{2}
D_{\mu} \Psi = \partial_{\mu} \Psi + \omega_{\mu} \Psi = \partial_{\mu} \Psi - i \frac{1}{4} \, \omega_{\mu}^{\; ab} \, \sigma_{ab} \, \Psi,
\end{equation}
where $\omega_{\mu}^{\; ab}$ is the spin connexion.  Notice that (1) is a real number, and it is zero for the "on-shell" spinor fields.  The energy-momentum components $T_{\mu \nu}$ are defined by this expression:
\begin{equation}\tag{3}
T_{\mu \nu} \, \delta g^{\mu \nu} = 2 \, \delta \mathscr{L} - g_{\mu \nu} \, \mathscr{L} \, \delta g^{\mu \nu}.
\end{equation}
The last term is 0 for the on-shell field.  So we need to compute the explicit variation $\delta\mathscr{L}$ of (1) under an arbitrary variation of the metric: $\delta g^{\mu \nu}$.  Since $g^{\mu \nu} = \eta^{ab} \, e_a^{\mu} \, e_b^{\nu}$ and $g_{\mu \nu} = \eta_{ab} \, e_{\mu}^a \, e_{\nu}^b$, we have the following variation terms (these exclude variations that don't come from the metric itself.  Infinitesimal local Lorentz transformations, which leave the local metric invariant, are explicitly excluded):
\begin{align}
\delta e_{\lambda}^a &= \frac{1}{2} \, e_{\mu}^a \, g^{\mu \nu} \, \delta g_{\lambda \nu} \equiv -\, \frac{1}{2} \, e_{\mu}^a \, g_{\lambda \nu} \, \delta g^{\mu \nu}, \tag{4} \\[2ex]
\delta e_a^{\mu} &= \frac{1}{2} \, e_a^{\lambda} \, g_{\lambda \nu} \, \delta g^{\mu \nu} \equiv -\, \frac{1}{2} \, e_a^{\lambda} \, g^{\mu \nu} \, \delta g_{\lambda \nu}, \tag{5} \\[2ex]
\delta \Gamma^{\mu} &= \gamma^a \, \delta e_a^{\mu} = \frac{1}{2} \, \Gamma_{\nu} \, \delta g^{\mu \nu}. \tag{6}
\end{align}
Using (2) and (6), the variation of the lagrangian density (1) is this:
\begin{equation}\tag{7}
\delta \mathscr{L} = i \frac{1}{4} \Bigl( \bar{\Psi} \, \Gamma_{\nu} \, D_{\mu} \, \Psi - (D_{\mu} \bar{\Psi}) \, \Gamma_{\nu} \, \Psi \Bigr) \, \delta g^{\mu \nu} + i \frac{1}{2} \bar{\Psi} \{ \, \Gamma^{\mu}, \, \delta\omega_{\mu} \} \Psi,
\end{equation}
where $\{ A, \, B \} = A B + B A$ (this anti-commutator comes from the covariant derivatives of $\Psi$ and $\bar{\Psi}$).  The first part is symmetric, since $\delta g^{\mu \nu}$ is symmetric.  As for the scalar field, the antisymmetric part is killed by the symmetric variation $\delta g^{\mu \nu}$.  So, from (3), we could write the following expression (remember that $\omega_{\mu} = -\, i \frac{1}{4} \, \sigma_{ab} \, \omega_{\mu}^{ab}$):
\begin{equation}
\begin{aligned}
T_{\mu \nu} \, \delta g^{\mu \nu} \equiv 2 \, \delta \mathscr{L} = i \frac{1}{4} \Bigl( \bar{\Psi} \, \Gamma_{\mu} \, D_{\nu} \, \Psi + \bar{\Psi} \, \Gamma_{\nu} \, D_{\mu} \, \Psi &- (D_{\mu} \bar{\Psi}) \, \Gamma_{\nu} \, \Psi - (D_{\nu} \bar{\Psi}) \, \Gamma_{\mu} \, \Psi \Bigr) \, \delta g^{\mu \nu} \\
&+ \frac{1}{4} \bar{\Psi} \{ \, \gamma^c, \, \sigma_{ab} \} \Psi \, e_c^{\mu} \, \delta \omega_{\mu}^{\; ab}.
\end{aligned}
\end{equation}
The last term can be expanded with some efforts (this is the most laborious part!).  Using the Dirac algebra, we find the spin density of the spinor field:
\begin{equation}\tag{8}
\frac{1}{4} \bar{\Psi} \{ \, \gamma^c, \, \sigma_{ab} \} \Psi \, e_c^{\mu} = \frac{1}{2} \, \varepsilon_{ab}^{\; \; cd} \, \bar{\Psi} \, \gamma_d \, \gamma^5 \, \Psi \, e_c^{\mu} \equiv S^c_{\; ab} \, e_c^{\mu} = S^{\mu}_{\; ab}.
\end{equation}
This expression is fully anti-symmetric in all the indices (this is crucial for the contraction on $\delta \omega_{\mu}^{\; ab}$).
Now, we have to find the variation of the spin-connexion, and this is where things are getting a bit tricky.  We assume the minimal coupling to gravity, i.e no torsion.  The connection is symmetric, as in classical general relativity (Levi-Civita connexion, from the symmetric Christoffel symbols):
\begin{equation}\tag{9}
{\omega_{\mu}}^a_{\; b} = e_{\lambda}^a \, \nabla_{\mu} \, e_b^{\lambda}.
\end{equation}
This thing is the simplest expression that defines the covariant derivative (2) under the local Lorentz transformations.  This is the place where torsion could enter the scene, but it isn't a necessity (contrary to what many authors may be saying).  The variation of (9) above is messy to do, but you could check the following result (using (4) and (5), and using the definition of the Christoffel symbols):
\begin{equation}\tag{10}
S^{\mu}_{\; ab} \, \delta \omega_{\mu}^{\; ab} \equiv {S^{\mu}}_{\lambda}^{\; \nu} \, \delta \Gamma_{\mu \nu}^{\lambda} - \frac{1}{2} \, S^{\mu \kappa \nu} \, \nabla_{\mu} \, \delta g_{\nu \kappa} = 0.
\end{equation}
This cancellation comes directly from the contraction of a fully anti-symmetric expression on two symmetric parts (${S^{\mu}}_{\lambda}^{\; \nu}$ on $\delta \Gamma_{\mu \nu}^{\lambda}$, and $S^{\mu \kappa \nu}$ on $\delta g_{\nu \kappa}$).  So the net result is that (3) gives our desired energy-momentum components:
\begin{equation}\tag{11}
T_{\mu \nu} \, \delta g^{\mu \nu} = i \frac{1}{4} \Bigl( \bar{\Psi} \, \Gamma_{\mu} \, D_{\nu} \, \Psi + \bar{\Psi} \, \Gamma_{\nu} \, D_{\mu} \, \Psi - (D_{\mu} \bar{\Psi}) \, \Gamma_{\nu} \, \Psi - (D_{\nu} \bar{\Psi}) \, \Gamma_{\mu} \, \Psi \Bigr) \, \delta g^{\mu \nu},
\end{equation}
so finally:
\begin{equation}\tag{12}
T_{\mu \nu} = i \frac{1}{4} \Bigl( \bar{\Psi} \, \Gamma_{\mu} \, D_{\nu} \, \Psi + \bar{\Psi} \, \Gamma_{\nu} \, D_{\mu} \, \Psi - (D_{\mu} \bar{\Psi}) \, \Gamma_{\nu} \, \Psi - (D_{\nu} \bar{\Psi}) \, \Gamma_{\mu} \, \Psi \Bigr).
\end{equation}
This is the well-known symmetric energy-momentum of the Dirac field, as defined from the canonical approach with the fancy Belinfante-Rosenfeld correction.
Take note that, except for the vanishing of $\mathscr{L}$ in (3) (for on-shell fields), we don't need to assume that the spinor field obeys its equation of motion.  The energy-momentum only gets a correction from the symmetric $\mathscr{L} \, g_{\mu \nu}$ if the field is off-shell (maybe is this related to the cosmological constant ?).  Finally, $T_{\mu \nu}$ is more complicated if we let torsion to kick in the game: the covariant derivatives get a contorsion tensor, but the full $T_{\mu \nu}$ stays symmetric.
