First I'll try to answer the last question you mentioned regarding commutativity in the path integral. One benefit of the path integral formulation is that the integrand behaves like an ordinary complex or Grassman number---if you're computing the trace of a product of matrices in Einstein notation, you can rearrange terms in the product as long as you keep track of index contractions and whether the numbers anticommute. In a sense, the operator character is encoded in the measure $\mathscr{D}\phi $ on the space of fields (the range of indices that you sum over).

Now your second question regarding counterterms. The reason quadratic counterterms are not included in the propagator ultimately has to do with the process of renormalization, which I will review briefly for clarity. In free theory, the propagator coincides with the two-point correlation function. In an interacting theory, we keep the propagator of the 'closest' free theory, and the two-point correlation function has perturbative corrections. The question is, what small parameter do we expand in when computing the perturbative corrections? Ideally, this parameter will be a small dimensionless ratio of two different energy scales. However, these physical quantities do not appear immediately as parameters in our model ('bare' coupling constants, masses, field normalization, etc.). We need to relate the model parameters to a small number of measured/fixed physical quantities first, like pole masses and known scattering amplitudes, and then use perturbation theory to relate the initial set of reference quantities to other `nearby' quantities (like scattering cross sections at energy scales comparable to the reference energy).

After you relate your model parameters $\mathbf m$ to some reference quantities $\mathbf x$ (fixed scattering amplitudes and pole masses), you can make new predictions $\mathbf x'(\mathbf x)$ once you measure $\mathbf x$ experimentally. The relationship between the reference amplitude $g_R$ and $g$ typically starts at first order in $g$: $g_R=g+\mathcal{O}(g^2)$. Any new physical prediction $\mathbf x'$ will be expressed as a power series (or asymptotic series) in $g_R$. In particular, the two point correlation function will equal the propagator plus some term of order $g_R$. Of course it will also include higher order corrections that generally change the locations of poles, and hence masses. But the zeroth order term is always just the propagator, the 2-point function of free theory. The terms $(Z_i-1)$ keep track of contributions that are only important at order $g_R$ and higher. In other words, you can think of $Z_i$ as `formal power series' in $g_R$, where the important information is in the form of the coefficients of the expansion, not in any function that could be associated with the series.

Finally, your first question. Constraint 2 on the normalization of single-particle states is subtle, because you must ensure that the masses of particles are fixed as well (your density of states remains on the same momentum-space hyperboloid). Imposing this constraint requires the $Z_2$ parameter. The $Z_4$ and $Z_3$ parameters are constrained by $\langle\phi\rangle_\Omega=0$, and the additional parameter can be absorbed into the interaction strength $g$. Typically $g$ is constrained by measuring a scattering amplitude.

First I'll address the last question you mentioned regarding commutativity in the path integral. One benefit of the path integral formulation is that the integrand behaves like an ordinary complex or Grassman number---if you're computing the trace of a product of matrices in Einstein notation, you can rearrange terms in the product as long as you keep track of index contractions and whether the numbers anticommute. In a sense, the operator character is encoded in the measure $\mathscr{D}\phi $ on the space of fields (the range of indices that you sum over).

Now your second question regarding counterterms. The reason quadratic counterterms are not included in the propagator ultimately has to do with the process of renormalization, which I will review briefly for clarity. In free theory, the propagator coincides with the two-point correlation function. In an interacting theory, we keep the propagator of the 'closest' free theory, and the two-point correlation function has perturbative corrections. The question is, what small parameter do we expand in when computing the perturbative corrections? Ideally, this parameter will be a small dimensionless ratio of two different energy scales. However, these physical quantities do not appear immediately as parameters in our model ('bare' coupling constants, masses, field normalization, etc.). We need to relate the model parameters to a small number of measured/fixed physical quantities first, like pole masses and known scattering amplitudes, and then use perturbation theory to relate the initial set of reference quantities to other `nearby' quantities (like scattering cross sections at energy scales comparable to the reference energy).

After you relate your model parameters $\mathbf m$ to some reference quantities $\mathbf x$ (e.g. fixed scattering amplitudes and pole masses), you can make new predictions $\mathbf x'(\mathbf x)$ once you measure $\mathbf x$ experimentally. The relationship between the reference amplitude $g_R$ and $g$ starts at linear order in $g$: $g_R=g+\mathcal{O}(g^2)$. After inverting this order-by-order to obtain $g(g_R)$, any new physical prediction $\mathbf x'$ will be expressed as a power series (or asymptotic series) in $g_R$. In particular, the two point correlation function will equal the propagator plus $\mathcal{O}(g_R)=\mathcal{O}(g)$. Of course the exact 2-point function has higher order corrections that in general can change the locations of poles, and hence masses, but the zeroth order term is always just the propagator, the 2-point function of free theory. The terms $(Z_i-1)$ keep track of contributions that are only important at order $g_R$ (equivalently order $g$) and higher, and hence do not appear in the propagator. [In other words, you can think of $Z_i$ as formal power series in $g_R$, with 0th order term equal to 1, where the important information is in the form of the coefficients of the expansion, not in any function that could be associated with the series.]

Finally, your first question. Constraint 2 on the normalization of single-particle states is subtle, because you must ensure that the masses of particles are fixed as well (your density of states remains on the same momentum-space hyperboloid). Imposing this constraint requires the $Z_2$ parameter. The $Z_4$ and $Z_3$ parameters are constrained by $\langle\phi\rangle_\Omega=0$, and the additional parameter can be absorbed into the interaction strength $g$. Typically $g$ is constrained by measuring a scattering amplitude.