Renormalization group theory, first-order phase transitions, perturbative calculations Consider a satistical field theory that is defined by the generic free energy (or action, in the case of dynamical theories)
\begin{equation}
   F[\phi] = \int_{k,t} a_k |\phi_k|^2 + \mathcal{V}(\phi),
\end{equation}
where $k$ is momentum and $\mathcal{V}$ contains all of the nonlinear interaction terms that respect the symmetries and constraints of the problem at hand. The probability of a certain configuration is given by $P[\phi] = Z^{-1} e^{-F[\phi]}$ where the partition function is given by $Z = \int \mathcal{D} [\phi] \, P[\phi]$.
In treatments of critical phenomena, e.g., using the Landau-Ginzburg free energy, one makes use of the scale-invariance caused by the lack of characteristic length/time scales in the system and looks for a theory for which $P[\phi]$ remains unchanged under a change of scale. Conceptually, this is the basis of renormalization group (RG) theory to find a theory with the scale symmetry, which is otherwise not a trivial symmetry to implement into the free energy.
Now imagine for a specific theory, all possible terms that can be included in $\mathcal{V}$ turn out to be irrelevant. This can, for example, be established through a naive (engineering) scaling analysis. Thus we know under RG transformation, the Gaussian theory (corresponding to the $|\phi_k|^2$ term in the free energy) is a stable fixed point. From what I understand, for this theory a second-order phase transition does not occur.
Questions
(1) Can we extract more information from the RG equation other than about the second-order phase transition? for instance, does the convergence of the RG flows to the Gaussian fixed point also discard the possibility of a first-order phase transition as well? (Bonus: can RG, in general, be used to detect first-order transitions?)
(2) Apart from phase transition, one might be interested in calculating various quantities/averages from this field theory. This can be done, for instance, by first obtaining the quantity of interest in the linearized theory, and then using the linear solution to obtain corrections coming from the nonlinear terms in an iterative manner.
My question is whether such a calculation should also be performed by keeping all the possible nonlinearities in $\mathcal{V}$, or only those that come directly from our microscopic theory?
From what I understand, we only add all possible terms to the free energy/Hamiltonian because under coarse-graining and RG they will be generated anyways. But now that we are just looking for perturbative expressions for various physical quantities (and not coarse-graining anymore), I am not sure if those added terms are actually required.
(3) For a quantity/averaged observable of interest, is there a way to identify whether they are universal (i.e., calculating them using the linear theory gives the correct long-distance asymptotic behaviour), or they are non-universal and depend on the microscopic details (including the irrelevant nonlinear terms in $\mathcal{V}$)?
PS: any hints or suggestions of resources are greatly appreciated.
 A: Regarding the premises: I am not sure you can claim with generality that there is no second order phase transition in a model solely from the fact that the Gaussian fixed point is a sink. You would need to know whether there are any other fixed points of the renormalization group transformation with unstable directions that take you from the critical manifold to the Gaussian sink. Keep in mind that 'relevant' and 'irrelevant' interactions are always with reference to a fixed point of the RG transformation, so in principle you could have a Gaussian fixed point at which all interactions are irrelevant but there is another fixed point at which some of the interactions are relevant. If such a critical point exists then you could have a second order transition even if the Gaussian fixed point is a sink. That said, off the top of my head I don't know of any examples in which the Gaussian fixed point is a sink in the first place.
For the other questions:
1) Yes, the RG can identify first order transitions. See, for example, Goldenfeld's "Lectures on Phase transitions and the Renormalization Group," Ch. 9 Sec. 9.6.6 (end) and 9.7. In particular, Goldenfeld states that if a fixed point has an eigenvalues of the RG transformation equal to $\ell^d$, where $\ell$ is the coarse-graining factor in real space and $d$ is the dimension, is indicative of a first order transition. Goldenfeld calls this the "Nienhuis-Nauenberg" criterion, referencing Niehuis and Nauenberg, Phys. Rev. Lett. 35, 477 (1975).
2) You do not need to use RG to calculate statistical moments, so you would calculate them using your microscopic model as given. For example, if you have a nearest-neighbor Ising with Hamiltonian $\mathcal H(\mathbf{s}) = - J_1 \sum_{\langle ij\rangle} s_i s_j + h \sum_i s_i$ you would calculate moments using exactly this form. If you include higher order interactions, like a next-nearest-neighbor term $J_2 \sum_{i,j \in {\rm n.n.n.}} s_i s_j$ then you will get a model with a different numerical set of means, correlations, etc. The next-nearest-neighbor terms being irrelevant at the ferromagnet phase transition is a statement about the interactions under rescaling of the model (the second part of the RG transformation, after coarse-graining), and shows that the ferromagnetic transition of the next-nearest-neighbor model is in the same universality class as the nearest-neighbor model, but this does not come into play in calculating the moments. For the statistical moments you are not rescaling anything, nor do you even need to be anywhere near a fixed point of the rescaled model in general, and so any terms in your Hamiltonian/action contribute to those moments. Thus, adding more interaction terms to your $F(\phi)$ would change the statistical moments of the model, though they may not alter the fixed point structure.
Edit: To elaborate a little bit to answer the follow-up question in the comments:

About the 2nd point: if we compute the statistical moments with the
original nonlinear terms as you mentioned, does this computation
correspond to a specific scale, or its applicability can be extended
to different (presumably larger) scales? I think I am still confused
about the relationship between computing averages in a theory with
some given nonlinearities, with the effect of RG on those nonlinear
terms which, possibly, will generate new nonlinearities.

Let's keep going with the Ising model as an example. Suppose you had an experimental system that were exactly described by a classical nearest neighbor 2d Ising model on a triangular lattice with lattice spacing $a$ and spin-spin coupling $J$ at temperature $T$. If you have an experimental probe with sufficient resolution that you can measure the orientation of individual spins separated only by a lattice spacing (equivalently momenta all the way from $0$ to $\Lambda = \pi/a$), then you have access to every scale of the problem. If you wanted to predict, e.g., the covariance between two spins located at lattice positions $\mathbf{x}_m$ and $\mathbf{x}_n$ then you would compute
$$\langle s_m s_n \rangle = Z^{-1} \sum_{\mathbf{s}} s_i s_j \exp\left(\frac{J}{k_B T} \sum_{\langle ij \rangle} s_i s_j \right).$$
There's no need to renormalize anything, no need to include any higher-order interactions here: the model is exact over all spatial scales (by assumption).
Now, let's suppose that your measurement device cannot resolve the orientation of spins all the way down to separations of $a$, but only down to separations of $R$, such that any spins separated by less than this spacing appear blurred together (essentially, coarse-grained together) into some local magnetization; let the local magnetization at location $\mathbf{x}$ be $M_R(\mathbf{x})$, where I have denoted dependence the resolution scale $R$. The underlying individual spins are still described by the nearest neighbor Ising model, but since the coarse-grained $M_R(\mathbf{x})$ values are the only thing you can measure, you might prefer to have a model for it. In principle, the way you would formally construct such a model would be to calculate
$$\exp(-\beta L[M_R(\mathbf{x}]) = Z^{-1}\sum_{\mathbf{s}} e^{-\beta \mathcal H(\mathbf{s})} \delta\left[\sum_{\mathbf{x}_m  \in \mathbf{x} + R} s_m - N_R(\mathbf{x}) M_R(\mathbf{x})  \right],$$
where $\delta$ is a delta function enforcing its argument, $\mathbf{x}_m  \in \mathbf{x} + R$ restricts the sum over spins whose position are within the radius $R$ of the position $\mathbf{x}$, $N_R(\mathbf{x})$ counts the number of spins within that radius, and $L[M_R(\mathbf{x})]$ is an effective Hamiltonian for the new coarse degrees of freedom. Why is it labeled $L$? Well, essentially this would be a formal derivation of the Landau free energy; this argument follows Ch. 5 in Goldenfeld's book referenced above. You may know that one does not really construct the Landau free energy this way in practice, but instead appeals to symmetry considerations to work out the terms that would be present in $L$. If we were to derive the Landau free energy this way, however, this is where the effective interactions come into play: those terms that we would include on symmetry considerations would be the terms generated by coarse-graining the Ising model as part of a renormalization group procedure. Now, typically when using the Landau free energy one is just interested in understanding phase transitions, so not all of those higher order interactions generated by this coarse graining procedure are relevant: we know that if we rescale the coarse-grained variables by appropriate powers of the coarse-graining length $R$ to put the model back onto the original scale then many of those higher order terms are unimportant as far as the critical properties of the model are concerned. However, if instead what you wanted to do is accurately calculate things like correlation functions of these coarse-grained spins, then you need those higher order terms if you want to be quantitatively accurate. This is because the following two calculations must agree:
$$\langle M_R(\mathbf{x}) M_R(\mathbf{y}) \rangle = \int \mathcal D M ~M_R(\mathbf{x}) M_R(\mathbf{y}) e^{-\beta L[M]},$$
where I've formally written this as a path integral, but one does not need to assume a continuum limit, and
$$\langle M_R(\mathbf{x}) M_R(\mathbf{y}) \rangle = \frac{1}{N_R(\mathbf{x}) N_R(\mathbf{y})}\sum_{m \in \mathbf{x} + R}\sum_{n \in \mathbf{y} + R} \left[\frac{1}{Z}\sum_{\mathbf{s}} s_m s_n e^{-\beta \mathcal{H}(\mathbf{s})} \right],$$
where I inserted the definition of $M$ in terms of the microscopic spins. If you drop the 'irrelevant' terms from $L$ these two calculations will not agree. As in the answer to another of your questions, the nonlinearities generated are necessary for these two calculations to agree. Dropping them won't change the critical properties of the model (around the critical point corresponding to the ferromagnetic transition), but it would affect the quantitative calculation of the statistical moments/cumulants. If you working with the rescaled variables, $M'_R(\mathbf{x}') \equiv R^{-[M]}M_R(\mathbf{x}/R)$, where $[M] = (2-d)/2$ is the scaling dimension of $M_R$, and you want to calculate, e.g., $\langle M'_R(\mathbf{x}') M'_R(\mathbf{y}')\rangle$, then those higher-order terms in $L$ will be suppressed by powers of $R^{-1}$, and can be neglected without affecting the value of $\langle M'_R(\mathbf{x}') M'_R(\mathbf{y}')\rangle$ too much (because those terms will affect the values of the moments on scales $|\mathbf{x}| \lesssim R$.
I think perhaps where some of the confusion may be coming from the situation in high energy physics, in which there are additional complications arising from the fact that one does not actually know the microscopic model, and so there are additional steps involved in safely taking the continuum limit. In this case, I believe the variables one is working with are more akin to the rescaled variables $M'_R(\mathbf{x}')$, though someone more familiar with the high-energy situation should chime in here.
3) The general rule of thumb is that critical exponents and scaling functions tend to be universal, as well as certain ratios of quantities sometimes, but I don't know of a comprehensive list (to the extent that one could even be compiled). Moreover, even critical exponents are not guaranteed to be universal: Baxter's "Exactly solved models in statistical mechanics" includes an example of the square lattice eight-vertex model (Ch. 10) in which most of the critical exponents vary continuously with the parameters of the model. (However, Baxter notes that others have suggested this might be an artifact of defining the critical exponents as powers of $T-T_c$, and that if one defines exponents as powers of the inverse correlation length $\xi^{-1}$, which essentially divides the normal exponents by the correlation length exponent $\nu$, then these modified exponents do appear to be universal). Another example in which microscopic details show up in critical exponents are networks: in models with phase transitions on networks with degree distributions decaying as $p(n) \sim n^{-\gamma}$ at large degree $n$, with $2 < \gamma < 3$ (such that the mean degree exists but the degree variance is infinite), then the critical exponents are functions of $\gamma$. See, e.g., this paper's Eq. 8 for one example of this result. (Although, in this case $\gamma$ is playing a role similar to the dimension of a lattice, so the exponents are in a sense still 'universal')
Edit: I struck out the last example about networks because thinking about it more the final parenthetical comment interpreting the degree distribution exponent as being similar to dimension is the better way to think about the situation, so it is not really an example of microscopic details affecting the critical exponents. A better example would be this paper by Duclut and Delamotte investigating erosion over intermediate spatial scales. They find that there is a line of fixed points in the renormalization group flow, such that the value of the anomalous exponent is somewhat dependent on the initial conditions of the RG flow (i.e., on details of the microscopic action).
