Question on how to actually use momentum space Feynman rules in $phi^4$-theory

1. Momentum-space Feynman rules are almost always used in calculations, and this entails transforming the propagators (from the field contractions) to momentum space. However, the integrals over spacetime that arise from the Dyson series expansion of the S-matrix still remain - interchanging the order of integration (as physicists do) and performing these spacetime integrals furnishes the delta functions at vertices that effect momentum conservation.

2. For simplicity, I'll talk about $phi^4$ theory here, but the argument here is easily generalisable. Every momentum-space propagator brings with it another $int mathrm d^4 k_{(i)}$ to the calculation - this means that a priori, the internal momenta are off-shell (to be expected from virtual propagators), and hence each such $k_{(i)}$ is an additional unconstrained degree of freedom. The delta functions produced by the integration over spacetime have the effect of constraining the momentum of the internal lines - killing the momentum integral, as it were. Obviously, there will be one such delta function for each vertex in the diagram. Exactly one of these forces overall momentum conservation, and each of the remaining ones constrain the momenta of one internal line (i.e. reduces the number of degrees of freedom by one). Consequently there will still remain $n_text{int. lines} + 1 - n_text{vertices}$ of these $int mathrm d^4 k_{(i)}$ integrals.
   
   Since a delta function like $delta(k_{(i)} - k_{(j)})$ can be used to kill either the $int mathrm d^4 k_{(i)}$ integral or the $int mathrm d^4 k_{(j)}$ integral, it really makes no difference which one you choose to call "unconstrained", since the "constrained" one can in principle be calculated from the momentum conservation equation at that vertex. Finally, an external line can never have undetermined momentum - these are essentially the "initial conditions" for your scattering process.

I'll use an explicit example here, take this diagram:

Let's work in one dimension and use "scalar momentum" - suppose the momentum going in is $10$. By overall momentum conservation, the momentum coming out should also be $10$. The Feynman rules dictate that the internal lines should take on all real values, and we should sum over all of these processes (so supposedly we have three undetermined momenta here). However, momentum conservation dictates that these three values can't be arbitrary real numbers, but should sum to $10$. So we would have to sum over processes with $k_1 = 5, k_2 = 4, k_3 = 1$ and $k_1 = -pi, k_2 = 10, k_3 = pi$ and what have you. The main thing to realise here, is that for any one of these processes, we can calculate any one of the momenta from the others - for example, if you know that the incoming momentum is $10$, $k_1 = 4$ and $k_2 = 1$, you know that $k_3$ is constrained to be $5$. The choice of line to be constrained is completely arbitrary (in general, any line connected to the same vertex can be chosen), so you are correct in that no one internal line is favoured over another. We finally see that there are two undetermined momenta in this diagram, which agrees with $n_text{int. lines} + 1 - n_text{vertices} = 3 + 1 - 2 = 2$

Explicitly (in 3+1D), this diagram is, up to symmetry factors, given by $$ F = -lambda^2iint mathrm d^4 x mathrm d^4 y e^{iqy} D_{xy} D_{xy} D_{xy} e^{-ipx} D_{xy} rightarrow int frac{mathrm d^4 k}{(2pi)^4} frac{i}{k^2-m^2+ivarepsilon} e^{-ikcdot(x-y)} = int frac{mathrm d^4 k}{(2pi)^4} D_k e^{-ikcdot(x-y)} F = -lambda^2iint mathrm d^4 x mathrm d^4 y e^{iqy} int frac{mathrm d^4 k_1}{(2pi)^4} D_{k_1} e^{-ik_1cdot(x-y)} int frac{mathrm d^4 k_2}{(2pi)^4} D_{k_2} e^{-ik_2cdot(x-y)} int frac{mathrm d^4 k_3}{(2pi)^4} D_{k_3} e^{-ik_3cdot(x-y)} e^{-ipx} = -lambda^2 iiint frac{mathrm d^4 k_1}{(2pi)^4} frac{mathrm d^4 k_2}{(2pi)^4} frac{mathrm d^4 k_3}{(2pi)^4} D_{k_1}D_{k_2}D_{k_3} intmathrm d^4 x e^{i(sum_i k_i - p)cdot x}underbrace{intmathrm d^4 y e^{-i(sum_i k_i - q)cdot y}}_text{Momentum Conservation at Vertex} = -lambda^2 iiint frac{mathrm d^4 k_1}{(2pi)^4} frac{mathrm d^4 k_2}{(2pi)^4} frac{mathrm d^4 k_3}{(2pi)^4} D_{k_1}D_{k_2}D_{k_3} intmathrm d^4 x e^{i(sum_i k_i - p)cdot x} (2pi)^4delta(k_1+k_2+k_3-q) = -lambda^2 iint frac{mathrm d^4 k_1}{(2pi)^4} frac{mathrm d^4 k_2}{(2pi)^4} D_{k_1}D_{k_2}D_{color{red}{q - k_1 - k_2}} underbrace{intmathrm d^4 x e^{i(q - p)cdot x}}_text{Overall Momentum Conservation} = -lambda^2 (2pi)^4 delta(p - q) iint frac{mathrm d^4 k_1}{(2pi)^4} frac{mathrm d^4 k_2} {(2pi)^4} D_{k_1}D_{k_2}D_{q - k_1 - k_2} $$

Normally this overall delta function is absorbed into the definition of the matrix element, see here. As you can see, $k_1$ and $k_2$ are the undetermined momenta (but the $delta(k_1+k_2+k_3-q)$ could have killed any one of the $k_{(i)}$ integrals).