Specific commutator calculation using the path integral Consider the path integral quantisation of a scalar field $\phi$ on flat spacetime. Let the Lagrangian be $\mathcal{L}$. I would like to prove the following equal time commutation relation:
\begin{equation}
\tag{1}
\left[\phi(\mathbf{x},t),\frac{\partial{\mathcal{L}}}{\partial (\partial_0\phi)}(\mathbf{y},t)\right]=i\hbar\delta^{(3)}(\mathbf{x}-\mathbf{y}).
\end{equation}
The motivation is hopefully clear: I am aiming to show that the path integral is able to reproduce the CCRs that we impose in canonical quantisation.
So far, I've tried writing commutator matrix elements as follows:
\begin{equation}
\tag{2}
\langle \phi_2|\left[\phi(\mathbf{x},t),\frac{\partial{\mathcal{L}}}{\partial (\partial_0\phi)}(\mathbf{y},t)\right]|\phi_1\rangle =\int D\phi \left(\phi(\mathbf{x},t+\epsilon)\frac{\partial{\mathcal{L}}}{\partial (\partial_0\phi)}(\mathbf{y},t)- \phi(\mathbf{x},t)\frac{\partial{\mathcal{L}}}{\partial (\partial_0\phi)}(\mathbf{y},t+\epsilon)\right)e^{iS[\phi]/\hbar},
\end{equation}
where $\epsilon> 0$ is small, and the path integral has some appropriate domain ($\mathbb{R}^3\times[0,\epsilon]$ is my guess) and BCs. Note that the introduction of $\epsilon$ ensures the path integral time-orders the products as desired. However I've been unable to proceed from here.
Am I on the right track, or is there a better way to prove (1)?
[Please note, before you tell me that (1) is simply the CCR, which is an assumption/axiom in QFT: that is indeed true in the canonical framework, but the point of this question is to take the path integral as our starting point. Given the path integral tells us everything about the QFT in principle, it should be able to derive (1).]
 A: Consider the Schwinger-Dyson equation
$$\langle\left[\partial_\mu\frac{\partial\mathcal{L}}{\partial\left(\partial_\mu\phi\right)}- \frac{\partial\mathcal{L}}{\partial\phi} \right](x) \,\,\phi(y)\rangle = -i\delta^{(4)}(x-y)$$
This is just the equations of motion inside a correlation function and the presence of the $\phi(y)$ operator leads to the contact term on the right side (I'm setting $\hbar=1$). It can be derived in the path integral approach by considering $\langle\phi(y)\rangle$ and doing a local change of variables $\phi\rightarrow \phi+\delta\phi$, see e.g. Polchinski's String Theory Vol 1 Sec 2.3.
Now integrate this equation over a cylinder where the spatial coordinates $x^i$ are integrated over a region $A$, and the time coordinate $x^0$ is integrated over a tiny interval from $y^0-\epsilon/2$ to $y^0+\epsilon/2$. If we have an ordinary Lagrangian which doesn't involve terms with products of both $\phi$ and $\dot{\phi}$ we would expect the `potential' term involving $\partial \mathcal{L}/\partial \phi$ to involve no extra divergences and after integrating over $x$ this term will vanish in the limit $\epsilon\rightarrow 0$ and we will drop it.
The other `kinetic' term is a total divergence so we may rewrite the integral over $x$ as a flux through the surface of the cylinder. The flux of $\partial \mathcal{L}/\partial \left(\partial_i\phi\right)$ through the sides of the cylinder should likewise vanish in the limit $\epsilon\rightarrow 0$. So we are left with
$$\lim_{\epsilon\rightarrow 0}\int_A d^3 x\left[ \langle\left.\frac{\partial\mathcal{L}}{\partial\left(\partial_0\phi\right)}(x)\,\,\phi(y)\rangle\right|_{x^0=y^0+\epsilon/2} - \langle\phi(y)\left.\frac{\partial\mathcal{L}}{\partial\left(\partial_0\phi\right)}(x) \rangle\right|_{x^0=y^0-\epsilon/2} \right] = -i\,1_A(y^i)$$
where $1_A$ is an indicator function that equals one if $y^i\in A$, and zero otherwise. So the quantity in brackets is clearly $-i\delta^{(3)}(x^i-y^i)$. And also from the relation of the path integral correlation functions to time ordered correlation functions in the operator approach it must be the commutator
$$\left[\frac{\partial\mathcal{L}}{\partial\left(\partial_0\phi\right)}(x),\phi(y)\right]=-i\delta^{(3)}(x^i-y^i)$$
