Your answer is not going to be the same as the answer in the paper, because you are calculating the canonical stress-energy tensor, which is conserved, but which is not symmetric, and which has a complicated relation to the angular momentum tensor. The issue is that there two different conserved quantities, the stress energy tensor, and the angular momentum tensor, and the information in the two partially overlaps for a theory with both rotational/Lorentz and translation invariance.
The standard easy fix is to calculate the stress tensor by differentiating with respect to $g_{\mu\nu}$. This is how you automatically get a symmetric stress tensor with the right properties that it makes the angular momentum tensor by just multiplying by x factors
$$ L_{\mu\nu\alpha} = x_{\alpha} T_{\mu\nu} - x_{\mu}T_{\alpha\nu} $$
In some convention, and assuming you are using a symmetric stress tensor. The derivative of the action with respect to $g_{\mu\nu}$ gives Harvey's result.
I should point out that you missed a term of $-\eta^{\mu\nu} L$ in your stress tensors, both in yours and in Harvey's. This doesn't affect the question.
If you don't impose power-counting renormalizability, there are a host of other possibilities, since higher order derivatives or higher order interactions can be introduced. For example, terms $(Tr(F^2)^m)^n$ and are gauge invariant but for $m>1$ or $n>1$ not renormalizable.
If you impose power-counting renormalizability, uniqueness is fairly straightforward up to trivial field transformations. To see this one first looks at monomials - products of fields and their derivatives. By renormalizability, the total degree is not allowed to be bigger than 4. Each partial derivative $d_j=\partial_j$ counts as degree 1, each Bose fields $A_j$ as degree 1, and each Fermion field $\psi_j$ as degree 3/2. Moreover, Fermions must appear an even number of times to yield a scalar Lagrangian. This leads to a quite short list of possibilities: Up to 4 $A$s and $d$s, or $\psi\psi, d\psi\psi, A\psi\psi$, all with all possible indices. The general renormalizable local Lagrangian density is a linear combination of these, at fixed $x$. Now impose Poincare invariance and gauge invariance, and the only linear combinations left are the ones that one sees everywhere. For a single Yang-Mills field and nothing else (i.e., your question in the narrow sense)
the only freedom left is rescaling the fields, which eliminates an arbitrary factor in front of the trace. In the presence of fermionic fields there is the additional freedom to take linear combinations of Fermionic fields as new fields, which may be used to reduce the associated bilinear forms to weighted sums of squares.
If one drops gauge invariance, there are lots of other possible Langangian densities, for example a mass term, products of it with the terms described, and even more.
Note that proving the renormalizability of nonabelian gauge theories with broken symmetries was a highly nontrivial achievement (around hundred pages of published argument) worthy of a nobel prize for Veltman nd 't Hooft. Thus it is unreasonable to explain in an answer the reasons for where precisely the borderline is bewteen renormalizable and nonrenormalizable.
The answer to your question, ''Maybe I can put my question in simpler terms: is there any room for modifications in the Standard Model without introducing new fields? can we add new interactions among the gauge bosons (W,Z,…) and/or the matter fields without spoiling unitarity, covariance, or renormalisability? (at the perturbative level at least; here I dont care about θ terms, etc.)'' related to the bounty (which will disappear in a few hours)
is no, essentially by an extension of the above reasoning (including the 100 pages of proof of renormalizability).
Best Answer
The current is just the LHS of your eq. $(6)$. In fact, the antisimmetry of $F^{\mu \nu}$ implies $\partial _\mu \partial _\nu F^{\mu \nu}=0$.
Observe that the gauge fields carry a non-zero charge, since they transform non-trivially under global transformations. In fact, they transform according to the adjoint representation of the group of simmetry. The part of the current which involves only $\psi$ is the covariant current $\mathscr J ^\mu$, which enters in:$$D_\mu F^{\mu \nu} =\mathscr J ^\nu.$$ This current (differently from the Noether current $J^\mu$) transforms like the Yang-Mills tensor under gauge transformations.