How to calculate the correlation function like $\langle \partial_i \phi(x) \partial_j\phi(0) \rangle$ by Gaussian path integral? From the standard text book about quantum field theory, we know that if we consider $$\mathcal{L}=\frac{1}{2}(\partial_{\mu} \phi)^2-\frac{m^2}{2}\phi^2,$$ the partition function of this Gaussian theory can be worked out
$$\frac{Z[J]}{Z[0]}=\frac{\int D\phi\,e^{\frac{i}{2}\int d^4x\,\phi[-(\partial_\mu^2+m^2)]\phi+i\int\,d^4x\,J\phi}}{\int D\phi\,e^{\frac{i}{2}\int d^4x\mathcal{L}}}=e^{-\frac{1}{2}\int d^4xd^4y J(x)G(x-y)J(y)},$$
where we use integration by parts $\int d^4x(\partial_\mu\phi)^2=-\int d^4x\phi\partial^2_{\mu}\phi$, and $iG(x,y)=(\partial_\mu^2+m^2)^{-1}$. On the other hand,
$$\langle T\phi(x)\phi(y)\rangle=-\frac{1}{Z[0]}\frac{\delta^2Z[J]}{\delta J(x)\delta J(y)}|_{J=0}=G(x-y).$$
But how about this kind of correlation function $\langle \partial_i \phi(x) \partial_j\phi(0) \rangle$ ? May be we can use some similar strategy,
$$\frac{Z'[J]}{Z[0]}=\frac{\int D\phi\,e^{\frac{i}{2}\int d^4x\,\phi[-(\partial_\mu^2+m^2)]\phi+i\int\,d^4x\,J^\mu\partial_\mu\phi}}{\int D\phi\,e^{\frac{i}{2}\int d^4x\mathcal{L}}}=e^{-\frac{1}{2}\int d^4xd^4y (\partial_\mu J^\mu(x))G(x-y)(\partial_\nu J^\nu(y))},$$
and $$\langle T\partial_i\phi(x)\partial_j\phi(y)\rangle=-\frac{1}{Z[0]}\frac{\delta^2Z'[J]}{\delta J^i(x)\delta J^j(y)}|_{J=0},$$ then we will meet some terms like $$\int d^4x_1d^4y_1\frac{\delta \partial_\mu J^\mu(x_1)}{\delta J^i(x)}G(x_1-y_1)(\partial_\nu J^\nu(y)),$$ can we exchange $\delta$ and $\partial$ here? Then finally we may find $$\langle \partial_i \phi(x) \partial_j\phi(0) \rangle\sim \partial_i\partial_j\,(\partial_\mu^2+m^2)^{-1},$$ and after Fourier transformation, we may have $\frac{k_ik_j}{\omega^2-k^2+m^2}$.
Actually, I find this correlation in some book like Altland and Simons' book, Condensed matter field theory, 2nd, page 393, exercise and Xiao-Gang Wen's book, quantum field theory of many body systems, problem 3.3.3. I want to know if there is any standard method to calculate this kind of correlation, since it seems not rigorous in above "derivation".
 A: Your idea of introducing the generating functional
$$ Z[J] := \int \mathrm{D}\phi \exp\!\left(\mathrm{i}\,S[\phi]+\mathrm{i}\int\mathrm{d}^d{x} \ J^\mu\partial_\mu\phi\right)$$
is correct. Indeed you see that1
$$-\frac{\delta^2\ Z[J]}{\delta J^\mu(x)\ \delta J^\nu(y)}\Bigg|_{J=0} = \Big\langle\partial_\mu\phi(x)\partial_\nu\phi(y)\Big\rangle.$$
In your Gaussian theory, you were worried about exchanging $\partial_\mu$ and $\delta$. You don't have to do that though. Note that in this case an integration by parts gives
\begin{align} Z[J] &= \exp\!\left(-\frac12\int\mathrm{d}^d x\,\mathrm{d}^d y\ \partial_\mu J^\mu(x)\;G(x-y)\;\partial_\nu J^\nu(y)\right)=\\
&=\exp\!\left(\frac12\int\mathrm{d}^d x\,\mathrm{d}^d y\ J^\mu(x)\;\partial_\mu \partial_\nu G(x-y)\; J^\nu(y)\right),
\end{align}
where both derivatives are wrt $x$ in the last line, since $\partial_y f(x-y) = - \partial_x f(x-y)$. From here on you are free to take your $J(x)$ variational derivatives as you usually would and find indeed
$$\Big\langle\partial_\mu\phi(x)\partial_\nu\phi(y)\Big\rangle = - \partial_\mu \partial_\nu\ G(x-y).$$
Finally, @QuantumEyedea's comment about first computing $\big\langle\phi(x)\phi(y)\big\rangle$ and then taking derivatives is absolutely correct. It does come with limitations, however. If you want to compute a correlation function of polynomials of the field and its derivatives:
$$ \Big\langle \mathcal{O}_1(x)\cdots \mathcal{O}_n(x)\Big\rangle,$$
where $\mathcal{O}_i(x) := \mathrm{P}_i[\phi,\partial\phi,\partial^2\phi,\cdots](x)$ are polynomials, with @QuantumEyedea's method you would have to compute each $m$-point function of $\phi$ and then take appropriate derivatives, which will be chaotic and daunting. With the generating functional method, the computation is much more organised and straightforward. This is exactly what happens in practice when people compute string theory amplitudes.

1 setting $Z[0]=1$ for simplicity
