A Supersymmetric Hierarchical Model for Weakly Disordered 3d Semimetals

In this paper, we study a hierarchical supersymmetric model for a class of gapless, three-dimensional, weakly disordered quantum systems, displaying pointlike Fermi surface and conical intersections of the energy bands in the absence of disorder. We use rigorous renormalization group methods and supersymmetry to compute the correlation functions of the system. We prove algebraic decay of the two-point correlation function, compatible with delocalization. A main technical ingredient is the multiscale analysis of massless bosonic Gaussian integrations with purely imaginary covariances, performed via iterative stationary phase expansions.


Introduction
An important conjecture in mathematical quantum mechanics is that disordered, noninteracting, 3d quantum systems display a localization/delocalization transition as a function of the disorder strength [1,7]. The simplest model that is expected to give rise to such transition is the Anderson model, described by a random Schrödinger operator with −Δ the lattice Laplacian and V ω a random potential, e.g., (V ω ψ)(x) = ω(x)ψ(x) with {ω(x)} x∈Z 3 i.i.d. random variables with variance O (1). From a mathematical viewpoint, a lot is known about this problem for strong disorder, |γ| 1. There, one expects wave packets not to spread in time, and transport to be suppressed (zero conductivity). This phenomenon has been rigorously understood for general d-dimensional models starting from the seminal work [41], where a KAM-type multiscale analysis approach to localization was developed, and later via the fractional moments method [3]. See [6] for a pedagogical review of mathematical results on Anderson localization.
Instead, for small disorder much less is known from a rigorous viewpoint. In three dimensions, one expects nontrivial transport and an emergent diffusive behavior of the quantum dynamics. Unfortunately, so far no fully satisfactory rigorous result is available on this problem. Results have been obtained for tree graphs and similar structures [4][5][6]13,14,40,51,52,69]. The analogous problem for random matrix models is much better understood; see [32] for a review of recent results. Concerning short-ranged lattice models, important progress has been obtained in [33][34][35], where diffusion for the Anderson model has been proven in the scaling limit, and in [27,28], where a localization/delocalization transition for a supersymmetric effective model has been established. (See also [24] for more recent extensions.) The starting point of [27,28] is a mapping of the disorder-averaged correlations of the Anderson model into those of an interacting supersymmetric quantum field theory model. This mapping was first introduced in physics in [30] (see also [75], for a related approach based on the replica trick) and allows to import field-theoretic methods to study random Schrödinger operators. Let us briefly describe it. Consider a general class of random Schrödinger operators, H ω = H + γV ω , with H a short-ranged lattice Schrödinger operator, on a finite sublattice Λ of Z 3 . Let G ω (x, y; μ − iε) be the Green's function: The parameter μ ∈ R plays the role of chemical potential, while ε > 0 introduces a regularization of the Green's function. The relevant setting for weakly disordered metals is μ ∈ σ(H) (as L → ∞). It is well known that the Green's function can be represented as the covariance of a Gaussian Grassmann field, as follows: with C −1 ω := −i(H ω − μ) + ε; the reason for the multiplication by the trivial factor i will be clear in a moment. The denominator is the determinant of the matrix C −1 ω , which is a random object; as a consequence, the expression (1.3) is not very useful for the purpose of computing the disorder average.
The key remark is that the reciprocal of a determinant can be written as a complex Gaussian integral: , the inverse covariance is C −1 = −i(H − μ) + ε and λ = γ 2 /2. The same trick can be applied to rewrite the average of the product of Green's functions, by introducing internal degrees of freedom for the fields, labeling different copies of G ω . Internal degrees of freedom for H (e.g., spin or sublattice labels) can also be taken into account in a similar way.
Equation (1.5) is an exact formula for the averaged Green's function of the model on a finite volume. It allows to recast the problem of computing the averaged Green's function for a random Schrödinger operator into a statistical mechanics/quantum field theory problem. The factor i allows to circumvent the fact that the operator H −μ need not be positive. Moreover, the parameter ε > 0 allows to avoid singularities in the determinant at the denominator and to make sense of the complex Gaussian integrals. The problem we now have to face is to construct this interacting quantum field theory model for λ small, uniformly in the volume of the system and as ε → 0 + .
A formal approach often adopted in the physics literature is to perform a saddle point analysis for the Gaussian superfield Φ ± = (φ ± , ψ ± ); see [31] for a review. As a result, one obtains remarkable predictions about the behavior of the systems, such as the emergence of random matrix statistics for the eigenvalue distribution of H ω . Making this strategy rigorous, however, presents very serious mathematical challenges, which so far have been rigorously tackled only for a class of effective supersymmetric models, [23,27,28], or in the context of random matrix models (mean field regime) [71,72].
Another possibility, less explored from a rigorous viewpoint, is to apply rigorous renormalization group (RG) methods to construct the Gibbs state of the interacting supersymmetric model for small λ, that is to evaluate the integral in Eq. (1.5) via a convergent multiscale analysis. Similar methods have been recently used in [59], for an analysis at all orders in renormalized perturbation theory of the correlation functions of an effective supersymmetric model of graphene in the presence of random gauge fields. See also [15,56,57,68] for earlier approaches to disordered systems via a combination of RG and random matrix techniques. For quantum systems with quasi-random disorder, rigorous RG techniques have been used to prove the existence of localization in the ground state of the interacting fermionic chains [62].
In the present context, the most serious difficulties one has to face in a nonperturbative application of RG methods are: (i) the large field problem, due to the unboundedness of the bosonic fields; (ii) the infrared problem, which arises whenever μ lies in the spectrum of H; (iii) the presence of a purely imaginary covariance for the bosonic integration. The goal of this paper is to present a rigorous solution to these problems, in a simple yet nontrivial case. As usual in condensed matter physics models, the geometry of the Fermi surface determines how severe are the infrared divergences appearing in the naive perturbative expansion. In particular, in the context of interacting fermionic systems, the rigorous study of the ground state of models with general extended Fermi surfaces is so far out of the reach of the existing rigorous RG methods. Important progress has been achieved in [25,26], for the low-temperature construction of jellium, in [17], for lowtemperature analysis of the 2d Hubbard model on the square lattice, and in [36], for the Fermi liquid construction of 2d models with asymmetric Fermi surface.
Here, we shall consider a class of 3d quantum systems with pointlike Fermi surface; these are models for Weyl semimetals; see [8] for a review. Weyl semimetals are a class of recently discovered condensed matter systems [50] that might be thought as a 3d generalization of graphene. In these models, the (translation invariant) Schrödinger operator H can be written as H = ⊕ T 3 dkĤ(k), withĤ(k) the Bloch Hamiltonian of the model, k the quasimomentum of the particle, and T 3 the Brillouin zone. The energy bands of H(k) display conical intersections at the Fermi level μ, at a finite number of Fermi points, also called Weyl nodes, k = k α F , α = 1, . . . , 2M . As a consequence of this fact, one has, up to oscillating prefactors: It turns out that the reduced dimensionality of the Fermi surface allows to use RG methods to construct the low/zero temperature interacting Gibbs state of the model, both in two (corresponding to graphene-like systems) and in three dimensions, and to prove universality results for transport coefficients; see [45][46][47][48][49]60,61]. We refer the reader to [58] for a review of recent applications of rigorous RG methods interacting condensed matter systems.
Here, we shall focus on three-dimensional disordered Weyl semimetals, in the presence of weak disorder and no interactions. Heuristic perturbative analysis suggests that disorder is irrelevant in the renormalization group sense; see [37,38,55] for early renormalization group approaches to disordered semimetals. Nevertheless, the emergence of localized states at the Weyl points for weak disorder has been recently proposed in [67] and then challenged in [21]; see also [22] and references therein for numerical simulations, showing that the delocalized phase in Weyl semimetals with well-separated nodes is robust against weak disorder. From a rigorous viewpoint, the stability of the Weyl phase against weak quasi-periodic external potentials has been recently proved in [63]. In this paper, we shall consider a SUSY hierarchical approximation for disordered Weyl semimetals: The connection between the hierarchical model and the original lattice model lies in the scaling of the superfield covariance, which will be chosen so to match the decay properties of the massless Green's function of the original model, Eq. (1.6). Let us also point out that our model describes an interacting SUSY field associated with one Weyl node: from the point of view of Weyl semimetals, it is relevant for the description of random potentials that do not couple different Fermi points.
Hierarchical models played an important role in the development of rigorous RG methods [16,19,29,42,43]. For instance, we mention the study of the hierarchical ϕ 4 4 theory [43], which paved the way to the construction of the full lattice ϕ 4 4 theory [44]. The connection between the two models is provided by a cluster expansion [44], technically similar to a high-temperature expansion in classical statistical mechanics. See also [10,11] for a recent extension of this result to SUSY ϕ 4 4 , relevant for the study of the weakly self-avoiding walk, and [12] for a detailed discussion of the hierarchical approximation of the model.
Hierarchical models have also been considered in the context of random Schrödinger operators [20,53,54,65,73,74]; see also [64,66] for discussions about the connection with the Anderson localization/delocalization transition. There, the model is defined on a one-dimensional lattice, and the range of the hierarchical hopping is tuned to fix the effective dimension of the system. The works [53,54,65,73,74] prove that, as long as the hopping is summable, the model is in the localized phase.
In this paper, we rigorously construct the SUSY hierarchical version of 3d Weyl semimetals with well-separated Weyl nodes, and we prove algebraic decay of correlations; the decay exponents are the same as those of the nondisordered model. Our RG analysis is inspired by the block spin transformation of [43,44]; in particular, the study of the bosonic sector of the theory is performed thanks to the careful control of the growth of the analyticity domain of the effective action as a function of the complex bosonic field, and to the iteration of suitable analyticity (Cauchy) estimates. With respect to [43], an important simplification in our case is due to the fact that the interaction (hence the disorder) is irrelevant in the renormalization group sense. However, in contrast to [43], the Gaussian covariances are purely imaginary, which means that the single step of RG has to be performed exploiting oscillations. Also, one has to deal with the extra presence of fermionic fields, which makes the analysis considerably more involved with respect to a purely bosonic theory.
A key role in our construction is played by supersymmetry, that allows to reduce the number of running coupling constants and to prove the equality (up to a sign) of fermionic and bosonic correlations. If combined with a suitable cluster expansion, we expect our result to extend to the full lattice model; we postpone this study to future work. The first application of cluster expansion techniques and RG methods to QFT models with complex covariances we are aware of is the work [9], on the construction of the ultraviolet sector of interacting three-dimensional lattice bosonic systems. Recently, a cluster expansion for supersymmetric lattice models with imaginary covariance, relevant for disordered quantum systems, has been developed by one of us in [39]. This expansion allows to study the localization regime, for strong disorder and for all energies, or for weak disorder and away from the unperturbed spectrum. We expect the combination of these tools with multiscale analysis to shed light on delocalized phases, at least in cases in which the disorder is irrelevant in the renormalization group sense.
The paper is organized as follows. In Sect. 2, we introduce the model we will study, and we will state our main result, Theorem 2.2. In Sect. 3, we develop the RG method that we will first apply to the construction of the effective potential of the theory. Then, in Sect. 4 we apply this strategy to the computation of the two-point function of the model, which allows to prove our main result. Finally, in Appendix A we discuss the flow of the counterterm fixing the choice of the interacting chemical potential; in Appendix B, we prove some key technical results, while in Appendix C we discuss the (super) symmetries of the model.

The Hierarchical Gaussian Superfield
Let N ∈ N, L ∈ 2N. Let Λ ⊂ N 3 be the set: Later, it will be convenient to look at Λ as being covered by disjoint blocks B (1) z of side L and labeled by z ∈ Λ (1) : More generally, for any 1 ≤ k ≤ N , we set Λ (k) := L −k Λ ∩ N 3 . Obviously, Λ (k+1) ⊂ Λ (k) . We set, for any z ∈ Λ (k) , with the understanding Λ (0) ≡ Λ: where a denotes the vector in N 3 which approximates a from below. Notice that the box B (2) L −2 x and so on. In this way, one defines a hierarchy of boxes, where a box at a given scale contains the lattice points of the previous scale. We are now ready to introduce the hierarchical Gaussian superfield. Roughly speaking, one associates with any point x ∈ Λ a sum of independent Gaussian variables, each of them corresponding to a box within the hierarchy described above. We will follow the definition of the hierarchical model of [43], which captures the main features of the multiscale decomposition of the full lattice Gaussian field [44]. We define a complex Gaussian field φ x,σ and a pair of Grassmann Gaussian fields ψ + x,σ , ψ − x,σ as follows, for all x ∈ Λ, and for all spins σ =↑↓: . The smallest squares constitute the blocks B (1) , which, if identified with their labels, form the lattice Λ (1) . The smallest blocks containing B (1) are the blocks B (2)  , whose covariance will be defined in Eq. (2.10). The function A is the scale-independent local version of the kernels that appear in the blockspin transformation [44]. It is chosen so that A y = ±1, for all y ∈ Λ (0) , with (1) , and it is invariant under translations by L in each direction.
Recall the defining properties of Grassmann variables 1 : Setting φ − := φ and φ + := φ, we shall collect both complex and Grassmann fields in a single Gaussian superfield Φ ± x,σ : ψ,x,σ ). In terms of those, we rewrite the superfield Φ ± as: Equation (2.7) defines the hierarchical Gaussian superfield. Notice that the term in the above sum labeled by a given h varies on the length scale L h , and it is of size L −h . In particular, replacing the above sum by a truncated sum starting from the scale k, one obtains a field that varies on length scale L k , and it is of size L −k . For later convenience, it is useful to rescale this truncated field, in such a way that it varies on scale 1, and it is of size 1, for all k ≥ 0. We introduce, for all x ∈ Λ (k) : satisfies the recursion relation: This decomposition has a clear meaning: the field Φ (≥k) x is written as the sum of a term which is constant in the block B (k+1) x/L , the average of the field in the block, plus a fluctuation with zero sum in the same block.
We shall choose the covariance of the independent single-scale superfields ζ (h) as, for = φ, ψ: (2.10) Thus, the covariance of the full superfield is: ;σ,σ (x, y) where d(x, y) is the hierarchical distance between x and y: (2.12) Notice that this covariance mimics the real space algebraic decay of the Green's function of the full lattice model, Eq. (1.6). In particular, the algebraic decay of the covariance implies that the Gaussian superfield is massless.

The Gibbs State of the Interacting Hierarchical Model
The goal of this paper is to study weak perturbations of the massless Gaussian superfield defined in the previous section. We define: The first term plays the role of the many-body interaction for the superfield, while the second term will fix the chemical potential of the system and will be suitably chosen later on. Given an analytic function P (Φ), with Φ as in Eq. (2.7), we define, for λ > 0 and |μ| ≤ C|λ|: (2.14) Fermionic integration, = ψ, is defined in Eq. (2.5), while for bosonic integration, = φ, we use the convention dζ φ,x,σ . The minus sign in front of the RHS of last equation in (2.14) is for does not provide any decay, the integral is well-defined for all λ > 0 by the presence of the quartic interaction.
The normalization factor Z N is the partition function of the model, as discussed in Appendix C, Z N = 1 by the localization theorem of supersymmetric integration, [18,70].
Remark 2.1. In the following, we shall use the symbols C, C, K,K for generic universal constants. Whenever a constant will depend on L, we shall denote it explicitly, i.e., with the symbol C L .

Main Result
In this paper, we shall construct the Gibbs state of the hierarchical SUSY model. In particular, we shall focus on the two-point correlation function; the same methods can be extended in a straightforward way to all higher-order correlations. The next theorem is our main result.
The proof of Theorem 2.2 is based on rigorous renormalization group methods. The parameter μ plays the role of bare chemical potential; technically, its choice allows to control the relevant direction of the RG flow and to prove the convergence to a Gaussian fixed point. In the original lattice model, the choice of μ would correspond to a shift of the Fermi level associated with the Weyl phase, induced by the disorder. We stated the theorem for μ ∈ C; we believe that with some extra effort one could actually prove that μ ∈ R, but we will not need this improvement in our analysis. Notice that for a general disordered lattice model μ has to be real, since otherwise the Green's function would decay exponentially fast by a Combes-Thomas estimate.
As mentioned in the Introduction, the main difficulties in the RG analysis are due to the massless covariance of the superfield, to the unboundedness of the bosonic field (large field problem) and to the fact that the covariance is purely imaginary. In order to perform the single-scale integration in the RG, one has to exploit oscillations in the Gaussian integration.
Relation with the full lattice model. Before discussing the proof, let us briefly comment on the relevance of this result for the understanding of the behavior of the full lattice model, beyond the hierarchical approximation. As discussed in the Introduction, the supersymmetric representation can be used to study the averaged resolvent of the lattice model, with Gaussian disorder: Notice that we will not be concerned with the average of the absolute value of the resolvent, which is expected to diverge as ε → 0. Our analysis focuses on a hierarchical approximation of the SUSY field-theoretical representation of the averaged Green's function. More precisely, our hierarchical field is related to the quasi-particle field associated with one conical intersection; from the point of view of the original lattice model, this would amount to introducing a momentum cutoff in the covariance of the disorder, which has the effect of keeping the Fermi points decoupled. The resolvent can be used to compute the Fermi projector in a finite volume, via functional calculus [2,6]: where C μ is the counterclockwise path in the complex plane, and where: Here, we are assuming that the point μ does not belong to the spectrum of H ω (which consists of eigenvalues, since we are in a finite volume); as discussed below, this is true with probability one. The u integration in Q 2 converges absolutely, uniformly in all disorder realizations. Moreover, a Combes-Thomas estimate [2] gives |Q 2 (H ω − μ)(x, y)| ≤ Ce −c x−y uniformly in the disorder and in the system size. Consider now the Q 1 term. Let γ be the disorder strength, recall Eq. (1.1). The following bound holds true, for a broad class of disordered lattice models in a finite volume [3] (including the finite volume version of (1.1)): This estimate implies that, uniformly in x, y and in the system size, and for some C γ,θ < ∞, see (2.12) of [3]: is absolutely integrable in η ∈ [−1; 1]. Therefore, we can interchange the η integration with the average over the disorder: Thus, one can deduce decay properties of the averaged Fermi projector starting from the decay properties of the averaged Green's function. This is exactly what we study in the present paper, in the hierarchical approximation, at the most singular point η = 0. We prove that the averaged Green's function at η = 0 decays as x−y −2 at large distances (as in the disorder-free case). More generally, for η = 0 a straightforward extension of the analysis performed in the present paper would give a decay at most as x − y −2 e −|η| x−y for the hierarchical approximation of the resolvent. Therefore, based on our result for the hierarchical approximation, the natural conjecture for the full lattice model is that: Notice that this is not what happens in the regime of strong disorder, where one has the fractional moment bound , an estimate which can be used to prove Anderson localization. (See [6] for a review.) In particular, the fractional moment bound implies that [2]: As discussed in the Introduction, in order to extend our main result, Theorem 2.2, to the full lattice model, one has to combine the RG analysis introduced in this paper with a cluster expansion; we defer this nontrivial extension to future work.

Renormalization Group Analysis: The Effective Potential Flow
In this section, we discuss the renormalization group analysis of our model. We shall start by computing the effective potential of the theory, at all scales. The main result is the expression (3.143), for the effective potential on an arbitrary scale. This will play an important role in the computation of the two-point correlation function, postponed to Sect. 4.

The Effective Potential
The Gibbs state of the model will be constructed in an iterative fashion, integrating the fluctuation fields starting from the scale h = 0 until the last scale, h = N . A key role in this iterative strategy will be played by the study of the following map (recall Eq. (2.9)): x/L ,σ and where the normalization factor N (h) is given by: The main simplification introduced in the hierarchical model with respect to the original lattice model is that, on any scale h, the argument of the integral factorizes: ). Therefore, for all x ∈ Λ (h) , one has the following local version of the map defined in Eq. (3.1): with normalization factor: The assumption A y = ±1, y∈B (1) x A y = 0 implies: (3.5) The main technical goal of this paper is the control of this map. We shall prove that, as N → ∞, the iteration (3.5) converges to a unique Gaussian fixed point, for a suitable choice of the bare chemical potential μ. In the following, we shall drop the dependence of the superfields on the scale and position labels, since they will play no role in the study of the single step of the iteration.

Integration of the Scale Zero
In this section, we will discuss the first iteration of the map defined in Eq. (3.3). The iteration of the map on later scales will be performed inductively; the correct inductive assumptions will be motivated by the discussion of this section.

Setting Up the Integration.
Given two superfields Φ ± = (φ ± , ψ ± ) and ζ ± = (ζ ± φ , ζ ± ψ ), we define: (3.6) T , the fermionic product can also be represented as: where iσ 2 = 0 1 −1 0 acts on the components of the Grassmann vector with a given spin. Also, it is easy to check that ψ · ψ = σ ψ + σ ψ − σ . Concerning the bosonic product φ · ζ φ in (3.6), notice that, setting φ ± σ = φ 1,σ ± iφ 2,σ , with φ i,σ ∈ R: (3.8) Hence, (3.8) coincides with the usual scalar product of the following vectors in R 4 : Later, we will be interested in considering the extension of Eq. (3.8) for complex φ 1,σ and φ 2,σ . In this general case, Eq. (3.8) does not define a scalar product on C 4 , and φ · φ = φ 2 = i,σ |φ i,σ | 2 : what is missing to define a scalar product is the complex conjugate on the first factor. Nevertheless, we still have |(φ · ζ φ )| ≤ φ ζ φ . The usual scalar product in C n will be denoted by ·, · , v, w = j v j w j . Setting we shall discuss the evaluation of: (3.11) Equation (3.11) defines the effective interaction of the hierarchical model on scale 1. In order to perform the integration, we explicitly rewrite: Therefore, where we introduced: Our goal will be to discuss the integration of the fluctuation superfield ζ. We shall start by integrating its fermionic component.

Integration of the Fermionic Fluctuation Field.
We rewrite the effective interaction on the scale zero as: (3.15) Explicitly: Therefore, we rewrite the integral in Eq. (3.13) as: The integration of the Grassmann field ζ ψ will be performed by writing the function exp{−V f (Φ, ζ)} as a linear combination of finitely many monomials in the ζ ψ variable, due to the fact that the Grassmann algebra is finite. Then, we evaluate the Gaussian integration using the rules of Grassmann calculus, Eq. (2.5). We get: where the functions C (0) n (φ/L, ζ φ ) are polynomials in φ/L, ζ φ , which we shall estimate. The (ψ · ψ)-dependence of the outcome of the integration follows by symmetry; see Corollary C.2. We shall proceed as follows. A general function of the supervectors Φ and ζ can be written as where the summation is over multi-indices a, b ∈ {0, 1} {±}×{↑,↓} and ψ a , ζ b ψ are the corresponding monomial (uniquely defined once an order for the set {±} × {↑, ↓} is chosen). For such multi-indices a = (a ε σ ) ε=± σ=↑↓ , we also set |a| := ε,σ a ε σ . In our setting, the coefficients of the expansion f a,b are functions of φ and ζ φ . (3.19) for some κ, N , M ≥ 0. (The convention 0 0 ≡ 1 is understood.) If the function only depends on ψ, f (ψ) = a f a ψ a , we shall say that it satisfies (κ, N )bounds if |f a | ≤ κ N |a| .
Remark 3.2. The interesting part of the definition is the power law, (i.e., analytic) structure of the bounds, which exhibits useful properties under algebraic manipulation and fermion integration, see below. Notice that many choices of κ, N , M can of course be made. In particular, since the number of coefficients to bound is finite, for any choice of N , M, a suitable κ can always be found.
Let us take the function f (Φ, ζ) to be the exponential e −V (0) f (Φ,ζ) . We would like to prove (κ, N , M)-bounds for this function. We write: where: given by a sum of products of an even number of Grassmann variables, we can write: where every term in the product can be expanded in a finite sum. To begin, we shall derive (κ, N , M)-bounds for every contribution to the product. Let us denote by (κ j , N j , M j ) the parameters of the (κ, N , M)-bound for the jth term in the product. A simple computation gives: As already pointed out in Remark 3.2, we stress that these choices are not unique. However, as it will be clear in the following, the choice N = λ 1/4 /L is consistent with the size of the "small field region" for the bosonic field, to be introduced later. Also, the value κ = 1 is natural, because it is equal to the value of all the entries in the product (3.22) for vanishing fermionic fields. Next, we would like to get a (κ, N , M)-bound for the product in Eq. (3.22). To this end, we shall use the following lemma.
The coefficients h a,b are given by the following expression: for a suitable sign sign(a 1 , b 1 ; a 2 , b 2 ) ∈ {±}, which we shall leave unspecified. By the triangle inequality and by the hypotheses on f a1,b1 and g a2,b2 : 25) and the claim is proven.
Let us go back to (3.22). By Eq. (3.23) together with Lemma 3.3, we get that exp{−V This is due to the fact that subtraction by 1 simply sets to zero the first coefficient of the Grassmann expansion, f 0,0 = 0. Also, the Grassmann expansion of the function exp{−V Next, we have to perform the Grassmann Gaussian integration with respect to the variable ζ ψ . The result is a function of the Grassmann variable ψ. The next lemma allows to get bounds on the coefficients of the new Grassmann polynomial, in terms of the (κ, N , M)-bounds of the integrand.
Furthermore, if f a,0 = 0 for all a, then: Proof. Set h(ψ) = a h a ψ a . Notice that: If |b| is even, the outcome of the integration is bounded by 2. Therefore: If furthermore f a,0 = 0 for all a, then the 1 in the last parenthesis is not present.

Integration of the Bosonic Fluctuation
Field. Next, we shall compute: Vol. 21 (2020)

A Supersymmetric Hierarchical Model 3517
We rewrite: that is: We then get: where: The bounds for these new coefficients are collected in the next proposition. n (φ/L, ζ φ ) satisfy the following bounds, for λ small enough and for all ζ φ ∈ C 4 , φ ∈ C 4 : Proof. The statement for n = 0 is trivial, since D 0 . Consider now n = 0. From Eq. (3.40), we get: Using the bounds for C n−k , Proposition 3.6, we get: which concludes the proof.
We are now ready to integrate the ζ φ variable. We shall perform the integration for φ ∈ C 4 : the reason being that the bounds on the kernels on the next scales will be obtained via Cauchy estimates. In order to integrate the field ζ φ , we shall exploit the oscillations of the complex Gaussian, via the next lemma.
The next lemma is the key technical tool that we will use to estimate the derivatives and the error terms arising from the stationary phase expansion.
Then, for all multi-indices α ∈ N 4 and for all z ∈ B R : Then, for all k ∈ N, there exists C k > 0 such that: In particular, let E m (f ) be the error term in the stationary phase expansion (3.44). Then, there exists a universal constant K m > 0 such that: As a test run, let us estimate the normalization factor N (0) in Eq. (3.11). Notice that, as a consequence of the localization theorem, see Theorem C.10, supersymmetry (in the sense of Definition C.4) implies that, see Remark C.12: Nevertheless, the simple procedure discussed below will be generalized to the computation of the effective potential and of the correlation functions, where the localization theorem cannot be applied because of the lack of supersymmetry. We rewrite: (3.52) let us estimate the various terms. By Proposition 3.7, we have |D To estimate the second term, we write (all derivatives correspond to the field ζ φ ): Moreover, Proposition 3.7 together with Cauchy estimates for the derivatives gives: The second estimate follows from the Cauchy estimate (3.45) and from the bound (3.41), after taking R = λ − 1 4 . Consider now the last term in Eq. (3.52). We claim that, for some Ldependent constant C L > 0: To prove this estimate, we proceed as follows.
. This function is trivially entire in ζ φ ∈ C 4 . Let us estimate it. We have: Also, Therefore, using that |μ| ≤ C|λ|, for some K > 0 and taking λ small enough: This concludes the check of Eq. (3.56). Next, we shall consider, for φ = 0: to do this, we shall discuss separately a small and a large field regime for the bosonic field φ ∈ C 4 . More precisely, we shall consider separately the following cases: φ/L ∈ S (0) , the small field set, and φ/L ∈ L (0) , the large field set:

Small Field Regime.
Let φ ∈ LS (0) . Let us define: As a preliminary remark, notice that E (0) n is analytic in φ ∈ LS (0) , since the integral is absolutely convergent in ζ φ uniformly for φ ∈ LS (0) and the integrand is entire in φ. (Analyticity follows from dominated convergence theorem and from Morera's theorem.) Let us now prove bounds for E We shall proceed as we did for the analysis of the normalization factor N (0) . Consider the first term in Eq. (3.65). By Proposition 3.7, we get, To estimate the second term in Eq. (3.65), we use a Cauchy estimate. To begin, notice that for φ ≤ λ − 1 4 L and for ζ φ ≤ λ − 1 4 , by Proposition 3.7: Therefore, by the Cauchy bound (3.45), with R = λ −1/4 and R = 0: On the other hand, recalling the definition of V Using that: we finally get, from the bounds (3.67), (3.68), (3.69), for λ small enough: Consider now the third term in Eq. (3.65), which is the remainder in the stationary phase expansion. We shall first extend the bound (3.60) to φ ∈ 1 4 , proceeding as in Eqs. (3.57)-(3.59): (3.73) In conclusion, the bounds (3.72), (3.73), together with the estimate (3.41) for D Hence, we are in the position to apply Lemma 3.9, with: We get, by Lemma 3.9: This bound together with (3.66), (3.71) implies, for λ small enough: Let us now consider E (0) n (φ) for n = 1, 2. By the stationary phase expansion: By Proposition 3.7: To estimate the second term, we proceed as follows. We write: the first term is estimated using the bound (3.69). We get: Concerning the second term, it is convenient to rewrite the derivative in terms of the C (0) n coefficients; recall the definition (3.40). We get: We estimate the right-hand side using the bounds for the C (0) n coefficients, (3.32), plus a Cauchy estimate with R = λ − 1 4 for the first term. We get: Therefore: (3.86) Therefore, by Lemma 3.9: In conclusion, for λ small enough, the expression (3.79) for E As it will be clear later on, the bounds (3.77), (3.88) are not enough to iterate the multiscale integration on higher scales. We shall isolate the dangerous contributions by introducing a localization operation, as follows.
, with a slight abuse of notation.
Localization and renormalization. We define, for φ ∈ LS (0) : (3.89) That is, the L operator extracts the first few orders in the Taylor expansion of E (0) n . To see this, notice that (φ · φ) is just the analytic continuation of φ 2 from φ ∈ R 4 to φ ∈ C 4 , recall Eq. (3.8). Also, notice that in the expansion in φ ∈ R 4 of E (0) n (φ) odd powers of φ are forbidden, due to the fact that they are not analytic in φ i,σ . Hence, Eq. (3.89) collects the first few orders in the Taylor expansion in φ, for φ ∈ LS (0) . The terms E (0) 0 (0) are relevant in the renormalization group terminology. As it will be clear later on, they correspond to an expanding direction in the RG flow. The other terms are irrelevant, thus strictly speaking that there should be no need to localize them. Nevertheless, the above procedure turns out to simplify the analysis of the large field regime.
Correspondingly, we define the renormalization operation R so that the function RE To estimate the derivatives, we will use again Cauchy bounds. More precisely, we will consider the functions E (0) n (φ) in a much smaller domain than the original analyticity domain LS (0) , for which the bounds (3.77), (3.88) hold true. By the general estimate (3.45), we will use that every derivative introduces a gain with respect to the L ∞ bound of E (0) n (φ) in LS (0) , proportional to the inverse of the distance between the smaller domain and LS (0) .

Large Field Regime.
Let φ ∈ LL (0) . With respect to the small field region, here we have to face the extra difficulty that the terms " might be large; recall Eq. (3.38). In the small field region, we could control these terms using the quartic term in ζ φ , and the smallness of φ. In the large field region, we shall exploit the sign of the real part of such terms. This is the content of the next inequality. Let 0 < ε < 1/4, and consider ζ φ such that Im ζ φ ≤ λ − 1 4 +ε . Then, for This holds as a special case of Proposition 3.15, proven in Appendix B. Therefore, we are in the position to apply Lemma 3.8. We rewrite: The functions g The bounds for the functions D (0) n in (3.41) imply the following (non-optimal) estimates, for all φ, ζ φ in C 4 , and for a universal constant C > 0 1 4 , the bound (3.99) together with the bound (3.97) implies: To estimate E (0) n (φ) efficiently, we shall perform a stationary phase expansion. We have: From the second of (3.99): instead, to estimate the remainder term from the stationary phase expansion, we use that, for W = λ − 1 4 +ε , thanks to the bound (3.100): Hence, by Lemma 3.9: Consider first n = 0. Using that φ > Lλ −1/4 we have, for any 1 2 ≤ δ < 1, for λ small enough and L large enough, from the bounds (3.102), (3.104): The first bound follows from the fact that the combination λL φ 2 − λ 8L φ 4 can be arbitrarily negative, for L large enough uniformly in φ ∈ LL (0) . The second bound follows from the observation that, for λ small enough and L large enough uniformly in φ ∈ LL (0) , the combination λL 4 (3.106) Consider now E n for n = 1, 2. By the above reasoning, for λ small enough and L large enough: Ann. Henri Poincaré Therefore, This concludes the discussion of the large field regime.
3.2.6. The Effective Potential on Scale h = 1. We obtained: where the functions E  To conclude the discussion of the scale zero, we have to renormalize the coupling constant and the chemical potential, by taking into account the terms extracted with the localization procedure in Sect. 3

.2.4.
Small field bounds. We rewrite: L . We define: We then rewrite: The function U (1) L (Φ) can be expanded in powers of (ψ ·ψ), in terms of suitable coefficients U (1) L;n (φ): (3.115) Notice that, by construction, thanks to the definitions (3.112), U L (0) = 1 and moreover the function U where U R (Φ) can also be expanded in powers of (ψ·ψ), for suitable coefficients U All in all: To conclude the integration of the scale zero, we shall estimate the coefficients R (1) n (φ). The functions R (1) n (φ) are analytic in S (1) , with: (3.120) L;n (φ). By construction, this function has a Taylor series in φ = 0 that starts from order 6 − 2n. By inspection, and using the bounds (3.113), we have: R;n (φ). By inspection: (3.122) The constant C takes into account the fact that, for φ ∈ S (1) and for λ small enough, there exists a universal constant K such that e |β4| φ 4 +|β2| φ 2 ≤ K.
Let us summarize the large field analysis. The R n (φ) functions are analytic for φ ∈ L (1) , and for those values of φ they satisfy the bounds: respectively, from darker to brighter gray. The dotted, dashed and thick circles represent, respectively, the boundary of the small field sets LS (h) , S (h+1) and S (h) 142) for e − 1 8 ≤ δ < 1. These bounds conclude the discussion of the integration of the scale zero.

General Integration Step
We are now ready to perform the integration of the general scale h ≥ 0. We shall show inductively that the effective potentials satisfy certain properties and bounds that allow to iterate the map. These properties and bounds are the content of the next theorem. Theorem 3.10 (Effective potential flow). Under the same assumptions of Theorem 2.2, the following is true. Let C > 0, 0 < ε < 1 4 and e − 1 8 ≤ δ < 1. Then, for λ small enough, there exists a unique choice μ = μ(λ) ∈ C, with |μ(λ)| ≤ Cλ, such that for any N ∈ N and for any 0 ≤ h ≤ N the effective potential U (h) (Φ) can be written as: where:

144)
and where R (h) n (φ) are analytic functions in φ ∈ S (h) ∪ L (h) ; the small field set S (h) and the large field set L (h) are defined as, see Fig. 2: , and R (h) 0 (0) = 1. They satisfy the following bounds, for a universal constant C > 0. Small field bounds. Let φ ∈ S (h) . Then: (3.146) Large field bounds. Let φ ∈ L (h) . Then: 147) with c 0 = 0 and, for a universal constant K > 0: All the statements in the theorem are trivially true on scale h = 0 and have been checked on scale h = 1 in Sect. 3.2.6. The goal of this section is to prove by induction that they propagate to scale h + 1. To be more precise, we shall prove all but the statement concerning the existence of μ(λ) and the bounds |μ h | ≤ C|λ h |: This part of the proof is postponed to Appendix A.

Setting up the Integration.
Recall the flow of the effective potentials, defined in Eq. (3.5): with normalization: 148). We introduce the notation: We have: We set: our first task will be to derive bounds for the functions B (h) a,b (φ/L, ζ φ ). To this end, the next lemma will be useful.

Proof. P roof of item(i). For simplicity, consider
ψ . The claim follows by noticing that and hence, |h a,b | = |f a+b | ≤ κN |a| N |b| .
Lemma 3.11 will be used to prove the following statements on the B for |a| + |b| > 0.
(ii) Let φ/L+ζ φ ∈ S (h) and φ/L−ζ φ ∈ L (h) . Then, for δ as in the assumptions of Theorem 3.10, see also Eq. (3.147): Proof. To begin, notice that and therefore R We also notice that for φ/L ∈ S (h) the inductive estimates ( On the other hand, for φ ∈ L (h) , the inductive estimates (3.147) imply that f (Φ/L) satisfies (κ, N ) bounds with:

Proof of item (i). By using the inductive assumption (3.146) on R
(h) 0 , we get: For the case n > 0, we proceed as follows. By Lemma 3.11 part (i) and by Eq.
We then write: We are now ready to estimate the coefficients of the function g(Φ/L, ζ) with n > 1. We have: This concludes the proof of item (i).

Proof of item (ii).
Here, we write: (3.172) From this, the proof of item (ii) easily follows.

Proof of item (iii).
We proceed as in the proof of item (ii), except that now all functions depend on large fields. Hence, we shall use the (κ, N , M)-bounds for g (h) generated by Eq. (3.162). We easily get:

173) which implies item (iii).
Recalling the expression (3.152), (3.153), we rewrite: where: we rewrite: In the next section, we shall discuss the integration of the fermionic fluctuation field ζ ψ .

Integration of the Fermionic Fluctuation.
The goal of this section is to compute: The integration is performed by expanding the exponential, and by using the fermionic Wick's rule to integrate the field ζ ψ . By Corollary C.2, the outcome of the integration depends on the field ψ only via the combination (ψ · ψ). That is: for suitable functions C

178)
Then, there exists a universal constant K > 0 such that the following bounds hold true. (3.180) Proof. The functions C  .26), with the only difference that now λ and μ are replaced, respectively, by λ h and μ h . We get: Recall the definition of the function g (h) (Φ/L, ζ), Eq. (3.160).

Consider first the term h
To deduce bounds for h (h) I (ψ; φ, ζ φ ), we use Lemma 3.5. We get that h (h) I (ψ; φ, ζ φ ) satisfies (κ , N )-bounds, with: Next, by Lemma 3.5, the outcome of the Grassmann integration satisfies (κ , N )-bounds with:  This proposition concludes the discussion of the integration of the fermionic fluctuation field. In the next section, we shall discuss the integration of the bosonic fluctuation field.

Integration of the Bosonic Fluctuation Field.
We now consider: To begin, we extract from V (h) b (Φ, ζ φ ) the contribution due to the fermionic external field ψ. We proceed as for the integration of the scale zero; see Eqs.
(3.36)-(3.40). We have: (3.192) where: . Then: . Then: Proof. The statement about analyticity follows immediately from the analyticity of the C Suppose now that n = 0. We write: (3.198) Plugging in the bound (3.179) for C (h) n−k , we get: where in the last step we chose L large enough. The final claim (3.196) follows from Proof of item (iii). The proof of item (iii) is identical to the proof of item (ii), the only difference being that one has to use the bound (3.180) for the functions C (h) n−k . We omit the details.
Proposition 3.14 implies that the functions D (h) n satisfy the following (non-optimal) bound, for all values of φ, ζ φ such that φ/L ± ζ φ ∈ S (h) ∪ L (h) : where we used that for λ small enough. This bound will be used to estimate the remainder terms of the stationary phase expansion.

Small Field Regime.
We define: (3.203) in terms of which the effective potential can be rewritten as: We proceed in a way analogous to the integration of scale zero; see Sect. 3.2.4.
To begin, we have to prove that the bounds we derived on the functions D (h) n are compatible with integrability. This is an immediate consequence of the following bound, whose proof is deferred to Appendix B.

Proposition 3.15.
Let ε > 0 as in the inductive assumptions (3.148). 1 4 . Suppose that c h is as in Eq. (3.148). Then: (3.205) Proposition 3.15 together with the bound (3.201) immediately implies the following estimate, for Im Next, recall that, by Proposition 3.14, the integrand is analytic in φ and ζ φ , provided φ/L ± ζ φ ∈ S (h) ∪ L (h) . Therefore, since in the integral ζ φ ∈ R 4 , the integrand is analytic for φ ∈ LS (h) . Moreover, thanks to Proposition 3.15, for these values of φ the integrand is absolutely integrable in ζ φ . In conclusion, by dominated convergence and by Morera's theorem, the function E Let us now derive estimates for E We write, by stationary phase expansion: By Eq. (3.194): Consider now the second term in Eq. (3.207). We have: the first term is estimated using the analogue of the bound (3.69). We have: Hence, using Proposition 3.14 part (i), together with the bound (3.210) we get: The second term in Eq. (3.209) is bounded using a Cauchy estimate, for ζ φ ≤ (1/2)|λ h | −1/4 , in order to make sure that φ/L±ζ φ ∈ S (h) . We have, proceeding as in Eqs. (3.83)-(3.85): Finally, consider the third term in Eq. (3.207). This term will be estimated using Lemma 3.9, which requires analyticity in ζ φ in R 4 W . We shall choose W = |λ h | − 1 4 +ε ; due to the restriction to the smaller set of fields φ ∈ (L/2)S (h) , the condition ζ φ ∈ R 4 W implies that φ/L±ζ φ ∈ S (h) ∪L (h) ; hence, the argument of the integral in Eq. (3.203) is analytic in φ ∈ (L/2)S (h) and in ζ φ ∈ R 4 W . Thanks to Lemma 3.9 and the bound (3.206), we get: In conclusion, for λ small enough uniformly in h, for all φ ∈ (L/2)S (h) : Notice that this bound is identical to the corresponding one obtained on scale zero, Eqs. (3.77), (3.88), except that now λ is replaced by |λ h |. Notice also that that, by supersymmetry, E Localization and renormalization. Next, we define a localization and a renormalization procedure, restricting the set of allowed values of φ. Given λ h+1 ∈ C (to be determined later) such that |Lλ h+1 − λ h | ≤ C|λ h | 3 2 , we define: Notice that, for some universal constant c > 0: (3.217) Arrived at this point, the discussion of the small field regime is identical to the one for scale zero, see Eqs. (3.93)-(3.96), except that in all estimates one has to replace λ with |λ h |. We get, by the bound (3.214), using Cauchy estimates for φ ∈ S (h+1) : (3.218) Concerning the contribution of the L-term, we set: By supersymmetry, see Corollary C.9, Appendix C: (3.220) Also, by Cauchy estimates: This concludes the small field analysis.

Large Field Regime.
To begin, notice that the function E . This follows from the analyticity in φ, ζ φ of the argument of the integral in Eq. (3.203), which holds provided φ/L ± ζ φ ∈ S (h) ∪ L (h) , and from the bound (3.206), which ensures integrability in ζ φ . As for the small field regime, analyticity in φ ∈ LL (h) follows from dominated convergence and from Morera's theorem. We shall now prove bounds for E (h) n (φ), in the domain: with λ h+1 as in the small field regime (to be determined later). Notice that φ ∈ L (h+1) implies Im φ ≤ (1/2)|λ h | − 1 4 provided that L is large enough.
We compute E h n (φ) performing one step of stationary phase expansion. We have: We shall discuss the relevant bounds separately in the regions Case φ ∈ L (h+1) ∩ LS (h) . We estimate the first term in Eq. (3.223) as in (3.209), since the bound holds for any φ ∈ LS (h) : Consider now the remainder term of the stationary phase expansion in Eq. (3.223). By the bound (3.206), together with Lemma 3.9, for W = |λ h | − 1 4 +ε : Therefore, proceeding as we did for the bound (3.139), for λ small enough: withδ < δ, sayδ = e − 1 6 . Similarly, for n = 1, 2: We estimate the first term in Eq. (3.223) by using Proposition 3.14 item (iii ). We have, for L large enough:  4 for a universal constant K, for λ small enough. Concerning the remainder term, it is bounded as before.
As we did after the integration of the scale zero, we now redefine the effective coupling constant and the chemical potential, taking into account the terms extracted in the renormalization procedure.
Small field bounds. Here, we proceed exactly as for the corresponding discussion on scale zero, Eqs. (3.110)-(3.124). The only difference is that now λ is replaced by λ h . We get, for φ ∈ S (h+1) : with: By supersymmetry, it will be enough to study the bosonic two-point function; see Eq. (4.7). Following Eq. (2.9), we rewrite the fields φ + x , φ − y as: Let L −1 x = L −1 y . For the sake of notation, in this section we shall drop the spin label, unless otherwise stated. We compute the two-point correlation function with equal spins; by spin symmetry, the two-point correlation function with different spins is trivially zero. Plugging the decomposition (4.1) in φ + x φ − y , and using that: 2) we get: This procedure can be iterated. Let k ∈ N be the first integer such that We are left with discussing the evaluation of φ where the functions U (k−1) (·) are the outcome of the construction of Sect. 3 (we again used that N (h) = 1 for all h, by SUSY). Again by SUSY, see Remark C.12: It therefore suffices to compute the bosonic two-point function, since Eq. (4.6) together with the previous discussion implies

Integration of the Nontrivial Scales
We write: plugging this identity in Eq. (4.5), and using again that the average of odd functions is zero:

Integration of the Scale k−1.
In this section, we discuss the integration of the first nontrivial scale, corresponding to the entry h = k − 1 in the sum in Eq. (4.9). To evaluate ζ , we proceed as follows. By definition of effective potential: which we rewrite as: where we introduced: the last step follows from the integration of the fermionic fields, discussed in Sects. 3.3.1, 3.3.2, and from the definition of the D (k−1) functions; recall Eq. (3.192). We now discuss the integration of the bosonic fluctuation field, which we shall perform via a stationary phase expansion. We have: where we defined: Plugging the expansion (4.12) into (4.12), we get: where we used that, by the definition of the D (k−1) functions: (4.15) and we defined: ·)) . (4.16) Let us now estimate the right-hand side of Eq. (4.16). We use the bound (3.206), together with Lemma 3.9. We get, recalling the definition (4.13), for Im ζ φ ≤ |λ k−1 | − 1 4 +ε , Im (φ/L ± ζ φ ) ≤ |λ k−1 | − 1 4 : (4.17) Vol. 21 (2020) The extra factor |λ k−1 | − 1 2 is due to the presence of ζ Hence, by Lemma 3.9:
(4.20) We shall compute the integral in a multiscale fashion. We have, for all such that k + ≤ N : where, denoting by x h (y) := L −h+k y ∈ Λ (h) for h > k − 1: x h+1 (y) ) := dμ(ζ) In particular, Eqs. (4.20), (4.21) imply: That is, Let us rewrite the map (4.22) in a more symmetric way. Recalling that A z = ± and that x∈B (h+1) x h+1 (y) /L + ζ) . (4.25) To prove this equality, we assumed that A x h (y) = 1. If not, we can reduce the discussion to this case by performing a ζ → −ζ change of variable. Let us now assume inductively that: , and such that: |G for some 0 < C h ≤ 2 K L . These assumptions are true for h = k; see Eq. (4.19) (there, C k = K L ). Our goal will be to show that these bounds propagate to scale h + 1.
To complete the evaluation of F (h+1) k−1 (Φ), we are now left with the integration of the bosonic fluctuation field in Eq. (4.28). By stationary phase expansion, Lemma 3.8: Consider the first term in Eq. (4.47). We have, thanks to Proposition 4.1, for φ/L ∈ S (h) ∪ L (h) : (4.48) where in the last step we used that: and we took into account the first exponential by (1 + 24|λ h | 1 16 )e 4CL|λ h | 1 4 ≤ (1 + 48|λ h | 1 16 ). Consider now the remainder term in the stationary phase expansion. In order to bound the error, we need the analogue of (3.206) for the argument of E 1 . We shall rely on Proposition 3.15, together with the estimate: where we used that, for λ small enough, ; recall the definition after Eq. (3.206). The factor −|λ h | 4ε will be useful later on. Hence, by Lemma 3.9: (4.52) All in all, putting together Eqs. (4.47), (4.48), (4.52), we get: (4.54) Notice that we used the factor L −2n in Eq. (4.48) to update the prefactor |λ h | n 8 to |λ h+1 | n 8 . Let us now update the running coupling constant appearing in the explicit exponential prefactor in Eq. (4.53). We write: where in the last step we expanded the overall exponential as a polynomial in (ψ · ψ), and collected terms of the same powers. We shall now prove bounds for the new functions G (4.58)

Conclusion.
We are now ready to compute the two-point correlation function; recall Eq. (4.24). We have: (4.59) where the last bound follows from F k−1;0 (0) together with (recall (4.27)): |G the various terms admit the following bounds: In order to improve the bounds we already obtained on the beta function, we furthermore notice that provided that |μ j | ≤ 2C|λ j |, for all j ≤ h, the following expression is attained: where the first inequality follows from Cauchy estimates in the ball of radius R = L|λ h | − 1 4 , andK L ≡K L (C). Therefore, for λ is small enough: for a universal constant K > 0, which does not depend on C. In particular, by taking C large enough, we have: Existence. We shall show that there exists a solution {μ h } ∞ h=0 of Eq. (A.1) with the desired properties. Later, we will comment on uniqueness. Our discussion closely follows [43]; see also [12]. Let μ ≡ μ 0 ∈ {z ∈ C | |z| < 2C|λ 0 |} =: I 0 . By construction, β (0) 2 is a continuous function of μ 0 ∈ I 0 , and |β (0) 2 | ≤ KL|λ 0 |. Thus, Eq. (A.1) implies that I 0 μ 0 → μ 1 (μ 0 ) is continuous and that: where the last step follows from Eq. (A.6). Thus, by continuity there exists I 1 ⊂ I 0 such that: This shows in particular that, for μ 0 ∈ I 1 , |μ 1 | ≤ 2C|λ 1 |. Now, suppose inductively that there exists I h ⊂ I 0 such that I h μ 0 → μ k (μ 0 ) is continuous for all k ≤ h and that: In particular, Eq. (A.9) implies that |μ k | ≤ 2C|λ k | for all k ≤ h. These assumptions are true for h = 1, as we just proved. Let us check the inductive assumptions on scale h + 1. By the RG construction |β is continuous in μ k , k ≤ h, and hence in μ 0 ∈ I h . Equation (A.1) implies μ h+1 (μ 0 ) is continuous in μ 0 ∈ I h , and that: where the last step follows from Eq. (A.6). Thus, by continuity there exists I h+1 ⊂ I h such that: This shows in particular that |μ k | ≤ 2C|λ k | for all k ≤ h + 1, which is what we wanted to prove.
Uniqueness. Here, we shall prove the uniqueness of μ(λ) : We shall show that the set I h shrinks to a point as h → ∞. To do this, we rely on the Lipschitz continuity of the beta function of the effective chemical potential, as a function of the effective chemical potentials on all the previous scales. This will be proven via a Cauchy estimate, which in turn relies on the analyticity properties of the effective potential U (h) as function on {μ j } h j=0 . In fact, in the next proposition, we shall regard the effective potentials U (h) as functions of the sequence {μ k } h k=0 : As it is clear from the induction of Sect. 3.3, the only information about the sequence of effective potentials that we required is that |μ k | ≤ C|λ k |. Also, as already pointed out at the beginning of the section, the constant C can be taken arbitrarily large, provided λ is small enough.
Proof. The proof is by induction. The case h = 0 is true by inspection.
Next, we assume that Theorem 3.10 and Lemma A2 are true on scales j < h, and we shall prove that they hold on scale h. Consider: is an analytic function in {μ j } h j=0 provided that |μ j | < C|λ j | for all j.
Proof. In Sect. 3.3, we proved that U (h) 0 (φ) is analytic in φ ∈ S (h) and that U (h+1) 0 (φ) depend on φ ∈ R 4 only via φ (see Appendix C). Therefore, the last identity follows from R (h+1) 0 (0) = 1, implied by SUSY (see Appendix C), and by the fact that R (h+1) 0 (φ) − 1 is at least of order φ 6 . By Lemma A2, we know that U (0); this together with the identity β We are now ready to prove uniqueness of μ(λ). We shall proceed by contradiction. Suppose that the function μ(λ) is not unique: There exist two sequences μ = {μ k } ∞ k=0 and μ = {μ k } ∞ k=0 such that |μ k | ≤ 2C|λ k | and |μ k | ≤ 2C|λ k |, solving Eq. (A.1), with C as in the statement of Proposition A1. We are denoting by {λ k } the effective quartic couplings associated with the sequence {μ k }. By the assumptions on the sequence {μ k }, it is true that L k for all k. We have: 2 (μ k , . . . , μ 0 ) ; (A.14) by Corollary A3, the beta functions are analytic in {μ j }, {μ j }, respectively, if |μ j | < C|λ j | and |μ j | < C|λ j |. Notice that we are free to assume that C > 4C, if λ is small enough. We rewrite Eq. (A.14) as: where we used that the sequences vanish at infinity. In order to estimate the difference of the beta functions, we proceed as follows. We use the telescopic representation: and we use that, in the integral, the distance between ν and the boundary of the analyticity domain of β (j) 2 in its rth variable is bounded below by CλL −r . Also, for these values of ν, the beta function satisfies the bound |β . . . , μ r+1 , ν, . . . , μ 0 )| ≤ KλL −j . Therefore, by a Cauchy estimate: with μ ∞ = sup k |μ k | and where K is proportional to K/C. Plugging this bound into Eq. (A.15), we get: For L −1 K small enough, Eq. (A.20) implies μ = μ . This concludes the proof of Proposition A1.

B Proof of Technical Lemmas
In this section, we collect the proofs of three key results, namely Lemmas 3.8, 3.9 and Proposition 3.15.

B.1 Proof of Lemma 3.8
Let S(R n ) be the Schwartz space, and let f ∈ S(R n ). Then: The proof of (B.1) goes as follows. Since f ∈ L 1 (R n ) we have, by dominated convergence: Then, since e −(i+ε) x 2 ∈ L 2 (R n ), f ∈ L 2 (R n ) we have, by Plancherel's theorem: where we used the assumption (3.47). Thus, from Eq. (B.9) we easily get that, for any p ∈ R n , m ∈ N, and for some C m > 0:
Supersymmetry of regular enough functions implies remarkable identities after integration over the superfields. Here, we shall only consider Schwartz functions of the superfield Φ, defined as follows.
Definition C.6. We say that a function f (ζ) The next lemma will be useful later on, to perform integration by parts over the superfields.
This lemma can be used to prove that the effective potential U (h) (Φ), restricted to φ i,σ ∈ R, is a supersymmetric function.
Proof. The proof is by induction. It is trivially true for h = 0. Suppose it is true for k < h, and let us prove it for k = h. We rewrite U (h) as: Let φ i,σ ∈ R. Being the explicit factor e − λ h L (Φ·Φ) 2 −iLμ h (Φ·Φ) supersymmetric, Eq. (C.11) implies that: Explicitly, Finally, we conclude the appendix by mentioning a well-known result on supersymmetric functions [18,70] (see also [12], Theorem 11.4.5).
Theorem C.10 (Localization theorem). Let the function f (ζ) be supersymmetric in the sense of Eq. (C.8) and of Schwartz type in the sense of Definition C. 6. Then: (C.24) Remark C.11. We are not formulating the localization theorem in its most general form. Furthermore, in the present context the theorem is a simple and elegant application of integration by parts [18].