A curious issue about Dyson-Schwinger equation (DSE): why does it work so well? This question comes out of my other question Time ordering and time derivative in path integral formalism and operator formalism, especially from the discussion with drake. The original post is somewhat badly composed because it contains too many questions, and till today I finally get energetic enough to compose another question that hopefully clarifies what was asked in that post.
I have no problem with all the textbook derivations of DSE, but after changing a perspective I found it very curious that DSE actually works, I'll take the equation of motion(EOM) of time-ordered Green's function of free Klein-Gordon(KG) field as an example and explain what I actually mean.
The EOM of free KG T-ordered Greens function is
$$(\partial^2+m^2)\langle T\{\phi(x)\phi(x')\}\rangle=-i\delta^4(x-x').\tag{1}$$
The delta function comes from the fact that $\partial^2$contains time derivatives and it doesn't commute with T-ordering symbol. In general for Bosonic operators
$$\partial_t \langle  T\{A(t)\,B(t')\}\rangle=\langle T \{ \dot A(t)B(t')\rangle+\delta (t-t')\,\langle  [A(t),B(t')]\rangle ,\tag{2}$$
$(1)$ can be derived from $(2)$ and the equal time canonical commutation relation of the fields.
However $(1)$ isn't very obvious from path integral approach:
$$(\partial^2+m^2)\langle T\{\phi(x)\phi(x')\}\rangle=(\partial^2+m^2)\int\mathcal{D}\phi e^{iS}\phi(x)\phi(x')\\\quad \quad \quad \qquad \qquad \qquad \qquad =\int\mathcal{D}\phi e^{iS}[(\partial^2+m^2)\phi(x)]\phi(x').\tag{3}$$
Now if we formally and naively think
$$\int\mathcal{D}\phi e^{iS}[(\partial^2+m^2)\phi(x)]\phi(x')=\langle T\{[(\partial^2+m^2)\phi(x)]\phi(x')\}\rangle,\tag{4}$$
with the notation $\langle\cdots\rangle$ always denoting expectation value in operator approach. Then the final result will be 0(due to the field equation) instead of $-i\delta^4(x-x')$, and a result like equation(2) cannot be obtained.
As drake has pointed out, this is because of the ambiguity in the definition of the equal time operator product when time derivative is present in the path integral, it's very important(e.g. CCR in path integral) to define clearly the time derivative on time lattice, that is, the discretization. There are 3 possible definitions of $\dot \phi$(omitting the spatial variables):

*

*(a) forward derivative $$\dot \phi(t)=\frac{\phi(t+\epsilon^+)-\phi(t)}{\epsilon^+};$$


*(b) backward derivative $$\dot \phi(t)=\frac{\phi(t)-\phi(t-\epsilon^+)}{\epsilon^+};$$


*(c) centered derivative $$\dot \phi(t)=\frac{\phi(t+\epsilon^+)-\phi(t-\epsilon^+)}{2\epsilon^+}=\frac{1}{2}(\text{forward}+\text{backward}).$$
These different time discretizations will lead to different equal-time operator orderings(see Ron Maimon's answer in Path integral formulation of quantum mechanics this post ), respectively they are:

*

*(a)$$\int\mathcal{D}\phi e^{iS}\dot \phi(t)\phi(t)=\langle \dot \phi(t)\phi(t)\rangle;$$


*(b)$$\int\mathcal{D}\phi e^{iS}\dot \phi(t)\phi(t)=\langle \phi(t)\dot \phi(t)\rangle;$$


*(c)$$\int\mathcal{D}\phi e^{iS}\dot \phi(t)\phi(t)=\frac{1}{2}[\langle \dot \phi(t)\phi(t)\rangle+\langle \phi(t)\dot \phi(t)\rangle].$$
With these in mind, we can now get equation $(1)$ from path integral(I'll just  show it for the point of equal time, because for $t\neq t'$ there isn't any inconsistency): first I take definition (c) for $\dot \phi$, but define $\ddot\phi$ using forward derivative(which I agree is contrived), so we have
$$\int\mathcal{D}\phi e^{iS}\ddot \phi(t)\phi(t)\equiv\int\mathcal{D}\phi e^{iS}\frac{1}{\epsilon^+}[\dot \phi(t+\epsilon^+)-\dot \phi(t)]\phi(t)
=\frac{1}{\epsilon^+}\{\langle\dot \phi(t+\epsilon^+)\phi(t)\rangle-\frac{1}{2}\langle \dot \phi(t)\phi(t)\rangle-\frac{1}{2}\langle \phi(t)\dot \phi(t)\rangle\}\\
=\frac{1}{\epsilon^+}\{\langle 
\dot \phi(t+\epsilon^+)\phi(t)\rangle-\langle \dot \phi(t)\phi(t)\rangle+\frac{1}{2}\langle [\dot \phi(t),\phi(t)]\rangle\}\\
=\langle \ddot \phi(t)\phi(t)\rangle+\frac{1}{2\epsilon^+ }\langle[\dot \phi(t),\phi(t)]\rangle=\langle \ddot \phi(t)\phi(t)\rangle+\frac{1}{2\epsilon^+ }\delta^3(\mathbf{x}-\mathbf{x'})\tag{5}$$
Now we can formally think $\lim_{\epsilon^+\to 0}\frac{1}{2\epsilon^+}=\delta(0)$ because $\delta (t)= \lim_{\epsilon^+\to 0}\,{1\over 2\epsilon^+}\,e^{-\pi\,t^2/(4\,{\epsilon^+}^2)}$. So $(5)$ becomes
$$\int\mathcal{D}\phi e^{iS}\ddot \phi(t)\phi(t)=\langle \ddot \phi(t)\phi(t)\rangle+\delta(0)\delta^3(\mathbf{x}-\mathbf{x'})\tag{6}$$
The rest is trivial, just apply the spatial derivatives, add it to $(6)$ and apply the field equation, then $(1)$ will be reproduced. The above derivation is mostly due to drake, in a more organized form.
Now it's clear that carefully defining time-derivative discretization is crucial to get the correct result, a wrong discretization won't give us $(1)$. However the derivation of DSE makes absolutely no reference to any discretization scheme, but it always gives a consistent result with the operator approach, why does it work so well?
Many thanks to who are patient enough to read the whole post!
UPDATE: Recently I had a lucky chance to communicate with Professor Dyson about this problem. His opinion is that neither these manipulations nor DSE is true math, because of the lack of mathematical rigor of the underlying theory, so there could be situations where DSE might just fail too, but unfortunately he couldn't immediately provide such an example. Although not very convinced (in the sense that even there's such an example, I still think the "degrees of naivety" of different approaches can be discerned, DSE is clearly more sophisticated and powerful than a direct application of $\partial^2+m^2$), I'd be partially satisfied if someone can provide a situation where DSE fails .
 A: The Schwinger-Dyson equation works so well because it conveniently keeps the ambiguity in time discretization. We have two ways of evaluating the path integral.
(1) The Schwinger-Dyson way:
$$\langle \mathcal T_* \{A(t)B(t')\}\rangle = \int \mathcal D\phi\, e^{i S} A(t)B(t')$$
where $\mathcal T_*$ denotes covariant time ordering, and $\mathcal T$ denotes the usual time ordering.
(2) Time discretized way: When $t=t'$, we must specify whether to use forward/backward/centered derivative for time derivatives.
Indeed, (2) involves making a choice, whereas (1) does not—it remains ambiguous. In particular, (1) will be of form
$$\langle \mathcal T_* \{A(t)B(t')\}\rangle=\langle \mathcal T\{A(t)B(t')\}\rangle + \delta(t-t')(\text{commutator terms})\,.$$
The time ordering term usually written in terms of step functions
$$\langle \mathcal T\{A(t)B(t')\}\rangle=\theta(t-t')A(t)B(t')+\theta(t'-t)B(t')A(t)\tag{A}\,.$$
Now, we will show that (A) is not correct if $A(t)$ or $B(t')$ contains time derivatives. This ambiguity for (A) is solved in the exact same way as in the path integral: by discretizing time by $\epsilon$, applying the operation ($\mathcal T$ or path integral), and taking the $\epsilon\to 0$ limit.
Explicitly, let $A(t)=\dot\phi(t),B(t')=\phi(t)$. Then
(a) Forward derivative
$$\left\langle \mathcal T\left\{\frac{\phi(t+\epsilon)-\phi(t)}{\epsilon}\phi(t)\right\}\right\rangle=\langle \dot \phi(t) \phi(t)\rangle$$
(b) Backward derivative
$$\left\langle \mathcal T\left\{\frac{\phi(t)-\phi(t-\epsilon)}{\epsilon}\phi(t)\right\}\right\rangle=\langle  \phi(t) \dot\phi(t)\rangle$$
(b) Centered derivative
$$\left\langle \mathcal T\left\{\frac{\phi(t+\epsilon)-\phi(t-\epsilon)}{2\epsilon}\phi(t)\right\}\right\rangle=\frac 1 2 (\langle  \dot\phi(t) \phi(t)\rangle+\langle  \phi(t) \dot\phi(t)\rangle)\,.$$
Still, we should note that (A) can be used safely in computations. The ambiguity arises only in a measure $0$ set ($t=t'$), therefore any smooth real world computation will be blind to the choice of the derivative.
A: TL;DR: The main reason why the naive path integral derivation of SD eqs. works well beyond what could be expected is that both concepts employ the same underlying notion of time ordering, namely the covariant time ordering $T_{\rm cov}$, cf. e.g. my Phys.SE answer here.
I) In the main part of this answer we would like to investigate in more detail some formal aspects of the correspondence between
$$ \text{Path integral formalism}\qquad \longleftrightarrow \qquad
\text{Operator formalism},\tag{1} $$
in particular the cohabitation of, on one hand the Schwinger-Dyson (SD) equations, and on another hand, the Heisenberg's equations of motion (eom).
The correspondence (1) is notoriously subtle, see e.g. Ref. 1. and this Phys.SE post. For a general interacting field theory, both sides of the path integral/operator correspondence (1) are typically not rigorously defined, cf. e.g. Ref. 1. Both sides of the correspondence may in principle receive quantum corrections due to operator ordering problems. So it is difficult to come with reliable statements at this formal level.
To simplify the discussion and gain intuition, we are going to make some assumptions.

*

*We are going to consider free (quadratic) theories only. OP's example is covered by this. Free theories have the advantage that we can present explicit formulas.


*It is most economically to argue via the Hamiltonian (as opposed to the Lagrangian) formulation, since we then are only going to have one (as opposed to two) time derivatives. So that's what we are going to do. (Also for simplicity we ignore cases with singular Legendre transformations. Then it is always possible to (Gaussian) integrate out the momenta to get to the corresponding Lagrangian formulation, if that's what one wants.)


*Also we treat field theory formally like point mechanics. All spatial coordinates are suppressed via DeWitt condensed notation. Only the time-variable $t$ is manifestly kept. Thus our variables are
$$ z^I, \qquad  I~\in~ \{1, \ldots, 2N\},\tag{2} $$
and they depend on time $t$, where $N$ could be infinite.


*Also to avoid annoying sign factors, we restrict attention to theories with only Grassmann-even variables.


*Furthermore, we assume for simplicity (or via Darboux's Theorem) that the equal-time Poisson bracket is constant and $z$-independent
$$\begin{align}\omega^{IK}
~=~& \{ z^I(t),z^K(t) \}_{PB}, \cr  
I,K~\in~&\{1, \ldots, 2N\}. \end{align}\tag{3}  $$
Much of the above assumptions can be relaxed, but we will not discuss that here.
II) Classically, with the above assumptions, the Hamiltonian action reads
$$ S_H[z;J]~=~ \int dt\ L_H(z(t);\dot{z}(t);J(t)),\tag{4} $$
where the Hamiltonian Lagrangian is
$$ L_H(z;\dot{z};J)~=~\frac{1}{2} z^I\omega_{IK} \dot{z}^K - H(z,J),\tag{5} $$
and the Hamiltonian is quadratic
$$ \begin{align} H(z,J)~=~&H_0(z) -J_I z^I, \cr 
H_0(z)~:=~& \frac{1}{2} z^I h_{IK} z^K. \end{align}\tag{6} $$
The Euler-Lagrange derivatives (with an index raised by the symplectic metric $\omega^{IL}$) correspond to Hamilton's eom
$$\begin{align} 0~\approx~&\omega^{IK}\frac{\delta S_H[z;J]}{\delta z^K(t)}\cr
~=~&\dot{z}^I(t) - \{z^I(t),H(z(t),J(t))\}_{PB} \cr ~=~&\left(\delta^I_K\frac{d}{dt}- h^I{}_K\right)z^K(t) + \omega^{IK} J_K(t). \end{align}\tag{7} $$
(Here the $\approx$ sign means equality modulo classical equations of motion.) Also we have defined the matrix
$$ h^I{}_L~:=~\omega^{IK} h_{KL}. \tag{8} $$
In the corresponding operator formalism, the Hamilton's eom (7) turns into Heisenberg's eom, which is an operator identity.
III) The Hessian reads
$$ \begin{align} {\cal H}_{IK}(t,t^{\prime})~:=~&\frac{\delta^2 S_H[z;J]}{\delta z^I(t)\delta z^K(t^{\prime})} \cr
~=~&\omega_{IK}\delta^{\prime}(t-t^{\prime})
-h_{IK}\delta(t-t^{\prime})\cr
~=~&\omega_{IL}\left(\delta^L_K\frac{d}{dt}- h^L{}_K\right)\delta(t-t^{\prime}).\end{align}\tag{9} $$
The corresponding Green's function $G={\cal H}^{-1}$ is the inverse of the Hessian (9):
$$  G^{IK}(t,t^{\prime})~=~\frac{1}{2} {\rm sgn}(t-t^{\prime})\left(e^{(t-t^{\prime})h}\right)^I{}_L ~\omega^{LK}, \tag{10} $$
$$ \left(\delta^I_L\frac{d}{dt}- h^I{}_L\right)G^{LK}(t,t^{\prime})~=~\omega^{IK}\delta(t-t^{\prime}). \tag{11} $$
IV) We define the partition function
$$\begin{align} 
&Z[J]\cr
~:=~ &\int\![dz] e^{\frac{i}{\hbar}S_H[z;J]} \cr
~=~& e^{\frac{i}{\hbar}W_c[J]} \cr
~=~& \int\![dz] \exp\left[\frac{i}{2\hbar}\iint \! dt~dt^{\prime}~z^I(t) ~{\cal H}_{IK}(t,t^{\prime})~ z^K(t^{\prime}) 
 + \frac{i}{\hbar}\int \! dt~J_I(t) z^I(t)\right] \cr
~=~& Z[0]\exp\left[-\frac{i}{2\hbar}
\iint \! dt~dt^{\prime}~J_I(t) ~G^{IK}(t,t^{\prime})~ J_K(t^{\prime})\right] ,\end{align}\tag{12} $$
and the generator of connected Feynman diagrams
$$\begin{align} W_c[J]~:=~&\frac{\hbar}{i}\ln Z[J]\cr 
~=~& W_c[0]- \frac{1}{2} \iint \! dt~dt^{\prime}~J_I(t) ~G^{IK}(t,t^{\prime})~ J_K(t^{\prime}) \end{align} ,\tag{13} $$
whose explicit form follows from Gaussian integration. The quantum-average/expectation-value in the Heisenberg picture is defined as
$$\begin{align}\left< T_{\rm cov}\{F[\widehat{z}]\} \right>_J
~=~& \frac{1}{Z[J]}\int\![dz]~F[z]~e^{\frac{i}{\hbar}S_H[z;J]}\cr
~=~&\frac{1}{Z[J]} F\left[ \frac{\hbar}{i} \frac{\delta}{\delta J} \right] Z[J] .\end{align}\tag{14} $$
Here $T_{\rm cov}$ is covariant time-ordering, i.e. time-differentiations inside its argument should be taken after/outside the usual time ordering $T$. Here time-ordering $T$ is defined as
$$\begin{align}T\left\{\widehat{A}(t)\widehat{B}(t^{\prime})\right\}~:=~& \theta(t-t^{\prime})\widehat{A}(t)\widehat{B}(t^{\prime})\cr
&+\theta(t^{\prime}-t)\widehat{B}(t^{\prime})\widehat{A}(t). \end{align}\tag{15} $$
It is crucial that time-differentiation and time-ordering do not commute:
$$\begin{align}T_{\rm cov}&\left\{\frac{d\widehat{A}(t)}{dt}\widehat{B}(t^{\prime})\right\} \cr
~:=~&\frac{d}{dt}T\left\{\widehat{A}(t)\widehat{B}(t^{\prime})\right\}\cr 
~=~& T\left\{\frac{d\widehat{A}(t)}{dt}\widehat{B}(t^{\prime})\right\}+\delta(t-t^{\prime})[\widehat{A}(t),\widehat{B}(t^{\prime})], \end{align}\tag{16}$$
but instead give rise to equal-time commutator (contact) terms.
It is important to realize that the time derivatives inside the Boltzmann factor $e^{\frac{i}{\hbar}S_H[z;J]}$ in the path integral should respect the underlying time slicing procedure. See e.g. this and this Phys.SE answer. This induces the covariant time-ordering $T_{\rm cov}$ in eq. (14).
The 1-point function reads
$$\begin{align} \left< \widehat{z}^I(t)\right>_J
~=~& \frac{\delta W_c[J]}{\delta J_I(t)}\cr
~=~& -\int \!dt^{\prime} ~G^{IK}(t,t^{\prime})~J_K(t^{\prime}), \end{align}\tag{17} $$
while the time-ordered 2-point function reads
$$\begin{align} &\left< T_{\rm cov}\{ \widehat{z}^I(t)\widehat{z}^K(t^{\prime}) 
\}\right>_J\cr
~=~& -\frac{\hbar^2}{Z[J]}
\frac{\delta^2 Z[J]}{\delta J_I(t)\delta J_K(t^{\prime})}\cr
~=~&\frac{\delta W_c[J]}{\delta J_I(t)}\frac{\delta W_c[J]}{\delta J_K(t^{\prime})}+\frac{\hbar}{i}\frac{\delta^2 W_c[J]}{\delta J_I(t)\delta J_K(t^{\prime})}\cr
 ~=~&\left< \widehat{z}^I(t)\right>_J 
\left< \widehat{z}^K(t^{\prime})\right>_J
+i\hbar G^{IK}(t,t^{\prime}). \end{align}\tag{18}$$
The corresponding 2-point function without time-ordering reads:
$$\begin{align} \left< \widehat{z}^I(t)\widehat{z}^K(t^{\prime}) \right>_J 
~=~&\left< \widehat{z}^I(t)\right>_J 
\left< \widehat{z}^K(t^{\prime})\right>_J\cr
&+\frac{i\hbar}{2} \left(e^{(t-t^{\prime})h}\right)^I{}_L ~\omega^{LK}.\end{align} \tag{19} $$
V) The Schwinger-Dyson (SD) equations read
$$ i\hbar\left< T_{\rm cov}\left\{\frac{\delta F[\widehat{z}]}{\delta z^I(t)} \right\}\right>_J 
~=~\left< T_{\rm cov}\left\{ F[\widehat{z}]\frac{\delta S_H[\widehat{z};J]}{\delta z^I(t)}\right\}\right>_J.\tag{20} $$
The SD eqs. (20) are here formally written in the operator language, but their justification is most easily argued via the path integral formalism, cf. e.g. this Phys.SE post. The SD eqs. (20) simply reflect the fact that a path integral of a total derivative vanishes if the boundary contributions are zero
$$  0~=~\int [dz]\frac{\delta}{\delta z^I(t)}\left\{F[z] e^{\frac{i}{\hbar}S_H[z;J]}\right\}, \tag{21} $$
cf. this Phys.SE post.
VI) Naively the rhs. of the SD-eqs. (20) is proportional to the Heisenberg eom (7), which is an operator expression that vanishes identically, so why then there is a non-zero quantum correction on the lhs. of the SD-eqs. (20)? The resolution to this apparent paradox is hidden in the fact that time-differentiation and time-ordering do not commute, cf. eq. (16). To see how this works, pick for simplicity the functional $F[z]=z^{K}(t^{\prime})$ of a single variable. (Then it is enough to use the Poisson bracket rather than the Moyal-Groenewold $\star$-product.)  After raising an index by the symplectic metric $\omega^{IL}$, the SD-eqs. (20) become
$$\begin{align} i\hbar\omega^{IK} &\delta(t-t^{\prime}) \cr
~\stackrel{(20)}{=}~&\left< T_{\rm cov}\left\{\omega^{IL}\frac{\delta S_H[\widehat{z};J]}{\delta z^L(t)}\widehat{z}^K(t^{\prime})\right\}  \right>_J\cr 
~\stackrel{(7)}{=}~&\left< \frac{d}{dt}T\left\{\widehat{z}^I(t)
\widehat{z}^K(t^{\prime})\right\}\right>_J \cr
&- \left< T\left\{\{\widehat{z}^I(t),H(\widehat{z}(t),J(t))\}_{PB} \widehat{z}^K(t^{\prime})\right\}\right>_J\cr
~=~&\left(\delta^I_L\frac{d}{dt}- h^I{}_L\right)\left< T\left\{\widehat{z}^L(t)
\widehat{z}^K(t^{\prime})\right\}\right>_J \cr
&+\omega^{IL} J_L(t)\left< \widehat{z}^K(t^{\prime})\right>_J  \cr
 ~\stackrel{(16)}{=}~& i\hbar\omega^{IK} \delta(t-t^{\prime}) \cr
&+\left< T\left\{ \underbrace{\left(\frac{d\widehat{z}^I(t)}{dt} - \{\widehat{z}^I(t),H(\widehat{z}(t),J(t))\}_{PB} \right)}_{=0}\widehat{z}^K(t^{\prime})\right\}\right>_J\cr
~\stackrel{(7)}{=}~&i\hbar\omega^{IK} \delta(t-t^{\prime}).\end{align}\tag{22} $$
Eq. (22) shows how SD-eqs. (20) and Heisenberg's eom (7) can co-exist.
References:

*

*F. Bastianelli and P. van Nieuwenhuizen, Path Integrals and Anomalies in Curved Space, 2006.

A: Just want to reenforce an important point: Path integral alone can NOT take care of the time ordering in operator formalism.
One has to supplement path integral with Feynman's $i\epsilon$ prescription (or a specific discretization scheme as mentioned in the OP) in the propagators to recover the proper time ordering in operator formalism. Feynman's $i\epsilon$ prescription is an add-on to the path integral formalism, it's not something automatically implied by the path integral.
The propagators are not well defined in the path integral formalism, due to the ambiguity of the position of the propagator poles in the complex plane. One has to specify on which side away from the real axis the propagator poles are located. Feynman's $i\epsilon$ prescription is just one option out of many. Other alternative prescriptions would imply a different propagator than the one in time ordered operator formalism.
One advantage of path integral is that one can use it to derive Dyson-Schwinger equation (DSE) rather easily, as opposed to the painstaking way in operator formalism. The magic is that the pole prescription (or some discretization scheme) could be totally ignored (or implicit) in the derivation of the DSE in path integral formalism. I think this is the part that the OP is curious about. This is more of a math question. But hey, path integral itself is not well defined mathematically after all.
