Symmetrizing the Canonical Energy-Momentum Tensor The Canonical energy momentum tensor is given by 
$$T_{\mu\nu} = \frac{\partial {\cal L}}{\partial (\partial^\mu \phi_s)} \partial_\nu \phi_s - g_{\mu\nu} {\cal L}. $$
A priori, there is no reason to believe that the EM tensor above is symmetric. To symmetrize it we do the following trick.
To any EM tensor we can add the following term without changing its divergence and the conserved charges:
$${\tilde T}_{\mu\nu} = T_{\mu\nu} + \partial^\beta \chi_{\beta\mu\nu}, $$
where $\chi_{\beta\mu\nu} = - \chi_{\mu\beta\nu}$. The antisymmetry of $\chi$ in its $\mu\beta$ indices implies that ${\tilde T}_{\mu\nu}$ is conserved. Also, all the conserved charges stay the same. 
Now even though $T_{\mu\nu}$ is not a symmetric tensor, it is possible to choose $\chi_{\beta\mu\nu}$ in such a way so as to make ${\tilde T}_{\mu\nu}$ symmetric. It can be shown that choosing
$$\chi_{\lambda\mu\nu} = - \frac{i}{2}\left[ \frac{\delta {\cal L}}{\delta (\partial^\mu \phi_r) } (I_{\nu\lambda})_{rs} \phi_s +   \frac{\delta {\cal L}}{\delta (\partial^\lambda \phi_r) } (I_{\mu\nu})_{rs} \phi_s + \frac{\delta {\cal L} }{\delta (\partial^\nu \phi_r) } (I_{\mu\lambda})_{rs} \phi_s  \right]$$
makes the new EM tensor symmetric. Here $(I_{\mu\nu})_{rs}$ is the representation of the Lorentz Algebra under which the fields $\phi_s$ transform. 
Here's my question - Is it possible to obtain the symmetric EM tensor directly from variational principles by adding a total derivative term to the Lagrangian. In other words, by shifting ${\cal L} \to {\cal L} + \partial_\mu X^\mu$, and choosing $X^\mu$ appropriately, can we exactly get the shift in the EM tensor required, in order to make the canonical EM tensor symmetric?
What I've done so far - It is possible to show that under a shift in the Lagrangian by a total derivative, one shifts the EM tensor by $T_{\mu\nu} \to T_{\mu\nu} + \partial^\lambda \chi_{\lambda\mu\nu}$ where
$$\chi_{\lambda\mu\nu} =   \frac{1}{2} \frac{\delta X_\lambda}{\delta (\partial^\mu \phi_r)}  \partial_\nu \phi_r   -  \frac{1}{2}   \frac{ \delta X_\mu }{\delta (\partial_\lambda \phi_r)}    \partial_\nu \phi_r   +  X_\mu g_{\lambda\nu} -    X_\lambda g_{\mu\nu} \,. $$
What I wish to do next - I now have a differential equation that I wish to solve:
\begin{align}
&\frac{1}{2} \frac{\delta X_\lambda}{\delta (\partial^\mu \phi_r)}  \partial_\nu \phi_r   -  \frac{1}{2}   \frac{ \delta X_\mu }{\delta (\partial_\lambda \phi_r)}    \partial_\nu \phi_r   +  X_\mu g_{\lambda\nu} -    X_\lambda g_{\mu\nu} \\
&~~~~~~= - \frac{i}{2}\left[ \frac{\delta {\cal L}}{\delta (\partial^\mu \phi_r) } (I_{\nu\lambda})_{rs} \phi_s +   \frac{\delta {\cal L}}{\delta (\partial^\lambda \phi_r) } (I_{\mu\nu})_{rs} \phi_s + \frac{\delta {\cal L} }{\delta (\partial^\nu \phi_r) } (I_{\mu\lambda})_{rs} \phi_s  \right] \,.
\end{align}
Any ideas on how to solve this?
 A: OP's question (v7) asks:

Is it possible to obtain a symmetric stress-energy-momentum (SEM) tensor directly from the canonical SEM tensor by adding a total derivative term to the Lagrangian? In other words, by shifting $\Delta{\cal L}=d_\mu X^\mu$, and choosing $X^\mu$ appropriately, can we exactly get the shift in the SEM tensor required, in order to make the canonical SEM tensor symmetric?

No, that project is doomed already for E&M with the Maxwell Lagrangian density
$$ {\cal L}_0~:=~ -\frac{1}{4}F_{\mu\nu}F^{\mu\nu}\tag{1} $$
with
$$ F_{\mu\nu}~=~A_{\nu,\mu}-A_{\mu,\nu}, \qquad 
\frac{\partial{\cal L}_0}{\partial A_{\mu,\nu}}~\stackrel{(1)}{=}~ F^{\mu\nu}.\tag{2}$$
The vacuum EL equations read
$$ 0~\approx~F^{\mu\nu}{}_{,\nu}~=~  d^{\mu}(A^{\nu}_{,\nu})-d_{\nu}d^{\nu}A^{\mu}\tag{3} $$
In E&M, the canonical SEM tensor is$^1$
$$\begin{align} \Theta^{\mu}{}_{\nu}~:=~&\delta^{\mu}_{\nu}{\cal L}_0+\left(-\frac{\partial{\cal L}_0}{\partial A_{\alpha,\mu}}+ d_{\beta}\frac{\partial{\cal L}_0}{\partial A_{\alpha,\mu\beta}}\right) A_{\alpha,\nu}
-\frac{\partial{\cal L}_0}{\partial A_{\alpha,\mu\beta}}A_{\alpha,\nu\beta}\cr 
~\stackrel{(1)}{=}~&\delta^{\mu}_{\nu}{\cal L}_0+F^{\mu\alpha}A_{\alpha,\nu}~,\end{align}\tag{4}$$
while the symmetric SEM tensor is
$$ T^{\mu}{}_{\nu}~=~\delta^{\mu}_{\nu}{\cal L}_0+F^{\mu\alpha}F_{\nu\alpha}.\tag{5}$$
So the difference is$^2$
$$\begin{align}  T^{\mu}{}_{\nu} -\Theta^{\mu}{}_{\nu}~\stackrel{(4)+(5)}{=}&~ F^{\mu\alpha}A_{\nu,\alpha}~=~ d_{\alpha}(F^{\mu\alpha}A_{\nu}) - \underbrace{F^{\mu\alpha}{}_{,\alpha}}_{~\approx~0}A_{\nu} \cr
~\stackrel{?}{\approx}~&\delta^{\mu}_{\nu}\Delta{\cal L}+\left(-\frac{\partial\Delta{\cal L}}{\partial A_{\alpha,\mu}}+ d_{\beta}\frac{\partial\Delta{\cal L}}{\partial A_{\alpha,\mu\beta}}\right) A_{\alpha,\nu}
-\frac{\partial\Delta{\cal L}}{\partial A_{\alpha,\mu\beta}}A_{\alpha,\nu\beta}\end{align}\tag{6} $$
for some total derivative term $\Delta{\cal L}=d_\mu X^\mu$, where $X^\mu$ depends on $A$ and $\partial A$. The question mark (?) in eq. (6) is OP's question. Note that the continuum equation is unaltered on-shell
$$ d_{\mu}T^{\mu}{}_{\nu} ~\approx~ d_{\mu}\Theta^{\mu}{}_{\nu}~\approx~0. \tag{7} $$
For dimensional reasons $X^\mu$ must be on the form$^3$
$$  X^{\mu}
~=~ a A^{\mu} A^{\nu}_{,\nu} + b  A^{\nu}  A^{\mu}_{,\nu} + c  A^{\nu} A_{\nu}^{,\mu}\tag{8} $$
for some constants $a,b,c$. Then
$$\begin{align} \Delta{\cal L}~&~=~d_\mu X^\mu
~\stackrel{(8)+(10)}{=}~\Delta{\cal L}_1+\Delta{\cal L}_2,\tag{9} \cr
\Delta{\cal L}_1~&:=~a (A^{\mu}_{,\mu})^2 
+ b A^{\nu}_{,\mu} A^{\mu}_{,\nu} 
+ c A^{\nu}_{,\mu} A_{\nu}^{,\mu},\tag{10} \cr
\Delta{\cal L}_2~&:=~
(a+b) A^{\mu} A^{\nu}_{,\nu\mu} 
+ c A^{\mu}A_{\mu,\nu}^{,\nu}~\stackrel{(3)}{\approx}~(a+b+c) A^{\mu} A^{\nu}_{,\nu\mu}.\tag{11} \end{align} $$
Consider the last term on the right-hand side of eq. (6):
$$\begin{align}\frac{\partial\Delta{\cal L}}{\partial A_{\alpha,\mu\beta}}A_{\alpha,\nu\beta}
&~=~\frac{\partial\Delta{\cal L}_2}{\partial A_{\alpha,\mu\beta}}A_{\alpha,\nu\beta}\cr
&~=~\frac{a+b}{2} \left(A^{\alpha} A^{\mu}_{,\alpha\nu}+A^{\mu} A^{\alpha}_{,\alpha\nu}\right) +c A^{\alpha} A^{,\mu}_{\alpha,\nu} \tag{12}\end{align}$$
Apart from the diagonal term $\delta^{\mu}_{\nu}\Delta{\cal L}_2$, the terms in eq. (12) are the only appearances of 2nd-derivatives on the right-hand side of eq. (6). We conclude that
$$ \Delta{\cal L}_2~=~0\qquad\Leftrightarrow\qquad a+b~=~0\quad\wedge\quad c~=~0.\tag{13}$$
Similar arguments shows that eq. (6) is not possible$^4$. $\Box$
--
$^1$ In eq. (4) we have indicated the canonical SEM tensor for a Lagrangian density with up to 2nd-order derivatives. Some references, e.g. Weinberg QFT, have the opposite notational conventions for $T\leftrightarrow\Theta$. Here we are using the $(-,+,+,\ldots,+)$ Minkowski sign convention.
$^2$ In formula (6) we have neglected terms in $\Delta{\cal L}$ that depends on  $\partial^3A$, $\partial^4A$, $\partial^5A$, $\ldots$, etc. Such terms are excluded for various reasons.
$^3$ In retrospect, this answer completely shares the premise/ideology/program/conclusion of this Phys.SE post.
$^4$ Interestingly, if we just take the trace of eq. (6), we get
$$\begin{align}  A^{\nu}_{,\mu} A^{\mu}_{,\nu} 
- A^{\nu}_{,\mu} A_{\nu}^{,\mu} &~=~F^{\mu\alpha}A_{\mu,\alpha}\cr
&~\stackrel{?}{\approx}~n \Delta{\cal L}+\left(-\frac{\partial\Delta{\cal L}}{\partial A_{\alpha,\mu}}+ d_{\beta}\frac{\partial\Delta{\cal L}}{\partial A_{\alpha,\mu\beta}}\right) A_{\alpha,\mu}
-\frac{\partial\Delta{\cal L}}{\partial A_{\alpha,\mu\beta}}A_{\alpha,\mu\beta}\cr
&~\stackrel{(9)}{=}~(n-2) \Delta{\cal L}_1+(n-1) \Delta{\cal L}_2+ A_{\alpha,\mu} d_{\beta}\frac{\partial\Delta{\cal L}_2}{\partial A_{\alpha,\mu\beta}}\cr
&~\stackrel{(11)}{=}~(n-2) \Delta{\cal L}_1+(n-1) \Delta{\cal L}_2+  \frac{a+b}{2}\left((A^{\mu}_{,\mu})^2+A^{\nu}_{,\mu} A^{\mu}_{,\nu} \right) +c A^{\nu}_{,\mu} A_{\nu}^{,\mu} ,\tag{14}\end{align} $$
which leads to the linear eq. system
$$ \begin{align}
0&~=~a+b+c, \tag{15}\cr
-1&~=~(n-1)c\qquad\qquad\qquad\Rightarrow\qquad c~=~-\frac{1}{n-1},\tag{16}\cr
0&~=~(n-2)a +\frac{a+b}{2}\qquad\Rightarrow\qquad a~=~-\frac{1}{2(n-1)(n-2)},\tag{17}\cr
1&~=~(n-2)b +\frac{a+b}{2}\qquad\Rightarrow\qquad b~=~\frac{2n-3}{2(n-1)(n-2)},\tag{18}\end{align} $$
which remarkably has a unique & consistent solution. So it is not enough to just take the trace of eq. (6). However together with eq. (13), we conclude that there is no solution. $\Box$
A: I will try to obtain the result using another way.
It is well known that the Lagrangian density determined up to divergence of some four-vector $\mathcal{L}(x)\to \mathcal{L}(x)+\partial_\mu\psi^\mu(x)$
Let's understand what contribution the second term gives in the energy-momentum tensor.
$$\hat T^\nu _\mu=\partial_\rho\Bigr(\frac{\delta {\cal\psi^\rho }}{\delta (\partial^\mu \phi_s)} \partial^\nu \phi_s - g_\mu^\nu {\cal \psi^\rho}\Bigr)=\partial_\rho \chi^{\rho\nu}_ \mu 
$$
$\psi^\rho$ is arbitrary four-vector, contains in $\phi_s$ and $\partial^\rho\phi_r$. Set that $\psi^\rho=f(\phi^2)\phi_r\partial^\rho\phi_r$.(If I require that the Lagrangian dependence only $\phi_r$ and first derivative of it. It will be general form) We obtain the following result
$$
\chi^{\rho\nu}_ \mu =g_\mu^\rho {\cal \psi^\nu}-g_\mu^\nu {\cal \psi^\rho}
$$
where $g^\mu_\rho=\delta^\mu_\rho$ is a Kronecker symbol. Thus we obtain that energy momentum tensor defined up to such term $T^{\mu\nu}\to T^{\mu\nu}+\partial_\rho\chi^{\rho\mu\nu}$ where $\chi^{\rho\mu\nu}=-\chi^{\mu\rho\nu}$. This fact is a consequence of Lagrangian feature(The Lagrangian density determined up to divergence of some four-vector $\mathcal{L}(x)\to \mathcal{L}(x)+\partial_\mu\psi^\mu(x)$).
Edit
Using previous formula it is easy to obtain that 
$$
\chi^{\mu\rho\nu} =g^{\mu\rho} {\cal \psi^\nu}-g^{\mu\nu} {\cal \psi^\rho}
$$
After contraction with $g_{\mu\rho}$ we obtan that
$$
\psi^\nu=\frac{1}{D-1}\chi^{\mu\rho\nu}g_{\mu\rho}
$$
where $D$-is dimensionality of space.
In spite of the Lagrangian contain second derivatives, all of it is true.Because it Lagrangian differ only for full derivative. If you interesting in this question, you should write General relativity. Because action of General relativity which contain Riemann curvature tensor(which contain second derivatives).
A: It is possible to choose a Lagrangian such that the Noether energy momentum tensor is symmetric, namely $${\cal L} = \frac{\epsilon_0}{2} \partial_\mu A_\nu \partial^\mu A^\nu \,.$$ This Lagrangian however differ from the standard one $${\cal L} = \frac{\epsilon_0}{4} F_{\mu\nu} F^{\mu\nu} \,.$$ by the term $${\cal L} = \frac{\epsilon_0}{2} \partial_\mu A_\nu \partial^\nu A^\mu \,.$$ which is not a total derivative.
