Consider an (active) infinitesimal Lorentz transformation:
$$ x^\mu \rightarrow x^\mu + {\omega^\mu}_\nu x^\nu, $$
so that any scalar field is transformed as
$$ \phi(x) \rightarrow \phi'(x) = \phi(x) - {\omega^\mu}_\nu x^\nu \partial_\mu \phi(x) + O(\omega^2). $$
Now consider a Lagrangian density function $\mathcal{L(\phi, \partial\phi)}$ (with no explicit spacetime dependence). Every scalar field is associated to a Lagrangian density field $\mathcal{L}[\phi](x) := \mathcal{L}(\phi(x), \partial\phi(x))$, which is itself a scalar field. Therefore, it transforms with variation:
$$ \delta \mathcal{L} = -{\omega^\mu}_\nu x^\nu \partial_\mu \mathcal{L}[\phi] = -\partial_\mu ({\omega^\mu}_\nu x^\nu \mathcal{L}[\phi]), \tag{1}$$
where the second equality arises because $\omega$ is antisymmetric. Since the Lagrangian only varies by a four-divergence, the action is unchanged. This makes perfect sense: all we've done is shift around spacetime by an orthogonal transformation, moving around the $d^4 x\ \mathcal{L}$ terms in the action integral, so the total action integrated over the whole of spacetime isn't going to change. So far so good.
The problem arises when I try to calculate $\delta \mathcal{L}$ a different way. I think we should be able to calculate the variation using:
\begin{align} \delta \mathcal{L} & = \mathcal{L}[\phi'] - \mathcal{L}[\phi] \\ & = \frac{\partial\mathcal{L}}{\partial\phi}\delta\phi + \frac{\partial\mathcal{L}}{\partial(\partial_\mu \phi)} \delta(\partial_\mu \phi) \\ & = - \left [ \frac{\partial\mathcal{L}}{\partial\phi} {\omega^\alpha}_\beta x^\beta \partial_\alpha \phi + \frac{\partial\mathcal{L}}{\partial(\partial_\mu \phi)} \partial_\mu ({\omega^\alpha}_\beta x^\beta \partial_\alpha \phi) \right ] \\ & = - {\omega^\alpha}_\beta x^\beta \left[ \frac{\partial\mathcal{L}}{\partial\phi}\partial_\alpha \phi + \frac{\partial\mathcal{L}}{\partial(\partial_\mu \phi)} \partial_\mu \partial_\alpha \phi \right] - \frac{\partial \mathcal{L}}{\partial (\partial_\mu \phi)}{\omega^\alpha}_\mu \partial_\alpha \phi \\ & = -\partial_\alpha({\omega^\alpha}_\beta x^\beta \mathcal{L}) - \frac{\partial \mathcal{L}}{\partial (\partial_\mu \phi)}{\omega^\alpha}_\mu \partial_\alpha \phi.\tag{2} \end{align}
As you can see, there's a second term that's appeared out of nowhere! Where have I gone astray?
I've tried double-checking that $\delta (\partial_\mu \phi) = \partial_\mu (\delta \phi)$ and it seems to work out, so I don't think that's the problem. I found this very old post but I find the argument that
$$ \frac{\partial \mathcal{L}}{\partial (\partial_\mu \phi)} \propto \partial^\mu \phi \tag{3}$$
unconvincing. Say, for example, you had a $\frac{1}{2} \partial_\alpha \partial_\beta A^{\alpha\beta}$ term - that's fine, because differentiating w.r.t. $\partial_\mu \phi$ symmetrises $A$, so it works out that the extra term is zero, but the proportionality doesn't hold. So is it possible, in complete generality, to prove that $\frac{\partial \mathcal{L}}{\partial(\partial_\mu \phi)} \partial_\nu \phi$ is always symmetric in $\mu$ and $\nu$?