(Note: Those who are already familiar with the notion of OPE's might prefer to skip directly to the Question section below)
Setup: Consider, for simplicity, the quantum theory of the Euclidean K-G field. That is, the "quantization" of the the classical field which obeys the equations of motion,
\begin{align} (\partial^2-m^2)\phi=0, \end{align} where $\partial^2\equiv \sum^4_{i=1} \partial^2/\partial (x^i)^2$ with $x^i$ denoting elements of 4-dimensional Cartesian coordinates.
Following the "canonical" quantization procedure for the free field a la Peskin and Schroeder, I find the 2-point (Schwinger?) function
\begin{align} \langle 0| \phi(x_1)\phi(x_2)|0\rangle=\Delta_E(x_1-x_2), \tag{1} \label{2point} \end{align} where $\Delta_E$ denotes the Euclidean propagator/Green's function,
\begin{align} \Delta_E (x_1-x_2)= \int \frac{d^4p}{(2\pi)^4}\frac{e^{i p\cdot (x_1-x_2)}}{p^2+m^2}, \end{align} where $p,x_1,x_2$ are 4-vectors and "$\cdot$" denotes the usual Euclidean inner product.
Now, the 2-point function \eqref{2point} clearly diverges in the coincidence limit $x_1\to x_2$, so it is customary to define the quadratic field operator via a "normal ordering" procedure. This procedure is described in Peskin and Schroeder as "moving all the annihilation operators, $a_\vec{p}$, to the right of the creation operators, $a^\dagger_\vec{p}$" by hand, but this is essentially equivalent to subtracting off the vacuum contribution as,
\begin{equation} :\phi(x_1)\phi(x_2):\equiv \phi(x_1)\phi(x_2)-\Delta_E(x_1-x_2), \tag{2}\label{normal} \end{equation} so that the vacuum expectation of the normal-ordered operators vanish: $\langle 0|:\phi(x_1)\phi(x_2):0\rangle=0$, by definition. Hence, one now has a well-defined quadratic operator,
\begin{align} :\phi^2(x_2):\equiv \lim_{x_1\to x_2} \left\{\phi(x_1)\phi(x_2)-\Delta_E (x_1-x_2)\right\}, \end{align}
in the sense that $\langle0|:\phi^2(x):|0\rangle$ is finite--and, in fact, zero for this choice of normal ordering. (Of course, one could get a different vacuum expectation value via a different choice of normal ordering)
Definition of OPE coefficients: By an "operator product expansion" (OPE), I have in mind a short-distance asymptotic relation (see Scholarpedia article) which holds in any arbitrary (well-behaved) state,
\begin{align} \langle \phi(x_1)\phi(x_2)\rangle \sim \sum_{A} C^{\phi\phi}_{A}(x_1,x_2) \langle \mathcal{O}^A(x_2) \rangle, \end{align}
where $C^{\phi\phi}_A$ are the OPE coefficients of this expansion around $x_1\to x_2$ and the summation runs over all monomials of the fields and its derivatives.
For the free theory, it would appear such an expansion can be obtained via Taylor expanding the terms inside the normal ordering sign of \eqref{normal} as
\begin{align} \phi(x_1)\phi(x_2)= \sum_{\alpha} \frac{(x_1-x_2)^\alpha}{\alpha!} :\partial_\alpha \phi(x_2)\phi(x_2):+\Delta_E (x_1-x_2), \label{expansion} \tag{3} \end{align} where $\alpha$ is a multi-index.
Clearly, the preceding expansion holds trivially in the vacuum state. Because the particle states of the free theory can be expressed in terms of the vacuum, I believe such an expansion will hold for an arbitrary particle state.
Reading off the nontrivial coefficients of the expansion \eqref{expansion}, I find
\begin{align} C^{\phi\phi}_{I}(x_1,x_2)&=\Delta_E (x_1-x_2)\\ C^{\phi\phi}_{\partial_\alpha\phi\phi}(x_1,x_2)&=\frac{(x_1-x_2)^\alpha}{\alpha!}. \tag{4} \\ \end{align}
I expect to obtain more complicated OPE's in the free theory simply from applications of Wick's theorem.
Question: In the "canonical" approach, I seem to obtain the OPE coefficients for the free theory via normal ordering and Wick's theorem. However, it's not clear to me how one defines/obtains the OPE coefficients working instead strictly in the path integral approach or even how to get started in this direction.