It has recently come to my attention, though reading the notes of a course on QFT that I've started, that there seems to be an "ambiguity" in, or at least two distinct ways of, calculating scattering amplitudes which include Feynman diagrams with propagators connected to the same interaction vertex on both ends. I'll use $\varphi^4$ theory as an example.
The first method is the one seen used in almost all answers to questions on $\varphi^4$ scattering amplitudes on this site, as well as every online lecture I've watched up until now and (I think) every QFT textbook I've read, which is where the following interaction Hamiltonian density is used:
$$\mathcal{H}_{int,I} = \frac{\lambda}{4!} \varphi^4.$$
This leads to the following $1 \rightarrow 1$ amplitude up to first order in $\lambda$:
$$\langle q|p \rangle = \frac{\langle 0|\mathcal{T}\{ \alpha_\textbf{q} \mathcal{S} \alpha_\textbf{p}^\dagger \}|0 \rangle}{\langle 0|\mathcal{S}|0 \rangle}$$
$$= \langle 0|\alpha_\textbf{q} \alpha_\textbf{p}^\dagger|0 \rangle - 12i \lambda \int d^4x \langle 0|\mathcal{T}\{ \alpha_\textbf{q} \varphi(x) \}|0 \rangle \langle 0|\mathcal{T}\{ \varphi(x) \varphi(x) \}|0 \rangle \langle 0|\mathcal{T}\{ \varphi(x) \alpha_\textbf{p}^\dagger \}|0 \rangle$$
$$= 2E_\textbf{p} (2\pi)^3 \delta^{(3)}(\textbf{p} - \textbf{q}) - 12 \lambda (2\pi)^4 \delta^{(4)}(p - q) \int \frac{d^4k}{(2\pi)^4} \frac{1}{k^2 + m^2 - i\varepsilon} + O(\lambda^2)$$
At this point, renormalisation is put forward as the method for dealing with the divergent piece.
The second method, which is the one in my lecture notes and has occasionally been asked about (such as here, here and here), apparently completely avoids this altogether by using a normal-ordered interaction Hamiltonian density:
$$\mathcal{H}_{int,I} = \frac{\lambda}{4!} \mathbf{:} \varphi^4 \mathbf{:}.$$
I think this is because, according to point (II) of this answer, "we should only include contractions between pairs of operators who belong to different normal order symbols", which immediately gets rid of the divergent piece due to the presence of the $\langle 0|\mathcal{T}\{ \varphi(x) \varphi(x) \}|0 \rangle$ factor.
The answer references a proof of this in the context of string theory and radial ordering. Is there a source for a similar proof in the context of QFT and time ordering?
Why is the second method not more widely used? I appreciate that renormalisation is needed anyway for other types of divergent diagrams, but it seems to me that it would allow us to save some time by not bothering writing some of the diagrams in the first place, and is also more in line with the fact that normal ordering is used in other contexts such as keeping the canonical energy density finite.
