Diagrammatics of a current-current correlation function $\langle 0| T\{J^{\mu}(x) J^{\nu}(0)\}|0\rangle$ Say $J^{\mu} = \bar{\psi} \gamma^{\mu} \psi$ is the QED current. While it is clear to me how to compute something like 
$$\langle 0 |T\{ \bar \psi(x) \psi(0)\} |0\rangle$$
using a Feynman diagram expansion, it is not clear to me how I would go about computing 
$$\langle 0| T\{J^{\mu}(x) J^{\nu}(0)\}|0\rangle.$$
How do I read off the Feynman rules? This is an important step for deep-inelastic scattering analysis, for instance (c.f. Peskin 18.5). And how about for a general local operator 
$$\langle 0| T\{\mathcal{O}_1^{\mu_i}(x) \mathcal{O}^{\nu_i}_2(0)\}|0\rangle?$$
Knowing how to compute this in perturbation theory would be useful, for instance, if I wanted to compute its OPE. If there is any physical interpretation, that would be helpful to know as well. A (sketch of a) derivation starting from fundamentals such as the path integral would be great.
 A: The diagrammatic rules are the standard ones of QED. Computing 
$
\langle T\{J^\mu(x)J^\nu(0)\}\rangle
$
amounts to compute 
$$
\langle T\{\bar\psi(x)\gamma^\mu\psi(x)\bar\psi(0)\gamma^\nu\psi(0)\}\rangle=\gamma^\mu_{ab}\gamma^\nu_{cd}\langle T\{\bar\psi_a(x)\psi_b(x)\bar\psi_c(0)\psi_d(0)\}\rangle
$$
in perturbation theory. One just needs to draw the diagrams corresponding to the four-point function on the rhs, and contract with the $\gamma$ matrices at the end. 
For example, the lowest order correspond to that of a free theory (no photon exchange), and one gets, thanks to Wick's theorem,
$$
\langle T\{\bar\psi_a(x)\psi_b(x)\bar\psi_c(0)\psi_d(0)\}\rangle=\langle T\{\bar\psi_a(x)\psi_b(x)\}\rangle\langle T\{\bar\psi_c(0)\psi_d(0)\}\rangle-\langle T\{\bar\psi_a(x)\psi_d(0)\}\rangle\langle T\{\bar\psi_c(0)\psi_b(x)\}\rangle.
$$
In the case of more complicated operators, one just needs to replace the operator by its definition in terms of its field content and proceed in the same way.
Alternatively, the operators can be seen as some kind of new vertices (here each $J$ is a two-leg vertex), and the diagrammatic game is to write all diagrams of a given order that connect the legs of one of these vertex to the legs of the other.
A: The way that this is usually done in QED is deeply related to the Ward Identities and QED's LSZ reduction formula. This discussion is taken from Chapter 67-8 in Srednicki's book "Ward Identities in QED I-II". 
In particular it is related to the amputated correlators of  $\langle 0 |T\{ A^\mu(x) A^\nu (0)\} |0\rangle$, which are the photon-photon scattering terms, which are related to the (momentum space) diagrams $\langle 0 |T\{ A^\mu(-p) A^\nu (p)\} |0\rangle_{amputated} $ of the form below (there could be higher loop contributions if you're masochistic enough, but this is the only diagram that contributes at 1 loop after amputation). We will consider amputated diagrams.

In particular, we'll have that 
$\langle 0 |T\{ J^\mu(x) J^\nu (0)\} |0\rangle = \int \frac{d^4 p}{(2 \pi)^4} e^{i p x} [\langle 0 |T\{ A^\mu(-p) A^\nu (p)\} |0\rangle_{amputated} + \text{contact terms} ]$
where the contact terms are delta-function like terms that do not affect the divergence structure of the amplitudes.
The way to see this is to use the Ward identities and (an analog of) the LSZ scattering formula. The LSZ formula for photons relates photon-photon scattering amplitudes functions to derivates of photon-photon correlators which amputate the external legs. It says 
$\langle f | i \rangle = i \epsilon_\mu^f \epsilon_\nu^i \int d^4x e^{-ipx} (-\partial_x^2) \int d^4x e^{+ipy} (-\partial_y^2) \langle 0 | A^\mu(x)  A^{\nu} (y) |0 \rangle)$
where $\langle f | i \rangle$ is the amplitude going from the initial to final photon state, and $\epsilon_\mu^j$ are the initial and final polarization vectors, and the $\partial^2$ terms get rid of the leg contributions. From here, we will 'amputate' further and only consider this amplitude before contracting with the polarization. 
The next ingredient is the Ward Identities for QED. The Ward identities state that if $\frac{\delta \mathcal{L}}{\delta A^\mu}$ are the equations of motion for your Lagrangian $\mathcal{L}$, then 
$\langle 0| \frac{\delta \mathcal{L}}{\delta A^\mu}(x) O_1(y) O_2(z)... |0 \rangle = 0 + \text{contact terms}$.
Here, the contact terms are zero if $x,y,z$ are distinct, and will vanish inside any scattering amplitude. In addition, if all of the $O_i$ fields are distinct from $A_\mu$, then they are identically zero. For QED, the equations of motion will be 
$\frac{\delta \mathcal{L}}{\delta A^\mu} = - \partial^2 A_\mu + J_\mu$ 
where $J_\mu = \bar \psi \gamma_\mu \psi$ is the coupling to the electromagnetic field. (I'm ignoring some counterterms that are needed actually. Srednicki includes these). This tells us that up to contact terms, we can replace $\partial^2 A_\mu$ with $\bar \psi \gamma_\mu \psi$. In particular, we'll have exactly what we stated earlier,
$\langle 0 |T\{ J^\mu(x) J^\nu (0)\} |0\rangle = \int \frac{d^4 p}{(2 \pi)^4} e^{i p x} [\langle 0 |T\{ A^\mu(-p) A^\nu (p)\} |0\rangle_{amputated} + \text{contact terms} ]$
where the contact terms are zero up to this loop order. 
