Dirac Current Spectral Representation

I'm reading Strocchi's book on The Non-Perturbative Foundations of Quantum Field Theory.

In the chapter concerning point-splitting regularization, where the free Dirac current is defined as follows $$j_\mu(x)\equiv \lim_{\epsilon\to0}\left[\bar\psi(x+\epsilon)\gamma_\mu\psi(x)-\langle\bar\psi(x+\epsilon)\gamma_\mu\psi(x)\rangle\right]$$ to give a precise meaning to the product of two distributions evaluated at the same point, there is the following claim about the spectral properties of the current commutator function: $$\langle[j_\mu(x),j_\nu(y)]\rangle = (\Box g_{\mu\nu}-\partial_\mu\partial_\nu)\int d\rho(\mu^2) i\Delta(x-y;\mu^2)$$ where the $i\Delta(x-y;\mu^2)$ is the commutator function of a free scalar field o mass $\mu$, and $$\rho(\mu^2)= \frac{1}{3(2\pi)^2}\sqrt{1-\frac{4m^2}{\mu^2}}\left(1+\frac{2m^2}{\mu^2}\right).$$ The hint is to calculate this spectrum by inserting a complete set of states of electron-positron pairs. How can I do this?

I tried insertig them in the 2-point function obtaining: (I keep the spin label implicit in the sum over the complete set) $$\langle j_\mu(x) j_\nu(y)\rangle= \int d\mu^2 \int \frac{d^3k}{(2\pi)^3}\frac{1}{2\omega_\mu(\mathbf{k})} \langle j_\mu(x)|\mathbf{k};\mu^2\rangle\langle\mathbf{k};\mu^2| j_\nu(y)\rangle=\\ =\int d\mu^2 \int \frac{d^3k}{(2\pi)^3}\frac{1}{2\omega_\mu(\mathbf{k})}e^{-ik\cdot(x-y)} \langle j_\mu(0)|\mathbf{k};\mu^2\rangle\langle\mathbf{k};\mu^2| j_\nu(0)\rangle$$ now, I can consider a Lorentz boost $\Lambda_\mathbf{k}$ such that it brings us in the rest frame of the single particle state $|\mathbf{k};\mu^2\rangle$: $U(\Lambda_\mathbf{k})|\mathbf{k};\mu^2\rangle=|\mathbf{0};\mu^2\rangle$, and by covariance of the (vector) current $j_\mu$ we have $$\langle j_\mu(x) j_\nu(y)\rangle= \int d\mu^2 \langle j_\rho(0)|\mathbf{0};\mu^2\rangle\langle\mathbf{0};\mu^2| j_\lambda(0)\rangle \int \frac{d^3k}{(2\pi)^3}\frac{1}{2\omega_\mu(\mathbf{k})}e^{-ik\cdot(x-y)}(\Lambda_\mathbf{k})^\rho_\mu(\Lambda_\mathbf{k})^\lambda_\nu.$$ Assuming that the above structure is correct, since the differential operator $\Box g_{\mu\nu}-\partial_\mu\partial_\nu$ is dictated by Lorentz covariance and the current conservation, we have indeed: $$\langle j_\mu(x) j_\nu(y)\rangle=(\Box g_{\mu\nu}-\partial_\mu\partial_\nu)\int d\rho(\mu^2) i\Delta^+(x-y;\mu^2)=\\ =(\Box g_{\mu\nu}-\partial_\mu\partial_\nu)\int d\rho(\mu^2) \int \frac{d^3k}{(2\pi)^3}\frac{1}{2\omega_\mu(\mathbf{k})}e^{-ik\cdot(x-y)}=\\ =\int d\rho(\mu^2) \int \frac{d^3k}{(2\pi)^3}\frac{-\mu^2 g_{\mu\nu}+k_\mu k_\nu}{2\omega_\mu(\mathbf{k})}e^{-ik\cdot(x-y)}$$ and tracing and equating both expressions we get, since $\Lambda_\mu^\rho g^{\mu\nu}\Lambda_\nu^\lambda=g^{\rho\lambda}$: $$\rho(\mu^2)=-\frac{g^{\rho\lambda}}{3\mu^2}\langle j_\rho(0)|\mathbf{0};\mu^2\rangle\langle\mathbf{0};\mu^2| j_\lambda(0)\rangle.$$ My problem is here: the matrix element $$\langle j_\rho(0)|\mathbf{0};\mu^2\rangle$$ seems to vanish identically since there are three annihilation/creation operator in a vacuum expectation value! What did I miss?

(EDIT)

Ok, I guess I must sum over particle-antiparticle intermediate states and not just single particle states, which lead the v.e.v. to vanish identically.

Indicating the invariant phase space element as $$\int d\Pi(\mathbf{k})\equiv\int \frac{d^3k}{(2\pi)^3}\frac{1}{2\omega_m(\mathbf{k})}$$ we get, letting $\mathbf{k}$ and $\mathbf{p}$ be the momentum respectively of electron and positron and $q^2$ the square of the center of mass energy $$\langle j_\mu(x) j_\nu(y)\rangle= \int dq^2 \int d\Pi(\mathbf{k})\int d\Pi(\mathbf{p}) \langle j_\mu(x)|\mathbf{k},\mathbf{p};q^2\rangle\langle\mathbf{k},\mathbf{p};q^2| j_\nu(y)\rangle$$ now, I can consider a Lorentz boost $\Lambda$ that brings us in the center of mass frame of the particle-antiparticle state and by covariance of the (vector) current $j_\mu$ we have $$\langle j_\mu(x) j_\nu(y)\rangle= \int dq^2 \langle j_\rho(0)|\mathbf{p},-\mathbf{p};q^2\rangle\langle\mathbf{p},-\mathbf{p};q^2| j_\lambda(0)\rangle \frac{|\mathbf{p}|}{(2\pi)^24E_{CM}} \int \frac{d^3P}{(2\pi)^3}\frac{1}{2\omega_m(\mathbf{P})}e^{-iP\cdot(x-y)}(\Lambda)^\rho_\mu(\Lambda)^\lambda_\nu,$$ where we have performed a simplification in the 2-body phase space integral (Peskin, page 107), even though I'm not sure on how to get rid of the integration over solid angle $d\Omega$, which usually enters in the definition of cross section.

By tracing with $g^{\rho\lambda}$ we get, comparing to the general expression above for the 2-point fuction: $$\rho(q^2)=-\frac{1}{3q^2}\frac{|\mathbf{p}|}{(2\pi)^24E_{CM}}\langle j_\rho(0)|\mathbf{p},-\mathbf{p};q^2\rangle\langle\mathbf{p},-\mathbf{p};q^2| j^\rho(0)\rangle;$$ with a bit of Dirac algebra: letting $p_1=(\omega(\mathbf{p}),\mathbf{p})$, $p_2=(\omega(\mathbf{p}),-\mathbf{p})$ $$\langle j_\rho(0)|\mathbf{p},-\mathbf{p};\mu^2\rangle\langle\mathbf{p},-\mathbf{p};\mu^2| j^\rho(0)\rangle=-\text{Tr}\left[(\hat p_1+m)\gamma_\rho(\hat p_2-m)\gamma^\rho\right]=-8(2m^2+p_1\cdot p_2)$$ and since $$|\mathbf{p}|=\sqrt{\omega(\mathbf{p})^2-m^2}=\frac{E_{CM}}{2}\sqrt{1-\frac{4m^2}{q^2}}$$ we get $$\rho(q^2)=-\frac{1}{3q^2}\frac{\frac{E_{CM}}{2}\sqrt{1-\frac{4m^2}{q^2}}}{(2\pi)^24E_{CM}}[-8(2m^2+p_1\cdot p_2)]=\frac{1}{6(2\pi)^2}\sqrt{1-\frac{4m^2}{q^2}}\left(\frac{2m^2}{q^2}+1\right)$$ which looks good a part from the extra factor 2 in the denominator.