Does there exist a square root of Euler-Lagrange equations of a field? (Factorization) Does there exist a square root of Euler-Lagrange equations $\partial_{\mu}\frac{\partial \mathcal{L}}{\partial(\partial_\mu \phi)}-\frac{\partial \mathcal{L}}{\partial \phi} = 0$ in the sense that $(x+iy)(x-iy) = x^2+y^2$? By which I mean, does there exist another partial differential equation that squares to $\partial_{\mu}\frac{\partial \mathcal{L}}{\partial(\partial_\mu \phi)}-\frac{\partial \mathcal{L}}{\partial \phi} = 0$? Just like how Dirac operator $(i\gamma^\mu\partial_\mu-m)\Psi = 0$ squares to Klein-Gordon operator $(-\eta^{\mu\nu}\partial_\mu\partial_\nu + m)\psi = 0$. Can it be used to derive Dirac equation from Klein-Gordon Lagrangian?
 A: The "square root" in the sense of $-(a^2 +b^2)=(ia-b)(ia+b)\,\,$
of the Klein-Gordon's differential operator  is the Dirac differential operator. That is the main message:
$$(i\gamma^\mu\partial_\mu-m)(i\gamma^\nu\partial_\nu +m) = -\gamma^\mu\gamma^\nu\partial_\mu\partial_\nu -m^2 = -(\gamma^{\left[\mu\right.}\gamma^{\left.\nu\right]} +    \gamma^{\left(\mu\right.}\gamma^{\left.\nu\right)})\partial_\mu\partial_\nu -m^2 = -\gamma^{\left(\mu\right.}\gamma^{\left.\nu\right)} \partial_\mu\partial_\nu -m^2 = -\eta^{\mu\nu}\partial_\mu\partial_\nu -m^2 =  -(\Box +m^2)$$
In the algebra we threw away $\gamma^{\left[\mu\right.}\gamma^{\left.\nu\right]}\partial_\mu\partial_\nu$ because $\gamma^{\left[\mu\right.}\gamma^{\left.\nu\right]}$ is antisymmetric whereas the sucession of two partial derivatives is symmetric. It is also essential for this decomposition that the "mass operator" commutes with $i\gamma^\mu\partial_\mu$.
The analogy of a square root can be still extended a little bit further. Instead of writing the first term of the KG-equation (short Klein-Gordon: KG) $\partial^\mu \phi \partial_\mu\phi $ we can apply partial integration and write:
$$\partial^\mu \phi \partial_\mu\phi = -\phi \partial^\mu\partial_\mu \phi  +  \partial^\mu (\phi \partial_\mu \phi)$$
and just neglect the total partial derivative $\partial^\mu (\phi \partial_\mu \phi)$. Such a Lagrangian also leads to the KG-equation.
In that case the KG-Lagrangian looks like this:
$${\cal{L}} = -\phi (\Box + m^2)\phi$$
In order to reach as close as possible to the Dirac Lagrangian we actually better start off from the Lagrangian of the complex KG-field:
$${\cal{L}} = -\phi^\dagger (\Box + m^2)\phi$$
We replace  $\phi \rightarrow \psi$  and $\phi^\dagger \rightarrow \bar{\psi}$ as a first step and then choose of the 2 factors $i\gamma^\mu\partial_\mu\pm m$ the one with the minus sign :
$${\cal{L}} = \bar{\psi}(i\gamma^\mu\partial_\mu-m)\psi$$
Therefore, it is not the Lagrangian that is square-rooted, but the differential operator inside.
Needless to say that we cannot just keep the complex scalar field sandwiching the differential operator as the $\gamma^\mu$ matrices have not only vector character but are also "bispinor tensors"(the indices are usually suppressed in order to avoid clutter) which requires spinors as sandwich in order to make the Lagrangian Lorentz-covariant.
It was the genius of P.Dirac to recognize this thereby giving up the scalar character of the involved fields.
A: The Euler Lagrange equations are the equations used to minimise the action. So, they do not exist in versions (i.e. square root of EL eqs etc). What you should be after if you would want (for some reason) a Dirac equation for a scalar field is a Lagrangian whose EL EoM is the Dirac equation.
This is however troublesome, as in the Dirac equation you have matrices and the solutions are spinor fields, which is not in agreement with having scalar fields which is your case now
