Inteprertation of renormalisation of composite operators

The fact that renormalization mixes a composite operator with other operators can be understood using lattice QFT, which is manifestly finite from the beginning. Lattice QFT is messy, but for this question, the only important thing is that the lattice provides a non-perturbative cutoff. Renormalization is about changing the scale of the cutoff, and the goal is to understand what this does to a composite operator like $phi^2$.

Work in euclidean signature for simplicity. In Wilson's picture, renormalization amounts to lowering the cutoff from $Lambda_H$ to $Lambda_L$ by integrating out modes with wavenumbers between those two scales. The subscripts $H$ and $L$ stand for high and low. Write $$ phi(x) = phi_H(x)+phi_L(x), tag{1} $$ where $phi_{H}(x)$ is the part involving only wavenumbers between the two cutoffs, and $phi_L(x)$ is the part involving only wavenumbers below $Lambda_L$. The correlation functions of interest are generated by $$ Z[J,K]propto int [dphi] expleft(-S[phi] +int phi(x) J(x) +int phi^2(x) K(x)right). tag{2} $$ The quadratic source term $intphi^2 K$ is redundant in equation (2), because we can generate insertions of $phi^2(x)$ either by taking one derivative with respect to $K(x)$ or by taking two derivatives with respect to $J(x)$.

When we integrate out modes with wavenumbers between $Lambda_H$ and $Lambda_L$, we're agreeing to consider only correlation functions of operators that have been smeared enough to avoid resolving anything above the scale $Lambda_L$. The key point is that smearing does not commute with squaring, so inserting $phi^2(x)$ in a correlation function and then smearing it is not the same as inserting $phi_L^2(x)$ in a correlation function: $$ big(phi^2(x)big)_Lneq phi_L^2(x). tag{3a} $$ More explicitly, $$ int dx' f(x')phi^2(x+x')neq left(int dx' f(x')phi(x+x')right)^2 tag{3b} $$ where the integral with kernel $f$ is designed to eliminate the high-wavenumber modes. The generating function for correlation functions of sufficiently-smeared operators is obtained from (2) by setting $J_H=0$ and $K_H=0$, which prevents us from inserting any high-resolution operators. This leaves $$ Z[J_L,K_L]propto int [dphi] expleft(-S[phi] +int phi(x) J_L(x) +int phi^2(x) K_L(x)right). tag{4} $$ Modes with different wavenumbers are orthogonal to each other, so we have $$ Z[J_L,K_L]propto int [dphi] expleft(-S[phi] +int phi_L(x) J_L(x) +int big(phi^2(x)big)_L K_L(x)right). tag{5} $$ Thanks to the inequality (3), the quadratic source term $int phi^2 K_L$ is no longer redundant.

When we integrate over the high-wavenumber modes $phi_H$, the $J_L$-source term is not affected, because it doesn't involve $phi_H$. In contrast, the $K_L$ source term is affected, because the inequality (3) says that the composite operator $big(phi^2(x)big)_L$ does involve $phi_H$. The details depend on the action $S[phi]$, but even without working through those details, the inequality (3) already suggests that renormalizing a composite operator will mix it with other operators.