Time ordering and time derivative in path integral formalism and operator formalism In operator formalism, for example a 2-point time-ordered Green's function is defined as
$$\langle\mathcal{T}\phi(x_1)\phi(x_2)\rangle_{op}=\theta(x_1-x_2)\phi(x_1)\phi(x_2)+\theta(x_2-x_1)\phi(x_2)\phi(x_1),$$
where the subscript "op" and "pi" refer to operator and path integral formalism, respectively. Now if one is to take a time derivative of it, the result will be
$$\frac{\partial}{\partial x_1^0}\langle\mathcal{T}\phi(x_1)\phi(x_2)\rangle_{op}=\langle\mathcal{T}{\frac{\partial \phi(x_1)}{\partial x_1^0}}\phi(x_2)\rangle_{op}+\delta(x_1^0-x_2^0)[\phi(x_1),\phi(x_2)],$$
the delta function comes from differentiating the theta functions. This means time derivative does not commute with time ordering.
If we consider path integral formalism, the time-ordered Green's function is defined as
$$\langle\mathcal{T}\phi(x_1)\phi(x_2)\rangle_{pi}=\int\mathcal{D}\phi\phi(x_1)\phi(x_2)e^{iS(\phi)}.$$
Of course
$$\langle\mathcal{T}\phi(x_1)\phi(x_2)\rangle_{op}=\langle\mathcal{T}\phi(x_1)\phi(x_2)\rangle_{pi},$$
as is proved in any QFT textbook. However in path integral case time derivative commutes with time ordering, because we don't have anything like a theta function thus
$$\frac{\partial}{\partial x_1^0}\int\mathcal{D}\phi\phi(x_1)\phi(x_2)e^{iS(\phi)}=\int\mathcal{D}\phi\frac{\partial}{\partial x_1^0}\phi(x_1)\phi(x_2)e^{iS(\phi)}.$$
I did a bit googling and found out that for the path integral case the time-ordered product is called "$\mathcal{T^*}$ product" and operator case just "$\mathcal{T}$ product".
I am not that interested in what is causing the difference(still explanations on this are welcomed), because I can already vaguely see it's due to some sort of ambiguity in defining the product of fields at equal time. The question that interests me is, which is the right one to use when calculating Feynman diagrams?
I did find a case where both give the same result, i.e. scalar QED (c.f. Itzykson & Zuber, section 6-1-4 and this post), but is it always the case? If these two formulations are not effectively equivalent, then it seems every time we write down something like $\langle\partial_0\phi\cdots\rangle$, we have to specify whether it's in the sense of the path integral definition or operator definition.
EDIT: As much as I enjoy user1504's answer, after thinking and reading a bit more I don't think analytic continuation is all the mystery. In Peskin & Schroeder chap 9.6 they manage to use path integral to get a result equivalent to operator approach, without any reference to analytic continuation. It goes like this: Consider a $T$-product for free KG field
$$\langle T\{\phi(x)\phi(x_1)\}\rangle=\int\mathcal{D}\phi\phi(x)\phi(x_1)e^{iS(\phi)}.$$ Apply Dyson-Schwinger equation, we get
$$\int\mathcal{D}\phi(\partial^2+m^2)\phi(x)\phi(x_1)e^{iS}=-i\delta^4(x-x_1),$$
then they just assume the $\partial^2$ commute with path integration (which is already weird according to our discussion) and they conclude
$$(\partial^2+m^2)\int\mathcal{D}\phi\phi(x)\phi(x_1)e^{iS}=(\partial^2+m^2)\langle T\{\phi(x)\phi(x_1)\}\rangle=-i\delta^4(x-x_1).$$
This is just the right result given by operator approach, in which $\delta(x^0-x_1^0)$ comes from $\theta$ function. Given my limited knowledge on the issue, this consistency looks almost a miracle to me. What is so wicked behind these maths?
Response to  @drake: If $a$ is a positive infinitesimal, then $$\int \dot A(t) B(t) \,e^{iS}\equiv\int D\phi\, {A(t+a)-A(t)\over a}B(t)\,e^{iS}=\frac{1}{a}\langle T\{A(t+a)B(t)\}\rangle-\frac{1}{a}\langle A(t)B(t)\rangle,$$
notice the second term has an ordering ambiguity from path integral (say $A=\dot{\phi},B=\phi$), and we can make it in any order we want by choosing an appropriate time discretization, c.f. Ron Maimon's post cited by drake. Keeping this in mind we proceed:
$$\frac{1}{a}\langle T\{A(t+a)B(t)\}\rangle-\frac{1}{a}\langle A(t)B(t)\rangle\\=\frac{1}{a}\theta(a)\langle A(t+a)B(t)\rangle+\frac{1}{a}\theta(-a)\langle B(t)A(t+a)\rangle-\frac{1}{a}\langle A(t)B(t)\rangle\\=\frac{1}{a}\theta(a)\langle A(t+a)B(t)\rangle+\frac{1}{a}[1-\theta(a)]\langle B(t)A(t+a)\rangle-\frac{1}{a}\langle A(t)B(t)\rangle\\=\frac{\theta(a)}{a}\langle [A(t+a),B(t)]\rangle+\frac{1}{a}[\langle B(t)A(t+a)\rangle-\langle A(t)B(t)\rangle].$$
Now taking advantage of ordering ambiguity of the last term to make it $\langle B(t)A(t)\rangle$ (this amounts to defining $A$ using backward discretization, say $A=\dot{\phi}(t)=\frac{\phi(t+\epsilon^-)-\phi(t)}{\epsilon^-}$), then the finally:
$$\frac{\theta(a)}{a}\langle [A(t+a),B(t)]\rangle+\frac{1}{a}\langle B(t)[A(t+a)-A(t)\rangle]\to \frac{1}{2a}\langle [A(t),B(t)]\rangle+\langle B(t)\dot{A}(t)\rangle.$$
(Here again a very dubious step, to get $\frac{1}{2a}$ we need to assume $\theta(a\to 0^+)=\theta(0)=\frac{1}{2}$, but this is really not true because $\theta$ is discontinuous.)
However on the other hand, since $a$ was defined to be a positive infinitesimal, at the very beginning we could've written
$$\frac{1}{a}\langle T\{A(t+a)B(t)\}\rangle-\frac{1}{a}\langle A(t)B(t)\rangle=\frac{1}{a}\langle A(t+a)B(t)\rangle-\frac{1}{a}\langle A(t)B(t)\rangle,$$
then all the above derivation doesn't work. I'm sure there are more paradoxes if we keep doing these manipulations.
 A: EDIT:  I'm leaving this up as background reading to @drake's answer.  (The point of the following is that the path integral does indeed give the correct time ordering, so it is producing the correct $\theta$-function weighted, time-ordered sums, which must be accounted for when differentiating its output.)
The two formalisms are equivalent; if they don't give the same result, something is wrong in the calculation.  To see this you have to understand a subtlety which is not usually well-explained in textbooks, namely that the path integral is not defined merely by taking the limit of a bunch of integrals of the form $\int_{\mbox{lattice fields}} e^{iS(\phi)} d\phi$.
The problem is that these finite-dimensional integrals are not absolutely convergent, because $|e^{iS(\phi)}| = 1$.  To define even the lattice path integral in Minkowski signature, you have to specify some additional information, to say exactly what is meant by the integral.  
In QFT, the additional information you want is that the path integral should be calculating the kernel of the time evolution operator $e^{iH\delta t}$, which is an analytic function of $\delta t$.  This fact is usually expressed by saying that the Minkowski signature path integral is the analytic continuation of a Euclidean signature path integral:   The Euclidean $n$-point functions $E(y_1,...,y_n)$ defined by
$E(y_1,...,y_n) = \int \phi(y_1)...\phi(y_n) e^{-S_E(\phi)} d\phi$
are analytic functions of the Euclidean points $y_i \in \mathbb{R}^d$.  This function $E$ can be continued to a function $A(z_1,...,z_n)$ of $n$ complex variables $z_i \in \mathbb{C}^d$.   This analytic function $A$ does not extend to the entire plane; it has singularities, and several different branches.  Each branch corresponds to a different choice of time-ordering.  One branch is the correct choice, another choice is the 'wrong sign' time-ordering.  Other choices have wrong signs on only some subsets of the points.  If you restrict $A$ to the set $B$ of boundary points of the correct branch, you'll get the Minkowski-signature $n$-point functions $A|_B =  M$, where $M(x_1,...,x_n) = \langle \hat{\phi}(x_1)...\hat{\phi}(x_n)\rangle_{op}$ and the $x_i$ are points in Minkowski space.
In perturbation theory, most of this detail is hidden, and the only thing you need to remember is that the $+i\epsilon$ prescription selects out the correct time-ordering.
A: Good question! It made me think. Here it's the answer: 


*

*Strictly speaking, it is not possible to compute $\theta(t-t')\langle \dot\phi(t)\phi(t')\rangle+\theta (t'-t)\langle \phi (t')\dot\phi (t)\rangle$ ( shorthand notation: $\langle\,\equiv \langle 0 |\,$. Also note that I'm omitting the spatial arguments of the fields ) using the Lagrangian version of the path integral, because to derive this, we are to assume that the insertions (factors multiplying $e^{iS}$ in the integrand of the path integral) are functionals  of the fields at a given time (i.e., they are independent of momenta $\pi$). Thus,


$$\partial_t \,\left[\theta(t-t')\langle \phi(t)\phi(t')\rangle+\theta (t'-t)\langle \phi (t')\phi (t)\rangle\right]=\lim_{a\to 0}{1\over a}\left[\theta(t+a-t')\langle \phi(t+a)\phi(t')\rangle+\theta (t'-t-a)\langle \phi (t')\phi (t+a)\rangle - \theta(t-t')\langle \phi(t)\phi(t')\rangle+\theta (t'-t)\langle \phi (t')\phi (t)\rangle\right]=\lim_{a\to 0}{1\over a}\int D\phi\, (\phi(t+a)\phi(t')-\phi(t)\phi(t'))\,e^{iS}=\int D\phi\, \dot\phi(t)\phi(t')\,e^{iS}$$
which agrees with  $\theta(t-t')\langle \dot\phi(t)\phi(t')\rangle+\theta (t'-t)\langle \phi (t')\dot\phi (t)\rangle$ only if $t\neq t'$.$^1$ 


*

*Example. Let $A(t)$ and $B(t)$ two functional at a given time, then you can check that


$$\int \dot A(t) B(t) \,e^{iS}=\langle\dot A(t)B(t)\rangle+\lim_{a\to 0}{1\over 2a}\langle [A(t),B(t)]\rangle$$
The last term represents the Dirac delta you found after deriving the step function.


*

*While your question is very interesting, I think that your example is unfortunate since $\delta (t-t')[\phi(x),\phi(x')]=\delta (t-t')[\phi(t,\vec x),\phi(t,\vec x')]=0$. Let me choose an example in which this last commutator is different from zero to show how the time derivative of a correlation function splits in the two terms you mention: $\partial_t \langle \, T\, \phi(t)\,A(t')\,\rangle$, where $A$ is a functional of fields and momenta at a given time. Let's make a general variation (which doesn't modified the path integral measure) of the Hamiltonian or phase-space path integral ($S=\int dt\, \dot \phi\,\pi-H \, $ ). Since the momentum is an integration variable, the integral may not change:


$$0={\delta\over \delta \pi (t)}\int D\phi D\pi \, A(t')\,e^{iS}\,\delta \pi=\\
\int D\phi D\pi\, \left( {\delta A(t')\over \delta \pi (t)}+A(t')i\dot\phi (t)-A(t')(-i){\delta H (t)\over \delta \pi (t)}\right)\,e^{iS}\,\delta \pi$$
As  $\delta \pi$ is a general variation and:
$${\delta A(t')\over \delta \pi (t)}=\delta(t-t')(-i)\,[\phi(t),A(t')]$$
$${\delta H (t)\over \delta \pi (t)}=-i\,[\phi(t), H]$$
we obtain:
$$\partial_t \langle \, T\, \phi(t)\,A(t')\,\rangle=i\langle \, T\, [\phi(t),H]\,A(t')\,\rangle+\delta (t-t')\,\langle \, [\phi(t),A(t')]\,\rangle \\ =\theta(t-t')\langle \, \dot\phi(t)\,A(t')\,\rangle+\theta(t'-t)\langle \, A(t')\,\dot\phi(t)\,\rangle+\delta (t-t')\,\langle \, [\phi(t),A(t')]\,\rangle$$
If, for example, $A(t')=\pi (t',\vec x')$, the last term gives $i\,\delta ^4 (x-x')$.
$\\$
Just in case it is not clear enough, let me remark that derivatives do commute with the path integral measure. The key point is that 
$$\partial_t \,\left[\theta(t-t')\langle \phi(t)\phi(t')\rangle+\theta (t'-t)\langle \phi (t')\phi (t)\rangle\right]=\int D\phi\, {\phi(t+\epsilon^+)-\phi (t)\over \epsilon ^+}\phi(t')\,e^{iS}\neq\theta(t-t')\langle \dot\phi(t)\phi(t')\rangle+\theta (t'-t)\langle \phi (t')\dot\phi (t)\rangle$$
In addition, I would like to emphasize that ${\delta H (t)\over \delta \pi (t)}$ is a functional at a given time, while $\dot\phi (t)$ in the integrand of a path integral must be interpreted as ${\phi(t+\epsilon^+)-\phi (t)\over \epsilon ^+}$, that is, as a difference between fields evaluated at different times.

$^1$ Note that the derivative can be defined in different ways giving rise to different operator orderings, see Maimon's excellent answer Path integral formulation of quantum mechanics, be careful with some typo in the expressions: where says $x(t)p(t)$ should say $p(t)x(t)$.

EDIT: To derive some of the results above, one needs to take $\theta (0)=1/2$. However, one can proceed in a slightly different manner to avoid such choice (whi, in my opinion, is totally right). For example,
$$\int \dot A(t) B(t) \,e^{iS}=\int {A(t+a)-A(t)\over a } B(t+a/2) \,e^{iS}=\langle\dot A(t)B(t)\rangle+\lim_{a\to 0}{1\over a}\langle [A(t),B(t)]\rangle$$
