We can write the Fourier transform of
$\langle 0|\mathcal{T}A_{\nu}(x)\psi(x_1)\bar\psi(x_2)|0\rangle$
as $$S(p) D_{\nu\alpha}(q) \ e\,\Gamma^{\alpha}(p,q,p+q)S(p+q)$$ where $S(p)$ is the full fermion propagator, $D_{\nu\alpha}(q)$ is the full photon propagator, $\Gamma^{\alpha}(p,q,p+q)$ is the proper vertex function, and an overall momentum conservation delta function has been dropped. Similarly, we can write the Fourier transform of $\langle 0|\mathcal{T}j^{\mu}(x)\psi(x_1)\bar\psi(x_2)|0\rangle
$ as
$$S(p)V^{\mu}(p,q,p+q)S(p+q)$$ where $V^{\mu}(p,q,p+q)$ is a vertex function that we want to relate to $\Gamma^{\mu}(p,q,p+q).$ The vertex function $V^{\mu}(p,q,p+q)$ enters into the derivation of the Ward-Takahashi identity in Peskin and Schroeder on page 311, but the Ward-Takahashi identity is normally stated in terms of $\Gamma^{\mu}(p,q,p+q)$. Your conundrum (as I understand it) is that according to your analysis of the Schwinger-Dyson equation, $V^{\mu}(p,q,p+q)$ and $\Gamma^{\mu}(p,q,p+q)$ ought to differ by a factor of $Z_3$, but this contradicts the usual statement of the Ward-Takahashi identity where no such factor of $Z_3$ appears. I will argue from the Schwinger-Dyson equation that the longitudinal parts (in $q^{\mu}$) of $V^{\mu}(p,q,p+q)$ and $\Gamma^{\mu}(p,q,p+q)$ are equal, but that the transverse parts differ by the factor of $Z_3$ that you have found. Since only the longitudinal part enters into the Ward-Takahashi identity, the factor of $Z_3$ does not inter into that identity. You may want to review page 246 of Peskin and Schroeder. There they show that only the transverse part of the photon propagator is modified by the self-energy, but that in calculating Feynman diagrams we can simplify the analysis by including the self-energy in the longitudinal part as well because the longitudinal part does not contribute to the Feynman diagrams due to the Ward identity. However, the Schwinger-Dyson equation involves an inverse propagator which does not arise in Feynman diagrams and we need to reevaluate where the self-energy does and does not enter.
Specializing the Schwinger-Dyson equation to the case of $\langle 0|\mathcal{T}A_{\nu}(x)\psi(x_1)\bar\psi(x_2)|0\rangle$ and Fourier transforming, we have
$$\tag{1} (D^{(0)\mu\nu}(q))^{-1} D_{\nu\alpha}(q) S(p) \ e\,\Gamma^{\alpha}(p,q,p+q)S(p+q) = \\ e\, S(p)V^{\mu}(p,q,p+q)S(p+q)$$
where $(D^{(0)\mu\nu}(q))^{-1}$ is the inverse of the non-interacting photon propagator. The Dyson equation for the photon propagator is
$$\tag{2} D_{\nu\alpha}(q) = D^{(0)}_{\nu\alpha}(q) + D^{(0)}_{\nu\beta}(q) i\Pi^{\beta\gamma}(q)D_{\gamma\alpha}(q) ,$$
so
$$\tag{3} (D^{(0)\mu\nu}(q))^{-1} D_{\nu\alpha}(q) = \delta^{\mu}_{\alpha} + i\Pi^{\mu\gamma}(q)D_{\gamma\alpha}(q).$$
Equation (1) then implies
$$\tag{4}\Bigl(\delta^{\mu}_{\alpha} + i \Pi^{\mu\gamma}(q)D_{\gamma\alpha}(q)\Bigr) \Gamma^{\alpha}(p,q,p+q) = V^{\mu}(p,q,p+q).$$
The Ward identity forces the longitudinal part of $\Pi^{\mu\gamma}(q)$ to vanish; that is, $q_{\mu}\Pi^{\mu\gamma}(q) = 0.$
Contracting equation (4) with $q_{\mu}$,
we therefore have
$$\tag{5} q_{\alpha} \Gamma^{\alpha}(p,q,p+q) = q_{\mu}V^{\mu}(p,q,p+q)$$
so no factor of $Z_3$ appears between the longitudinal parts of $\Gamma^{\alpha}(p,q,p+q)$ and $V^{\mu}(p,q,p+q)$ and therefore no factor of $Z_3$ appears in the Ward-Takahashi identity.
The transverse component does not enter the Ward identity for the vertex function but it is useful to consider the transverse component to illustrate where the factor of $Z_3$ does arise. Define $\Pi(q^2)$ by the equation
$\Pi^{\mu\nu}(q) = q^2(g^{\mu\nu} - q^{\mu}q^{\nu}/q^2)\Pi(q^2).$ The quantity $(g^{\mu\nu} - q^{\mu}q^{\nu}/q^2)$ can be described as a projection operator that projects out the transverse part of a vector.
Contracting equation (4) with $(g_{\nu\mu} - q_{\nu}q_{\mu}/q^2)$ and using the fact that $\Pi^{\mu\gamma}(q)$ is already transverse, we have
$$\tag{7} \Bigl(g_{\nu\alpha} - q_{\nu}q_{\alpha}/q^2 + i q^2\Pi(q^2)D^T_{\nu\alpha}(q)\Bigr) \Gamma^{\alpha}(p,q,p+q) =\\ \bigl(g_{\nu\mu} - q_{\nu}q_{\mu}/q^2\bigr)V^{\mu}(p,q,p+q),$$
where $$D^T_{\nu\alpha}(q) = \frac{-i}{q^2 (1-\Pi(q^2))} \bigl(g_{\nu\alpha} - q_{\nu}q_{\alpha}/q^2\bigr)$$ is the transverse part of the photon propagator (see page 246, Peskin and Schroeder). Equation (7) can then be written
$$\tag{8} \bigl(g_{\nu\alpha} - q_{\nu}q_{\alpha}/q^2\bigr) \Bigl(1/\bigl(1-\Pi(q^2)\bigr)\Bigr) \Gamma^{\alpha}(p,q,p+q) =\\ \bigl(g_{\nu\mu} - q_{\nu}q_{\mu}/q^2\bigr)V^{\mu}(p,q,p+q).$$
Now consider $q^2$ small enough that $\Pi(q^2)\approx \Pi(0)$ and use the relation (Peskin and Schroeder, page 246)
$$Z_3 = \Bigl(1/\bigl(1-\Pi(0)\bigr)\Bigr).$$ We have
$$\tag{9} \bigl(g_{\nu\alpha} - q_{\nu}q_{\alpha}/q^2\bigr) Z_3 \Gamma^{\alpha}(p,q,p+q) = \bigl(g_{\nu\mu} - q_{\nu}q_{\mu}/q^2\bigr)V^{\mu}(p,q,p+q).$$ So we see that the transverse parts of $V^{\mu}(p,q,p+q)$ and $\Gamma^{\mu}(p,q,p+q)$ differ by a factor of $Z_3.$
The idea is exactly the same as in single particle QM. Take the expectation value $\langle 0 | \hat\phi(x) | 0 \rangle$ and insert a complete set of field eigenstates $|\varphi\rangle$; these are the states that satisfy $\hat\phi(x) |\varphi\rangle = \varphi(x)|\varphi\rangle$. We do this with a path integral:
$$\langle 0 | \hat\phi(x) | 0 \rangle = \int \mathcal{D}\varphi\ \langle 0 | \hat\phi(x)|\varphi\rangle \langle\varphi | 0 \rangle = \int \mathcal{D}\varphi\ \langle 0 |\varphi\rangle \langle\varphi | 0 \rangle \varphi(x).$$
So your probability distribution is $p[\varphi] = |\langle\varphi | 0 \rangle|^2$, just like in regular QM. I'm not sure yet how to calculate this, but by analogy with the harmonic oscillator it would probably be along the lines of $\exp(-\varphi^2)$.
Edit: I think I've managed to calculate it, though I make no guarantees of this being correct. In particular I'm pretty sure there should be some $2\pi$'s around.
We proceed by analogy with the harmonic oscillator. You can find the ground state wavefunction by using the fact that $\langle x | \hat a | 0 \rangle = 0$; this gives you a differential equation for $\langle x | 0 \rangle = 0$. In much the same way, we write $\langle \varphi | \hat a_\mathbf{p} | 0 \rangle = 0$. We then express $\hat a_\mathbf{p}$ in terms of $\hat \phi$ and $\hat \pi$, and use the fact that acting on a wavefunctional $\hat \pi$ is $-i \delta / \delta \varphi(x)$. I did it following Peskin and Schroeder's conventions for normalizations and working in the Schrödinger picture, and arrived at the following equation (kx is the dot product $\mathbf{k} \cdot \mathbf{x}$):
$$ \int \mathrm{d}^3x\ e^{-ikx} \left( \omega_k \varphi(x) \langle\varphi | 0 \rangle + \frac{\delta}{\delta \varphi(x)}\langle\varphi | 0 \rangle \right) = 0$$
We will eventually express everything in terms of the Fourier transform $\hat \varphi$ of $\varphi$ (who knew, right?), so we will need to use the fact that
$$ \int \mathrm{d}^3x\ e^{-ikx} \frac{\delta}{\delta \varphi(x)}\langle\varphi | 0 \rangle = \frac{\delta}{\delta \hat \varphi} \langle\varphi | 0 \rangle \bigg|_{-k}$$
so our equation becomes
$$\frac{\delta}{\delta \hat \varphi (k)} \langle\varphi | 0 \rangle = - \omega_k \hat \varphi (-k) \langle\varphi | 0 \rangle$$
and we can guess the solution
$$\langle\varphi | 0 \rangle = \exp \left( -\frac12 \int \mathrm{d}^3k\ \omega_k \hat \varphi(k) \hat \varphi(-k) \right).$$
You could rewrite this in position space but the Fourier transform of $\omega_k$ involves Bessel functions so I don't know if it's worth it.
Best Answer
You just need to show that $1_\mathrm{1p}$ acts as the identity on one-particle states: $$ 1_\mathrm{1p}|\boldsymbol q\rangle=\left[\int\frac{d^3\textbf{p}}{(2\pi)^{3}}|\textbf{p}\rangle\frac{1}{2E_\textbf{p}}\langle\textbf{p}|\right]|\boldsymbol q\rangle=\int\frac{d^3\textbf{p}}{(2\pi)^{3}}|\textbf{p}\rangle\frac{1}{2E_\textbf{p}}\overbrace{\langle\textbf{p}|\boldsymbol q\rangle}^{(2\pi)^32E_p\delta(\boldsymbol p-\boldsymbol q)} $$ which indeed equals $|\boldsymbol q\rangle$, as one would expect for the identity operator. Moreover, using $\langle \boldsymbol p|\boldsymbol q_1,\boldsymbol q_2\cdots,\boldsymbol q_n\rangle=0$, we see that $(1)_\mathrm{1p}$ is orthogonal to multiparticle states.
In other words, $1_\mathrm{1p}$ is the identity when acting on one-parcle states, and it annihilates states with two or more particles. Therefore, it is the projector into the one-particle subspace.
In general, if your basis is normalised to $\langle \varphi|\varphi'\rangle=f(\varphi)\delta(\varphi-\varphi')$, then the identity is $$ 1=\int\frac{\mathrm d\varphi}{f(\varphi)}|\varphi\rangle\langle\varphi| $$ as can be seen by acting with $1$ on a general state $|\varphi'\rangle$, and checking that one indeed obtains $1|\varphi'\rangle=|\varphi'\rangle$.
The operator $1_\mathrm{1p}$ projects into the one-particle subspace. You can say that $1\approx 1_\mathrm{1p}$ if and only if you can neglect the multiparticle contribution. For an example, see Källén–Lehmann spectral representation: here, the one-particle subspace generates a pole at $p^2=m^2$. Therefore, as long as $p^2\sim m^2$, the multiparticle contribution is negligible as compared to the one-particle contribution: the latter generates a pole, while the former is regular. In this case, you can approximate $1\approx 1_\mathrm{1p}$. For other cases, there is no general rule: you can neglect the multiparticle contribution if and only if it the one-particle part dominates, which you can only conclude if you explicitly calculate both parts and compare them.
The projector into the $n$ particle subspace is straightforward: $$ 1_n=\frac{1}{n!}\int\frac{d^3\textbf{p}_1}{(2\pi)^{3}2E_\textbf{p}^1}\cdots \frac{d^3\textbf{p}_n}{(2\pi)^{3}2E_\textbf{p}^n}|\textbf{p}_1,\cdots,\boldsymbol p_n\rangle\langle\textbf{p}_1,\cdots\boldsymbol p_n| $$ as can be seen by acting on a general $n$ particle state $|\boldsymbol q_1,\cdots,\boldsymbol q_n\rangle$.
With this, the identity over the full Hilbert space is $$ 1=|0\rangle\langle 0|+\sum_{n=1}^\infty 1_n $$