The Feynman propagator for the free electron field is the Fourier transform w.r.t. $y$ of the time-ordered 2-point VEV $\left<0\right|\mathcal{T}[\hat\psi(x)\hat\psi(x+y)]\left|0\right>$, taking $\hbar=c=1$,
$$\mathcal{F}[\left<0\right|\mathcal{T}[\hat\psi(x)\hat\psi(x+y)]\left|0\right>](k)=
\frac{k\!\cdot\!\gamma+m_e}{k^2-m_e^2+\mathrm{i}\epsilon}.$$
In QED, after renormalization, we obtain, valid in the infrared, when $k^2-m_e^2$ is small (given in this neat form, representing the sum of many perturbative integrals, in, for example, Appelquist & Carrazone, Phys.Rev.D 11, 2856 (1975)),
$$\mathcal{F}[\left<0\right|\mathcal{T}[\hat\psi(x)\hat\psi(x+y)]\left|0\right>](k)=
\frac{k\!\cdot\!\gamma+m_e}{(k^2-m_e^2+\mathrm{i}\epsilon)^{1-{\alpha_{EM}}/{\pi}}},$$
where $\alpha_{EM}\approx 1/137$.
This is often called an infraparticle propagator or a dressed particle propagator (I'm rehearsing this stuff, please bear with me, and feel free to comment on anything that seems to need some change).
When we move to renormalization group methods, the coupling constant becomes a function of the renormalization scale. In this case, if we measure the Fourier transform of the 2-point VEV, using the operator $\mathcal{F}[\mathcal{T}[\hat\psi(x)\hat\psi(x+y)]](k)$, the measurement scale is determined by the wave-number $k$, so I take it we can write the 2-point VEV as
$$\mathcal{F}[\left<0\right|\mathcal{T}[\hat\psi(x)\hat\psi(x+y)]\left|0\right>](k)
{{?\atop =}\atop\ }
\frac{k\!\cdot\!\gamma+m_e}{(k^2-m_e^2+\mathrm{i}\epsilon)^{1-{\alpha_R}(k)/{\pi}}},$$
where the running coupling constant $\alpha_R(k)$ is $\alpha_R(m_e)\approx 1/137$ when $k=m_e$, at about 0.5 MeV and $\alpha_R(m_Z)\approx 1/127$ when $k=m_Z$, the mass of the $Z$ particle, at about 90 GeV (see here).
My understanding is that in QED the function $\alpha_R(k)$ is an increasing function of $k$, even to the extent that there is a Landau pole at finite inverse length $m_L$, $\alpha_R(m_L)=\infty$, but that the high-energy behavior of QED has only been calculated perturbatively, so that an analytic form for $\alpha_R(k)$ is not known. I would like, however, to look at a good article, perhaps a review, that gives as closed a form for $\alpha_R(k)$ as is currently known, in as neat a form as possible (this is an implicit question that may be too much to ask of the literature, given that almost everyone has moved on to supersymmetry, noncommutative geometry, string theory, etc., etc., and QED is obviously not empirically useful at high energy).
I would also like better to understand the relationship between the running coupling constant formalism and the Källén–Lehmann formalism. Can we equate the two, at least approximately,
$$\frac{k\!\cdot\!\gamma+m_e}{(k^2-m_e^2+\mathrm{i}\epsilon)^{1-{\alpha_R}(k)/{\pi}}}
{{?\atop \approx}\atop\ }
\int\limits_{m_e}^\infty \frac{(k\!\cdot\!\gamma+m)f(m^2)\mathrm{d}m^2}{k^2-m^2+\mathrm{i}\epsilon}.$$
Here, $f(m^2)\ge 0$ is undefined for $m>m_L$ if there is a Landau pole in QED, but I suppose $f(m^2)$ would not have any poles and would approach zero fast enough for at least the integral $\int\limits_0^\infty f(m^2)\mathrm{d}m^2$ to exist in a well-defined theory? We could equate the two exactly if we were talking about a scalar Feynman propagator, because we could solve for $\alpha_R(k)$, but the change from $k\!\cdot\!\gamma+m_e$ to $k\!\cdot\!\gamma+m$ doesn't look good because the two expressions have different effects for different components of the Dirac spinor.
Better understanding the 2-point VEVs would of course leave the $n$-point VEVs to think about. Uggh.
Finally, there are three explicit Questions here (marked by question marks!), but I would also be interested in any ruminations from this starting point that include references.
Best Answer
The propagator as defined by the time-ordered 2-point VEV ⟨0|T[ψ^(x)ψ^(x+y)]|0⟩ is clearly not invariant under the local gauge transformation
$$A_{\mu}(x)\to A_{\mu}(x)+\partial_\mu\Lambda(x)\ \ \ ; \ \ \ \psi(x) \to \exp(ie\Lambda(x)\psi(x)$$
since the VEV acquires an extra phase factor $\exp(ie\Lambda(y))$ under the transformation. Schwinger, among many others, noted that you can restore gauge invariance by inserting a factor
$$ \Phi(x) = \exp\left(-ie\int_0^y A_{\mu}dx^{\mu}\right)$$
between the two fermion operators in the time-ordered product.
The non-integer exponent in the first equation with a ? is actually gauge dependent (depends on the coefficient of the $k_\mu k_\nu/k^2$ term in the photon propagator). It was already studied by Soviet physicists in the 1950s -- as I recall, it is discussed in the field theory textbook by Bogolyubov and Shirkov. I believe there has been a calculation including the phase factor that leads to the electron propagator having a simple pole in $k^2$, but I do not have a reference handy. In any event, the running of the fine structure constant is not relevant to the gauge invariant propagator.