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".