How to integrate out the $W$-boson fields? What does it mathematically mean to 'integrate out' the $W$-boson fields to obtain the Fermi Lagrangian from the electroweak theory? How does one achieve this mathematically? It will be helpful if someone can explain this both in the path integral formalism and the operator formalism of quantum field theory.
In particular, I know that the $W$-boson propagator, $\frac{1}{p^2-m_w^2}$ in the limit $p^2\ll m_w^2$ becomes $-\frac{1}{m_w^2}$. But why is such an approximation called 'integrating out' the fields? What are we really doing to the EOM when we are making this approximation?
 A: In classical physics, equations of motion are obtained by varying an action $S$. When we quantise, we consider a path integral $Z\left[\varphi\right]=\int\mathcal{D}\varphi \text{e}^{iS\left[\varphi\right]}$, a function integral over a collection $\varphi$ of all fields, so that $S\left[\varphi\right]$ depends only on $\varphi$ and its derivatives. Then operators have means$$\left\langle\mathcal{O}\right\rangle=\frac{\int\mathcal{D}\varphi\mathcal{O} \text{e}^{iS\left[\varphi\right]}}{\int\mathcal{D}\varphi \text{e}^{iS\left[\varphi\right]}}.$$As the definition of functional integration is in its infancy, the starting point is to generalise the famous identity $$\int_{\mathbb{R}^N} d^Nx\text{e}^{-\frac{1}{2}x^TAx+J^Tx}=\frac{\left(2\pi\right)^{N/2}e^{\frac{1}{2}J^2}}{\sqrt{\det A}}$$for $N\times N$ symmetric invertible real matrices $A$. Taking the $N\to\infty$ limit of an integration measure designed with a normalisation that cancels the power of $2\pi$, we obtain a functional integral I'll get to in a moment. First, let's note that the dot-product $U^TV=\sum_i U_iV_i$ is generalised with an integral, so the result is$$\int\mathcal{D}\varphi\text{e}^{\int\left(-\frac{1}{2} \varphi A\varphi+J\varphi\right)d^Dx}=\text{e}^{\frac{1}{2}J^2}\sqrt{\det D}$$with $D:=A^{-1}$ for an invertible operator $A$, whatever an operator "determinant" may be! (This assumes a $D$-dimensional spacetime in which $\varphi$ lives. The operator $D$ is, of course, a propagator.) A similar result is obtainable with the replacement $A\to iA$, and we can use this to consider many integrals of the form $\int\mathcal{D}\varphi\text{e}^{i\int d^Dx\mathcal{L}\left[\varphi,\,J\right]}$, with $\mathcal{L\left[\varphi,\,J\right]}$ a Lagrangian density dependent on the fields and (this is where we get another physical difference) a current $J$. When I say we can consider these integrals, I mean that taking ratios to compute operator means is now possible. As an example, if $\mathcal{O}=e^{-aW^2}$ for a constant $a$ and field $W$ included in $\varphi$, $\left\langle\mathcal{O}\right\rangle$ is computed as a ratio of two integrals, and under the right circumstances the determinants helpfully cancel. More generally, a ratio of the determinants of two operators survives, which we can treat by generalising the matrix result$$\frac{\det \left(A+\epsilon\right)}{\det A}=\det \left(1+D\epsilon\right)\approx\exp\text{tr}D\epsilon$$with $\epsilon$ "small" to a functional result,$$\frac{\det \left(A+\epsilon\right)}{\det A}=\exp{\int d^D x D\epsilon}.$$(For your question, the "square" $W^2:=W_\mu W^\mu$.) Note also that the parts of $\phi$ on which $\mathcal{O}$ doesn't depend needn't be integrated over in either the numerator or denominator. Similarly, we can "throw away" $W$ instead when computing the mean  of a variable on which it does not depend.
A: The tree level 4-fermion amplitude you get when you collapse the propagator to a point is a $J_\mu^+ J^{\mu ~-}$ contact term, the 4-Fermi interaction of the 30s.
With the advent of functional integrals to the SM, it became easier to understand the UV origin of this low energy effective theory. The relevant part of the SM action contributing to this charged current amp is 
$$
{\cal L}_{eff}= m^2_W W_{\mu}^+ W^{\mu ~-}+ \frac{g}{\sqrt 2}( W_{\mu}^+ J^{\mu -} +W^{\mu ~-} J_\mu^+) +O(p^2/M_W^2).\tag{1}
$$
Since the Ws do not have derivatives in this part of the action, they are superfluous fields with algebraic equations of motion and may be eliminated, by outright performing the functional integration involving them as variables, as follows.
Complete the complex square,
$$
{\cal L}_{eff} = m^2_W \left (W_{\mu}^+  +\frac{g}{\sqrt {2} m_W^2} J^{\mu~+}\right) \left(W_{\mu}^-  +\frac{g}{\sqrt {2} m_W^2}J^{\mu~-} \right ) -\frac{g^2}{2m_W^2} J_{\mu}^{+} J^{\mu~-} .
$$
Now observe the first term represents a shift in the definition of the Ws, that is, they may be redefined to absorb the current pieces;
integrated over spacetime and stuck in the exponent of the functional integral, the first term amounts to two Gaussians, when you resolve it into the original "new", shifted  W1, W2; functional integration of these Gaussians w.r.t. the shifted Ws leaves no trace of the Ws in this low energy part of the path integral. They have been "integrated out", as per your question. 
The sole usable residue of their presence, is the "constant" (as far as W degrees of freedom are concerned) second term, the current-current interaction, $-\frac{2G_F}{\sqrt 2}  J_{\mu}^{+} J^{\mu~-} $, where one defined 
$G_F\sqrt 2 \equiv g^2/4m_w^2=2/v^2$. Note you would get the very same answer from the merely algebraic equations of motion of ${\cal L}_{eff}$, namely $W_\mu^{\pm}=-g J_\mu^{\pm}/\sqrt{2} m_w^2$; using these to eliminate the Ws would result in the same current-current residual interaction.
(By the way, in 1933, this was essentially the first application of QFT: its linchpin  feature of creation and annihilation of fermion species.) 
A very analogous procedure, naturally, obtains for the neutral current amplitudes involving Z exchange, mutatis mutandis. 
A: Heuristically, there are two ways. 


*

*In the electroweak Lagrangian, you substitute for W by its classical solution. That is, you find the Euler-Lagrange equation for W, solve for it, and plug the solution in Lagrangian. In doing this, the W fields are "integrated out". 

*When an internal line of the W boson is present in a Feynman diagram, you erase it. For instance, the tree level diagram of electron scattering with neutrino involves a W propagator. Erasing it, you get a four-point vertex of fermions. In terms of Lagrangian, you have an operator of four spinors, which is a dimension 6 operator. The numerical prefactor of it can be fixed by matching the scattering amplitude. 
