# The Renormalization Group for Large-Scale Structure: Primordial non-
Gaussianities

Charalampos Nikolis  Henrique Rubira  Fabian Schmidt

###### Abstract

The renormalization group for large-scale structure (RG-LSS) describes the
evolution of galaxy bias and stochastic parameters as a function of the cutoff
$\Lambda$. In this work, we introduce interaction vertices that describe
primordial non-Gaussianity into the Wilson-Polchinski framework, thereby
extending the free theory to the interacting case. The presence of these
interactions forces us to include new operators and bias coefficients to the
bias expansion to ensure closure under renormalization. We recover the
previously-derived “scale-dependent bias” contributions, as well as a new
(subdominant) stochastic contribution. We derive the renormalization group
equations governing the RG-LSS for a large class of interactions which account
for vertices at linear order in $f_{\rm NL}$ that parametrize interacting
scalar and massive spinning fields during inflation. Solving the RG equations,
we show the evolution of the non-Gaussian contributions to galaxy clustering
as a function of scale.

##  1 Introduction

Primordial non-Gaussianities (PNG) in the initial conditions of the Universe
are a direct probe of distinct inflationary scenarios [1, 2, 3, 4, 5, 6, 7].
Observations from the cosmic microwave background (CMB) have put strong
constraints on the magnitude of the non-Gaussianities [8]. Interest in PNG has
been renewed due to recent progress in calculating cosmological correlators
via the “cosmological bootstrap” program [9] (see [10] and references
therein).

The large-scale structure (LSS) of the Universe provides an important probe to
search for PNG on complementary scales to those probed by the CMB. Based on
the framework of the effective field theory of LSS (EFTofLSS) [11, 12, 13, 14,
15, 16] (for applications of PNG in the context of EFTofLSS see [17, 18, 19,
20, 21, 22, 23, 24, 25, 26, 27]), we can expand the galaxy density as a sum of
all operators $O$ allowed by the symmetries of non-relativistic General
Relativity, and assign to each of them a bias (or Wilson) coefficient $b_{O}$
[28]

$\delta_{g}(\bm{k},\tau)=\sum_{O}b_{O}^{\Lambda}(\tau)O[\delta^{(1)}_{\Lambda}](\bm{k},\tau)\,.$
(1.1)

This relation implies that physical processes inherently non-local in time
such as galaxy formation [29, 30] can be factorized in the time-dependent bias
coefficients. This is made possible by properties of the perturbation theory
for matter and the gravitational potential at each order [31, 32].

A crucial point of the bias expansion Eq. (1.1) is that it contains composite
operators, meaning that it suffers from small-scale Fourier modes back-
reacting onto large scales [33, 34, 35]. Therefore, one has to smooth it with
a cutoff $\Lambda$ that regularizes the theory. A central piece of the
EFTofLSS and the bias expansion above is that it allows us to treat the dark
matter and galaxy density perturbatively, while allowing for a consistent way
to parametrize non-linear non-perturbative modes in the theory by imposing
this cutoff $\Lambda$ that smooths out small-scale modes.

In the presence of PNG, the bias expansion needs a non-trivial generalisation
[36, 37, 38, 39]. The reason is that the horizon-scale dynamics during
inflation that generated PNG correspond to large-scale correlations in the
late universe that cannot be generated by subhorizon gravitational or non-
gravitational forces. Since the statistics of galaxies which we want to
describe via the EFT involve much larger scales than the scales that govern
local galaxy formation, the additional galaxy bias terms depend on specific
kinematic regimes of primordial interactions, in particular those where one of
the wavelengths is soft compared to the others, commonly called the squeezed
limit. The unique new bias operators offer the opportunity to robustly detect
signatures from different models of primordial interactions. However, these
signatures are multiplied by additional Wilson coefficients $b_{O_{\rm NG}}$.
Therefore, there is strong motivation to put independent constraints on the
values of these coefficients, as they are largely degenerate with the
amplitude of PNG [38, 40]. Moreover, the values of these bias coefficients are
found to depend sensitively on the tracer population considerered [41, 42,
43], at least for the leading bias parameter introduced by local-type PNG. Our
goal here in particular is to derive how the $b_{O_{\rm NG}}$ run from one
scale $\Lambda$ to another. While this does not tell us about the fully non-
perturbative physics setting the absolute scale of the bias coefficients, it
can be compared with measurements of the bias coefficients as a function of
scale (for example in simulations), and thus provides another validation test
of the theory.

Different renormalization schemes have been proposed for the bias expansion.
The most used option so far is to renormalize the bias operators at the
$n$-point function level [34, 35], defining the (renormalized) bias operators
deep in the IR via the corresponding $n$-point function for which they appear
at tree level. In this work, we follow the path-integral approach of [13, 44,
45], keeping the cutoff $\Lambda$ explicitly. We generalize the
renormalization of bias parameters of [44, 45] to non-Gaussian initial
conditions. We use the Wilson-Polchinski approach to derive the
renormalization group (RG) equations for the bias parameters including PNG.
The RG flow approach in principle allows for the resummation of Feynman
diagrams, and therefore there are some grounds to expect that additional
information can be extracted mainly when adding PNG interactions.

Our Renormalization Group Equations (RGE) describe a large class of cubic
interactions during inflation. The equations are general enough to depend only
on the characteristic squeezed limit of the non-Gaussian interaction $\sim
k^{\Delta}$. This already enables us to calculate the RG flow for different
known isotropic (spin zero) and non-isotropic (spin 2) shapes of PNG. More
importantly, novel calculations for bispectra (e.g from the bootstrap
programme) could lead to shapes of non-Gaussianity with different
characteristic soft scalings than previously known. These shapes can directly
be implemented in the RG-flow we present in this work. Thus, the effect of new
PNG shapes in the galaxy statistics can easily be retrieved from our general
RG equations.111For other work on the different regularization schemes for PNG
see [46].

The outline of the paper is the following. In the rest of this section we set
the notation. In Sec. 2 we introduce PNG, generalize the operator basis to
include them and derive the non-Gaussian partition function. In Sec. 3 we
derive the running of the bias parameters via the Wilson-Polchinski framework.
We present the results in Sec. 4 and conclude in Sec. 5. We have four
appendices: App. A is dedicated to an alternative formulation of PNG in the
partition function, coupling them directly to the current $J$; App. B for
examples on how to derive $n$-point functions from the partition function;
App. C discussing the interacting PNG vertex and App. D for explicit
evaluations of the bias renormalization.

For numerical results we use the transfer function as well as the linear
matter power spectrum computed by CLASS [47] with Planck 2018 Euclidian
$\Lambda$CDM cosmology [48].

### Notation and conventions

We review now the notation adopted in [44].

#### Momentum conventions.

We use $\bm{k},\bm{p}$ for momenta variables where with bold letters we denote
Euclidean vectors. We adopt a shorthand notation $\bm{p}_{1\dots
n}=\bm{p}_{1}+\dots+\bm{p}_{n}$ for the sum of vectors. The momentum-space
integrals are written as

$\int_{\bm{p}_{1},\dots,\bm{p}_{n}}=\int\frac{d^{3}p_{1}}{(2\pi)^{3}}\dots\int\frac{d^{3}p_{n}}{(2\pi)^{3}}\,.$
(1.2)

We keep the variable $\bm{q}$ for the Lagrangian position.

#### Cutoff filters.

We can coarse-grain a field $f$ by applying a sharp filter $W$ on a length
scale $1/\Lambda$. In Fourier space this is simply a product

$f_{\Lambda}(\bm{k})=\,W_{\Lambda}(\bm{k})f(\bm{k})\,.$ (1.3)

#### The kernels.

We consider three types of kernels and vertex diagrams (beyond the non-
Gaussian kernel presented below). First, we write the operators $O$ of the
bias expansion Eq. (1.1) in terms of the dark matter overdensity field
$\delta(\bm{p}_{i})$ as

$O[\delta^{(1)}_{\Lambda}](\bm{k})=\int_{\bm{p}_{1},...,\bm{p}_{n}}\hat{\delta}_{\rm
D}(\bm{k}-\bm{p}_{12\dots
n})S_{O}(\bm{p}_{1},...,\bm{p}_{n})\delta(\bm{p}_{1})...\delta(\bm{p}_{n})\,,$
(1.4)

with $\hat{\delta}_{\rm D}=(2\pi)^{3}\delta_{\rm D}$ for the Dirac delta
function. These operators can be classified according to their perturbative
order, so in this case, we have $O=O^{[n]}$. Second, we expand the matter
overdensity with the usual Eulerian kernels $F_{n}$ [49]

$\displaystyle\delta[\delta^{(1)}_{\Lambda}](\bm{k})=\sum_{n=1}^{\infty}\int_{\bm{p}_{1},...,\bm{p}_{n}}\hat{\delta}_{\rm
D}(\bm{k}-\bm{p}_{12..n})F_{n}(\bm{p}_{1},...,\bm{p}_{n})\delta^{(1)}_{\Lambda}(\bm{p}_{1})...\delta^{(1)}_{\Lambda}(\bm{p}_{n})\;,$
(1.5)

where $\delta^{(1)}_{\Lambda}$ is the smoothed linear density field. The third
type of vertex $K$ represents a union of the former two

$O(\bm{k})=\sum_{\ell={\rm
order}(O)}^{\infty}\int_{\bm{p}_{1},\ldots,\bm{p}_{\ell}}\hat{\delta}_{\rm
D}(\bm{k}-\bm{p}_{1\ldots\ell})K_{O}^{(\ell)}(\bm{p}_{1},\ldots\bm{p}_{\ell})\delta^{(1)}(\bm{p}_{1})\cdots\delta^{(1)}(\bm{p}_{\ell})\,,$
(1.6)

in which we expanded the operators in terms of linear legs. Notice it is
possible to derive the $K^{\ell}_{O}$ by inserting Eq. (1.5) into Eq. (1.4).
Moreover, we define $K_{O}=\sum_{\ell={\rm order}(O)}^{\infty}K^{(\ell)}_{O}$
and we assume $K^{(\ell)}_{O}$ to be symmetrized in its arguments. The three
types of vertices are represented by

$\raisebox{0.0pt}{\parbox[c]{0.0pt}{\includegraphics[scale={1}]{figs/diags/diag_opvertex.pdf}}}\raisebox{0.0pt}{\parbox[c]{0.0pt}{\includegraphics[scale={1}]{figs/diags/diag_Fvertex.pdf}}}\raisebox{0.0pt}{\parbox[c]{0.0pt}{\includegraphics[scale={1}]{figs/diags/diag_Kvertex.pdf}}}\,,\quad$
(1.7)

with the differentiation between large boxes and the small boxes, and also the
dashed lines representing linear legs.

#### Power spectra and field variance.

From the linear matter overdensity, we can define its power spectrum (i.e.,
the linear propagator)

$P_{\rm
L}^{\Lambda}(k)=\left<\delta^{(1)}_{\Lambda}(\bm{k})\delta^{(1)}_{\Lambda}(\bm{k}^{\prime})\right>^{\prime}=\raisebox{0.0pt}{\parbox[c]{0.0pt}{\includegraphics[scale={1.0}]{figs/diags/diag_Plin_cut}}}\,,$
(1.8)

where the prime notation in $\langle\dots\rangle^{\prime}$ indicates that the
momentum-conservation factor $\hat{\delta}_{\rm D}$ has been dropped. Notice
that we imposed momentum restrictions on the propagators, cutting them off at
$\Lambda$. If we shift the cutoff to $\Lambda^{\prime}=\Lambda+\lambda$, we
can separate the contributions into one with support up to wavenumber
$\Lambda$ and one with support in the momentum shell
$k\in[\Lambda,\Lambda+\lambda]$,

$\displaystyle P_{\rm shell}(k)$ $\displaystyle=\left<\delta^{(1)}_{\rm
shell}(\bm{k})\delta^{(1)}_{\rm
shell}(\bm{k}^{\prime})\right>^{\prime}=\raisebox{0.0pt}{\parbox[c]{0.0pt}{\includegraphics[scale={1.0}]{figs/diags/diag_Plin_shell}}}\,.$
(1.9)

Furthermore, from the matter linear propagator, we can calculate the matter
overdensity variance

$\displaystyle\sigma^{2}_{\Lambda}=\int_{\bm{p}}^{\Lambda}\left<\delta^{(1)}_{\Lambda}(\bm{p})\delta^{(1)}_{\Lambda}(\bm{p}^{\prime})\right>^{\prime}=\int_{\bm{p}}^{\Lambda}P_{\rm
L}^{\Lambda}(p)\,.$ (1.10)

#### Gaussian partition function.

For the case of Gaussian initial conditions we can generate all $n$-point
functions of matter ($m$) and galaxies ($g$) from a partition function

$\displaystyle\mathcal{Z}_{\rm
G}[J_{m,\Lambda},J_{g,\Lambda}]=\int\mathcal{D}\delta^{(1)}_{\Lambda}\mathcal{P}[\delta^{(1)}_{\Lambda}]\exp\left\\{{\int_{\bm{k}}J_{m,\Lambda}(\bm{k})\delta[\delta^{(1)}_{\Lambda}](-\bm{k})+\int_{\bm{k}}J_{g,\Lambda}(\bm{k})\delta_{g}[\delta^{(1)}_{\Lambda}](-\bm{k})}\right\\}\;,$
(1.11)

with the full nonlinear solutions for the dark matter density field and the
galaxy density fields coupled to their respective currents
$J_{m,\Lambda},\;J_{g,\Lambda}$. This partition function generates the matter,
galaxy and cross $n$-point correlations (see discussion in App. B). Here we
have a path-integral over all initial configurations of the linear overdensity
$\delta^{(1)}_{\Lambda}$ with support over wavenumbers up to $\Lambda$.
$\mathcal{P}[\delta^{(1)}_{\Lambda}]$ refers to the Gaussian measure:

$\mathcal{P}[\delta^{(1)}_{\Lambda}]=\left(\prod_{\bm{k}}^{\Lambda}2\pi P_{\rm
L,\Lambda}(k)\right)^{-1/2}\exp{\left[-\frac{1}{2}\int_{\bm{k}}^{\Lambda}\frac{|\delta^{(1)}_{\Lambda}|^{2}}{P_{\rm
L,\Lambda}(k)}\right]}\;,$ (1.12)

and $\delta[\delta^{(1)}_{\Lambda}]$ is the perturbative expansion of the
dark-matter Eq. (1.5). For other constructions of the large-scale structure
partition function, see [50, 51, 52, 53, 54].

This enables us to reproduce all the $n$-point statistics for the matter and
the galaxies by taking derivatives with respect to the matter or galaxy
currents $J_{m}\;,J_{g}$ (see e.g. [54]). For example, the Fourier-space
$m$-point correlation function of a given biased tracer is given by

$\displaystyle\langle\delta_{g}(\bm{k}_{1})\dots\delta_{g}(\bm{k}_{m})\rangle$
$\displaystyle=\int\mathcal{D}\delta^{(1)}_{\Lambda}\mathcal{P}[\delta^{(1)}_{\Lambda}]\,\delta_{g}(\bm{k}_{1})\dots\delta_{g}(\bm{k}_{m})\,,$
(1.13)

which can be obtained by taking derivatives of the partition function w.r.t.
$J_{g,\Lambda}$

$\langle\delta_{g}(\bm{k}_{1})\dots\delta_{g}(\bm{k}_{m})\rangle=\frac{1}{\mathcal{Z}[0]}\frac{\delta^{m}\mathcal{Z}}{\delta
J_{g,\Lambda}(\bm{k}_{1})\dots\delta
J_{g,\Lambda}(\bm{k}_{m})}\Big{|}_{J_{m,\Lambda}=0,J_{g,\Lambda}=0}\;,$ (1.14)

and we can also have cross-correlations, like the matter-galaxy cross power
spectrum

$\langle\delta(\bm{k}_{1})\delta_{g}(\bm{k}_{2})\rangle=\frac{1}{\mathcal{Z}[0]}\frac{\delta^{2}\mathcal{Z}}{\delta
J_{m,\Lambda}(\bm{k}_{1})\delta
J_{g,\Lambda}(\bm{k}_{2})}\Big{|}_{J_{m,\Lambda}=0,J_{g,\Lambda}=0}\,.$ (1.15)

We discuss some basic examples in App. B.

#### List of Gaussian operators.

The basis of of Gaussian operators is given by

First order: $\displaystyle\quad\delta;$ Second order:
$\displaystyle\quad\delta^{2},\,\mathcal{G}_{2};\;$ Third order:
$\displaystyle\quad\delta^{3},\;\delta\mathcal{G}_{2},\;\Gamma_{3},\;\mathcal{G}_{3}\,,$
(1.16)

with the Galileon operators defined as

$\displaystyle\mathcal{G}_{2}(\Phi_{g})$ $\displaystyle\equiv$
$\displaystyle(\nabla_{i}\nabla_{j}\Phi_{g})^{2}-(\nabla^{2}\Phi_{g})^{2}\,,$
(1.17) $\displaystyle\mathcal{G}_{3}(\Phi_{g})$ $\displaystyle\equiv$
$\displaystyle-\frac{1}{2}\left[2\nabla_{i}\nabla_{j}\Phi_{g}\nabla^{j}\nabla_{k}\Phi_{g}\nabla^{k}\nabla^{i}\Phi_{g}+(\nabla^{2}\Phi_{g})^{3}-3(\nabla_{i}\nabla_{j}\Phi_{g})^{2}\nabla^{2}\Phi_{g}\right]\,,$
(1.18)

with $\Phi_{g}\equiv\nabla^{-2}\delta$ being the scaled gravitational
potential. We also have

$\displaystyle\Gamma_{3}(\Phi_{g},\Phi_{v})\equiv\mathcal{G}_{2}(\Phi_{g})-\mathcal{G}_{2}(\Phi_{v})\,,$
(1.19)

with $\Phi_{v}\equiv\nabla^{-2}\bm{\nabla}\cdot\bm{v}$, the velocity
potential. Moreover, we define as

$\sigma^{2}_{\bm{k}_{1},\bm{k}_{2}}=\left(\bm{k}_{1}\cdot\bm{k}_{2}/k_{1}k_{2}\right)^{2}-1\,,$
(1.20)

the Fourier-transformed coefficient of $\mathcal{G}_{2}$, such that in the
notation of Eq. (1.4), we can clearly see that
$\sigma^{2}_{\bm{k}_{1},\bm{k}_{2}}\equiv S_{\mathcal{G}_{2}}$. Finally, we
normalize the operators as

$O(\bm{k})\to O(\bm{k})-\langle O(\bm{k}=0)\rangle\,,$ (1.21)

to avoid tadpole contributions.

##  2 The generalized galaxy-bias partition function with non-Gaussianity

In this section we start by discussing the connection of the partition
function formalism for the bias expansion to the usual textbook interactions
in quantum field theory (QFT). Next, we present how to introduce non-
Gaussianities in cosmology and also fix the notation for the new interacting
vertex. Later we describe how to modify the galaxy bias expansion to include
PNG. Finally, we derive the partition function for the galaxy bias including
non-Gaussianities.

###  2.1 Parallelism with QFT

We begin with a QFT example, that will serve as basis for our later analysis.
Consider a real scalar field with a cubic interaction of the type

$\displaystyle\mathcal{L}=-\frac{1}{2}\phi(\square+m^{2})\phi-\frac{g}{3!}\phi^{3}=\mathcal{L}_{0}+\mathcal{L}_{\rm
int}\,.$ (2.1)

We define the free partition function as

$\displaystyle\mathcal{Z}_{0}[J]$
$\displaystyle=\int\mathcal{D}\phi\,e^{iS_{0}[\phi]+i\int
d^{4}xJ\phi}\equiv\int\mathcal{D}\phi\mathcal{P}[\phi]\,e^{i\int
d^{4}xJ\phi}\,,$ (2.2)

where $\mathcal{P}[\phi]$ denotes the free (Gaussian) measure of the path-
integral. The textbook approach of handling QFT interactions is to set up a
perturbative expansion of the interaction vertices, and find quantum
corrections to the free theory, that will induce deviations from the free-
Gaussian case

$\displaystyle\mathcal{Z}[J]$
$\displaystyle=\int\mathcal{D}\phi\,e^{iS_{0}[\phi]+iS_{\rm int}[\phi]+i\int
d^{4}xJ\phi}=\int\mathcal{D}\phi\mathcal{P}[\phi]\,e^{+i\int d^{4}xJ\phi}$
$\displaystyle\quad\times\left(1-\frac{ig}{3!}\int_{\bm{x}_{1}}\phi^{3}(\bm{x}_{1})+\frac{1}{2}\left(\frac{-ig}{3!}\right)^{2}\int_{\bm{x}_{1},\bm{x}_{2}}\phi^{3}(\bm{x}_{1})\phi^{3}(\bm{x}_{2})+\dots\right)\,.$
(2.3)

An important difference between the (usual) QFT case and the EFTofLSS, to be
explored below, is the canonical linear coupling of the current to the field.

The $n$-point function for the interacting theory are calculated as the
perturbative expansion of the interaction exponential

$\displaystyle\langle\phi(\bm{k}_{1})\dots\phi(\bm{k}_{n})\rangle=\frac{1}{\mathcal{Z}[0]}\frac{\delta^{n}\mathcal{Z}}{\delta
J_{i}\dots\delta
J_{n}}\Big{|}_{J=0}=\int\mathcal{D}\phi\mathcal{P}[\phi]\,\phi(\bm{k}_{1})\dots\phi(\bm{k}_{n})\,e^{iS_{\rm
int}[\phi]}$
$\displaystyle=\int\mathcal{D}\phi\mathcal{P}[\phi]\,\phi(\bm{k}_{1})\dots\phi(\bm{k}_{n})\left(1+iS_{\rm
int}[\phi]+\frac{i^{2}}{2}S_{\rm int}[\phi]^{2}\dots\right)\,.$ (2.4)

Corrections from the Gaussian free field are sourced by the interactions in
the Lagrangian, or in other words by the $J^{0}$ term in the partition
function. Moreover, a linear coupling of the current with the field is added
as a tool to calculate $n$-point interactions.

An alternative way of formulating interactions is to couple the current to
nonlinear operators of the (free or interacting) field. These are termed
composite operators. As we show now, in general one can switch between the two
formulations via a nonlinear field redefinition, where we restrict ourselves
to linear order in $g$. The classical equation of motion

$\displaystyle\left(\square+m^{2}\right)\phi=-\frac{g}{2}\phi^{2}\,,$ (2.5)

admits a formal solution of the form

$\displaystyle\phi_{\rm cl}=\phi-i\frac{g}{2}\int
d^{4}x^{\prime}G\left(x,x^{\prime}\right)\phi_{\rm
cl}^{2}(x^{\prime})=\phi-i\frac{g}{2}\int
d^{4}x^{\prime}G\left(x,x^{\prime}\right)\phi^{2}(x^{\prime})+\mathcal{O}(g^{2})\,,$
(2.6)

with $G(x,x^{\prime})$ the Green’s function of the free equation of motion,

$\displaystyle\left(\square+m^{2}\right)G(x,x^{\prime})=-i\delta_{\rm
D}(x-x^{\prime})\,,$ (2.7)

and $\phi$ the free-field solution. This motivates us to do the following
field redefinition

$\displaystyle\tilde{\phi}=\phi+i\frac{g}{3!}\int
d^{4}x^{\prime}G\left(x,x^{\prime}\right)\phi^{2}(x^{\prime})\,,$ (2.8)

then

$\displaystyle\phi=\tilde{\phi}-i\frac{g}{3!}\int
d^{4}x^{\prime}G\left(x,x^{\prime}\right)\tilde{\phi}^{2}(x^{\prime})+\mathcal{O}(g^{2})\,.$
(2.9)

Plugging it back to the Lagrangian we get

$\displaystyle\mathcal{L}=-\frac{1}{2}\tilde{\phi}\left(\square+m^{2}\right)\tilde{\phi}+J\tilde{\phi}\,,$
(2.10)

where the interaction term has been cancelled and the current is coupled to
the full solution of the field. The composite-operator formulation is useful
when no fundamental interaction Lagrangian is known, but one can still write
down the equations of motion.

For the case of LSS, the tracer field $\delta_{g}[\delta^{(1)}]$ in Eq. (1.1)
is the main object from which we want to calculate correlators and therefore
have a partition function for. We differentiate between two types of
interactions that appear for matter tracers when constructing the partition
function: the interactions present in the initial conditions due to inflation
[imprinted in $\delta^{(1)}$ as e.g. Eq. (2.14) below], and the non-
linearities introduced by the gravitational evolution [both in the operator
definition Eq. (1.4) and in the $F$ kernels of Eq. (1.5)].

For the first type of interactions, we can directly formulate the inflationary
initial conditions as nonlocal interactions of $\delta^{(1)}$ [see e.g. Eq.
(2.14) below]. Therefore, it is convenient to include those interactions as
new vertices in the $J$-independent (interaction) term, following the
canonical steps outlined above leading to Eq. (2.1).

For the interactions from nonlinear evolution, on the other hand, the
composite-operator formulation is useful since the forward evolution of the
density field can be described at the level of equations of motion via the
EFTofLSS framework (whereas no fundamental Lagrangian is known for this
problem). We then couple the current to the (classical) solution of the
equations of motion.222In general, the renormalization of composite operators
requires contributions to the effective action that are higher-order in the
current [55, 56]; for a general treatment of composite operators and
renormalization methods see [57], and for the application to LSS see [45]. If
these have a local form, such as $\int_{\bm{x}}J^{2}(\bm{x})\cdots$, they are
known as contact terms. In the case of the partition function for LSS, terms
coupled to higher powers of the current are in one-to-one correspondence with
stochastic (noise) contributions to tracer $n$-point functions [54, 58, 45].

In summary, for constructing the partition function we add the non-linearities
from non-Gaussian initial condition in the interacting $J$-independent term,
while the non-linearities from the later evolution couple to interacting
$J$-dependent terms.

As a final note, one can also consider a formulation where both types of
interactions are represented in terms of composite operators; that is, moving
the non-Gaussian initial conditions from the Lagrangian to the current
coupling via a field redefinition such as Eq. (2.8). We provide one example
for this in App. A and App. C.3. In fact, such reformulation in general offers
independent physical insights.

###  2.2 Primordial non-Gaussianities

We now introduce non-Gaussianities in the context of cosmology. In the usual
approach [39, 35], weakly non-Gaussian initial conditions are defined by
expanding around a Gaussian field $\varphi_{G}$,

$\displaystyle\varphi(\bm{k})$ $\displaystyle={\varphi_{\rm G}}(\bm{k})+f_{\rm
NL}\int_{\bm{p}_{1},\bm{p}_{2}}\hat{\delta}_{\rm D}(\bm{k}-\bm{p}_{12})K_{\rm
NL}(\bm{p}_{1},\bm{p}_{2}){\varphi_{\rm G}}(\bm{p}_{1}){\varphi_{\rm
G}}(\bm{p}_{2})+\mathcal{O}[f_{\rm NL}^{2},g_{\rm NL}]\,,$ (2.11)

with $K_{\rm NL}$ the PNG kernel that depends on the structure of interactions
during inflation. Here, $\varphi$ is the superhorizon gravitational potential
during matter domination, which is directly related to the curvature
perturbation in comoving gauge $\mathcal{R}$ by $\varphi=(3/5)\mathcal{R}$.
This leads to a non-vanishing bispectrum for the primordial fluctuations,

$\displaystyle\left<\varphi(\bm{k}_{1})\varphi(\bm{k}_{2})\varphi(\bm{k}_{3})\right>^{\prime}\equiv
B_{\varphi}\left(k_{1},k_{2},k_{3}\right)=2f_{\rm NL}K_{\rm
NL}(\bm{k}_{1},\bm{k}_{2})P_{\varphi}(k_{1})P_{\varphi}(k_{2})+2\,\text{perms}\,,$
(2.12)

where two cyclic permutations of the three momenta $k_{1},k_{2},k_{3}$ shall
be considered. An important limit for galaxy clustering is the squeezed limit
defined as

$\displaystyle
B_{\varphi}\left(k_{\ell},|\bm{k}_{s}-\frac{1}{2}\bm{k}_{\ell}|,|\bm{k}_{s}+\frac{1}{2}\bm{k}_{\ell}|\right)\xrightarrow{k_{\ell}\ll
k_{s}}$ $\displaystyle\hskip 56.9055pt2f_{\rm NL}\Big{[}K_{\rm
NL}(\bm{k}_{s},\bm{k}_{\ell})P_{\varphi}(k_{s})P_{\varphi}(k_{\ell})+K_{\rm
NL}(-\bm{k}_{s},\bm{k}_{\ell})P_{\varphi}(k_{s})P_{\varphi}(k_{\ell})\Big{]}$
$\displaystyle\hskip 56.9055pt+2f_{\rm NL}K_{\rm
NL}(-\bm{k}_{s},\bm{k}_{s})\left[P_{\varphi}(k_{s})\right]^{2}+\mathcal{O}\left(\frac{k^{2}_{\ell}}{k^{2}_{s}}\right)\,.$
(2.13)

The first two terms in the RHS are what we usually refer to as the squeezed
limit of the kernel, which can be uniquely determined by the bispectrum. We
refer to those terms as the scattering interaction, for reasons that will
become clear in Sec. 3.2. The last term is what we refer to as the
annihilation interaction which will contribute to the stochastic bias
operators.

We can immediately translate the relation Eq. (2.11) to one involving the
linearly evolved matter overdensity

$\boxed{\begin{aligned} \delta^{(1)}(\bm{k})&=\delta^{(1)}_{\rm
G}(\bm{k})+f_{\rm NL}\int_{\bm{p}_{1},\bm{p}_{2}}\hat{\delta}_{\rm
D}(\bm{k}-\bm{p}_{12})K_{\rm
NL}(\bm{p}_{1},\bm{p}_{2})\frac{M(|\bm{p}_{1}+\bm{p}_{2}|)}{M(p_{1})M(p_{2})}\delta^{(1)}_{\rm
G}(\bm{p}_{1})\delta^{(1)}_{\rm G}(\bm{p}_{2})\\\ &\hskip
142.26378pt+\mathcal{O}[f_{\rm NL}^{2},g_{\rm NL}],\end{aligned}}$ (2.14)

with $g_{\rm NL}$ the quartic interaction and

$\delta^{(1)}(\bm{k},\tau)=M(k,\tau)\varphi(\bm{k})\,,$ (2.15)

and

$M(k,\tau)=\frac{2}{3}\frac{k^{2}T(k)D(\tau)}{H_{0}^{2}\Omega_{m,0}}\,,$
(2.16)

where $T(k)$ is the transfer function describing linear evolution through
horizon crossing, normalized such that $T(k)\to 1$ for $k\ll k_{\rm eq}$,
$H_{0}$ is the value of the Hubble constant today, $\Omega_{m,0}$ the
fractional matter density and $D$ the growth factor. For simplicity, we have
dropped in Eq. (2.14) the time argument $\tau$. For the rest of the text it
should be kept in mind that the transfer function should be evaluated at the
time of observation $\tau$. Finally, $\delta^{(1)}_{\rm G}$ is the Gaussian
overdensity acquired from acting with $M(k,\tau)$ on the Gaussian counterpart
of the primordial fluctuation $\varphi_{\rm G}$.

Assuming that the initial conditions are homogeneous and isotropic, we can
deduce that $K_{\rm NL}(\bm{k}_{1},\bm{k}_{2})$ depends only on
$k_{1}=|\bm{k}_{1}|$, $k_{2}=|\bm{k}_{2}|$ and the angle between the two
vectors $\hat{\bm{k}}_{1}\cdot\hat{\bm{k}}_{2}$ [39]. This allows us to
express the squeezed limit333A point worth mentioning is that the scaling of
the kernel $K_{\rm NL}$ in a regime different than the squeezed limit (such as
for the annihilation interaction) is ambiguous and depends on the complete
form of the kernel. Since we cannot completely define this form from the
bispectrum alone, studies investigating the kernel have reversed-engineered an
expression [39, 59] and find that this term contributes with a non trivial
scaling for momentum regions different than the squeezed limit even for the
equilateral case, while preserving the squeezed limit behaviour as expected.
($k_{l}\ll k_{s}$) of the kernel with Legendre polynomials of even order
$\mathcal{P}_{L}$. If we make the further assumption of a scale invariant
bispectrum, we expect only dimensionless quantities in the expansion, and we
have [35]

$K_{\rm NL}(\bm{k}_{s},\bm{k}_{l})\xrightarrow{k_{l}\ll
k_{s}}\sum_{L,i}a_{L,i}\left(\frac{k_{l}}{k_{s}}\right)^{\Delta_{i}}\mathcal{P}_{L}\left(\hat{\bm{k}}_{s}\cdot\hat{\bm{k}}_{l}\right)\,,$
(2.17)

where $i$ refers to the different non-Gaussianities, $L$ to the rank (spin) of
the multipole expansion, and $a_{L,i}$ are dimensionless constants. In the
following, we always focus on one fixed $\Delta_{i}$ and hence keep a single
coefficient $a_{L}$ for the amplitude, which moreover we fix to $a_{L}=1$ in
our numerical results. This is a way to express a class of different micro-
physical models444Even though it cannot account for recent non-perturbative
models of non-Gaussianity such as [60]. such as:

  * •

Spin-0 non-Gaussianity ($L=0$): In that case, we have $\mathcal{P}_{0}=1$ and
we can recover the well known local $\Delta=0$ [61], and the equilateral
($\Delta=2$) cases [62, 6]. In different inflationary scenarios we can have
intermediate values of $\Delta$ [63, 64].

  * •

Higher-spin non-Gaussianity ($L=2,4\dots$): This type of non-Gaussianity can
be a smoking gun of massive spinning fields during inflation (see e.g [65, 66,
67]).

###  2.3 The non-Gaussian bias basis

The non-Gaussianity in the initial conditions will lead to new terms appearing
in the bias expansion, as discussed in Sec. 3 and was originally pointed out
for local non-Gaussianity in [36, 68], and for non-local shapes in [39]. In
this part we discuss modification in the basis of the operators to account for
those new terms. While we appear to introduce these terms ad-hoc here, we will
later show that it is precisely these terms that are generated by the RG flow
in the presence of PNG.

#### Spin-0.

Assuming we have only one source of non-Gaussianity ($i=1$), for the spin-0
($L=0$) case we add to the bias expansion a linear operator $\Psi$ with the
following form

$a_{0}f_{\rm NL}\Psi_{\Lambda}(\bm{k})=a_{0}f_{\rm NL}\left(\frac{k}{k_{\rm
NG}}\right)^{\Delta}\frac{\delta^{(1)}_{\Lambda}(\bm{k})}{M(k)}\equiv
a_{0}f_{\rm NL}\left(\frac{k}{k_{\rm
NG}}\right)^{\Delta}\varphi_{\Lambda}(\bm{k})\;,$ (2.18)

with $k_{\rm NG}$ an arbitrary momentum scale to make the operator
dimensionless, so we can keep the bias parameters as dimensionless constants.
If we choose $k_{\rm NG}\sim k_{*}=R_{*}^{-1}$ to be of the order of the
characteristic spatial length scale governing the formation of the tracer
considered, then we expect the associated bias coefficients to be of order 1
[28] at the scale $\Lambda\to 0$. Different choices of $k_{\rm NG}$ can always
be absorbed in a rescaling of the bias coefficients. For the local
($\Delta=0$) spin-zero type of non-Gaussianities, we then find

$\Psi_{\Lambda}(\bm{k})=\varphi_{\Lambda}(\bm{k})\;,$ (2.19)

which is usually referred in the literature with the notation $b_{\Psi}\to
b_{\varphi}$ (or $b_{\phi}$). For any $\Delta$, the new basis will consist of
operators of the form $\left\\{O_{\rm G},O_{\rm G}\Psi\right\\}$, where
$O_{\rm G}$ refers to the Gaussian basis Eq. (1) and $O_{\rm NG}=O_{\rm
G}\Psi$ are the new non-Gaussian operators added, with respective bias
coefficients $b_{O_{\rm NG}}$ (recall that we are working to linear order in
$f_{\rm NL}$). We have now an expanded list of operators due to non-
Gaussianity in the bias expansion:

First order: $\displaystyle\quad\delta,\,\Psi;$ Second order:
$\displaystyle\quad\delta^{2},\,\mathcal{G}_{2},\,\delta\Psi;\;$ Third order:
$\displaystyle\quad\delta^{3},\,\delta^{2}\Psi,\;\delta\mathcal{G}_{2},\;\Psi\mathcal{G}_{2},\;\Gamma_{3},\;\mathcal{G}_{3}\,.$
(2.20)

For example, the first few terms can be written as [28]

$\displaystyle\delta_{g}[\bm{k}]$
$\displaystyle=b_{\delta}^{\Lambda}\delta[\delta^{(1)}_{\Lambda}](\bm{k})+a_{0}f_{\rm
NL}b_{\Psi}^{\Lambda}\Psi_{\Lambda}[\delta^{(1)}_{\Lambda}](\bm{k})$
$\displaystyle\quad+b_{\delta^{2}}^{\Lambda}\delta^{2}[\delta^{(1)}_{\Lambda}](\bm{k})+b_{\mathcal{G}_{2}}^{\Lambda}\mathcal{G}_{2}[\delta^{(1)}_{\Lambda}](\bm{k})+a_{0}f_{\rm
NL}b_{\Psi\delta}^{\Lambda}\Psi_{\Lambda}[\delta^{(1)}_{\Lambda}](\bm{k})\delta[\delta^{(1)}_{\Lambda}](\bm{k})+\dots\;.$
(2.21)

Notice that throughout we have assumed that the initial conditions can be
written in terms of a single degree of freedom, highlighted by explicitly
writing $\Psi_{\Lambda}[\delta^{(1)}_{\Lambda}]$. In multifield inflationary
scenarios, the potential perturbations $\varphi$ are in general given by a sum
over several independent degrees of freedom. To capture this, our formalism
can be straightforwardly generalized to integrate out these fields as well,
using additional propagators and interaction terms appropriate to these
fields. In that case, there would be additional contributions to Eq. (2.21)
that are independent of $\delta^{(1)}_{\Lambda}$ [69, 70].

Beyond the linear matter propagators defined in Eq. (1.8) and Eq. (1.9), we
now also have propagators defined for the new non-Gaussian (linear) operator
$\Psi$

$\displaystyle P^{\Lambda}_{\rm L,\Psi}(k)$
$\displaystyle=\left<\delta^{(1)}_{\Lambda}(\bm{k})\Psi_{\Lambda}(\bm{k}^{\prime})\right>^{\prime}=\left(\frac{k}{k_{\rm
NG}}\right)^{\Delta}\frac{1}{M(k)}P^{\Lambda}_{\rm L}(k)\,,$ $\displaystyle
P_{\rm shell}^{\Psi}(k)$ $\displaystyle=\left<\delta^{(1)}_{\rm
shell}(\bm{k})\Psi_{\rm
shell}(\bm{k}^{\prime})\right>^{\prime}=\left(\frac{k}{k_{\rm
NG}}\right)^{\Delta}\frac{1}{M(k)}P_{\rm shell}(k)\,.$ (2.22)

The propagators including two instances of $\Psi$, such as
$\langle\Psi_{\Lambda}(\bm{k})\Psi_{\Lambda}(\bm{k}^{\prime})\rangle^{\prime}$,
contribute at $\mathcal{O}(f_{\rm NL}^{2})$ and are not included in this work.
Including those terms would also force us to include new terms in the
interaction of order $g_{\rm NL}$ (quartic interactions). We expect both
contributions to be strongly suppressed by the fact that $\varphi\sim 10^{-5}$
[35], although this of course depends on the value of the coefficients $f_{\rm
NL},g_{\rm NL}$. However, even though terms proportional to $(f_{\rm NL})^{n}$
are ignored in the evaluation of the RG flow, they are generated at the
$n$-point function level even from a linear in $f_{\rm NL}$ partition function
and can be significant at sufficiently large scales if one of the modes is of
order $aH$ [28]. Similar to Eq. (1.10), we also have the variance due to the
mixed propagator

$\displaystyle\sigma^{2}_{\Psi,\Lambda}=\int_{\bm{p}}^{\Lambda}\left<\delta^{(1)}_{\Lambda}(\bm{k})\Psi_{\Lambda}(\bm{k}^{\prime})\right>^{\prime}=\int_{\bm{p}}^{\Lambda}P_{{\rm
L},\Psi}^{\Lambda}\,.$ (2.23)

#### Spin-2.

For spin-2 ($L=2$) non-Gaussianity, the squeezed limit of the kernel motivates
the introduction of a trace-free two-tensor

$\Psi^{\Lambda}_{ij}(\bm{k})=\frac{3}{2}\left(\frac{k_{i}k_{j}}{k^{2}}-\frac{\delta_{ij}}{3}\right)\left(\frac{k}{k_{\rm
NG}}\right)^{\Delta}\varphi_{\Lambda}(\bm{k})\,.$ (2.24)

The non-Gaussian bias operator at leading order will be
$\Psi^{ij}\Pi_{ij}^{[1]}=\operatorname{Tr}\left[\Psi\Pi^{[1]}\right]$, with
[28]

$\Pi_{ij}^{[1]}(\bm{k})=\frac{k_{i}k_{j}}{k^{2}}\delta\,.$ (2.25)

For spin-2 non-Gaussianity, the operator basis becomes:

First order: $\displaystyle\quad\delta;$ Second order:
$\displaystyle\quad\delta^{2},\;\mathcal{G}_{2},\;\operatorname{Tr}\left[\Psi\Pi^{[1]}\right];\;$
Third order:
$\displaystyle\quad\delta^{3},\;\delta\mathcal{G}_{2},\;\delta\;\operatorname{Tr}\left[\Psi\Pi^{[1]}\right],\;\Gamma_{3},\;\mathcal{G}_{3},\;\operatorname{Tr}\left[\Psi\Pi^{[2]}\right],$
(2.26)

with

$\Pi^{ij}_{[2]}=\partial^{i}\partial_{k}\Phi_{g}\partial^{k}\partial^{j}\Phi_{g}-\frac{5}{7}\frac{\partial^{i}\partial^{j}}{\nabla^{2}}\mathcal{G}_{2}\,.$
(2.27)

Notice that spin-2 PNG contributes to the galaxy bias expansion only at second
order in perturbation theory. For non-scalar tracers such as galaxy shapes,
this changes however: in that case there is a leading-order contribution
$\propto\Psi_{ij}$ [71, 72]. Scalar contractions of order $(\Psi_{ij})^{n}$
with $n>1$ should in principle be included, but those terms are
$\mathcal{O}(f_{\rm NL}^{n})$ which are beyond the scope of this work.

One final important point to note is that the fields $\Psi,\Psi_{ij}$ are to
be evaluated at the Lagrangian position $\bm{q}$ corresponding to the Eulerian
position of observation $\bm{x}$ (see App. D for a more detailed
justification).

###  2.4 The non-Gaussian partition function for the galaxy bias

We now generalize the Gaussian partition function to include PNG as
interactions which will then be translated into non-Gaussian contributions to
$n$-point functions. The guiding point for this is to construct a partition
function that reproduces the primordial curvature bispectrum at first order of
the nonlinear kernels. We of course want to keep the integrals to be over the
Gaussian initial conditions after expanding the interactions, such as Eq.
(2.1), such that we can still use QFT tools. We can then modify the Gaussian
measure to account for slightly non-Gaussian fields. Inserting Eq. (2.14) into
the kinetic term of the action, we get

$\displaystyle-\frac{1}{2}\int_{\bm{k}}^{\Lambda}\frac{|\delta^{(1)}_{\Lambda}|^{2}}{P_{\rm
L,\Lambda}(k)}\rightarrow-\frac{1}{2}\int_{\bm{k}}^{\Lambda}\frac{|\delta^{(1)}_{\Lambda}|^{2}}{P_{\rm
L,\Lambda}(k)}-f_{\rm
NL}\int_{\bm{p}_{1},\bm{p}_{2},\bm{p}_{3}}\left[\hat{\delta}_{\rm
D}(\bm{p}_{123})\right.$ (2.28) $\displaystyle\hskip
142.26378pt\times\left.K_{\rm
NL}(\bm{p}_{2},\bm{p}_{3})\frac{M(p_{1})}{M(p_{2})M(p_{3})}\frac{\delta^{(1)}_{\Lambda}(\bm{p}_{1})}{P_{\rm
L}^{\Lambda}(p_{1})}\delta^{(1)}_{\Lambda}(\bm{p}_{2})\delta^{(1)}_{\Lambda}(\bm{p}_{3})\right]\;.$

Notice that the new cubic kernel includes a $(P_{\rm L}^{\Lambda})^{-1}$
factor, which has an important role when reproducing $n$-point statistics. We
hereafter drop the subscript “G” in $\delta^{(1)}_{\rm\Lambda,G}$. The
equation above implies that we can add an interaction term to the Lagrangian
which is equivalent (at linear order in $f_{\rm NL}$) to modifying the
Gaussian measure in the partition function, therefore it will induce non-
Gaussian statistics (i.e non-vanishing bispectrum).

Our modified partition function is therefore given by

$\mathcal{Z}[J_{m,\Lambda},J_{g,\Lambda}]=\int\mathcal{D}\delta^{(1)}_{\Lambda}\mathcal{P}[\delta^{(1)}_{\Lambda}]\exp\left\\{{S_{\rm
int}+\int_{\bm{k}}J_{m,\Lambda}(\bm{k})\delta[\delta^{(1)}_{\Lambda}](-\bm{k})+\int_{\bm{k}}J_{g,\Lambda}(\bm{k})\delta_{g}[\delta^{(1)}_{\Lambda}]}(-\bm{k})\right\\}\,,$
(2.29)

with $S_{\rm int}$ a cubic interaction of the form

$\boxed{\begin{aligned} S_{\rm int}[\delta^{(1)}_{\Lambda}]&=f_{\rm
NL}\int_{\bm{p}_{1},\bm{p}_{2},\bm{p}_{3}}\hat{\delta}_{\rm
D}(\bm{p}_{123})K_{\rm
NL}(\bm{p}_{2},\bm{p}_{3})\frac{M(p_{1})}{M(p_{2})M(p_{3})}\frac{\delta^{(1)}_{\Lambda}(\bm{p}_{1})}{P_{\rm
L,\Lambda}(p_{1})}\delta^{(1)}_{\Lambda}(\bm{p}_{2})\delta^{(1)}_{\Lambda}(\bm{p}_{3})\\\
&\hskip
113.81102pt=\raisebox{0.0pt}{\parbox[c]{0.0pt}{\includegraphics[scale={.7}]{figs/diags/diag_vertex0.pdf}}}\,,\end{aligned}}$
(2.30)

where in the diagrammatic representation of the $f_{\rm NL}$ vertex the dashed
lines represent the $\delta^{(1)}_{\Lambda}$ legs and the doubled-dashed line
corresponds to the leg that is divided by a propagator,
$\delta^{(1)}_{\Lambda}(\bm{p})/P_{\rm L,\Lambda}(p)$, while the momenta are
all incoming due to the Dirac-delta function $\hat{\delta}_{\rm
D}(\bm{p}_{123})$. We stress again that the terms coupled to the current obey
the Gaussian equations of motion, since all non-Gaussianites are induced from
the interaction term. From this we can reproduce all the non-Gaussian matter
and bias $n$-point statistics, by taking derivatives with respect to the
different currents of the partition function. For example, for the matter
bispectrum

$\displaystyle\langle\delta(\bm{k}_{1})\delta(\bm{k}_{2})\delta(\bm{k}_{3})\rangle$
$\displaystyle=\frac{1}{\mathcal{Z}_{\Lambda}[0]}\frac{\delta^{3}\mathcal{Z}}{\delta
J_{m,\Lambda}\delta J_{m,\Lambda}\delta
J_{m,\Lambda}}\Big{|}_{J_{m,\Lambda}=0,J_{g,\Lambda}=0}$
$\displaystyle=\int\mathcal{D}\delta^{(1)}_{\Lambda}\mathcal{P}[\delta^{(1)}_{\Lambda}]e^{S_{\rm
int}}\;\delta[\delta^{(1)}_{\Lambda}](\bm{k}_{1})\delta[\delta^{(1)}_{\Lambda}](\bm{k}_{2})\delta[\delta^{(1)}_{\Lambda}](\bm{k}_{3})\;.$
(2.31)

We see that by taking three derivatives, we attach a propagator to each of the
three incoming lines. The double dashed line, which represents the
$\delta^{(1)}_{\Lambda}(\bm{p})/P_{\rm L,\Lambda}(p)$, signals that this line
will be cut, resulting in the usual form of the bispectrum with the two
propagators (power spectra) attached to it [see also Eq. (2.12)]

$\raisebox{0.0pt}{\parbox[c]{0.0pt}{\includegraphics[scale={.7}]{figs/diags/bispec1.pdf}}}\,.$
(2.32)

We discuss some specific examples in App. B, showing that this partition
function reproduces the known results for the matter and galaxy bispectra.

Finally, we can see that when calculating the expectation value of an operator

$\displaystyle\langle
O\rangle\equiv\int\mathcal{D}\delta^{(1)}\mathcal{P}[\delta^{(1)}]\;O\;e^{S_{\rm
int}}$
$\displaystyle=\int\mathcal{D}\delta^{(1)}\mathcal{P}[\delta^{(1)}]\;\left[O+O\times
S_{\rm int}+\mathcal{O}\left(S_{\rm int}^{2}\right)\right]$
$\displaystyle\equiv\langle O\rangle_{\rm free}+\langle O\rangle_{\rm int}\;,$
(2.33)

where we have split the terms into a free and an interacting part, with

$\displaystyle\langle O\rangle_{\rm free}$
$\displaystyle\equiv\int\mathcal{D}\delta^{(1)}\mathcal{P}[\delta^{(1)}]\;O=\langle
O\rangle_{\rm G}\,,$ $\displaystyle\langle O\rangle_{\rm int}$
$\displaystyle\equiv\langle O\times S_{\rm int}\rangle_{\rm
G}+\mathcal{O}\left(S_{\rm int}^{2}\right)\,,$ (2.34)

where we use the notation $\langle\dots\rangle_{\rm G}$ to represent
expectation values w.r.t the Gaussian measure. Hereafter, to not clutter the
text, we drop the contributions quadratic in $f_{\rm NL}$ from
$\mathcal{O}(S_{\rm int}^{2})$. In summary, we can separate expectation values
into the usual expectation value for Gaussian initial conditions and a new
term induced purely by the non-Gaussian interaction $S_{\rm int}$.

##  3 RG flow via the partition function

Now that we have the non-Gaussian partition function Eq. (2.29) for the galaxy
bias, we proceed to integrate out the shell modes similarly to [44]. We start
with the Wilson-Polchinski approach for the shell integrals in the presence of
PNG. Later we describe the contribution of the new interaction vertex to the
running of bias terms and then we derive the RGE for the bias parameters.

###  3.1 The Wilson-Polchinski shell integrals

We want to consider the running over an infinitesimal momentum difference
$\Lambda=\Lambda^{\prime}-\lambda$ integrating out field modes that have
support in the shell with infinitesimal width $\lambda$. This means that the
modes can be split in

$\delta^{(1)}_{\Lambda^{\prime}}(\bm{k})=\delta^{(1)}_{\Lambda}(\bm{k})+\delta^{(1)}_{\rm
shell}(\bm{k})\,.$ (3.1)

We then start from the partition function for the bias expansion making
explicit the split between shell and $\Lambda$ modes

$\displaystyle\mathcal{Z}[J_{\Lambda}]=\int\mathcal{D}\delta^{(1)}_{\Lambda}\mathcal{P}[\delta^{(1)}_{\Lambda}]\int\mathcal{D}\delta^{(1)}_{\rm
shell}\mathcal{P}[\delta^{(1)}_{\rm shell}]\exp\left(S_{\rm
int}[\delta^{(1)}_{\Lambda}+\delta^{(1)}_{\rm shell}]]\right)$ (3.2)
$\displaystyle\hskip
113.81102pt\times\exp\left(\int_{\bm{k}}J_{\Lambda}(\bm{k})[\sum_{O}b_{O}^{\Lambda^{\prime}}O[\delta^{(1)}_{\Lambda}+\delta^{(1)}_{\rm
shell}](-\bm{k})]+\mathcal{O}\left[J^{2}\right]\right)\,,$

where the current $J$ has support only up to $\Lambda$, in order to guarantee
the orthogonality condition $\int_{\bm{p}}J_{\Lambda}(\bm{p})\delta^{(1)}_{\rm
shell}(\bm{p})=0$. We can then expand the operators evaluated at $n^{\rm
th}$-order $O^{(n)}[\delta^{(1)}_{\Lambda}+\delta^{(1)}_{\rm shell}]$ in
powers of $\delta^{(1)}_{\rm shell}$ (see Appendix A.1 of [44] for examples)

$\displaystyle O^{(n)}[\delta^{(1)}_{\Lambda}+\delta^{(1)}_{\rm
shell}]=O^{(n)}[\delta^{(1)}_{\Lambda}]+O^{(n),{(1)_{\rm
shell}}}[\delta^{(1)}_{\Lambda},\delta^{(1)}_{\rm shell}]+O^{(n),{(2)_{\rm
shell}}}[\delta^{(1)}_{\Lambda},\delta^{(1)}_{\rm shell}]$ (3.3)
$\displaystyle\hskip 142.26378pt+\ldots+O^{(n),{(n-1)_{\rm
shell}}}[\delta^{(1)}_{\Lambda},\delta^{(1)}_{\rm
shell}]+O^{(n)}[\delta^{(1)}_{\rm shell}]\,.$

In addition, the partition function with PNG includes the interaction term Eq.
(2.30)

$\displaystyle S_{\rm int}[\delta^{(1)}_{\Lambda}+\delta^{(1)}_{\rm shell}]=$
(3.4) $\displaystyle\hskip 28.45274ptf_{\rm
NL}\int_{\bm{p}_{1},\bm{p}_{2},\bm{p}_{3}}\hat{\delta}_{\rm
D}(\bm{p}_{123})K_{\rm
NL}\left(\bm{p}_{2},\bm{p}_{3}\right)\frac{M(\bm{p}_{1})}{M\left(\bm{p}_{3}\right)M(\bm{p}_{2})}\frac{\delta^{(1)}_{\Lambda^{\prime}}\left(\bm{p}_{1}\right)}{P_{{\rm
L},\Lambda^{\prime}}\left(\bm{p}_{1}\right)}\delta^{(1)}_{\Lambda^{\prime}}\left(\bm{p}_{2}\right)\delta^{(1)}_{\Lambda^{\prime}}\left(\bm{p}_{3}\right)\,.$

Hence, on top of the expansion Eq. (3.3) for the bias operators, this
interaction term also has to be expanded in powers of $\delta^{(1)}_{\rm
shell}$, as we discuss below in Sec. 3.2 and with more details in App. C.

After integrating out the shell modes, we find the effective partition
function

$\displaystyle\mathcal{Z}[J_{\Lambda}]$
$\displaystyle=\int\mathcal{D}\delta^{(1)}_{\Lambda}\mathcal{P}[\delta^{(1)}_{\Lambda}]\exp\left(S_{\rm
int}[\delta^{(1)}_{\Lambda}]\right)$
$\displaystyle\times\exp\left(\int_{\bm{k}}J_{\Lambda}(\bm{k})\sum_{O}b_{O}^{\Lambda^{\prime}}\left[O[\delta^{(1)}_{\Lambda}](-\bm{k})+\mathcal{S}_{O}[\delta^{(1)}_{\Lambda}](-\bm{k})\right]+\mathcal{O}(J_{\Lambda}^{2})\right)\;,$
(3.5)

where $\mathcal{S}_{O}[\delta^{(1)}_{\Lambda}]$ is the shell correction, which
we can write as

$\displaystyle\mathcal{S}_{O}(\bm{k})$
$\displaystyle\equiv\left[\mathcal{S}_{O}^{2}\right]_{\rm
free}(\bm{k})+\left[\mathcal{S}_{O}\right]_{\rm
int}(\bm{k})+\mathcal{O}(\lambda^{2})\,,$ (3.6)

where we have dropped higher-loops contributions proportional to
$\lambda^{2}$; as argued in [44], we can always choose $\lambda$ sufficiently
small so that $\mathcal{O}(\lambda^{2})$ corrections can be neglected. We have
split into a free and an interacting part similar to Eq. (2.4), with the free
part being defined as in [44]

$\displaystyle\left[\mathcal{S}_{O}^{2}\right]_{\rm free}(\bm{k})$
$\displaystyle\equiv\int\mathcal{D}\delta^{(1)}_{\rm
shell}\mathcal{P}[\delta^{(1)}_{\rm shell}]\sum_{n\geq 2}O^{(n),(2)_{\rm
shell}}[\delta^{(1)}_{\Lambda},\delta^{(1)}_{\rm
shell}](\bm{k})=\sum_{\ell}\raisebox{0.0pt}{\parbox[c]{0.0pt}{\includegraphics[scale={.9}]{figs/diags/diag_S2_ell_legs.pdf}}}\,,$
(3.7)

where we see that the contribution from two shells is still the relevant term
for the free component, since $\delta^{(1)}_{\rm shell}$ is still a Gaussian
field. The new interaction arising due to the inclusion of PNG can be written
as

$\displaystyle\left[\mathcal{S}_{O}\right]_{\rm int}(\bm{k})$
$\displaystyle\equiv\int\mathcal{D}\delta^{(1)}_{\rm
shell}\mathcal{P}[\delta^{(1)}_{\rm
shell}]O[\delta^{(1)}_{\Lambda},\delta^{(1)}_{\rm shell}](\bm{k})\left(S_{\rm
int}[\delta^{(1)}_{\Lambda}+\delta^{(1)}_{\rm shell}]\right)\,.$ (3.8)

The structure of this interaction term is the topic of next section.

###  3.2 The cubic interaction

We now discuss the shell expectation values in the presence of cubic
interaction PNG, focusing only on the central points. The main conclusion is
summarized in Eq. (3.15) and we leave a detailed discussion for App. C.

We start by splitting the interaction term evaluated for
$\delta^{(1)}_{\Lambda^{\prime}}=\delta^{(1)}_{\Lambda}+\delta^{(1)}_{\rm
shell}$ into four components

$\displaystyle S_{\rm int}[\delta^{(1)}_{\Lambda}+\delta^{(1)}_{\rm
shell}]=S_{\rm int}[\delta^{(1)}_{\Lambda}]+S_{\rm int}^{(1)_{\rm
shell}}[\delta^{(1)}_{\Lambda},\delta^{(1)}_{\rm shell}]$ (3.9)
$\displaystyle\hskip 113.81102pt+S_{\rm int}^{(2)_{\rm
shell}}[\delta^{(1)}_{\Lambda},\delta^{(1)}_{\rm shell}]+S_{\rm int}^{(3)_{\rm
shell}}[\delta^{(1)}_{\rm shell}]\,,$

where $S_{\rm int}^{(i)_{\rm shell}}$ refers to the number of (upstairs) legs
in Eq. (3.4) that have support in the shell.555Notice however that the net
shell counting can change by a factor -2 depending whether
$\delta^{(1)}/P_{{\rm L},\Lambda^{\prime}}$ has support on the shell or not,
due the $1/P_{\rm shell}$ factor. Moreover, the support of the terms up and
downstairs of $\delta^{(1)}_{\Lambda^{\prime}}\left(\bm{p}\right)/P_{{\rm
L},\Lambda^{\prime}}\left(\bm{p}\right)$ are always related since
$\frac{\delta^{(1)}_{\Lambda^{\prime}}\left(\bm{p}\right)}{P_{{\rm
L},\Lambda^{\prime}}\left(\bm{p}\right)}=\frac{\delta^{(1)}_{\Lambda}\left(\bm{p}\right)}{P_{{\rm
L},\Lambda}\left(\bm{p}\right)}+\frac{\delta^{(1)}_{\rm
shell}\left(\bm{p}\right)}{P_{\rm shell}\left(\bm{p}\right)}\,.$ (3.10) We
find for each one of those terms:

  * •

$S_{\rm int}[\delta^{(1)}_{\Lambda}]$ captures the interactions of modes that
are not integrated out, and therefore does not act as source for the RGE.

  * •

The contribution with one shell is given by

$\displaystyle S_{\rm int}^{(1)_{\rm
shell}}[\delta^{(1)}_{\Lambda},\delta^{(1)}_{\rm shell}]$
$\displaystyle=f_{\rm
NL}\int_{\bm{p}_{1},\bm{p}_{2},\bm{p}_{3}}\hat{\delta}_{\rm
D}(\bm{p}_{123})K_{\rm
NL}\left(\bm{p}_{2},\bm{p}_{3}\right)\frac{M(p_{1})}{M(p_{2})M(p_{3})}$ (3.11)
$\displaystyle\times\left[\frac{\delta^{(1)}_{\rm shell}(\bm{p}_{1})}{P_{\rm
shell}(p_{1})}\delta^{(1)}_{\Lambda}\left(\bm{p}_{2}\right)\delta^{(1)}_{\Lambda}\left(\bm{p}_{3}\right)+2\,\frac{\delta^{(1)}_{\Lambda}\left(\bm{p}_{1}\right)}{P_{\rm
L,\Lambda}\left(p_{1}\right)}\delta^{(1)}_{\rm
shell}\left(\bm{p}_{2}\right)\delta^{(1)}_{\Lambda}\left(\bm{p}_{3}\right)\right]$
$\displaystyle=\raisebox{0.0pt}{\parbox[c]{0.0pt}{\includegraphics[scale={.7}]{figs/diags/diag_vertex1cut.pdf}}}+\raisebox{0.0pt}{\parbox[c]{0.0pt}{\includegraphics[scale={.7}]{figs/diags/diag_vertex1.pdf}}}\,,$

in which, again, the double lines represent the $\delta^{(1)}_{\rm
shell}(\bm{p}_{1})/P_{\rm shell}(p_{1})$ (in the left diagram) or
$\delta^{(1)}_{\Lambda}\left(\bm{p}_{1}\right)/P_{\rm
L,\Lambda}\left(\bm{p}_{1}\right)$ legs (in the right diagram). The red wiggly
lines represent shell modes, whereas dashed lines represent lines with support
in $\Lambda$. Diagrams of this type are not allowed by the momentum constrain
$\hat{\delta}_{\rm D}(\bm{p}_{123})$ when taking two of the legs with momentum
much smaller than the third one that is on the shell.

  * •

The contributions with two shell modes demand more attention. They are given
by

$\displaystyle S_{\rm int}^{(2)_{\rm
shell}}[\delta^{(1)}_{\Lambda},\delta^{(1)}_{\rm shell}]$
$\displaystyle=f_{\rm
NL}\int_{\bm{p}_{1},\bm{p}_{2},\bm{p}_{3}}\hat{\delta}_{\rm
D}(\bm{p}_{123})K_{\rm
NL}\left(\bm{p}_{2},\bm{p}_{3}\right)\frac{M(p_{1})}{M(p_{2})M(p_{3})}$
$\displaystyle\times\left[2\,\frac{\delta^{(1)}_{\rm
shell}(\bm{p}_{1})}{P_{\rm shell}(p_{1})}\delta^{(1)}_{\rm
shell}\left(\bm{p}_{2}\right)\delta^{(1)}_{\Lambda}\left(\bm{p}_{3}\right)+\frac{\delta^{(1)}_{\Lambda}\left(\bm{p}_{1}\right)}{P_{\rm
L,\Lambda}\left(p_{1}\right)}\delta^{(1)}_{\rm
shell}\left(\bm{p}_{2}\right)\delta^{(1)}_{\rm
shell}\left(\bm{p}_{3}\right)\right]$ $\displaystyle\equiv\left[S_{\rm
int}^{(3),(2)_{\rm shell}}\right]_{\rm scatter}+\left[S_{\rm
int}^{(3),(2)_{\rm shell}}\right]_{\rm annihilation}$ (3.12)
$\displaystyle=\raisebox{0.0pt}{\parbox[c]{0.0pt}{\includegraphics[scale={.7}]{figs/diags/diag_vertex2.pdf}}}+\raisebox{0.0pt}{\parbox[c]{0.0pt}{\includegraphics[scale={.7}]{figs/diags/diag_vertex2cut.pdf}}}\,,$

where we have split into a scattering part, in which a shell UV mode interacts
with a $\Lambda$ IR mode leading to another shell mode, and an annihilation
part, with two shell UV modes interacting to generate a large-scale mode below
$\Lambda$. We show (up to our knowledge for the first time in the literature)
in App. C.3 that the annihilation modes lead to a stochastic contribution that
starts to contribute at $\mathcal{O}(J^{2})$ . We leave a full treatment of
stochastic contributions in the presence of PNG for a future work (for
Gaussian initial conditions, see [45]). After contracting the scattering
interaction with two other shell modes we find (see also Appendix A.2 of [35])

$\boxed{\begin{aligned} \left[\mathcal{S}_{O}\right]_{\rm
int}(\bm{k})&=\int\mathcal{D}\delta^{(1)}_{\rm
shell}\mathcal{P}[\delta^{(1)}_{\rm shell}]\sum_{n\geq 2}O^{(n),(2)_{\rm
shell}}(\bm{k})\left[S_{\rm int}^{(2)_{\rm shell}}\right]_{\rm scatter}\\\
&=\sum_{n\geq 2}4\int_{\bm{p}_{1},\dots,\bm{p}_{n}}\hat{\delta}_{\rm
D}(\bm{k}-\bm{p}_{1\dots
n})K_{O}^{(n)}(\bm{p}_{1},\bm{p}_{2},\dots,\bm{p}_{n})\delta^{(1)}_{\Lambda}\left(\bm{p}_{1}\right)\dots\delta^{(1)}_{\Lambda}\left(\bm{p}_{n-2}\right)\\\
&\hskip 56.9055ptK_{\rm
NL}(\bm{p}_{n-1}+\bm{p}_{n},-\bm{p}_{n})\frac{\delta^{(1)}_{\Lambda}\left(\bm{p}_{n-1}+\bm{p}_{n}\right)}{M\left(|\bm{p}_{n-1}+\bm{p}_{n}|\right)}P_{\rm
shell}(p_{n})\\\ &\hskip
85.35826pt=\sum_{\ell}\raisebox{0.0pt}{\parbox[c]{0.0pt}{\includegraphics[scale={.7}]{figs/diags/NG_inter.pdf}}}\,.\end{aligned}}$
(3.13)

We notice two things about the interaction term. First, that it contains an
extra external leg compared to the free term Eq. (3.7). Second, since in Eq.
(3.13) $\bm{p}_{n}$ is constrained to be in the shell and
$\bm{p}_{n-1}+\bm{p}_{n}$ is constrained to have support in $\Lambda$, we can
use the squeezed limit of the kernel $K_{\rm
NL}(\bm{p}_{n-1}+\bm{p}_{n},-\bm{p}_{n})$, Eq. (2.17).

  * •

The contribution with three shell modes

$\displaystyle S_{\rm int}^{(3)_{\rm shell}}[\delta^{(1)}_{\rm shell}]$
$\displaystyle=f_{\rm
NL}\int_{\bm{p}_{1},\bm{p}_{2},\bm{p}_{3}}\hat{\delta}_{\rm
D}(\bm{p}_{123})K_{\rm
NL}\left(\bm{p}_{2},\bm{p}_{3}\right)\frac{M(p_{1})}{M(p_{2})M(p_{3})}$ (3.14)
$\displaystyle\hskip 113.81102pt\times\frac{\delta^{(1)}_{\rm
shell}(\bm{p}_{1})}{P_{\rm shell}(p_{1})}\delta^{(1)}_{\rm
shell}\left(\bm{p}_{2}\right)\delta^{(1)}_{\rm
shell}\left(\bm{p}_{3}\right)\,.$
$\displaystyle=\raisebox{0.0pt}{\parbox[c]{0.0pt}{\includegraphics[scale={.7}]{figs/diags/diag_vertex3.pdf}}}\,,$

only leads to higher-loop or tadpole contributions.

In summary, we find that the shell corrections when including both free and
interacting terms are represented by two diagrams

$\boxed{\begin{aligned}
\mathcal{S}_{O}(\bm{k})&\equiv\left[\mathcal{S}_{O}^{2}\right]_{\rm
free}(\bm{k})+\left[\mathcal{S}_{O}\right]_{\rm int}(\bm{k})\\\
&=\sum_{\ell}\raisebox{0.0pt}{\parbox[c]{0.0pt}{\includegraphics[scale={.9}]{figs/diags/diag_S2_ell_legs.pdf}}}+\sum_{\ell}\raisebox{0.0pt}{\parbox[c]{0.0pt}{\includegraphics[scale={.9}]{figs/diags/NG_inter.pdf}}}\;.\end{aligned}}$
(3.15)

Via the free part, Gaussian operators source other Gaussian operators (as
calculated in [44]). Non-Gaussian operators in the new basis Eq. (2.3) or Eq.
(2.3) can source either Gaussian and non-Gaussian terms, as we calculate
below. The new interacting part, originating from the cubic PNG vertex, are
only due to the scattering term in Eq. (3.12). The interaction also allows for
the sourcing of non-Gaussian operators by Gaussian terms, while the
contributions from non-Gaussian to Gaussian terms that $S_{\rm int}$ sources
are suppressed by $f_{\rm NL}^{2}$.

###  3.3 On the $\mathcal{S}_{O^{\prime}}$ corrections to $O$ and the RG
equations

We can now derive the RG equations describing the evolution of the bias
parameters following the steps of [44, 45]. We compare the terms linear in the
current $J_{\Lambda}$ in Eq. (3.5) finding

$\displaystyle\sum_{O}b_{O}^{\Lambda}O[\delta^{(1)}_{\Lambda}](\bm{k})=\sum_{O}b_{O}^{\Lambda^{\prime}}\left(O[\delta^{(1)}_{\Lambda}](\bm{k})+\mathcal{S}_{O}[\delta^{(1)}_{\Lambda}](\bm{k})\right)\,.$
(3.16)

For the Gaussian initial condition case in [44], we could use that the shell
correlators $\mathcal{S}$ can be written as an operator expansion with
coefficients $s_{O^{\prime}}^{O}$

$\mathcal{S}_{O^{\prime}}[\delta^{(1)}_{\Lambda}]=\mathcal{S}_{O^{\prime}}^{2}[\delta^{(1)}_{\Lambda}]=\lambda\frac{d\sigma^{2}_{\Lambda}}{d\Lambda}\Big{|}_{\Lambda}\sum_{O}s_{O^{\prime}}^{O}O[\delta^{(1)}_{\Lambda}],$
(3.17)

where we pulled out the linear dependence on the shell width $\lambda$ through

$\int_{\bm{p}}P_{\rm
shell}(p)=\int_{\Lambda}^{\Lambda+\lambda}\frac{p^{2}dp}{2\pi^{2}}P_{\rm
L}(p)=\frac{d\sigma^{2}_{\Lambda}}{d\Lambda}\Big{|}_{\Lambda}\lambda+\mathcal{O}(\lambda^{2})\,.$
(3.18)

See the first block of Tab. 1 for a summary of the $s_{O^{\prime}}^{O}$
coefficients for the Gaussian case.

When including non-Gaussian interactions, on top of the shell integral above,
we will also find terms proportional to

$\displaystyle\int_{\bm{p}}P_{\rm
shell}^{\Psi}(p)=\int_{\Lambda}^{\Lambda+\lambda}\frac{dp\;p^{2}}{2\pi^{2}}P_{L,\Psi}(p)=\frac{d\sigma^{2}_{\Psi}}{d\Lambda}\lambda+\mathcal{O}(\lambda^{2})\;.$
(3.19)

We can relate this mixed variance to the variance $\sigma_{\Lambda}^{2}$ by
writing

$\displaystyle\int_{\Lambda}^{\Lambda+\lambda}\frac{dp\;p^{2}}{2\pi^{2}}P^{\Psi}_{\rm
shell}(p)$ $\displaystyle\simeq\left(\frac{\Lambda}{k_{\rm
NG}}\right)^{\Delta}\left(\frac{H_{0}}{\Lambda}\right)^{2}\frac{3\;\Omega_{m}}{2\;T(\Lambda)}\frac{d\sigma_{\Lambda}^{2}}{d\Lambda}\Big{|}_{\Lambda}\lambda\;,$
(3.20)

using Eq. (2.16). From the squeezed limit of the non-Gaussian kernel, we also
encounter terms of the form

$\displaystyle\int_{\Lambda}^{\Lambda+\lambda}\frac{dp\;p^{2}}{2\pi^{2}}\left(\frac{k_{\rm
NG}}{p}\right)^{\Delta}P_{\rm shell}(p)$
$\displaystyle\simeq\left(\frac{k_{\rm
NG}}{\Lambda}\right)^{\Delta}\frac{d\sigma_{\Lambda}^{2}}{d\Lambda}\Big{|}_{\Lambda}\lambda\;.$
(3.21)

Eq. (3.20) appears when considering contributions from non-Gaussian to
Gaussian operators, whereas Eq. (3.21) appears when considering the opposite
case. Notice that $(H_{0}/\Lambda)^{2}\ll 1$ acts as a suppression factor.
Also since $(\Lambda/k_{\rm NG})\lesssim 1$ in order to keep the perturbative
expansion under control, and $\Delta\geq 0$, we find that this factor acts as
a suppression in Eq. (3.20) but as an enhancing factor in Eq. (3.21).

The appearence of new types of shell integrals from PNG forces us to slightly
modify Eq. (3.17) to accommodate these new $\Lambda$-dependent factors via
functions $f_{O^{\prime}}^{O}(\Lambda)$, writing

$\mathcal{S}_{O^{\prime}}[\delta^{(1)}_{\Lambda}]=\lambda\frac{d\sigma^{2}_{\Lambda}}{d\Lambda}\Big{|}_{\Lambda}\sum_{O}s_{O^{\prime}}^{O}f_{O^{\prime}}^{O}(\Lambda)O[\delta^{(1)}_{\Lambda}]\,.$
(3.22)

Therefore, we can always parametrize the running using
$d\sigma_{\Lambda}^{2}/d\Lambda$: despite having more than one propagator,
they relate to each other in the shell region. We then find the ODE describing
the running of the bias parameters

$\boxed{\frac{d}{d\Lambda}b_{O}(\Lambda)=-\frac{d\sigma^{2}_{\Lambda}}{d\Lambda}\sum_{O^{\prime}}s_{O^{\prime}}^{O}\,f_{O^{\prime}}^{O}(\Lambda)b_{O^{\prime}}(\Lambda)\,,}$
(3.23)

similar to [44] but now with the extra factor $f_{O^{\prime}}^{O}$ that
satisfies (at linear order in $f_{\rm NL}$)

$\displaystyle f_{O^{\prime}}^{O}=1\,,\quad$ if $O$ and $O^{\prime}$ are both
Gaussian, $\displaystyle f_{O^{\prime}}^{O}=a_{L}f_{\rm NL}\,,\quad$ if $O$
and $O^{\prime}$ are both non-Gaussian, $\displaystyle
f_{O^{\prime}}^{O}(\Lambda)=a_{L}f_{\rm NL}\left(\frac{\Lambda}{k_{\rm
NG}}\right)^{\Delta}\left(\frac{H_{0}}{\Lambda}\right)^{2}\frac{3\;\Omega_{m}}{2\;T(\Lambda)}\,,\quad$
if $O$ is Gaussian and $O^{\prime}$ non-Gaussian, $\displaystyle
f_{O^{\prime}}^{O}(\Lambda)=a_{L}f_{\rm NL}\left(\frac{k_{\rm
NG}}{\Lambda}\right)^{\Delta}\,,\quad$ if $O$ is non-Gaussian and $O^{\prime}$
Gaussian. (3.24)

Eq. (3.23) fully describes the structure of the bias RGE with any second-order
PNG at $\mathcal{O}(f_{\rm NL})$, generalizing the previously derived Gaussian
case which corresponds to the first line of Eq. (3.24). In the following, we
explicitly derive the RGE for two types of PNG, the spin-0 and spin-2 cases.
We focus in this work on corrections to first and second-order operators by
operators that are up to third-order.666For spins higher than two, the
procedure to derive the RGE is similar, by adding the appropriate operators to
the respective bias expansion. The operators associated to $L>2$ are order
higher than 2 and therefore are not considered in this work. Odd ranked spins
(e.g., $L=1$) are zero at leading order due to the invariance of the squeezed
bispectrum in $\bm{k}_{s}\rightarrow-\bm{k}_{s}$ as mentioned in [35] (for
more details see Appendix A of [22] and [73, 74]).

#### Spin-0 non-Gaussianity.

We start by calculating the shell correlators (see App. D.2 for a more
complete explanation) and deriving the RGE for the spin-0 case. The assumption
that $k\ll p_{\rm shell}$ allow us to use the squeezed limit Eq. (2.17) of the
$K_{\rm NL}$ kernels, writing them in terms of indices $\Delta$ and projecting
onto the Legendre polynomial of order zero $\mathcal{P}_{0}=1$

$K_{\rm NL}(\bm{p}_{\rm shell},\bm{k})\xrightarrow{k\ll p_{\rm
shell}}a_{0}\left(\frac{k}{p_{\rm shell}}\right)^{\Delta}\,.$ (3.25)

Corrections to this approximation lead to higher-derivative non-Gaussian
operators such as $\nabla^{2}\Psi$, which we will neglect here; see [75, 76]
for a discussion and attempts to resum these beyond-squeezed-limit
corrections.

For the Gaussian operators, those that we already include in the bias
expansion without having non-Gaussianities, we find that their new
contributions due to the interaction vertex are

$\displaystyle\left[\mathcal{S}_{\delta^{2}}\right]_{\rm int}(\bm{k})$
$\displaystyle=\left(4f_{\rm NL}a_{0}\Psi(\bm{k})+\frac{272}{21}f_{\rm
NL}a_{0}[\Psi\delta](\bm{k})\right)\int_{\bm{p}}\left(\frac{k_{\rm
NG}}{p}\right)^{\Delta}P_{\rm shell}(p)+\text{h.d.}\;,$ (3.26)
$\displaystyle\left[\mathcal{S}_{\delta^{3}}\right]_{\rm int}(\bm{k})$
$\displaystyle=12f_{\rm
NL}a_{0}[\Psi\delta](\bm{k})\int_{\bm{p}}\left(\frac{k_{\rm
NG}}{p}\right)^{\Delta}P_{\rm shell}(p)+\text{h.d.}\;,$ (3.27)
$\displaystyle\left[\mathcal{S}_{\mathcal{G}_{2}\delta}\right]_{\rm
int}(\bm{k})$ $\displaystyle=-\frac{8}{3}f_{\rm
NL}a_{0}\left[\Psi\delta\right](\bm{k})\int_{\bm{p}}\left(\frac{k_{\rm
NG}}{p}\right)^{\Delta}P_{\rm shell}(p)+\text{h.d.}\;,$ (3.28)

in which we omitted higher-derivative (h.d.) terms. We can notice a few
things. First, non-Gaussian operators are sourced by Gaussian operators via
the interaction and one therefore has to include operators with $\Psi$ in the
basis Eq. (2.3) in order to be closed under normalization, as already pointed
out by [35]. Second, as expected, at linear-order in $f_{\rm NL}$ only non-
Gaussian operators are sourced by Gaussian operators via the interaction term.
Third, also expected, the amplitude of the RG running is proportional to
$f_{\rm NL}a_{0}$. Finally, we notice the factor $\left(k_{\rm
NG}/p\right)^{\Delta}$ in the integral, which is due to the squeezed limit Eq.
(3.25) and will lead to developments we comment on below.

For the new non-Gaussian operators, those are already linear order in $f_{\rm
NL}$. Therefore, their (linear-in-$f_{\rm NL}$) contribution cannot contain
any interaction vertex but has to come from the free term in Eq. (3.15). We
find

$\displaystyle\left[\mathcal{S}^{2}_{\Psi}\right]_{\rm free}(\bm{k})$
$\displaystyle=a_{0}f_{\rm
NL}\left(-\frac{13}{21}\delta(\bm{k})+\frac{43}{135}\left(\delta^{2}\right)(\bm{k})-\frac{1699}{13230}\left(\mathcal{G}_{2}\right)(\bm{k})\right)\int_{\bm{p}}P_{\rm
shell}^{\Psi}(p)+\text{h.d.}\,,$ (3.29)
$\displaystyle\left[\mathcal{S}^{2}_{\Psi\delta}\right]_{\rm free}(\bm{k})$
$\displaystyle=\left[\frac{13}{21}\delta(\bm{k})+\frac{478}{135}\delta^{2}(\bm{k})+\frac{79}{2205}\mathcal{G}_{2}(\bm{k})\right]\int_{\bm{p}}P_{\rm
shell}^{\Psi}(p)+\text{h.d.}\,,$ (3.30)
$\displaystyle\left[\mathcal{S}_{\psi\delta^{2}}^{2}\right]_{\rm
free}(\bm{k})$ $\displaystyle=a_{0}f_{\rm NL}\Psi(\bm{k})\int_{\bm{p}}P_{\rm
shell}(p)+2a_{0}f_{\rm NL}\delta(\bm{k})\int_{\bm{p}}P_{\rm shell}^{\Psi}$
(3.31) $\displaystyle+\frac{68}{21}a_{0}f_{\rm
NL}\left[\Psi\delta\right](\bm{k})\int_{\bm{p}}P_{\rm
shell}(p)+\frac{47}{21}a_{0}f_{\rm NL}\delta^{2}(\bm{k})\int_{\bm{p}}P_{\rm
shell}^{\Psi}(p)+\text{h.d.}\,,$
$\displaystyle\left[\mathcal{S}_{\rm\Psi\mathcal{G}_{2}}^{2}\right]_{\rm
free}(\bm{k})$ $\displaystyle=a_{0}f_{\rm
NL}\left[-\frac{4}{3}\delta(\bm{k})-\frac{31}{21}\delta^{2}(\bm{k})-\frac{1}{21}\mathcal{G}_{2}(\bm{k})\right]\int_{\bm{p}}P_{\rm
shell}^{\Psi}(p)+\text{h.d.}\,.$ (3.32)

We notice they contribute both to Gaussian [via $P^{\Lambda}_{\rm L,\Psi}$]
and also non-Gaussian operators [via $P_{\rm L,\Lambda}(p)$]. We stress again
that the operator $\Psi$ is calculated at the Lagrangian point
$\Psi(\bm{q})=\Psi(\bm{x}-\mathbf{s})$ and expanded up to third-order in
displacements (total perturbative order 4; see App. D.1). In addition, we
highlight that the first-order operator $\Psi$ contributes to other (leading-
in-derivative) terms due to the Lagrangian displacements, differently than
$\delta$ that only contributes to higher-derivative terms due to mass and
momentum conservation. We present the details of this calculation in App. D.
We emphasize that terms that break the Galilean invariance are cancelled, as
required by the equivalence principle, and the coefficients of the operator
$\delta^{(n)}$ are the same at first ($n=1$) and second order ($n=2$) in
perturbation theory. Specifically, the operators in Eqs. (3.29)–(3.32) are
calculated up to second order; in particular, $\delta$ here corresponds to
$\delta^{(1)}+\delta^{(2)}$.

After having calculated the shell contributions, we can derive the RGE for
spin-0 non-Gaussianities and a generic scaling index $\Delta$ using Eq.
(3.23), leading to

$\displaystyle\frac{db_{\delta}}{d\Lambda}$
$\displaystyle=\frac{db_{\delta}}{d\Lambda}\Big{|}_{\rm free}-a_{0}f_{\rm
NL}\left[-\frac{13}{21}b_{\Psi}+\frac{13}{21}b_{\Psi\delta}+2b_{\Psi\delta^{2}}-\frac{4}{3}b_{\Psi\mathcal{G}_{2}}\right]\left(\frac{\Lambda}{k_{\rm
NG}}\right)^{\Delta}\left(\frac{H_{0}}{\Lambda}\right)^{2}\frac{3\;\Omega_{m}}{2\;T(\Lambda)}\frac{d\sigma_{\Lambda}^{2}}{d\Lambda}\,,$
(3.33) $\displaystyle\frac{db_{\delta^{2}}}{d\Lambda}$
$\displaystyle=\frac{db_{\delta^{2}}}{d\Lambda}\Big{|}_{\rm free}-a_{0}f_{\rm
NL}\left[\frac{43}{135}b_{\Psi}+\frac{478}{135}b_{\Psi\delta}+\frac{47}{21}b_{\Psi\delta^{2}}-\frac{31}{21}b_{\Psi\mathcal{G}_{2}}+b^{\\{\delta^{2}\\}_{\rm
NG}}_{n=4}\right]$ $\displaystyle\hskip
199.16928pt\times\left(\frac{\Lambda}{k_{\rm
NG}}\right)^{\Delta}\left(\frac{H_{0}}{\Lambda}\right)^{2}\frac{3\;\Omega_{m}}{2\;T(\Lambda)}\frac{d\sigma_{\Lambda}^{2}}{d\Lambda}\,,$
(3.34) $\displaystyle\frac{db_{\mathcal{G}_{2}}}{d\Lambda}$
$\displaystyle=\frac{db_{\mathcal{G}_{2}}}{d\Lambda}\Big{|}_{\rm
free}-a_{0}f_{\rm
NL}\left[-\frac{1699}{13230}b_{\Psi}+\frac{79}{2205}b_{\Psi\delta}-\frac{1}{21}b_{\Psi\mathcal{G}_{2}}+b^{\\{\mathcal{G}_{2}\\}_{\rm
NG}}_{n=4}\right]$ $\displaystyle\hskip
199.16928pt\times\left(\frac{\Lambda}{k_{\rm
NG}}\right)^{\Delta}\left(\frac{H_{0}}{\Lambda}\right)^{2}\frac{3\;\Omega_{m}}{2\;T(\Lambda)}\frac{d\sigma_{\Lambda}^{2}}{d\Lambda}\,,$
(3.35)

where $\frac{db}{d\Lambda}\Big{|}_{\rm free}$ refer to the running of Gaussian
coefficients without any interactions described in [44]. We also have included
$b^{\\{O\\}_{\rm NG}}_{n=4}$, the contributions of fourth-order non-Gaussian
operators to the $b_{O}$, which are beyond the scope of this work.777Notice
that $\frac{db}{d\Lambda}\Big{|}_{\rm free}$ includes $b^{\\{O\\}_{\rm
G}}_{n=4}$, the contribution of fourth-order Gaussian operators to $O$. In
general an operator of order $n$ contributes to all operators that are order
$n-2$ [45]. We find that the running of Gaussian coefficients due to non-
Gaussian operators is doubly suppressed by $(H_{0}/\Lambda)^{2}$ and
$(\Lambda/k_{\rm NG})^{\Delta}$.

For the non-Gaussian operators, we find, including both the free and
interaction terms,

$\displaystyle\frac{db_{\Psi}}{d\Lambda}$ $\displaystyle=-a_{0}f_{\rm
NL}b_{\Psi\delta^{2}}\frac{d\sigma^{2}_{\Lambda}}{d\Lambda}-4a_{0}f_{\rm
NL}b_{\delta^{2}}\left(\frac{k_{\rm
NG}}{\Lambda}\right)^{\Delta}\frac{d\sigma_{\Lambda}^{2}}{d\Lambda}\,,$ (3.36)
$\displaystyle\frac{db_{\Psi\delta}}{d\Lambda}$ $\displaystyle=-a_{0}f_{\rm
NL}\left[\frac{68}{21}b_{\Psi\delta^{2}}+b^{\\{\Psi\delta\\}_{\rm
NG}}_{n=4}\right]\frac{d\sigma^{2}_{\Lambda}}{d\Lambda}$ $\displaystyle\quad-
a_{0}f_{\rm
NL}\left[\frac{272}{21}b_{\delta^{2}}+12b_{\delta^{3}}-\frac{8}{3}b_{\mathcal{G}_{2}\delta}+b^{\\{\Psi\delta\\}_{\rm
G}}_{n=4}\right]\left(\frac{k_{\rm
NG}}{\Lambda}\right)^{\Delta}\frac{d\sigma_{\Lambda}^{2}}{d\Lambda}\,,$ (3.37)

where $b^{\\{O\\}_{\rm G}}_{n=4}$ is the contribution of fourth-order Gaussian
operators to $b_{O}$. We see that the running of non-Gaussian coefficients due
to Gaussian operators is enhanced by $(k_{\rm NG}/\Lambda)^{\Delta}$. We
summarize those findings in Tab. 1, with the coefficients $s_{O^{\prime}}^{O}$
of Eq. (3.23). We further discuss solutions for those equations in Sec. 4.

$s_{O^{\prime}}^{O}$ | $\delta^{2}$ | $\delta^{3}$ | $\delta\mathcal{G}_{2}$ | $\Psi$ | $\Psi\delta$ | $\Psi\delta^{2}$ | $\Psi\mathcal{G}_{2}$ | $\operatorname{Tr}\Psi\Pi^{[1]}$ | $\delta\operatorname{Tr}\Psi\Pi^{[1]}$ | $\operatorname{Tr}\Psi\Pi^{[2]}$
---|---|---|---|---|---|---|---|---|---|---
$\delta$ | $68/21$ | 3 | $-4/3$ | $-13/21$ | $13/21$ | $2$ | $-4/3$ | $34/21$ | $1$ | $34/21$
$\delta^{2}$ | $8126/2205$ | $68/7$ | $-376/105$ | $43/135$ | $478/135$ | $47/21$ | $-31/21$ | $124/315$ | $178/105$ | $14347/6027$
$\mathcal{G}_{2}$ | $254/2205$ | - | $116/105$ | $-1699/13230$ | $79/2205$ | - | $-1/21$ | $-661/4410$ | $4/35$ | $-241/735$
$\Psi$ | $4$ | - | - | - | - | 1 | - | - | - | -
$\delta\Psi$ | $272/21$ | $12$ | $-8/3$ | - | - | $68/21$ | - | - | - | -
$\operatorname{Tr}\Psi\Pi^{[1]}$ | $64/105$ | - | $16/15$ | - | - | - | - | - | 8/105 | 58/305

Table 1: The $s_{O^{\prime}}^{O}$ coefficients from Eq. (3.23), that describe
how the $O^{\prime}$ operators (columns, up to third-order, including only
terms that source leading-in-derivative operators) contribute to $O$ operators
(rows, up to second order). This generalizes Table 1 of [45] to include non-
Gaussian contributions, for both spin-0 and spin-2 types.

#### Spin-2 non-Gaussianity.

We proceed to spin-2 non-Gaussianity using again the squeezed limit Eq. (2.17)
of the $K_{\rm NL}$ kernels and now projecting onto the Legendre polynomial of
order two,

$K_{\rm NL}(\bm{p}_{\rm shell},\bm{k})\xrightarrow{k\ll p_{\rm
shell}}a_{2}\left(\frac{k}{p_{\rm
shell}}\right)^{\Delta}\mathcal{P}_{2}\left(\hat{\bm{k}}\cdot\hat{\bm{p}}_{\rm
shell}\right)\,.$ (3.38)

The evaluation of the shell integrals for spin-2 field is included in more
details in App. D.3, but we summarize them here

$\displaystyle\left[\mathcal{S}_{\delta^{2}}\right]_{\rm int}(\bm{k})$
$\displaystyle=a_{2}f_{\rm
NL}\frac{64}{105}\operatorname{Tr}\left[\Psi\Pi^{[1]}\right](\bm{k})\int_{\bm{p}}\left(\frac{k_{\rm
NG}}{p}\right)^{\Delta}P_{\rm shell}(p)+\text{h.d.}\,,$ (3.39)
$\displaystyle\left[\mathcal{S}_{\delta\mathcal{G}_{2}}\right]_{\rm
int}(\bm{k})$ $\displaystyle=a_{2}f_{\rm
NL}\frac{16}{15}\operatorname{Tr}\left[\Psi\Pi^{[1]}\right](\bm{k})\int_{\bm{p}}\left(\frac{k_{\rm
NG}}{p}\right)^{\Delta}P_{\rm shell}(p)+\text{h.d.}\,,$ (3.40)

for the contribution of the Gaussian terms via the new interaction and

$\displaystyle\left[\mathcal{S}^{2}_{\rm Tr[\Psi\Pi^{[1]}]}\right]_{\rm
free}(\bm{k})$ $\displaystyle=a_{2}f_{\rm
NL}\left[\frac{34}{21}\delta(\bm{k})+\frac{124}{315}\delta^{2}(\bm{k})-\frac{661}{4410}\mathcal{G}_{2}(\bm{k})\right]\int_{\bm{p}}P_{\rm
shell}^{\Psi}(p)+\text{h.d.}\,,$ (3.41)
$\displaystyle\left[\mathcal{S}^{2}_{\rm Tr[\delta\Psi\Pi^{[1]}]}\right]_{\rm
free}(\bm{k})$ $\displaystyle=a_{2}f_{\rm
NL}\left[\delta(\bm{k})+\frac{178}{105}\delta^{2}(\bm{k})+\frac{4}{35}\mathcal{G}_{2}(\bm{k})\right]\int_{\bm{p}}P_{\rm
shell}^{\Psi}(p)\,,$ $\displaystyle\quad+a_{2}f_{\rm
NL}\frac{8}{105}\operatorname{Tr}\left[\Psi\Pi^{[1]}\right](\bm{k})\int_{\bm{p}}P_{\rm
shell}(p)+\text{h.d.}$ (3.42) $\displaystyle\left[\mathcal{S}^{2}_{\rm
Tr[\Psi\Pi^{[2]}]}\right]_{\rm free}(\bm{k})$ $\displaystyle=a_{2}f_{\rm
NL}\left[\frac{34}{21}\delta(\bm{k})+\frac{14347}{6027}\delta^{2}(\bm{k})-\frac{241}{735}\mathcal{G}_{2}(\bm{k})\right]\int_{\bm{p}}P_{\rm
shell}^{\Psi}(p)$ $\displaystyle\quad+a_{2}f_{\rm
NL}\frac{58}{105}\operatorname{Tr}\left[\Psi\Pi^{[1]}\right](\bm{k})\int_{\bm{p}}P_{\rm
shell}(p)+\text{h.d.}\,,$ (3.43)

for the contribution of the non-Gaussian operators via the free term in Eq.
(3.15). We again note the same kind of features pointed out for the spin-0
scenario, such as the need to include non-Gaussian operators in Eq. (2.3) to
have a basis that is complete under renormalization. Moreover, differently
than for the spin-0 case, in which we have a first-order non-Gaussian
operator, we find that spin-2 non-Gaussianities start to contribute at second-
order via $\operatorname{Tr}\left[\Psi\Pi^{[1]}\right]$.

Putting everything together, we find the RGE for the Gaussian operators to be

$\displaystyle\frac{db_{\delta}}{d\Lambda}$
$\displaystyle=\frac{db_{\delta}}{d\Lambda}\Big{|}_{\rm free}-a_{2}f_{\rm
NL}\left[\frac{34}{21}b_{\Psi\Pi^{[1]}}+b_{\delta\Psi\Pi^{[1]}}+\frac{34}{21}b_{\Psi\Pi^{[2]}}\right]$
$\displaystyle\hskip 199.16928pt\times\left(\frac{\Lambda}{k_{\rm
NG}}\right)^{\Delta}\left(\frac{H_{0}}{\Lambda}\right)^{2}\frac{3\;\Omega_{m}}{2\;T(\Lambda)}\frac{d\sigma_{\Lambda}^{2}}{d\Lambda}\,,$
(3.44) $\displaystyle\frac{db_{\delta^{2}}}{d\Lambda}$
$\displaystyle=\frac{db_{\delta^{2}}}{d\Lambda}\Big{|}_{\rm free}-a_{2}f_{\rm
NL}\left[\frac{124}{315}b_{\Psi\Pi^{[1]}}+\frac{178}{105}b_{\delta\Psi\Pi^{[1]}}+\frac{14347}{6027}b_{\Psi\Pi^{[2]}}+b^{\\{\delta^{2}\\}_{\rm
NG}}_{n=4}\right]$ $\displaystyle\hskip
199.16928pt\times\left(\frac{\Lambda}{k_{\rm
NG}}\right)^{\Delta}\left(\frac{H_{0}}{\Lambda}\right)^{2}\frac{3\;\Omega_{m}}{2\;T(\Lambda)}\frac{d\sigma_{\Lambda}^{2}}{d\Lambda}\,,$
(3.45) $\displaystyle\frac{db_{\mathcal{G}_{2}}}{d\Lambda}$
$\displaystyle=\frac{db_{\mathcal{G}_{2}}}{d\Lambda}\Big{|}_{\rm
free}-a_{2}f_{\rm
NL}\left[-\frac{661}{4410}b_{\Psi\Pi^{[1]}}+\frac{4}{35}b_{\delta\Psi\Pi^{[1]}}-\frac{241}{735}b_{\Psi\Pi^{[2]}}+b^{\\{\mathcal{G}_{2}\\}_{\rm
NG}}_{n=4}\right]$ $\displaystyle\hskip
199.16928pt\times\left(\frac{\Lambda}{k_{\rm
NG}}\right)^{\Delta}\left(\frac{H_{0}}{\Lambda}\right)^{2}\frac{3\;\Omega_{m}}{2\;T(\Lambda)}\frac{d\sigma_{\Lambda}^{2}}{d\Lambda}\,,$
(3.46)

where again $\frac{db}{d\Lambda}\Big{|}_{\rm free}$ represents the running of
Gaussian coefficients due to the free theory described in [44] and for the
non-Gaussian operator

$\displaystyle\frac{db_{\Psi\Pi^{[1]}}}{d\Lambda}$ $\displaystyle=-a_{2}f_{\rm
NL}\left[\frac{8}{105}b_{\delta\Psi\Pi^{[1]}}+\frac{58}{305}b_{\Psi\Pi^{[2]}}+b^{\\{\Psi\Pi^{[1]}\\}_{\rm
NG}}_{n=4}\right]\frac{d\sigma^{2}_{\Lambda}}{d\Lambda}$ $\displaystyle\quad-
a_{2}f_{\rm
NL}\left[\frac{64}{105}b_{\delta^{2}}+\frac{16}{15}b_{\delta\mathcal{G}_{2}}+b^{\\{\Psi\Pi^{[1]}\\}_{\rm
G}}_{n=4}\right]\left(\frac{k_{\rm
NG}}{\Lambda}\right)^{\Delta}\frac{d\sigma_{\Lambda}^{2}}{d\Lambda}\,.$ (3.47)

Those RGE are again summarized by Eq. (3.23), with the coefficients in Tab. 1.
Results for those equations will be discussed in the following section.

##  4 Results

After having derived the RGE for the spin-0 and spin-2 scenarios, we discuss
the solutions of these equations in this section. We start with the evolution
of the Gaussian bias parameters. We then focus on the evolution of the non-
Gaussian parameters, splitting the discussion into local and non-local PNG.

###  4.1 Gaussian coefficients

The evolution of Gaussian coefficients (and the truncation of the RGE) was
first developed for the free theory in [44]. Here we discuss how this is
affected by the interaction term and the new non-Gaussian bias operators,
focussing on spin-0 local PNG, since similar conclusions also hold for other
PNG types. We can write Eqs. (3.33)–(3.35) as

$\displaystyle\frac{db_{\delta}}{d\Lambda}$
$\displaystyle=-\left[\frac{68}{21}b_{\delta^{2}}(\Lambda)+b_{n=3}^{\ast\\{\delta\\}_{\rm
G}}\right]\frac{d\sigma^{2}_{\Lambda}}{d\Lambda}$ $\displaystyle\hskip
56.9055pt-a_{0}f_{\rm
NL}\left[-\frac{13}{21}b_{\Psi}+\frac{13}{21}b_{\Psi\delta}+b_{n=3}^{\ast\\{\delta\\}_{\rm
NG}}\right]\left(\frac{H_{0}}{\Lambda}\right)^{2}\frac{3\;\Omega_{m}}{2\;T(\Lambda)}\frac{d\sigma_{\Lambda}^{2}}{d\Lambda}\,;$
(4.1) $\displaystyle\frac{db_{\delta^{2}}}{d\Lambda}$
$\displaystyle=-\left[\frac{8126}{2205}b_{\delta^{2}}(\Lambda)+b_{n=3+4}^{\ast\\{\delta^{2}\\}_{\rm
G}}\right]\frac{d\sigma^{2}_{\Lambda}}{d\Lambda}$ $\displaystyle\hskip
56.9055pt-a_{0}f_{\rm
NL}\left[\frac{43}{135}b_{\Psi}+\frac{478}{135}b_{\Psi\delta}+b^{*\\{\delta^{2}\\}_{\rm
NG}}_{n=3+4}\right]\left(\frac{H_{0}}{\Lambda}\right)^{2}\frac{3\;\Omega_{m}}{2\;T(\Lambda)}\frac{d\sigma_{\Lambda}^{2}}{d\Lambda}\,;$
(4.2) $\displaystyle\frac{db_{\mathcal{G}_{2}}}{d\Lambda}$
$\displaystyle=-\left[\frac{254}{2205}b_{\delta^{2}}(\Lambda)+b_{n=3+4}^{\ast\\{\mathcal{G}_{2}\\}_{\rm
G}}\right]\frac{d\sigma^{2}_{\Lambda}}{d\Lambda}$ $\displaystyle\hskip
56.9055pt-a_{0}f_{\rm
NL}\left[-\frac{1699}{13230}b_{\Psi}+\frac{79}{2205}b_{\Psi\delta}+b^{*\\{\mathcal{G}_{2}\\}_{\rm
NG}}_{n=3+4}\right]\left(\frac{H_{0}}{\Lambda}\right)^{2}\frac{3\;\Omega_{m}}{2\;T(\Lambda)}\frac{d\sigma_{\Lambda}^{2}}{d\Lambda}\,,$
(4.3)

noticing again the different $\Lambda$ dependence of non-Gaussian sources
relative to the Gaussian sources discussed in [44]. We also use
$b^{\ast}=b(\Lambda_{\ast})$ to denote the parameters evaluated at a fixed
renormalization scale $\Lambda_{\ast}$ with

$\displaystyle b_{n=3}^{\\{\delta\\}_{\rm
G}}=3b_{\delta^{3}}-\frac{4}{3}b_{\mathcal{G}_{2}\delta}\quad$
$\displaystyle\textrm{and}\quad b_{n=3}^{\\{\delta\\}_{\rm
NG}}=2b_{\Psi\delta^{2}}-\frac{4}{3}b_{\Psi\mathcal{G}_{2}}\,,$ (4.4)
$\displaystyle b_{n=3+4}^{\\{\delta^{2}\\}_{\rm
G}}=\frac{{68}}{7}b_{\delta^{3}}-\frac{376}{105}b_{\mathcal{G}_{2}\delta}+b_{n=4}^{\\{\delta^{2}\\}}\quad$
$\displaystyle\textrm{and}\quad b_{n=3+4}^{\\{\delta^{2}\\}_{\rm
NG}}=\frac{68}{21}b_{\Psi\delta^{2}}-\frac{31}{21}b_{\Psi\mathcal{G}_{2}}+b^{\\{\delta^{2}\\}_{\rm
NG}}_{n=4}\,,$ (4.5) $\displaystyle b_{n=3+4}^{\\{\mathcal{G}_{2}\\}_{\rm
G}}=\frac{116}{105}b_{\mathcal{G}_{2}\delta}+b_{n=4}^{\\{\mathcal{G}_{2}\\}}\quad$
$\displaystyle\textrm{and}\quad b_{n=3+4}^{\\{\mathcal{G}_{2}\\}_{\rm
NG}}=\frac{20}{21}b_{\Psi\mathcal{G}_{2}}+b^{\\{\mathcal{G}_{2}\\}_{\rm
NG}}_{n=4}\,,$ (4.6)

with $b_{n=3}^{\\{O\\}}$ ($b_{n=4}^{\\{O\\}}$) representing contributions from
third(fourth)-order operators to $O$, and
$b_{n=3+4}^{\\{O\\}}=b_{n=3}^{\\{O\\}}+b_{n=4}^{\\{O\\}}$ from either Gaussian
or non-Gaussian operators. We generically refer to all those higher-order bias
parameters as $b_{\rm H.O.}$ when discussing solutions to the ODE. As pointed
out by [44], and verified below, these higher-operator terms constitute an
important source for the running of lower-order operators, but one can safely
neglect their running in the RGE.

Figure 1: We show the running of the Gaussian bias parameters $b_{\delta}$ and
$b_{\delta^{2}}$ and $b_{\mathcal{G}_{2}}$ for spin-0 and local ($\Delta=0$)
PNG. We present results for $f_{\rm NL}=10$ (dotted) and $f_{\rm NL}=50$
(dashed) and for two different initial conditions in the left and right
panels. Furthermore we assume all higher-order bias terms to be zero
($b^{*}_{\rm H.O.}=0$).

We display in Fig. 1 the evolution of the Gaussian bias operators for two
different initial conditions (left and right panels) and for different values
of $f_{\rm NL}$ for spin-0 local PNG ($\Delta=0$, discussed in more detail in
the next section). We neglect higher-order operators fixing $b^{*}_{\rm
H.O.}=0$. Notice the nonlinear dependence on $f_{\rm NL}$: while the influence
of $f_{\rm NL}=10$ on the running of the Gaussian bias parameters is modest,
it becomes significant for $f_{\rm NL}=50$. This is due to a back-reaction
effect: the Gaussian parameters sensitively affect the non-Gaussian bias
coefficients and those, in turn, contribute to the running of the Gaussian
coefficients. The overall effect of interactions therefore involves higher
powers of $f_{\rm NL}$. In other words, despite the fact that the RGE are
linear in $f_{\rm NL}$ and in the biases, their solutions are not.888A linear
analysis based on the RHS of Eqs. (3.33)–(3.35) and Eqs. (3.44)–(3.46) would
lead to an estimate for the running of the Gaussian parameters from the
renormalization scale $\Lambda_{\ast}$ to another scale $\Lambda_{\rm min}$
due to non-Gaussian parameters to be (setting $a_{L}=1$)  $\displaystyle\Delta
b_{\rm G}\Big{|}_{\rm lin}\propto$ $\displaystyle\,f_{\rm
NL}\int_{\Lambda_{\rm
min}}^{\Lambda_{\ast}}\frac{\Lambda^{2}\,d\Lambda}{2\pi^{2}}\left[\left(\frac{H_{0}}{\Lambda}\right)^{2}\frac{3\;\Omega_{m}}{2\;T(\Lambda)}\frac{d\sigma_{\Lambda}^{2}}{d\Lambda}\right]=f_{\rm
NL}\frac{3}{2}\Omega_{m}H_{0}^{2}\int_{\Lambda_{\rm
min}}^{\Lambda_{\ast}}\frac{\Lambda^{2}\,d\Lambda}{2\pi^{2}}\frac{P_{\rm
L}(\Lambda)}{\Lambda^{2}}\frac{1}{T(\Lambda)}\simeq 10^{-5}\,f_{\rm NL}\,.$
(4.7) In the last equality we assumed $\Lambda_{\ast}=0.25\,h/$Mpc and we
integrate from $\Lambda_{\rm min}=0$. This $10^{-5}$ suppression can be
misleading since it is only valid at linear-order in $f_{\rm NL}$ and neglects
the coupled ODE system. Therefore, when considering $f_{\rm NL}\gtrsim 10$,
one has to take into account the contribution of PNG to the running of
Gaussian coefficients and solve the joint set of RGE as we do throughout this
section (unless explicitly stated otherwise).

Furthermore, another point evident from Fig. 1 is that the scales in which the
PNG effect for the running is relevant are different than the scales where the
Gaussian running is important. While the latter are determined by the scale at
which $P_{\rm shell}$ varies, the dynamical range of PNGs is proportional to
$P^{\Psi}_{\rm shell}$ [see Eq. (3.20)]. For the case of $\Delta=0$ shown in
Fig. 1, we find PNG to be relevant on larger scales.

###  4.2 Local (spin-0) PNG

Figure 2: Top row: the running of the bias parameters $b_{\Psi}$ (blue) and
$b_{\Psi\delta}$ (black) for local-type spin-0 PNG for different values of
$f_{\rm NL}^{\rm loc}$ (solid and dashed) and assuming all higher-order bias
terms $b^{*}_{\rm H.O.}=0$ at $\Lambda_{\ast}$. Middle and bottom rows: the RG
flow of $b_{\Psi}$ and $b_{\Psi\delta}$ as a function of $b_{\delta^{2}}$ for
$f_{\rm NL}=1$ and $f_{\rm NL}=50$ respectively. Each line displays different
initial conditions (with all other bias parameters fixed to zero at
$\Lambda\to 0$). The arrows point towards $\Lambda\to 0$.

We now turn to the running of the non-Gaussian coefficients. For local
($\Delta=0$) spin-0 non-Gaussianity, Eqs. (3.36)–(3.37) reduce to

$\displaystyle\frac{db_{\Psi}}{d\Lambda}$ $\displaystyle=-a_{0}f_{\rm
NL}b^{*\\{\Psi\\}_{\rm
NG}}_{n=3}\frac{d\sigma^{2}_{\Lambda}}{d\Lambda}-4a_{0}f_{\rm
NL}b_{\delta^{2}}\frac{d\sigma^{2}_{\Lambda}}{d\Lambda}\,,$ (4.8)
$\displaystyle\frac{db_{\Psi\delta}}{d\Lambda}$ $\displaystyle=-a_{0}f_{\rm
NL}\left[\frac{272}{21}b_{\delta^{2}}+b^{*\\{\Psi\delta\\}_{\rm
G}}_{n=3+4}+b^{*\\{\Psi\delta\\}_{\rm
NG}}_{n=3+4}\right]\frac{d\sigma^{2}_{\Lambda}}{d\Lambda}\,,$ (4.9)

in which we fixed the higher-order operators at $\Lambda_{\ast}$ with

$\displaystyle b^{*\\{\Psi\\}_{\rm NG}}_{n=3}$
$\displaystyle=b_{\Psi\delta^{2}}^{\ast}\,,$ $\displaystyle
b^{*\\{\Psi\delta\\}_{\rm G}}_{n=3+4}$
$\displaystyle=12b_{\delta^{3}}^{\ast}-\frac{8}{3}b_{\mathcal{G}_{2}\delta}^{\ast}+b^{*\\{\Psi\delta\\}_{\rm
G}}_{n=4}\,,$ $\displaystyle b^{*\\{\Psi\delta\\}_{\rm NG}}_{n=3+4}$
$\displaystyle=\frac{68}{21}b_{\Psi\delta^{2}}^{*}+b^{*\\{\Psi\delta\\}_{\rm
NG}}_{n=4}\,.$ (4.10)

We discuss this approximation later in this section.999We highlight that these
equations can be solved analytically after neglecting the running of Gaussian
coefficients due to non-Gaussian terms. In that case, we can change variables
in Eqs. (4.8)–(4.9) writing the equations in terms of the variance
$\sigma^{2}(\Lambda)$ $\displaystyle\frac{db_{\Psi}}{d\sigma^{2}}$
$\displaystyle=-a_{0}f_{\rm NL}b^{*\\{\Psi\\}_{\rm NG}}_{n=3}-4a_{0}f_{\rm
NL}b_{\delta^{2}}\,,$ (4.11)
$\displaystyle\frac{db_{\Psi\delta}}{d\sigma^{2}}$ $\displaystyle=-a_{0}f_{\rm
NL}\left[\frac{272}{21}b_{\delta^{2}}+b^{*\\{\Psi\delta\\}_{\rm
G}}_{n=3+4}+b^{*\\{\Psi\delta\\}_{\rm NG}}_{n=3+4}\right]\,,$ (4.12) which can
be directly integrated to find a solution similar to those provided in Sec.
2.4 of [44]. However, as mentioned in Sec. 4.1, the coupling between the
Gaussian and non-Gaussian coefficients is relevant and Eqs. (4.1)–(4.3) do not
admit the same variable transformation. Therefore, we numerically solve the
entire coupled ODE system.

Figure 3: Similar to Fig. 2, but with non-zero higher-order bias operators. On
the left, we fix all higher-order bias operators $b^{*}_{\rm H.O.}=1$. On the
right we set all initial conditions $b^{*}=0$ and all higher order bias
parameters to zero, except $b^{*\\{\Psi\\}_{\rm
NG}}_{n=3}=b^{*\\{\Psi\delta\\}_{\rm G}}_{n=3+4}=b^{*\\{\Psi\delta\\}_{\rm
NG}}_{n=3+4}=10$. We notice in the right panels that even when setting the
initial values of $b_{\Psi},b_{\Psi\delta}$ to zero, a nontrivial running is
induced by non-zero higher order terms.

We display in the top panels of Fig. 2 some examples for the evolution of the
local spin-0 bias parameters with $f_{\rm NL}=1$ (solid) and $f_{\rm NL}=50$
(dashed), which are still in agreement with structure formation probes. We
fixed the initial conditions at $\Lambda_{\ast}$ and neglect all higher-order
operators $b_{\rm H.O.}=0$ ($b^{*\\{\Psi\\}_{\rm
NG}}_{n=3}=b^{*\\{\Psi\delta\\}_{\rm G}}_{n=3+4}=b^{*\\{\Psi\delta\\}_{\rm
NG}}_{n=3+4}=0$). Again, we fix $a_{0}=1$. The middle (bottom) panels show the
RG flow of the those bias parameters as a function of $b_{\delta^{2}}$ for
$f_{\rm NL}=1$ (50). Lower values of $f_{\rm NL}$ lead to a rather simple RG
structure, similarly to [44] while higher values of $f_{\rm NL}$ (bottom
panels) lead to a more complex structure in the RG flow due to the back-
reaction of the Gaussian-non-Gaussian coupling for the parameters running
discussed in Sec. 4.1.

We next turn to the impact of higher-order operators, whose coefficients we
have set to zero so far. We consider in the panels of Fig. 3 scenarios with
non-zero higher-order biases. In the left, we take all higher-order bias
parameters $b^{*}_{\rm H.O.}=1$ at $\Lambda_{\ast}\simeq 0.2h/$Mpc. We see
that the RG flow of the non-Gaussian operators present a more complex
structure (bottom panels) compared to the bottom panels of Fig. 2. The right
panels consider the scenario in which all first- and second-order bias
parameteres are zero and only $b^{*\\{\Psi\delta\\}_{\rm
G}}_{n=3+4}+b^{*\\{\Psi\delta\\}_{\rm NG}}_{n=3+4},\;b^{*\\{\Psi\\}_{\rm
NG}}_{n=3}$ is non-zero in Eq. (4.12). We find in that case that despite all
parameters starting at zero, the higher-order operators source $\Psi$ and
$\Psi\delta$ via Eq. (4.11) and Eq. (4.12). We find then substantial
differences in the bias evolution when compared to the scenario in which
higher-order bias parameters are zero (Fig. 2).

Finally, we discuss the truncation of the RGE hierarchy by neglecting the
running of higher-order bias terms. We display in Fig. 4 the solutions for
local PNG bias operators when neglecting higher-order operators (solid) and
when taking them as constant-in-$\Lambda$ factors (dashed), now setting the
initial conditions at $\Lambda_{\ast}=10^{-5}\,h/$Mpc. We notice a substantial
difference when considering those terms, especially for $b_{\Psi}$. In fact,
this is in agreement with the results of [44] and there is no reason to expect
that those contributions are suppressed. Similar to Eq. (2.40) of [44], we
also include as dotted lines in Fig. 4 the results for considering a heuristic
ansatz

$\displaystyle
b_{n=3+4}^{\\{O\\}}(\sigma^{2})=b^{*\\{O\\}}_{n=3+4}\,e^{-s^{\\{O\\}}(\sigma^{2}-\sigma^{2}_{\ast})}\,,$
(4.13)

for the evolution of higher-order bias parameters, with $s^{\\{O_{G}\Psi\\}}$
the equivalent of the shell prefactors for higher order non-Gaussian bias
parameters and $b^{*\\{O\\}}_{n=3+4}$ the higher order bias. They are fixed at
$s^{\\{O_{G}\Psi\\}}=b^{*\\{O\\}}_{n=3+4}=1$. We can see that both for
$b_{\Psi}$ and $b_{\Psi\delta}$ the difference of the dynamical ansatz Eq.
(4.13) for a constant higher-order bias operator is minimal, especially on
large scales. In summary, we find similar results to those reported in [44],
in which third and fourth-order operators are found to be relevant for the
evolution of first and second-order operators, but it is sufficient (on large
scales) to approximate those terms as $\Lambda$-independent constants. Thus,
this offers a self-consistent way to truncate the bias hierarchy within the RG
flow.

Figure 4: Running of the non-Gaussian bias parameters $b_{\Psi}$ (left panel)
and $b_{\Psi\delta}$ (right panel) neglecting the higher-order bias parameters
(solid), assuming them as constants (dashed) or considering the exponential
ansatz Eq. (4.13) (dotted). We set the initial conditions of the RGE at
$\Lambda_{\ast}=10^{-5}\,h/$Mpc. We neglected here the non-Gaussian
counterpart related to the running of the Gaussian parameters, since for
$f_{\rm NL}=1$ it is suppressed.

###  4.3 Non-local (spin-0 and spin-2) PNG

We now focus on non-local ($\Delta\neq 0$) PNG both for spin-0 and spin-2
fields. For spin-0, we have for the non-Gaussian operators

$\displaystyle\frac{db_{\Psi}}{d\Lambda}$ $\displaystyle=-a_{0}f_{\rm
NL}b^{*}_{\Psi\delta^{2}}\frac{d\sigma^{2}_{\Lambda}}{d\Lambda}-4a_{0}f_{\rm
NL}b_{\delta^{2}}\left(\frac{k_{\rm
NG}}{\Lambda}\right)^{\Delta}\frac{d\sigma_{\Lambda}^{2}}{d\Lambda}\,,$ (4.14)
$\displaystyle\frac{db_{\Psi\delta}}{d\Lambda}$ $\displaystyle=-a_{0}f_{\rm
NL}\left[\frac{272}{21}b_{\delta^{2}}+b^{*\\{\Psi\delta\\}_{\rm
G}}_{n=3+4}\right]\left(\frac{k_{\rm
NG}}{\Lambda}\right)^{\Delta}\frac{d\sigma_{\Lambda}^{2}}{d\Lambda}-a_{0}f_{\rm
NL}b^{*\\{\Psi\delta\\}_{\rm
NG}}_{n=3}\frac{d\sigma^{2}_{\Lambda}}{d\Lambda}\,,$ (4.15)

and for the spin-2 case, the running of
$\operatorname{Tr}\left[\Psi\Pi^{[1]}\right]$ is given by

$\displaystyle\frac{db_{\Psi\Pi^{[1]}}}{d\Lambda}=-a_{2}f_{\rm
NL}b^{*\\{\Psi\Pi^{[1]}\\}_{\rm
NG}}_{n=4}\frac{d\sigma^{2}_{\Lambda}}{d\Lambda}-a_{2}f_{\rm
NL}\left[\frac{64}{105}b_{\delta^{2}}+b^{*\\{\Psi\Pi^{[1]}\\}_{\rm
G}}_{n=3+4}\right]\left(\frac{k_{\rm
NG}}{\Lambda}\right)^{\Delta}\frac{d\sigma_{\Lambda}^{2}}{d\Lambda}\,.$ (4.16)

Note that now $\left(k_{\rm NG}/\Lambda\right)^{\Delta}$ enhances the running
of non-Gaussian coefficients due to Gaussian contributions for $\Delta>0$. The
term proportional solely to $d\sigma^{2}_{\Lambda}/d\Lambda$, arising due to
the free contribution in the action is therefore suppressed compared to the
enhanced term. Therefore, the running of non-Gaussian $\Delta>0$ PNG is mainly
driven by Gaussian terms due to the interaction term in the action.
Furthermore, note the extra dependence with $\Lambda$ compared to the local
type described in the former section.

Figure 5: Left: running of all the non-Gaussian parameters for $\Delta=2$ and
$f_{\rm NL}=10\,(50)$ in solid (dashed) lines for spin $L=0$ and $L=2$
operators. Right: same, but for $\Delta=3/2$. We present both spins together
for convenience, but the respective ODE are solved separately since they do
not mix.

We display the numerical solution for the $\Delta=2$, $L=0,2$ cases in the
left panel of Fig. 5, with the evolution of $b_{\Psi}$ (first-order spin-0),
$b_{\Psi\delta}$ (second-order spin-0) and $b_{\Psi\delta}$ (second-order
spin-2) bias parameters as a function of $\Lambda$. We fix $f_{\rm NL}=10$
(solid) and $f_{\rm NL}=50$ (dashed), values that are still allowed by CMB
[77] and LSS data. Note that the running of the $\Delta=2$ bias is enhanced
due to the extra $\left(k_{\rm NG}/\Lambda\right)^{\Delta}$ factor, as
compared to the local scenario described in Sec. 4.2 and Fig. 2. The RG flow
of those parameters as a function of $b_{\delta^{2}}$ are displayed in Fig. 6
for $f_{\rm NL}=10,\,50$ and $200$ for some representative initial conditions.
We find that increasing $f_{\rm NL}$ leads to a more rich RG flow structure,
since the coupling with Gaussian parameters start to be significant: we change
from a weakly coupled scenario (in which we could in principle solve
separately for Gaussian and non-Gaussian terms) to a strongly coupled scenario
in which both ODE systems have to be solved together. For this last case, we
caution that the $f_{\rm NL}^{2}$ contributions neglected in this work could
start to play a significant role. Notice that most of the trajectories in Fig.
6 converge toward an attractor trajectory as we evolve towards the IR
($\Lambda\to 0$).

Moreover, to exemplify other PNG cases, we also display in the right panel of
Fig. 5 the evolution of the bias parameters as a function of $\Lambda$ for the
$\Delta=3/2$ case for different $f_{\rm NL}$, as we did for $\Delta=2$. The
RG-flow for $\Delta=3/2$ is displayed in Fig. 7. In the $\Delta=3/2$ case, we
included an enlarged set of initial conditions, leading to trajectories
through different regions of the phase space. We see a very similar behaviour
to the $\Delta=2$ scenario. For both $\Delta=2$ and $\Delta=3/2$ we can see a
change in the main force driving the RG flow as we increase $f_{\rm NL}$: for
low $f_{\rm NL}$ the $b_{\delta^{2}}$ term on the RHS of Eq. (4.2) seems to
drive the evolution of $b_{\delta^{2}}$, but as $f_{\rm NL}$ increases the
$b_{\Psi}$ and $b_{\Psi\delta}$ contributions become relevant. Moreover, we
notice that the running of the spin-2 $b_{\Psi\Pi^{[1]}}$ parameter is
substantially slower compared to the other spin-0 parameters, since
$s_{\Psi\Pi^{[1]},\delta^{2}}=64/105$ is smaller compared to the other
coefficients in Tab. 1. In summary, we find a richer structure for the RG
evolution when compared to the solely Gaussian case of [44]. This could lead
to relevant changes in the fixed point structure of the system, a subject that
we plan to address in a future project.

Figure 6: Running of non-Gaussian parameters for $\Delta=2$ and $f_{\rm
NL}=10$ (top), $50$ (middle) and 200 (bottom). Note the change in the regime
as $f_{\rm NL}$ increases: for low $f_{\rm NL}$ the $b_{\delta^{2}}$ term on
the RHS of Eq. (4.2) drives the evolution of $b_{\delta^{2}}$, but as $f_{\rm
NL}$ increases, $b_{\Psi}$ and $b_{\Psi\delta}$ become relevant. Each line
displays different initial conditions with all other bias parameters fixed to
zero at $\Lambda\to 0$, the direction indicated by the arrows.

Figure 7: Same as Fig. 6 but for $\Delta=3/2$. Here we have added more
trajectories starting at different values of $b_{\delta^{2}}$ of the RG flow
to better illustrate the flow in the RG phase space.

It is evident that for both $\Delta=2$ and $\Delta=3/2$ cases the bias
parameters grow a few orders of magnitude due to the enhancement mentioned
above. One may then wonder about the observational consequences of such high
values for the bias parameters. Notice that the ultimate observable in the
case of the $2$-point function at tree-level is (see App. B for a complete
derivation)

$\displaystyle\left<\delta_{g}(\bm{k})\delta_{g}(\bm{k}^{\prime})\right>_{\rm
tree}^{\prime}$
$\displaystyle=(b_{\delta}^{\Lambda})^{2}\left<\delta^{(1)}_{\Lambda}(\bm{k})\delta^{(1)}_{\Lambda}(\bm{k}^{\prime})\right>^{\prime}+2a_{0}f_{\rm
NL}b_{\delta}^{\Lambda}b_{\Psi}^{\Lambda}\left<\delta^{(1)}_{\Lambda}(\bm{k})\Psi_{\Lambda}(\bm{k}^{\prime})\right>^{\prime}$
$\displaystyle\quad+\left(a_{0}f_{\rm
NL}\right)^{2}(b_{\Psi}^{\Lambda})^{2}\left<\Psi_{\Lambda}(\bm{k})\Psi_{\Lambda}(\bm{k}^{\prime})\right>^{\prime}\,.$
(4.17)

For $\Delta\gtrsim 1$, and as long as $f_{\rm NL}$ is not extremely large, the
last term is numerically suppressed. The main PNG contribution is therefore

$a_{0}f_{\rm
NL}b_{\delta}^{\Lambda}b_{\psi}^{\Lambda}\left<\delta^{(1)}_{\Lambda}(\bm{k})\Psi_{\Lambda}(\bm{k}^{\prime})\right>^{\prime}=a_{0}f_{\rm
NL}b_{\delta}^{\Lambda}b_{\Psi}^{\Lambda}\left(\frac{k}{k_{\rm
NG}}\right)^{\Delta}\frac{1}{M(k)}P_{\rm L,\Lambda}(k)\,.$ (4.18)

We aim to compare the term to the linear Gaussian term
$(b_{\delta}^{\Lambda})^{2}\left<\delta^{(1)}(\bm{k})\delta^{(1)}(\bm{k}^{\prime})\right>^{\prime}$
as a function of $\Lambda$. We focus on $\Delta=2$, for which the only
$k$-dependent part is the transfer function $T(k)$ inside $M(k)$. In that
case, if we focus on wavenumbers smaller than the equality scale, $k<k_{\rm
eq}$, we have $T(k)\sim 1$ and Eq. (4.18) becomes $k$ independent. We display
in Fig. 8 the ratios101010Notice that in that case the ratio of operators is
meant to display the same $(k/k_{\rm NG})^{\Delta}\;M^{-1}(k)$ dependence as
the one in Eq. (4.18)

$b_{\Psi}^{\Lambda}\frac{\Psi_{\Lambda}}{\delta^{(1)}_{\Lambda}}\,,\quad\textrm{and}\quad
b_{\Psi\delta}^{\Lambda}\frac{\Psi_{\Lambda}}{\delta^{(1)}_{\Lambda}}\,,$
(4.19)

to estimate the relative size of this contribution. We find that, despite the
running in the $b^{\Lambda}_{O}$ varying by a few orders of magnitude, the
relative contribution of this PNG $\Lambda$-running term only changes by a
small factor relative to the linear bias parameter. Ultimately, the $n$-point
functions are independent of $\Lambda$ by construction, and the
$\Lambda$-dependence part of $b_{\Psi}^{\Lambda}$ and
$b_{\Psi\delta}^{\Lambda}$ described in Fig. 8 is absorbed by the running of
other parameters, including $b_{\delta}^{\Lambda}$. Thus, Fig. 8 shows that
the apparently large running of $b_{\Psi}^{\Lambda}$, for example, can in fact
be absorbed by a modest running of $b_{\delta}^{\Lambda}$.

Figure 8: The running of the non-Gaussian contributions to the $n$-point
galaxy statistics $f_{\rm
NL}a_{0}b_{\Psi}\frac{\Psi_{\Lambda}}{\delta^{(1)}_{\Lambda}}$ and $f_{\rm
NL}a_{0}b_{\Psi\delta}\frac{\Psi_{\Lambda}}{\delta^{(1)}_{\Lambda}}$ evaluated
at different $\Lambda$, such that we can compare the contribution of the RG
running for those terms relative to the linear power spectra. We assume all
quantities are calculated at a wavelength $k\lesssim k_{\rm eq}$ such that
$T(k)=1$. We observe a moderate running in contrast to the “bare” bias
operators $b_{\Psi}$ and $b_{\Psi\delta}$. Notice that this running is to be
absorbed by other parameters, such that the $n$-point function becomes
$\Lambda$-independent.

Finally, we discuss the truncation of the RG equations for non-local PNG,
similarly to Sec. 4.2. We display in Fig. 9 the same scenarios considered in
Fig. 4, but now for $\Delta=2$ $L=0,2$ cases. We find a similar conclusion:
the inclusion of higher-order terms substantially changes the running of the
lower-order parameters. However, the ansatz Eq. (4.13) for the running of
those higher-order parameters does not significantly contribute to lower-order
operators. Taking those contributions as constant is therefore sufficient and
the truncation in terms of the order of the coefficients still works.

Figure 9: Similar to Fig. 4 but for the $\Delta=2$ scenario with $f_{\rm
NL}=50$. In order to solve the equations with the analytical ansatz Eq.
(4.13), we approximate the running of non-Gaussian bias as dominated by the
part that involves the enhancing factor $\propto\left(k_{\rm
NG}/\Lambda\right)^{\Delta}$ and we neglect the contribution to the evolution
controlled by just $d\sigma^{2}_{\Lambda}/d\Lambda$.

##  5 Discussion and conclusions

In this work, we have extended the RG-LSS framework of [44, 45] to include
interaction terms from primordial non-Gaussianities. We have considered a
cubic vertex Eq. (2.14), and work to linear order in $f_{\rm NL}$ throughout.
We have found interaction terms Eq. (3.9) that can be added directly to the
Lagrangian. The new interaction can be easily incorporated using the Wilson-
Polchinski framework to understand the running of bias parameters when
including PNG, providing a master formula Eq. (3.23) valid at linear order in
$f_{\rm NL}$. We identify the main diagrams [see Eq. (3.15) and App. C] from
the interaction term that contribute to the running of the bias. In addition
to the well-known squeezed limit of the interaction kernel (scattering
contribution), which has been the basis for the literature on the effect of
PNG on LSS so far, we identify a new annihilation contribution which is in a
different kinematic regime than the squeezed limit. We argue that it will be
absorbed by a new stochastic contribution to the partition function, but defer
a detailed investigation of this new contribution to future work.

Using the squeezed limit of the interaction kernel, parametrized as $\propto
f_{\rm NL}(k_{\ell}/k_{s})^{\Delta}$, we derive the RGE describing the running
of the bias parameters for operators induced by spin-0 and 2 interactions
during inflation, and show that non-Gaussian bias terms are in fact generated
by the RG flow. The RGE for the bias evolution are still linear, as in the
case for Gaussian initial conditions. We provide solutions for the bias
running for some specific scenarios (spin-0 local and equilateral and spin-2
equilateral). These solutions in fact scale non-linearly with the interaction
strength $f_{\rm NL}$. This also applies to the Gaussian bias coefficients,
i.e. those present in the bias expansion without PNG. Moreover, the running of
non-Gaussian coefficients also depends on the Gaussian terms and is enhanced
for larger values of $\Delta$.

We find that the solutions for the evolution of the coefficients present a
rich structure (see e.g. Fig. 6), much richer than the case for bias evolution
with Gaussian initial conditions, which could substantially change the
discussion around fixed points.

The RG solutions predict the value of the PNG bias coefficients, once fixed at
a given scale, at any other scale. This evolution can be compared with
measurements of the PNG bias coefficients as a function of scale (e.g. in
simulations at the field level [78]), which provides an independent validation
of the EFT description of galaxy statistics. Should PNG be detected in actual
data, this validation can also be applied there. Moreover, the evolution of
the non-Gaussian bias coefficients with scale is important to put priors on
these coefficients in field-level inference analyses of PNG [79], which
necessarily work at a finite scale [80]. Whether the RGE is able to ressum
part of diagrams improving the theoretical prediction is something we leave
for a future project.

In the future it would also be interesting to generalize the interaction terms
to include higher-order vertices, such as $g_{\rm NL}$ and higher-in-$f_{\rm
NL}$ contributions. Going to higher-order in PNG including new types of
vertices could eventually induce ($\Lambda$-dependent) corrections to $f_{\rm
NL}$ itself, similarly to what happens in the context of QFT. Also, the
structure of stochasticity when including non-Gaussianities changes, as we
point out in App. C. A related issue is the ambiguity of the kernel $K_{\rm
NL}$ beyond the squeezed limit when given only the bispectrum as constraint
[39, 59]. We leave a more in-depth study of these issues for a future work.

## Acknowledgments

HR is supported by the Deutsche Forschungsgemeinschaft under Germany’s
Excellence Strategy EXC 2094 ‘ORIGINS’. (No. 390783311). CN acknowledges
support from DAAD through the scholarship programme “Master studies for all
academic disciplines.” We thank Mathias Garny, Julia Stadler, Beatriz Tucci
and Leander Thiele for discussions.

##  Appendix A Non-Gaussianities via direct coupling with the current

In this section we will show that at linear order in $f_{\rm NL}$ our
description including an interaction term in the partition function is
equivalent to using non-Gaussian statistics in the matter density directly.
From the partition function Eq. (2.29), and writing the Gaussian measure
explicitly, we have

$\displaystyle\mathcal{Z}[J_{\rm m,\Lambda},J_{\rm
g,\Lambda}]=\int\mathcal{D}\delta^{(1)}_{\Lambda}\exp\left\\{-\frac{1}{2}\int_{\bm{k}}^{\Lambda}\frac{|\delta^{(1)}_{\Lambda}|^{2}}{P_{\rm
L,\Lambda}(k)}+S_{\rm int}\right\\}$ $\displaystyle\hskip
85.35826pt\times\exp\left\\{\int_{\bm{k}}J_{m,\Lambda}(\bm{k})\delta[\delta^{(1)}_{\Lambda}](-\bm{k})+\int_{\bm{k}}J_{g,\Lambda}(\bm{k})\delta_{g}[\delta^{(1)}_{\Lambda}](-\bm{k})\right\\}\,.$
(A.1)

Note that the $J$-independent term can be written as

$\displaystyle-\frac{1}{2}\int_{\bm{k}}^{\Lambda}\frac{|\delta^{(1)}_{\Lambda}|^{2}}{P_{\rm
L,\Lambda}(k)}+S_{\rm
int}=-\frac{1}{2}\int_{\bm{k}}^{\Lambda}\frac{\delta^{(1)}_{\Lambda}(\bm{k})}{P_{\rm
L,\Lambda}(k)}\left(\delta^{(1)}_{\Lambda}(-\bm{k})\right.$ (A.2)
$\displaystyle\hskip 85.35826pt\left.-2f_{\rm
NL}\int_{\bm{p}_{1},\bm{p}_{2}}\hat{\delta}_{\rm D}(\bm{k}-\bm{p}_{12})K_{\rm
NL}(\bm{p}_{1},\bm{p}_{2})\frac{M(k)}{M(p_{1})M(p_{1})}\delta^{(1)}_{\Lambda}(\bm{p}_{1})\delta^{(1)}_{\Lambda}(\bm{p}_{2})\right)\,.$

After applying the transformation

$\tilde{\delta}_{\Lambda}^{(1)}(\bm{k})=\delta^{(1)}_{\Lambda}(\bm{k})-f_{\rm
NL}\int_{\bm{p}_{1},\bm{p}_{2}}\hat{\delta}_{\rm D}(\bm{k}-\bm{p}_{12})K_{\rm
NL}(\bm{p}_{1},\bm{p}_{2})\frac{M(k)}{M(p_{1})M(p_{1})}\delta^{(1)}_{\Lambda}(\bm{p}_{1})\delta^{(1)}_{\Lambda}(\bm{p}_{2})\,,$
(A.3)

the inverse transformation can be written as

$\delta^{(1)}_{\Lambda}(\bm{k})=\tilde{\delta}_{\Lambda}^{(1)}(\bm{k})+f_{\rm
NL}\int_{\bm{p}_{1},\bm{p}_{2}}\hat{\delta}_{\rm D}(\bm{k}-\bm{p}_{12})K_{\rm
NL}(\bm{p}_{1},\bm{p}_{2})\frac{M(k)}{M(p_{1})M(p_{1})}\tilde{\delta}_{\Lambda}^{(1)}(\bm{p}_{1})\tilde{\delta}_{\Lambda}^{(1)}(\bm{p}_{2})+\mathcal{O}(f_{\rm
NL}^{2})\,,$ (A.4)

and we get

$-\frac{1}{2}\int_{\bm{k}}^{\Lambda}\frac{|\delta^{(1)}_{\Lambda}|^{2}}{P_{\rm
L,\Lambda}(k)}+S_{\rm
int}=-\frac{1}{2}\int_{\bm{k}}^{\Lambda}\frac{|\tilde{\delta}_{\Lambda}^{(1)}|^{2}}{P_{\rm
L,\Lambda}(k)}.$ (A.5)

Furthermore, at linear order in $f_{\rm NL}$ we see that111111The Jacobian of
this transformation does not contribute at $\mathcal{O}(f_{\rm NL})$.

$\int\mathcal{D}\tilde{\delta}_{\Lambda}^{(1)}\Big{|}\frac{\mathcal{D}\delta^{(1)}_{\Lambda}}{\mathcal{D}\tilde{\delta}_{\Lambda}^{(1)}}\Big{|}\mathcal{P}[{\tilde{\delta}_{\Lambda}^{(1)}}]=\int\mathcal{D}\tilde{\delta}_{\Lambda}^{(1)}\mathcal{P}[{\tilde{\delta}_{\Lambda}^{(1)}}]+\mathcal{O}(f_{\rm
NL}^{2})\,.$ (A.6)

In conclusion, the partition function for the field
$\tilde{\delta}_{\Lambda}^{(1)}$ will have an interaction-less $J$-dependent
form

$\mathcal{Z}[J_{\rm m,\Lambda},J_{\rm
g,\Lambda}]=\int\mathcal{D}\tilde{\delta}_{\Lambda}^{(1)}\mathcal{P}[\tilde{\delta}_{\Lambda}^{(1)}]\exp\left\\{\int_{\bm{k}}J_{\rm
m,\Lambda}(\bm{k})\delta[\tilde{\delta}_{\Lambda}^{(1)}](-\bm{k})+\int_{\bm{k}}J_{\rm
g,\Lambda}(\bm{k})\delta_{\rm
g}[\tilde{\delta}_{\Lambda}^{(1)}](-\bm{k})\right\\}\,,$ (A.7)

and the non-Gaussianities are coupled to the current directly. An example of
how an operator constructed from $\tilde{\delta}_{\Lambda}^{(1)}$ will look
like is

$\displaystyle\delta^{2}[\tilde{\delta}_{\Lambda}^{(1)}](\bm{k})=(\tilde{\delta}^{2})(\bm{k})+2f_{\rm
NL}\int_{\bm{p}_{1},\bm{p}_{2},\bm{p}_{3}}\hat{\delta}_{\rm
D}(\bm{k}-\bm{p}_{123})\left[K_{\rm
NL}(\bm{p}_{2},\bm{p}_{3})\frac{M(p_{1})}{M(p_{2})Mp_{3})}\right.$ (A.8)
$\displaystyle\hskip
199.16928pt\times\left.\tilde{\delta}_{\Lambda}^{(1)}(\bm{p}_{1})\tilde{\delta}_{\Lambda}^{(1)}(\bm{p}_{2})\tilde{\delta}_{\Lambda}(\bm{p}_{3})\right]+\mathcal{O}(f_{\rm
NL}^{2})\,.$

##  Appendix B The $n$-point statistics from the non-Gaussian partition
function

In this section we show power spectra and bispectra generated from the non-
Gaussian partition function Eq. (2.29). We start evaluating the power spectrum
taking derivatives of the partition function w.r.t the currents (either
$J_{m,\Lambda}$ or $J_{g,\Lambda}$). For the galaxy power spectrum

$\displaystyle\langle\delta_{g}(\bm{k})\delta_{g}(\bm{k}^{\prime})\rangle$
$\displaystyle=\frac{1}{\mathcal{Z}_{\Lambda}[0]}\frac{\delta^{2}\mathcal{Z}}{\delta
J_{g,\Lambda}\delta
J_{g,\Lambda}}=\int\mathcal{D}\delta^{(1)}_{\Lambda}\mathcal{P}[\delta^{(1)}_{\Lambda}]\;\delta_{g}(\bm{k})\delta_{g}(\bm{k}^{\prime})e^{S_{\rm
int}}\,.$ (B.1)

We see then that by using the bias expansion Eq. (2.21) we get

$\displaystyle\left<\delta_{g}(\bm{k})\delta_{g}(\bm{k}^{\prime})\right>_{\rm
tree}^{\prime}$
$\displaystyle=(b_{\delta}^{\Lambda})^{2}\left<\delta^{(1)}(\bm{k})\delta^{(1)}(\bm{k}^{\prime})\right>^{\prime}+2a_{0}f_{\rm
NL}b_{\delta}^{\Lambda}b_{\Psi}^{\Lambda}\left<\delta^{(1)}(\bm{k})\Psi_{\Lambda}(\bm{k}^{\prime})\right>^{\prime}$
$\displaystyle\quad+\left(a_{0}f_{\rm
NL}\right)^{2}(b_{\Psi}^{\Lambda})^{2}\left<\Psi_{\Lambda}(\bm{k})\Psi_{\Lambda}(\bm{k}^{\prime})\right>^{\prime}\,,$
(B.2)

where we notice that non-Gaussianities appear already at tree-level power
spectrum for matter tracers. We can similarly calculate the (tree-level)
matter-galaxy cross power spectrum

$\displaystyle\left<\delta_{m}(\bm{k})\delta_{g}(\bm{k}^{\prime})\right>_{\rm
tree}^{\prime}=\int\mathcal{D}\delta^{(1)}_{\Lambda}\mathcal{P}[\delta^{(1)}_{\Lambda}]\,\,\delta_{m}(\bm{k})\delta_{g}(\bm{k}^{\prime})\left[1+S_{\rm
int}+\mathcal{O}(f_{\rm NL}^{2})\right]$
$\displaystyle\quad=\int\mathcal{D}\delta^{(1)}_{\Lambda}\mathcal{P}[\delta^{(1)}_{\Lambda}]\delta^{(1)}_{\Lambda}(\bm{k})\left(b_{\delta}^{\Lambda}\delta^{(1)}_{\Lambda}(\bm{k})+f_{\rm
NL}b_{\Psi}^{\Lambda}\Psi[\delta^{(1)}_{\Lambda}](\bm{k})\right)$
$\displaystyle=\quad
b_{\delta}\left<\delta^{(1)}(\bm{k})\delta^{(1)}(\bm{k}^{\prime})\right>^{\prime}+f_{\rm
NL}b_{\Psi}\left<\delta^{(1)}(\bm{k})\Psi_{\Lambda}(\bm{k}^{\prime})\right>^{\prime}\,.$
(B.3)

We now move to the bispectra analysis, defined via three derivatives of the
partition function. For the matter bispectrum we have

$\displaystyle\langle\delta(\bm{k}_{1})\delta(\bm{k}_{2})\delta(\bm{k}_{3})\rangle$
$\displaystyle=\frac{1}{\mathcal{Z}_{\Lambda}[0]}\frac{\delta^{3}\mathcal{Z}}{\delta
J_{m,\Lambda}\delta J_{m,\Lambda}\delta J_{m,\Lambda}}$
$\displaystyle=\int\mathcal{D}\delta^{(1)}\mathcal{P}[\delta^{(1)}_{\Lambda}]\;\delta(\bm{k}_{1})\delta(\bm{k}_{2})\delta(\bm{k}_{3})e^{S_{\rm
int}}\,,$ (B.4)

and by expanding the exponential we get

$\displaystyle\langle\delta(\bm{k}_{1})\delta(\bm{k}_{2})\delta(\bm{k}_{3})\rangle=\int\mathcal{D}\delta^{(1)}_{\Lambda}\mathcal{P}[\delta^{(1)}_{\Lambda}]\;\delta^{(1)}_{\Lambda}(\bm{k}_{1})\delta^{(1)}_{\Lambda}(\bm{k}_{2})\delta^{(1)}_{\Lambda}(\bm{k}_{3})\left(1+S_{\rm
int}+\dots\right)$ $\displaystyle=f_{\rm
NL}\int_{\bm{p}_{1},\bm{p}_{2},\bm{p}_{3}}\hat{\delta}_{\rm
D}(\bm{p}_{123})P_{\rm L,\Lambda}^{-1}\left(\bm{p}_{1}\right)K_{\rm
NL}\left(\bm{p}_{3},\bm{p}_{2}\right)\frac{M(\bm{p}_{1})}{M\left(\bm{p}_{3}\right)M(\bm{p}_{2})}$
$\displaystyle\times\Big{\langle}\delta^{(1)}_{\Lambda}(\bm{k}_{1})\delta^{(1)}_{\Lambda}(\bm{k}_{2})\delta^{(1)}_{\Lambda}(\bm{k}_{3})\delta^{(1)}_{\Lambda}\left(\bm{p}_{1}\right)\delta^{(1)}_{\Lambda}\left(\bm{p}_{2}\right)\delta^{(1)}_{\Lambda}\left(\bm{p}_{3}\right)\Big{\rangle}\dots$
$\displaystyle=0+2f_{\rm NL}\frac{M(k_{1})}{M(k_{2})M(k_{3})}P_{\rm
L,\Lambda}(k_{2})P_{\rm L,\Lambda}(k_{3})\hat{\delta}_{\rm
D}(\bm{k}_{123})+\text{perms}+\mathcal{O}\left[f_{\rm
NL}^{2},\left(\delta^{(1)}_{\Lambda}\right)^{4}\right]\,,$ (B.5)

where we imply two cyclic permutations on all 3-momenta. We again see that for
the Gaussian case ($f_{\rm NL}=0$), the bispectrum vanishes at linear order of
all fields.

The galaxy bispectrum is defined as

$\displaystyle\langle\delta_{g}(\bm{k}_{1})\delta_{g}(\bm{k}_{2})\delta_{g}(\bm{k}_{3})\rangle$
$\displaystyle=\frac{1}{\mathcal{Z}_{\Lambda}[0]}\frac{\delta^{3}\mathcal{Z}}{\delta
J_{g,\Lambda}\delta J_{g,\Lambda}\delta J_{g,\Lambda}}$
$\displaystyle=\int\mathcal{D}\delta^{(1)}\mathcal{P}[\delta^{(1)}_{\Lambda}]\;\delta_{g}(\bm{k}_{1})\delta_{g}(\bm{k}_{2})\delta_{g}(\bm{k}_{3})e^{S_{\rm
int}}\,,$ (B.6)

and following a similar procedure, we find at tree-level

$\displaystyle\left<\delta_{g}(\bm{k}_{1})\delta_{g}(\bm{k}_{2})\delta_{g}(\bm{k}_{3})\right>_{\rm
tree}=\int\mathcal{D}\delta^{(1)}_{\Lambda}\mathcal{P}[\delta^{(1)}_{\Lambda}]\delta_{g}(\bm{k}_{1})\delta_{g}(\bm{k}_{2})\delta_{g}(\bm{k}_{3})(1+S_{\rm
int})$ $\displaystyle=b_{\delta}^{3}B_{111}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$
$\displaystyle+\left(b_{\delta}+a_{0}f_{\rm
NL}b_{\Psi}^{\Lambda}\left(\frac{k_{1}}{k_{\rm
NG}}\right)^{\Delta}\frac{1}{M(k_{1})}\right)\left[b_{\delta}+a_{0}f_{\rm
NL}b_{\Psi}^{\Lambda}\left(\frac{k_{2}}{k_{\rm
NG}}\right)^{\Delta}\frac{1}{M(k_{2})}\right]$
$\displaystyle\left[\frac{34}{21}b_{\delta}+2b_{\delta^{2}}+a_{0}f_{\rm
NL}b_{\Psi\delta}\left(\left(\frac{k_{1}}{k_{\rm
NG}}\right)^{\Delta}\frac{1}{M(k_{1})}+\left(\frac{k_{2}}{k_{\rm
NG}}\right)^{\Delta}\frac{1}{M(k_{2})}\right)\right]P_{\rm
L,\Lambda}(k_{1})P_{\rm L,\Lambda}(k_{2})\,.$ (B.7)

One can similar calculate cross terms for matter and galaxy.

##  Appendix C The cubic vertex

In this appendix, we study the structure of the cubic interaction defined in
Eq. (2.30), which we repeat here for convenience

$\displaystyle S_{\rm int}[\delta^{(1)}_{\Lambda}]$ $\displaystyle=f_{\rm
NL}\int_{\bm{p}_{1},\bm{p}_{2},\bm{p}_{3}}\hat{\delta}_{\rm
D}(\bm{p}_{123})K_{\rm
NL}\left(\bm{p}_{2},\bm{p}_{3}\right)\frac{M(p_{1})}{M(p_{2})M(p_{3})}\frac{\delta^{(1)}_{\Lambda}\left(\bm{p}_{1}\right)}{P_{\rm
L,\Lambda}\left(p_{1}\right)}\delta^{(1)}_{\Lambda}\left(\bm{p}_{2}\right)\delta^{(1)}_{\Lambda}\left(\bm{p}_{3}\right)$
$\displaystyle=\raisebox{0.0pt}{\parbox[c]{0.0pt}{\includegraphics[scale={.7}]{figs/diags/diag_vertex0.pdf}}}\,,$
(C.1)

in which the doubled lines represent the
$\delta^{(1)}_{\Lambda}\left(\bm{p}_{1}\right)/P_{\rm
L,\Lambda}\left(\bm{p}_{1}\right)$ leg. When evaluating the interaction at
$S_{\rm int}[\delta^{(1)}_{\Lambda}+\delta^{(1)}_{\rm shell}]$, we can use Eq.
(3.9)

$\displaystyle S_{\rm int}[\delta^{(1)}_{\Lambda}+\delta^{(1)}_{\rm
shell}]=S_{\rm int}[\delta^{(1)}_{\Lambda}]+S_{\rm int}^{(1)_{\rm
shell}}[\delta^{(1)}_{\Lambda},\delta^{(1)}_{\rm shell}]$ (C.2)
$\displaystyle\hskip 113.81102pt+S_{\rm int}^{(2)_{\rm
shell}}[\delta^{(1)}_{\Lambda},\delta^{(1)}_{\rm shell}]+S_{\rm int}^{(3)_{\rm
shell}}[\delta^{(1)}_{\rm shell}]\,,$

in which $S_{\rm int}^{(\ell)_{\rm shell}}$ refer to the number of external
shell legs $\ell$ in the interaction. We now describe each one of those terms
in more details.

###  C.1 The $S_{\rm int}^{(1)_{\rm shell}}$ term

We start from the contribution with one single leg in the shell

$\displaystyle S_{\rm int}^{(1)_{\rm
shell}}[\delta^{(1)}_{\Lambda},\delta^{(1)}_{\rm shell}]$
$\displaystyle=f_{\rm
NL}\int_{\bm{p}_{1},\bm{p}_{2},\bm{p}_{3}}\hat{\delta}_{\rm
D}(\bm{p}_{123})K_{\rm
NL}\left(\bm{p}_{2},\bm{p}_{3}\right)\frac{M(p_{1})}{M(p_{2})M(p_{3})}$ (C.3)
$\displaystyle\hskip 28.45274pt\times\left[\frac{\delta^{(1)}_{\rm
shell}(\bm{p}_{1})}{P_{\rm
shell}(p_{1})}\delta^{(1)}_{\Lambda}\left(\bm{p}_{2}\right)\delta^{(1)}_{\Lambda}\left(\bm{p}_{3}\right)+2\,\frac{\delta^{(1)}_{\Lambda}\left(\bm{p}_{1}\right)}{P_{\rm
L,\Lambda}\left(p_{1}\right)}\delta^{(1)}_{\rm
shell}\left(\bm{p}_{2}\right)\delta^{(1)}_{\Lambda}\left(\bm{p}_{3}\right)\right]\,,$
$\displaystyle=\raisebox{0.0pt}{\parbox[c]{0.0pt}{\includegraphics[scale={.7}]{figs/diags/diag_vertex1.pdf}}}+\raisebox{0.0pt}{\parbox[c]{0.0pt}{\includegraphics[scale={.7}]{figs/diags/diag_vertex1cut.pdf}}}\,.$

These diagrams vanish by the momentum constrain $\hat{\delta}_{\rm
D}(\bm{p}_{123})$ when taking two of the legs with momentum much smaller than
the third one that is on the shell.

###  C.2 The $S_{\rm int}^{(3)_{\rm shell}}$ term

We skip for now the $S_{\rm int}^{(2)_{\rm shell}}$ term, focusing on the
contribution with three external legs in the shell

$\displaystyle S_{\rm int}^{(3)_{\rm shell}}[\delta^{(1)}_{\rm shell}]$
$\displaystyle=f_{\rm
NL}\int_{\bm{p}_{1},\bm{p}_{2},\bm{p}_{3}}\hat{\delta}_{\rm
D}(\bm{p}_{123})K_{\rm
NL}\left(\bm{p}_{2},\bm{p}_{3}\right)\frac{M(p_{1})}{M(p_{2})M(p_{3})}$ (C.4)
$\displaystyle\hskip 113.81102pt\times\frac{\delta^{(1)}_{\rm
shell}(\bm{p}_{1})}{P_{\rm shell}(p_{1})}\delta^{(1)}_{\rm
shell}\left(\bm{p}_{2}\right)\delta^{(1)}_{\rm shell}\left(\bm{p}_{3}\right)$
$\displaystyle=\raisebox{0.0pt}{\parbox[c]{0.0pt}{\includegraphics[scale={.7}]{figs/diags/diag_vertex3.pdf}}}\,.$

For this diagram contractions, we have two options. First, if it is contracted
with another operator $O^{(n),(3)_{\rm shell}}$, leading to higher-order
contributions of the form $\left[\int_{\bm{p}}P_{\rm shell}(p)\right]^{2}$,
that is a two-loop term and is not relevant to the order of expansion we are
looking at in this work. The second option is to contracted with
$O^{(n),(1)_{\rm shell}}$, which forces the other two remaining shell lines to
contract among themselves internally in the vertex. This contribution
contracted with the shell from another operator gives a tadpole that evaluates
to

$\displaystyle\langle\delta^{(1)}_{\rm shell}(\bm{k})\times S_{\rm
int}^{(3)_{\rm shell}}\rangle_{\rm shell}$ $\displaystyle=$
$\displaystyle\hat{\delta}_{\rm D}(\bm{k})f_{\rm
NL}\left[M(k)\int_{\bm{p}}K_{\rm NL}(\bm{p},-\bm{p})\frac{P_{\rm
shell}(p)}{M^{2}(p)}+2\frac{P_{\rm shell}(k)}{M(k)}\int_{\bm{p}}K_{\rm
NL}(-\bm{k},\bm{p})\right]$ (C.5)
$\displaystyle=\raisebox{0.0pt}{\parbox[c]{0.0pt}{\includegraphics[scale={1.0}]{figs/diags/blob1.pdf}}}+\raisebox{0.0pt}{\parbox[c]{0.0pt}{\includegraphics[scale={1.0}]{figs/diags/blob2.pdf}}}\,,$

which is cancelled by noticing that the expected value of the single linear
density field is no longer zero

$\displaystyle\langle\delta^{(1)}(\bm{k})\rangle$
$\displaystyle=\frac{1}{\mathcal{Z}[0]}\frac{\delta\mathcal{Z}}{\delta J_{\rm
m}(\bm{k})}\Big{|}_{J_{\rm
m}=0}=\int\mathcal{D}\delta^{(1)}_{\Lambda}\mathcal{P}[\delta^{(1)}_{\Lambda}]\;\delta^{(1)}(\bm{k})e^{S_{\rm
int}}$ $\displaystyle=\hat{\delta}_{\rm D}(\bm{k})f_{\rm
NL}\left[M(k)\int_{\bm{p}}\frac{K_{\rm NL}(\bm{p},-\bm{p})}{M^{2}(p)}P_{\rm
L}(p)+2\frac{P_{\rm L}(k)}{M(k)}\int_{\bm{p}}K_{\rm
NL}(-\bm{k},\bm{p})\right]\;,$ (C.6)

where we neglected terms $\mathcal{O}(f_{\rm NL}^{2})$. We can then change the
definition of the linear density to subtract this part as

$\delta^{(1)}(\bm{k})\to\delta^{(1)}(\bm{k})-\langle\delta^{(1)}(\bm{k})\rangle\,.$
(C.7)

After this subtraction, $S_{\rm int}^{(3)_{\rm shell}}$ only provides higher-
order in $\lambda$ contributions.

###  C.3 The $S_{\rm int}^{(2)_{\rm shell}}$ term

Therefore, we are left with the $S_{\rm int}^{(2)_{\rm shell}}$ diagrams,
which as we will see leads to the single relevant contribution to the RG flow.
There are of two types of $S_{\rm int}^{(2)_{\rm shell}}$ diagrams, either
with a cut in a shell propagator or with a cut in a linear external leg:

$\displaystyle S_{\rm int}^{(2)_{\rm
shell}}[\delta^{(1)}_{\Lambda},\delta^{(1)}_{\rm shell}]$
$\displaystyle=f_{\rm
NL}\int_{\bm{p}_{1},\bm{p}_{2},\bm{p}_{3}}\hat{\delta}_{\rm
D}(\bm{p}_{123})K_{\rm
NL}\left(\bm{p}_{2},\bm{p}_{3}\right)\frac{M(p_{1})}{M(p_{2})M(p_{3})}$
$\displaystyle\hskip 28.45274pt\times\left[2\,\frac{\delta^{(1)}_{\rm
shell}(\bm{p}_{1})}{P_{\rm shell}(p_{1})}\delta^{(1)}_{\rm
shell}\left(\bm{p}_{2}\right)\delta^{(1)}_{\Lambda}\left(\bm{p}_{3}\right)+\frac{\delta^{(1)}_{\Lambda}\left(\bm{p}_{1}\right)}{P_{\rm
L,\Lambda}\left(p_{1}\right)}\delta^{(1)}_{\rm
shell}\left(\bm{p}_{2}\right)\delta^{(1)}_{\rm
shell}\left(\bm{p}_{3}\right)\right]$ $\displaystyle\equiv\left[S_{\rm
int}^{(2)_{\rm shell}}\right]_{\rm scatter}+\left[S_{\rm int}^{(2)_{\rm
shell}}\right]_{\rm annihilation}$ (C.8)
$\displaystyle=\raisebox{0.0pt}{\parbox[c]{0.0pt}{\includegraphics[scale={.7}]{figs/diags/diag_vertex2.pdf}}}+\raisebox{0.0pt}{\parbox[c]{0.0pt}{\includegraphics[scale={.7}]{figs/diags/diag_vertex2cut.pdf}}}\,.$

The $\left[S_{\rm int}^{(2)_{\rm shell}}\right]_{\rm scatter}$ diagram with a
cut in the shell propagator has the arguments of the $K_{\rm NL}$ kernel in
the squeezed limit. It therefore represents the scattering of a UV mode by a
IR mode. This term will be the main contribution to the running of non-
Gaussian terms.

The $\left[S_{\rm int}^{(2)_{\rm shell}}\right]_{\rm annihilation}$ diagram,
with a double line in the external propagator, was already mentioned in [22]
where it was shown that it couples to the stochastic part of the energy
momentum tensor. When dealing with the bias, the literature has so far focused
on the squeezed limit of the interaction kernel (e.g., [35]) as the relevant
limit for the galaxy statistics and neglected this term. We formally show here
(up to our knowledge) for the first time that this term leads to a stochastic
contribution to galaxy clustering. This is to some extent expected since this
contribution comes from the annihilation of two shell (small scale) modes into
a large-scale mode.121212For the $\left[S_{\rm int}^{(2)_{\rm
shell}}\right]_{\rm scatter}$ diagram, only one of the shell (UV) modes couple
through the kernel. Here the arguments of $K_{\rm NL}$ have the same orders of
magnitude, which does not allow us to use the squeezed limit expansion of this
kernel.

Given that it is sufficient to consider the $S_{\rm int}^{(2)_{\rm shell}}$
terms, we then start from the partition function

$\displaystyle\mathcal{Z}[J_{\Lambda}]$
$\displaystyle=\int\mathcal{D}\delta^{(1)}_{\Lambda}\mathcal{P}[\delta^{(1)}_{\Lambda}]\int\mathcal{D}\delta^{(1)}_{\rm
shell}\mathcal{P}[\delta^{(1)}_{\rm shell}]\exp\left(\left[S_{\rm
int}^{(2)_{\rm shell}}\right]_{\rm scatter}+\left[S_{\rm int}^{(2)_{\rm
shell}}\right]_{\rm annihilation}\right)$ $\displaystyle\hskip
76.82234pt\times\exp\left(\int_{\bm{k}}J_{\Lambda}(\bm{k})\sum_{O}b_{O}^{\Lambda^{\prime}}O[\delta^{(1)}_{\Lambda}+\delta^{(1)}_{\rm
shell}](-\bm{k})+\mathcal{O}\left[J^{2}\right]\right)$
$\displaystyle=\int\mathcal{D}\delta^{(1)}_{\Lambda}\int\mathcal{D}\delta^{(1)}_{\rm
shell}$
$\displaystyle\quad\exp\left(-\frac{1}{2}\int_{\bm{k}}\frac{\left(\delta^{(1)}_{\Lambda}(\bm{k})\right)^{2}}{P_{\rm
L,\Lambda}(k)}-\frac{1}{2}\int_{\bm{k}}\frac{\left(\delta^{(1)}_{\rm
shell}(\bm{k})\right)^{2}}{P_{\rm shell}(k)}+\left[S_{\rm int}^{(2)_{\rm
shell}}\right]_{\rm scatter}+\left[S_{\rm int}^{(2)_{\rm shell}}\right]_{\rm
annihilation}\right)$
$\displaystyle\quad\times\exp\left(\int_{\bm{k}}J_{\Lambda}(\bm{k})\sum_{O}b_{O}^{\Lambda^{\prime}}O[\delta^{(1)}_{\Lambda}+\delta^{(1)}_{\rm
shell}](-\bm{k})+\mathcal{O}\left[J^{2}\right]\right)\,,$ (C.9)

where we wrote the kinetic terms explicitly. The $\left[S_{\rm int}^{(2)_{\rm
shell}}\right]_{\rm scatter}$ integral, as said before, contributes to the
squeezed limit corrections and we calculate them in App. D. Our goal is to
understand the effect of the $\left[S_{\rm int}^{(2)_{\rm shell}}\right]_{\rm
annihilation}$ in the RG-flow. In order to gain more physical insights about
this term we will apply a similar transformation as the one described in App.
A:

$\displaystyle\tilde{\delta}_{\Lambda}^{(1)}(\bm{k})$
$\displaystyle=\delta^{(1)}_{\Lambda}(\bm{k})-f_{\rm
NL}M(k)\int_{\bm{p}_{2},\bm{p}_{3}}\hat{\delta}_{\rm
D}(\bm{k}-\bm{p}_{23})\frac{K_{\rm
NL}(\bm{p}_{2},\bm{p}_{3})}{M(p_{2})M(p_{3})}\delta^{(1)}_{\rm
shell}(\bm{p}_{2})\delta^{(1)}_{\rm shell}(\bm{p}_{3})\,,$
$\displaystyle\tilde{\delta}_{\rm shell}^{(1)}(\bm{k})$
$\displaystyle=\delta^{(1)}_{\rm shell}(\bm{k})-f_{\rm
NL}M(k)\int_{\bm{p}_{2},\bm{p}_{3}}\hat{\delta}_{\rm
D}(\bm{k}-\bm{p}_{23})\frac{K_{\rm
NL}(\bm{p}_{2},\bm{p}_{3})}{M(p_{2})M(p_{3})}\delta^{(1)}_{\rm
shell}(\bm{p}_{2})\delta^{(1)}_{\rm shell}(\bm{p}_{3})\,,$ (C.10)

with the integration measure of the path-integral [see Eq. (A.6)]

$\displaystyle\mathcal{D}\delta^{(1)}_{\Lambda}$
$\displaystyle=\mathcal{D}\tilde{\delta}^{(1)}_{\Lambda}+\mathcal{O}(f_{\rm
NL}^{2})\,,$ $\displaystyle\mathcal{D}\delta^{(1)}_{\rm shell}$
$\displaystyle=\mathcal{D}\tilde{\delta}_{\rm shell}^{(1)}+\mathcal{O}(f_{\rm
NL}^{2})\,.$ (C.11)

We notice that with this transformation the annihilation part is absorbed in
the field redefinition

$\displaystyle-\frac{1}{2}\int_{\bm{k}}\frac{\left(\delta^{(1)}_{\Lambda}(\bm{k})\right)^{2}}{P_{\rm
L,\Lambda}(k)}+\left[S_{\rm int}^{(2)_{\rm shell}}\right]_{\rm
annihilation}=-\frac{1}{2}\int_{\bm{k}}\frac{\left(\tilde{\delta}^{(1)}_{\Lambda}\right)^{2}(\bm{k})}{P_{\rm
L,\Lambda}(k)}\,,$ (C.12)

while the scattering term will remain unchanged at linear order in $f_{\rm
NL}$

$\left[S_{\rm int}^{(2)_{\rm shell}}\right]_{\rm
scatter}\left(\delta^{(1)}_{\Lambda},\delta^{(1)}_{\rm
shell}\right)=\left[S_{\rm int}^{(2)_{\rm shell}}\right]_{\rm
scatter}\left(\tilde{\delta}_{\Lambda}^{(1)},\tilde{\delta}_{\rm
shell}^{(1)}\right)\,.$ (C.13)

The kinetic term with the shell transforms as

$\displaystyle-\frac{1}{2}\int_{\bm{k}}\frac{\left(\delta^{(1)}_{\rm
shell}(\bm{k})\right)^{2}}{P_{\rm
shell}(k)}=-\frac{1}{2}\int_{\bm{k}}\frac{\left(\tilde{\delta}_{\rm
shell}^{(1)}\right)^{2}(\bm{k})}{P_{\rm shell}(k)}-S_{\rm int}^{(3)_{\rm
shell}}(\tilde{\delta}_{\rm shell}^{(1)})\,,$ (C.14)

and the remaining three-shell interaction will exactly cancel the interaction
vertex Eq. (C.4). After the field redefinition and grouping the terms together
we have

$\displaystyle\mathcal{Z}[J_{\Lambda}]$
$\displaystyle=\int\mathcal{D}\tilde{\delta}^{(1)}_{\Lambda}\mathcal{P}[\tilde{\delta}^{(1)}_{\Lambda}]\mathcal{D}\tilde{\delta}^{(1)}_{\rm
shell}\mathcal{P}[\tilde{\delta}^{(1)}_{\rm shell}]\exp\left(\left[S_{\rm
int}^{(2)_{\rm shell}}\right]_{\rm scatter}\right)$ $\displaystyle\hskip
28.45274pt\times\exp\left(\int_{\bm{k}}J_{\Lambda}(\bm{k})\sum_{O}b_{O}^{\Lambda^{\prime}}O[\tilde{\delta}^{(1)}_{\Lambda}+\tilde{\delta}^{(1)}_{\rm
shell}](-\bm{k})+\mathcal{O}(J_{\Lambda}^{2})\right)\,,$ (C.15)

in which we observe that the annihilation interaction is absorbed by the field
transformation leaving only the scattering vertex. Furthermore, the terms
coupled to the current now contain a non-Gaussian contribution: we have to
solve the inverse transformation of Eq. (C.3) to get $\delta^{(1)}$ in terms
of $\tilde{\delta}^{(1)}$ at linear order in $f_{\rm NL}$ and the resulting
operators coupled to the current will then be of the form
$O[\delta^{(1)}]=O[\delta^{(1)}[\tilde{\delta}^{(1)}]]\equiv
O[\tilde{\delta}^{(1)}]$. Since we work only at linear order in $f_{\rm NL}$
the effect of the scattering vertex is the same as before the transformation
because we only consider coupling with the Gaussian counterpart of the
$\tilde{\delta}_{\Lambda}^{(1)}$ [couplings with the non-Gaussian will give
terms of $\mathcal{O}(f_{\rm NL}^{2})$].

We can immediately see that the field redefinition Eq. (C.3) will not change
the shell corrections $[\mathcal{S}_{O}^{2}]_{\rm free}$ defined in Eq. (3.15)
since

$\displaystyle[\mathcal{S}_{\delta^{2}}^{2}]_{\rm free}\supset
M(0)\int_{\bm{p}}\frac{K_{\rm NL}(\bm{p},-\bm{p})}{M^{2}(p)}P_{\rm
shell}(p)=0\,,$ (C.16)

and therefore the annihilation term does not source any $J^{1}$ corrections to
the bias parameters. Despite not changing the result for the $J^{1}$ terms
calculated in this work, the annihilation interactions induce the running of
stochastic coefficients (see [45] for a description of those terms for the
Gaussian case). As seen in [45], the main source of stochasticity (or $J^{m}$
terms with $m>1$) come from shell correlators such as

$\displaystyle\mathcal{S}_{O_{1},O_{2}}^{22}[\delta^{(1)}_{\Lambda}](\bm{k}_{1},\bm{k}_{2})=\sum_{n_{1}\geq
2,n_{2}\geq 2}\left\langle O_{1}^{(n_{1}),(2)_{\rm
shell}}[\delta^{(1)}_{\Lambda},\delta^{(1)}_{\rm
shell}](\bm{k}_{1})O_{2}^{(n_{2}),(2)_{\rm
shell}}[\delta^{(1)}_{\Lambda},\delta^{(1)}_{\rm
shell}](\bm{k}_{2})\right\rangle_{\rm shell}\,.$ (C.17)

We have then, for instance, a new non-Gaussian contribution

$\displaystyle\mathcal{S}_{\delta\delta^{2}}^{22}(\bm{k},\bm{k}^{\prime})=M(k)\hat{\delta}_{\rm
D}(\bm{k}+\bm{k}^{\prime})\int_{\bm{p}}\frac{K_{\rm