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.