In the canonical quantization of a massive Lagrangian, we obtain equations of motion via the Euler-Lagrange formalism. As I have read, these can become rather long but decouple into the Klein-Gordon (or Dirac) plus a set of lower-order equations. E.g., for the spin-1 Proca field: $$ (\Box + m^2) A^\mu - \partial^\mu (\partial^\nu A_\nu) = 0 \Rightarrow (\Box + m^2)A^\mu = 0,~ \partial^\mu A_\mu = 0. $$ As I see it, there is no obvious way of going back, i.e., deriving the coupled equations from the decoupled ones.
When constructing quantum fields of spin $s$ via irreducible representations of the Poincaré group (let's call it the Weinberg approach), we end up with completely symmetric tensors (or spinors, in the half-integer spin case) with $s$ Lorentz indices. As was shown in this answer, it follows that these fields fulfill a number of differential equations, namely:
$$ (\Box + m^2) \Phi_{\mu_1\cdots\mu_s} = 0 \text{ or } (i\not\partial + m) \Phi_{\mu_1\cdots\mu_s} = 0 ~,~~ \partial^{\mu_1} \Phi_{\mu_1\cdots\mu_s} = 0 \text{ or } \gamma^{\mu_1} \Phi_{\mu_1\cdots\mu_s} = 0~,~~ \eta^{\mu_i\mu_j} \Phi_{\mu_1\cdots\mu_s} = 0, $$
depending on whether it is a tensor or tensorial spinor. So here, we directly obtain the decoupled equations of motion.
My question now is: Are the "coupled" equation of motion an artifact of the Lagrange formalism or is there a way to derive them from the Weinberg approach?
I am curious about this because time-ordered products of the fields do have a special relation to the coupled equations of motion. They are the respective greens function. E.g., for the Proca field: $$ \Big(\eta^{\mu\nu} (\Box + m^2) - \partial^\mu \partial^\nu \Big) \langle T A_\nu(x) A_\rho(x') \rangle = -i \eta^\mu_\rho \delta(x-x'). $$ This does not necessarily follow from the decoupled equations, even though we have $D(\partial) \langle T A_\nu(x) A_\rho(x') \rangle \propto \delta(x-x')$ for $D = (\Box + m^2), \partial^\mu, ... $.
I am happy about any insight on how these properties are related!
[Edit:] With some help, I was able to narrow down what I am looking for. Let me try to (hopefully) clarify.
It is always possible to combine the "decoupled" equations such that the resulting equations are equivalent to the ones we started with. The equations of motion resulting from the Lagrangian are special because they allow to describe a deformation with an interaction Lagrangian $\mathcal{L}$ via a simple Dyson-series(?) (due to renormalization, this holds only when coupling linealy to an external potential, more general interactions will deviate), e.g. for spin-0 or spin-1: $$ (\Box + m^2)\phi_{int} = - \Big( \frac{\partial \mathcal{L}}{\partial \phi}\Big)_{int},~ \Big(\eta^\mu_\nu (\Box + m^2) - \partial^\mu \partial_\nu \Big) A^\nu = - \Big( \frac{\partial \mathcal{L}}{\partial A_\mu}\Big)_{int}. $$ The $_{int}$ here denotes a perturbative series in the propagators, here the time-ordered products (this corresponds to the equation for the propagator of $A^\mu$ I wrote above).
So my question would rather be:
Is there a generalization of this formalism to higher-spin fields? Are there still equations of motion and time-ordered products that act as their Green's functions?
Thanks again for any help!