Schwinger-Dyson equations in the Coulomb gauge Introduction
As far as I know (and please correct me if something is wrong!) the usual narrative to deal with perturbation theory in QED with the Coulomb gauge goes as follows:
First, the gauge field is fixed to
$$
\nabla\cdot\vec A(x)=0\tag{1}
$$
so that the equations of motion read
$$
\begin{aligned}
\nabla^2 A^0(x)&=J^0(x)\\
\partial^2\vec A(x)&=P(\nabla)\vec J(x)
\end{aligned} \tag{2}
$$
where $\partial^2=\partial_0^2-\nabla^2$ and
$$P^{ij}(\vec k)=\delta^{ij}-\frac{k^ik^j}{\vec k^2}\tag{3}$$
is the projector into the transverse part of the current.
Now, the equation of motion for $A^0$ is not dynamical, and so we integrate that field out by setting
$$
A^0=\frac{1}{\nabla^2}J^0 \tag{4}
$$
into the Lagrangian, generating the well-known non-local Coulomb term.
On the other hand, the true field is $\vec A$, and its propagator is
$$
\Delta^{ij}(k)=\frac{P^{ij}(\vec k)}{k^2+i\epsilon} \tag{5}
$$
The important point is the following: the non-local term is, for all practical purposes, equivalent to an extended propagator $\Delta^{\mu\nu}$ such that it coincides with $\Delta^{ij}$ on the spatial components, and
$$
\Delta^{00}(\vec k)=\frac{1}{\vec k^2} \tag{6}
$$
One can show that this extended propagator is the same as the covariant propagator in the Feynman gauge $\xi=1$, up to (non-covariant) terms proportional to $k^\mu$, and therefore the theory is, after all, covariant.
While it is still true that $\Delta^{ij}=\langle A^i A^j\rangle$, now there is no operator corresponding to $\Delta^{00}$.

My doubt: An example
So far so good (?). Now, the cancellation of the $k^\mu$ terms can be proven, on general terms, thanks to the Ward-Takahashi identities, which state that general correlation functions are zero when contracted with $k^\mu$ (kind of: they need not be zero, but their form is highly constrained).
But it seems to me that this identities cannot be naively used in the Coulomb gauge, because we don't have $A^0$ any more, and so correlation functions do not include the $\mu=0$ component. In other words, the true correlation functions are
$$
\langle 0|T\ A^i(x)A^j(y)\psi(z)\cdots|0\rangle \tag{7}
$$
and it is impossible to contract this with $\partial_\mu$. So I guess we can still write
$$
\partial_\mu\langle J^\mu \psi_1\psi_2\cdots\rangle=\text{contanct-term} \tag{8}
$$
but, in the Coulomb gauge, there is no simple relation bewtween $\langle J^\mu\cdots\rangle$ and $\langle A^\mu\cdots\rangle$. In generalised gauges we have
$$
\langle A^\mu \cdots\rangle=\frac{1}{\partial^2}\langle J^\mu\cdots\rangle+\text{contanct-term} \tag{9}
$$
but this is no longer true in the Coulomb gauge.
Or put it another way: we can write the WT identities with $J^\mu$ instead of $A^\mu$, but this doesn't take us much further: in order to show cancellation of $k^\mu$ terms we must write the correlation function in terms of $A^\mu$. In covariant gauges this is easy: using $-\partial^2 A^\mu=J^\mu$, we can replace any $A^\mu$ in a corralation function, up to a propagator $1/k^2$ and a contact term (due to the $T$ symbol). But this is no longer true in the Coulomb gauge: for one thing, $\langle A^\mu\cdots\rangle$ is undefined for $\mu=0$.
After all, the equations $(8)$ and $(9)$ are an essential ingredient for the proof of the fact that $k^\mu$ terms do not contribute to measurable predictions (for example, see Srednicki's book, chapters 67-68; in particular, eqs. [67.9] - [67.12]).

My question
In a more general context, what can be said about the Schwinger-Dyson equations in the Coulomb gauge? Are they still valid? can we use $A^\mu$, or must we restrict ourselves to $A^i$? How can we efficiently use the fact that the non-local interaction is equivalent to an extended propagator? How is this implemented on the level of the DS equations instead of on the level of Feynman diagrams?
As far as I know, for general covariant gauges the DS equations work just fine, and they are valid for the four components of $A^\mu$.

Further info
Weinberg, in his QFT book, chapter 9.6 introduces an auxiliary field $A^0$ and groups it together with the old, physical $A^i$ in such a way that $A^\mu$ can be formally thought of as a vector field, with $A^0$ behaving like it was its true temporal component (see page 415, in particular the discussion below 9.6.6). This is pretty much what I would like: I want to use the Coulomb gauge, but I want to use $A^\mu$ as well, in order to be able to use the WT identities or any DS equation in general. But here Weinberg uses path integrals instead of operators. Sigh.
If any one knows some good source where QED in the Coulomb gauge is thoroughly discussed (without path integrals!) it would be very very welcome.
 A: Although $A^0$ is not a dynamical field, it is still a perfectly well-defined quantum field and hence can be used to contract with other vector expressions.
The standard source for canonical QED in the Coulomb gauge is the classical (but now somewhat dated) textbook by Bjorken and Drell, Relativistic quantum field theory.
A: *

*OP is asking about different gauge-fixing choices in QED, such as, e.g. Coulomb gauge, Lorenz gauge in Feynman gauge $\xi=1$, etc. Physical observables are invariant under gauge symmetry, and do not depend on gauge-fixing choices.

*OP is specifically asking about the fate of non-gauge-invariant $n$-point correlators involving the $A^0$-field using the Coulomb gauge. This is perhaps easiest seen in the path integral formulation. One should keep in mind the gauge-fixing terms upstairs$^1$ in the gauge-fixed QED action. 


*

*If there are no $A^0$-fields downstairs in the $n$-point correlator, the $A^0$ path integration is a Gaussian integral, which produces the usual Coulomb potential term between charge sources upstairs.

*If there are $A^0$-fields downstairs in the $n$-point correlator, the path integral can still be performed via the usual bag of Feynman tricks. 


*The pertinent Schwinger-Dyson equations and Ward identities with appropriate gauge-fixing terms still hold in arbitrary gauge.
-- 
$^1$ A correlator function $\langle F \rangle$ in the path integral formulation is  schematically of the form $\langle F \rangle=\frac{1}{Z} \int F e^{\frac{i}{\hbar}S}$. The words downstairs and upstairs refer to $F$ and $S$, respectively, for hopefully obvious reasons.
