Functional renormalization group for multilinear disordered Langevin dynamics I: Formalism and first numerical investigations at equilibrium

This paper aims at using the functional renormalization group formalism to study the equilibrium states of a stochastic process described by a quench--disordered multilinear Langevin equation. Such an equation characterizes the evolution of a time-dependent $N$-vector $q(t)=\{q_1(t),\cdots q_N(t)\}$ and is traditionally encountered in the dynamical description of glassy systems at and out of equilibrium through the so-called Glauber model. From the connection between Langevin dynamics and quantum mechanics in imaginary time, we are able to coarse-grain the path integral of the problem in the Fourier modes, and to construct a renormalization group flow for effective Euclidean action. In the large $N$-limit we are able to solve the flow equations for both matrix and tensor disorder. The numerical solutions of the resulting exact flow equations are then investigated using standard local potential approximation, taking into account the quench disorder. In the case where the interaction is taken to be matricial, for finite $N$ the flow equations are also solved. However, the case of finite $N$ and taking into account the non-equilibrum process will be considered in a companion investigation.


Introduction
Multilinear Langevin dynamics are characterized by a specific notion of disorder, said to be quenched, and referring to the existence of two fluctuations times: the time τ of some random variables {q i } and the time T of the couplings between these random variables. For a glassy system, we have T τ , meaning that the couplings are essentially constant during the typical fluctuation time of the random variables {q i }. Initially introduced as a mathematical description of the magnetic field, it allows to obtain strange thermodynamics properties through for instance the so-called Glauber model [24,29,38,46,61]. Glassy system have acquired an independent existence, essentially due to their mathematical complexity, and are encountered in a large area of research fields, among which we find the condensed matter physics, chemistry, genetics, computer science (as a description of neural networks behavior), collective processes in markets, economy, tensorial principal component analysis, etc. (see [2,13,14,28,41,42,79] for more details).
The existence of a freezing temperature T f , below which the system may be trapped into an equilibrium non-ergodic region of the phase space, and where the system is enforced to fluctuate around some states but cannot visit every other possible equilibrium states, is the fundamental behavior of the disordered system. The physical meaning of this phenomenon leads to the existence of many deep minima of the energy landscape, with large energy barriers between the different minima, such that the system requires a very long time to escape from them. The existence of such an ergodicity breaking is one of the standard signatures of the glassy transition, occurring generally in dynamical investigations based on a stochastic, Langevin equation. Other standard approaches for glass transition use replica method, the cavity method through a generalized mean-field approach, or the so-called "TAP" (Thouless-Anderson-Palmer [85]) approach focusing on complexity (i.e. configuration entropy, counting the number of energy minima) [24,26,27,38,67,61]. All of them provide the characterization of the transition with different signatures, replica symmetry breaking for the first one and the divergence of the magnetic susceptibility for the second one and then the divergence of the complexity for the third one, which highlights their complementarity.
The Langevin equation formulation has the advantage to be able to describe systems from a dynamical point of view both near and far from equilibrium, therefore we particularly focus on this formulation in the following work because it corresponds well to the system that we want to study. This kind of equation has been widely studied recently, and we can in particular mention some relevant results [5,17,30,54,60,80]. A large part of these investigations combine numerical and analytic methods; but analytic solutions were still found only in some limited cases, in particular the references [54,60] studied the large N limit where self-averaging property holds. Another specific case is matrices disorder where an exact analytic solution can be found due to the well-known large-N properties of the random matrices, whose spectrum converges toward the Wigner semi-circle law [9,52,69,74,91]. However, such a simplification does not occur for the multilinear case, which is the main topic of this paper, essentially because no exact analogue of the Wigner law exists for tensors 4 . In this paper, we propose an alternative approach using coarse-graning in frequencies modes to construct effective states of the system for large times. More precisely, the degrees of freedom can be labelled with frequency modes ω, and we are aiming to use renormalization group (RG) and nonperturbative techniques to construct an RG flow for effective quantities, integrating out large frequencies modes (i.e. rapid modes effects) 1/ω T , for some running time-scale T .
In particular, we deal with the large-time behavior of the system, describing equilibrium states from initial, out-of-equilibrium conditions. To this end, we consider a coarsegraining approach in time, through an RG description allowing us to investigate physics at different temporal scales integrating out large frequencies modes. Let us recall that RG realizes an interpolation between two physical descriptions of the same system, through a progressive integration of degrees of freedom; which provides an effective description valid at a large scale from a microscopic description [33,73,77,78,98,99,100]. This coarse-graining procedure dilutes information from each step, and the effective description discards some irrelevant effects which decrease along with the RG flow. Hence, RG is usually considered as a powerful tool to discuss universality in quantum and classical field theories; and generally in systems involving intrinsic variability. There is a recent literature on the link between RG and information theory, see for instance [16,71,37]. In [71] for instance, authors propose a point of view where RG looks like a maximum entropy dynamics where coarse-graining is generated by a decrease of the relative entropy (with respect to a fixed background probability). For recent connections with data analysis see [15,55,56,59,57,58].
In the present manuscript we use the functional renormalization group (FRG) formalism due to Wetterich & al. [10,31,32,33,36,64,65]. The main interest of this method is that it avoids some difficulties occurring in other mathematical incarnation of the RG idea, especially concerning the strong coupling regimes; where instabilities may occur with methods based on other functional relation as Polchinski equation [36,70,81]. Thus, its ability to support crude truncation is especially relevant in a case where the coupling's magnitude fluctuates and can become arbitrarily large. Coarse-graining in time has been considered through the FRG formalism for quantum mechanical systems [47,81,96]; especially regarding supersymmetric aspects of quantum systems. Exploiting the link between Langevin equation and Schrödinger equation in imaginary time, the same techniques have been used for stochastic systems [35,75,92,93]; focusing on the effective Euclidean action. In the quantum description, stochastic systems exhibit elementary supersymmetry which simplifies the choice of the truncation; explicit supersymmetric truncation strongly improves the approximate RG flow for non-supersymmetric ones [18,19,75,93].
RG has been considered as a powerful tool of investigation for such a kind of problem, especially in regard to the behavior of spin glass [21,22,23,63,72] or the so-called randomfield Ising problem [6,7,26,48,82,83,84,86,87,88,89]. In this paper we propose to investigate dynamics of a zero-dimensional p-spin like system, through FRG formalism, mathematically incarnated as a Langevin equation with quenched disorder, describing dynamics of a N dimensional random vector q i (t). A large part of this paper is devoted to the multilinear case, with tensorial disorder. Our aim in this paper is essential to establish the future the formalism at equilibrium, meaning that we are only able to discuss vitreous transitions above T c . Beyond physical applications, such a formalism can be helpful for other areas of sciences, like data analysis, where such a kind of equations appears, and we plan to investigate these aspects in a future work.
The outline is the following. In section 2 we define the models and review some aspects about the standard connection between Langevin equation, Fokker-Planck, imaginary time Schrödinger equations and supersymmetry. In section 3 we introduce the FRG formalism able to deal with time-reversal and quenched disorder. In section 4, we use standard local potential approximation to construct solutions of the exact RG equation compatible with the supersymmetry. We introduce a graphical formalism allowing us to draw flow equations in a very convenient way, especially relevant for multilinear models where the number of terms in the truncations become large. Finally, we investigate a similar analysis for the same model with the spherical constraint. In section 5 we provide a numerical analysis of the flow equations, and compare them with numerical simulations for the Langevin equation, especially in the large N limit. In section 6 we consider the special case of matricial disorders. Finally, section 7 is devoted to the conclusions on which we provide some remarks and ingredients which can help to investigate the non-symmetric phase and the non-equilibrium process which will be addressed taking taken into account in forthcoming works.

Preliminaries
The theory of Brownian motion is perhaps the simplest approximate way to treat the dynamics of nonequilibrium systems. The fundamental equation that governs this dynamic is called the Langevin equation. This equation is a stochastic differential equation describing the time evolution of a subset of the degrees of freedom. These degrees of freedom typically are collective (macroscopic) variables changing only slowly in comparison to the other (microscopic) variables of the system. In this section, we provide some definitions and basics for the disordered Langevin equation, required for the rest of this paper. For more detail concerning this notion, the reader can consult the standard references as [98,100].

The model
We introduce the functional renormalization group methods to investigate the behavior of a disordered system described by a random dynamical vector Q(t) = {q 1 (t), · · · q N (t)} ∈ V N depending on time t ∈ R and where V N designates the configuration space. We assume that the dynamics obeys to the non-linear Langevin equation: where η(t) = (η 1 (t), · · · , η N (t)) ∈ R N is a random vector assumed to be Gaussian distributed, with zero mean value and the variance given by: Two specific cases are relevant: 1. The unconstrained model (or soft p-spin model ), for which V N ≡ R N , and the stochastic Hamiltonian H 0 (J, Q) is defined such that: This Hamiltonian content two parts. The first one is a stochastic expression linear on the (symmetric) tensor J. The second one is the deterministic part with M ∈ N and Q 2 : , the N dimensional sphere of radius √ N , for which the stochastic Hamiltonian must be modified from a Lagrange parameter µ to ensure the spherical constraint Q 2 (t) = N : with the initial condition which satisfies the spherical constraint.
Note that we assumed the Einstein convention for summation of repeated indices. The time independent symmetric rank-p tensor J ∈ R is assumed to be randomly distributed, following to a Gaussian distribution p(J) such that: where λ ∈ R + , and: the sum running over permutations π of the first p positives integers. Note that since we are aiming to take the average with respect to two sources of fluctuation, we have to fix the notation to designate the mean quantities. Thus, we use the notation x for averaging for the noise η, and we use the notation x for averaging for J. The power of N is chosen to ensure that the extensive quantities like H 0 are proportional to N . √ λ define the size of the fluctuating tensorial coupling J which interpolates between two regimes: • For λ = 0, the fluctuations are frozen and J goes toward a fixed value J 0 . When J 0 becomes large, the influence of the noise fades out, and the dynamics reach a deterministic regime. For J 0 → 0 on the other hand, the noise dominates, and the dynamics is essentially Brownian.
• For λ → +∞, the fluctuations of the coupling dominate, and the system reach a glassy regime. We note that the nature of the dynamics could be characterized using the Hurst Exponent [44]. Indeed, this statistical measure allows to investigate the variability of a time series and thus identify a random walk from deterministic dynamics.
The spherical model is the one generally considered in the spin-glass literature. The first model however is as well considered as a simplification of the constrained case, which it reduces in some limits [1,45,50,76,80]. For instance, with J 0 = 0, and for arbitrary size N , the spherical model can be recovered as a limit case, with quartic deterministic part such that h 1 → −∞, h 2 → +∞ but h 1 /h 2 = const.. In the large N limit, the selfaveraging of O(N ) invariant quantities as Q 2 := i q 2 i (stemming as for the central limit theorem from the assumption that random variables q i are weakly correlated) ensures that quantity as Q 2 Q 2 satisfies the cluster properties [66,97,98]: Thus, by projecting equation (1) along Q, we get for the quartic potential (h m = 0 ∀m > 2): In addition to Q 2 = 0, the potential admits another minimum when h 1 < 0, for the value: The value of the mean potential for this value being V (Q) = −h 2 1 N/4h 2 . Thus, in the large N limit, the stable minimum of the potential becomes infinitely deep, and all the trajectories are trapped in effective spherical dynamics. We also note that the nonlinearities present in the local potential could induce interesting non-gaussian effects such as described in [20,3,49].

Associated Euclidean path integral and supersymmetry
The stochastic Langevin equation can be suitably reformulated as an Euclidean quantum mechanical issue involving a real wave function that obeys an imaginary time Schrödinger equation. This representation admits a path integral formulation, expressing the transition probability between two configurations as a sum over histories bounding them, weighted by the exponential of the classical action e −S [Q] . This representation moreover admits interesting supersymmetry, which is of great interest to construct an approximation of the exact RG flow in section 4.
The randomness of the trajectory Q(t) induced by the Gaussian noise η(t) can be equivalently described in terms of a functional probability distribution P (Q f , t) giving the probability to reach the final point Q f in the lapse t from the initial condition Q(t = 0) := c. Fixing disorder the probability distribution, suitably normalized, can be read as: [dη(t)] denoting the standard path integral measure. This probability distribution satisfies a Fokker-Planck type equation. Moreover, defining Ψ(Q, t) =: e H 0 [Q]/2D P (Q, t), we find that Ψ(Q, t) satisfies the Euclidean Schrödinger equationΨ = −ĤΨ, with positive definite Hamiltonian:Ĥ The ground state Ψ 0 such thatĤΨ 0 = 0 is formally solved by Ψ 0 ∝ e −H 0 (q)/2D (if it is normalizable). This has an important consequence in regard to the large time behavior of the distribution. Indeed, we expect that the transition operatorÛ (t) := e −tĤ projects into the ground state for t → ∞, and thus P (Q, t) reach its equilibrium solution lim t→∞ P (q, t) → e −H 0 (q)/D . The quantum transition probability can be formally deduced from the standard Feynman prescription. Alternatively it can be defined from the integral (10) from the formal identity 5 : where ∇ Q := {∂/∂q i } and M is the matrix with entries: Inserting the identity (12) into the integral (10), we get: The integration over the noise η(t) can be straightforwardly performed. Moreover, the determinant can be computed from the identity det M = e Tr ln M . Taking into account the normalization of the Gaussian integral over η and expanding in power of the trace, we get: .
The standard Stratonovich prescription imposes θ(0) = 1/2 for the value of the Heaviside function at origin. We thus obtain: where we used the identity, valid only for the Stratonovich prescription, see [67] for more detail: Remark that the integral in the expression (16) must be identified with the path integral associated to the quantum states Ψ. The corresponding classical action S is defined as: Also let us specify that this dynamic action associated to a stochastic process admits a natural BRS (Becchi-Rouet-Stora) symmetry [8] that can be embedded in a quantum mechanical supersymmetric formulation. This supersymmetry is known to have important consequences on the RG flow, especially relevant in the construction of nonperturbative approximate solutions of the exact RG flow [18,19]. A simple way to reveal this supersymmetry is to rewrite the determinant det M in (10) as a Grassmann integral: Now by imposing the following redefinition of the fields as: and introducing the auxiliary fieldsφ i ∈ R in order to break the square 1 2 (∇ φ W 0 ) 2 as ϕ 2 /2 + iφ i ∂W 0 /∂φ i , the classical action (18) takes the form: The supersymmetry can be made even more explicit by introducing the superfield: where θ andθ are Grassmann valued, as well as the operators: which satisfies the anti-commutator relations {D, D} = {D,D} = 0, and {D,D} = −2i∂ t . Then the following relation holds: Moreover, expanding H 0 (Φ) to the power of θ andθ, and using the properties of Grassmann algebra, it is not hard to see that the superpath integral: is reduced to the off-shell classical action (18), ifK := −[D,D]/2. Therefore, the supersymmetry properties may be explicitly checked by introducing the nilpotent operators Q := i∂θ + θ∂ t andQ := i∂ θ +θ∂ t such that {Q,Q} = 2i∂ t , and then the infinitesimal supersymmetric transformations acts on the superfield Φ through the operator δ =¯ Q − Q as Φ → Φ + δ Φ. A few algebraic computation show that the action (25) is invariant under such a transformation at first order in . More details on supersymmetry can be found in standard references [53,81,98].

Generating functional
A generating functional of the correlation functions can be defined through a standard method known as Martin-Siggia-Rose-Janssen-de Dominicis formalism which is largely discussed in the literature (see [26,43], and references therein). From a solution of the Langevin equation by fixing J, φ i (t), we may define formally the generating functional: Because of the normalization of the Gaussian noise, the averaging over η conserves the normalization and Z[j = 0] = 1. This condition allows us to perform an averaging over disorder without introducing the replica method [26]. The computation of the averaging follows the same strategy as for the previous section. The end conditions may be enforced with a suitable delta function, and the trick (12) can be used in the same way. The determinant can be rewritten using Grassmann fields, and after integrating out the delta function δ Q + ∇ Q H 0 − η , we may break the square (Q + ∇ Q H 0 ) 2 introducing an auxiliary vectorial fieldφ i . We obtain: where we introduced the short notation j(t)φ(t) := N i=1 j i (t)φ i (t) and: with the normalization conditionZ[j = 0,j = 0] = 1. As explained above, this property allows us to perform an integral over quenched disorder J without introducing replica: where Ξ := (φ, ψ,ψ,φ) denotes collectively all the fields involved in the path integral.
3 Functional renormalization group formalism

Basics
The RG aims to interpolate between a microscopic description and a macroscopic description, taking into account fluctuations at all scales. RG provides a progressive description of effective physics at different scales, integrating out firstly the modes having a small wavelength and ending by the ones having a large wavelength. This coarse-graining is the most operational incarnation of the block-spin idea of Kadanoff. For stochastic and non-equilibrium systems, temporal fluctuations cannot be ignored and can be included in the coarse-graining through a progressive integration in the frequencies space. This has been considered in some recent works, among them [35]. In this paper, we follow the same strategy, and consider a coarse-graining in time that interpolate between two regimes: 1. The UV regime, where fluctuations are frozen and fields configurations are determined by stationary points of the classical actionS.
2. The IR regime, where fluctuations are all integrated out and field configurations described through the effective action Γ, Legendre transform of the Gibbs free energy.
The standard procedure consists in adding to the classical actionS a suitable non-local functional (for some energy scale k), assuming R k , the regulator, to be differentiable in k and t. It plays the role of a momentum dependent mass, such that high energy modes concerning k are integrated out whereas low energy modes are frozen. In such a way we are aiming to construct a smooth interpolation Γ k ofS and Γ i.e. between some UV scale k = Λ (for some Λ) where all fluctuations are frozen and the IR scale k = 0 where all fluctuations are integrated out. The last condition generally imposes that R k vanish for k = 0, and the symmetry must be formally restored in the deep IR for the exact flow equations. However, this equation in itself cannot be solved exactly except for very special cases; and the approximations used to solve which work generally into a reduced dimensional phase space may introduce a spurious dependence on the regulator [70]. To avoid these difficulties, we have to be able to construct an RG that preserve as many symmetries as possiblly allowed from equilibrium condition. Among them, we demand that the RG flow will be such that it preserves causality, time reversal symmetry and supersymmetry. Note that time-reversal symmetry and supersymmetry are not truly independents, in the sense that all the Ward-Takahashi identities coming from supersymmetry can be used to derive equilibrium relations as a fluctuation-dissipation theorem, which are the consequence of the time-reversal symmetry (as the Onsager relation). However, as pointed out in [4] the inverse is not true and supersymmetry fails to provide relations where time-reversal appears explicitly. The origin of this non-reciprocity can be traced back to the fact that time-reversal symmetry cannot be decomposed as a sum of supersymmetry generators. From this observation, one expects that time-reversal symmetry implies supersymmetry, but that the reverse is false. Thus, we essentially focus on the time-reversal, and we will check at the end that supersymmetry and the corresponding Ward-Takahashi identities hold. Finally, the last requirement concerns the averaging over the disorder. Indeed, if we aim to coarse-grain before averaging over the disorder, we show from (29) and (30) that such a program requiresZ[j = 0,j = 0] = 1 before averaging, and the introduction of the regulator must not change this normalization.

Physical constraints on the regulator at equilibrium
1. Coarse-grained partition function. We start our analysis by constructing a coarsegrained partition function, labeled with a scale index k, in a way compatible with the normalization condition. Following the analyze in [35] we add to the force ∂H 0 /∂q i on the right-hand side of equation (1) a non-local driving force f i [q(t)] such that: The additional driving force is constructed as a non-local functional of the random trajectory q(t), depending on the history of the trajectory: In addition, we modify the noise correlation in a non-local manner, adding a non-local function to the Dirac delta. In such a way, the inverse propagator for η becomes: which can be understood as the introduction of a short memory in the dynamics of the system. Note that the averaging over disorder also produces such an explicit memory effect which couples different times. The non-local functions R (1) k and R (2) k aims to cut-off small frequencies modes ω k, for some infrared cut-off k. In such a way, the formal derivation of the previous section leads to the one family parameter of models described by the generating functional: the additional term ∆S k [Ξ] being defined as: Finally, we can perform the integral over J in equation (34). The integration is easier using the supersymmetric formalism, the integral that we need to compute reads as: is defined from the one to one mapping (22). From the properties (5) defining the Gaussian distribution p(J), we get the effective field theory non-local in time [26,34]:Z with: the super-coordinate z = (t, θ,θ) having measure dz := dtdθdθ and the effective (non-local) stochastic potential W 0 being given by: where we introduced the notations: Taking derivative with respect to k of the partition function (37), we deduce an equation describing the evolution of the free energy W k := lnZ k , which can be translated as an equation describing the evolution of the effective averaged action Γ k , which interpolates between the microscopic action (38) at k = Λ (for some cut-off Λ in high frequency modes) and the full effective action Γ for k = 0. It is defined as: where we introduced the inner product between superfields: including external sources (χ i ,χ i ) for fermions; and the notation M α ≡ ∂W k /∂Ξ α denoting collectively the classical fields -the extra index α running over all fermionic and bosonic fields involved into the set Ξ. The resulting flow equation for Γ k , describes the move through the full theory space from UV scales (k ≈ Λ) to IR scales (k ≈ 0) is [31,32]: the traces Tr B and Tr F denoting respectively traces over all relevant bosonic and fermionic indices; and the matrix R k is defined as: and: Note that in the rest of this paper we will denote as Γ (n) k the nth functional derivative with respect to the classical superfield.
2. Time reversal symmetry and causality. We now move on to the main physical constraints, namely the time reversal symmetry and the causality, and following [35] we will investigate the conditions on R (1) k and R (2) k which preserve these physical requirements. Indeed, from the hypothesis that the system relaxes toward equilibrium for sufficiently large time, one infers time reversal invariance of the move equations, and thus of the classical action (28). The fermionic part of the action has to be invariant from time reversal by construction. Indeed, because all the prescriptions to express the determinant must be physically true, it is sufficient to check invariance from Stratonovich prescription (15). The remaining part of the action is invariant in the transformation, up to a total derivative. These transformations acting on the bosonic sector of ∆S k [Ξ] leads to: or, from the symmetries of R Thus, integrating by part, the requirement that the global action is invariant from time reversal thus writes as, up to a total derivative: This constraint must be satisfied by all fieldsφ i and φ i . Moreover, an integration by part of the term dtdt iφ i (t)R (1) k and R (2) k is the following: Equation (45) is the first physical requirement that we impose on the regulators. The second requirement comes from causality. Indeed, because the additional driving force where θ(t) denotes the standard Heaviside function.
3. Supersymmetry and Ward-Takahashi identities. The original model (30) exhibits a non-trivial supersymmetry which can be easily checked from the discussion of the previous sections. Indeed, the kinetic part of the classical action (38) can be rewritten explicitly in a supersymmetric form [98]. Introducing the differential operatorsD and D as:D satisfying D 2 =D 2 = 0 and {D,D} = −2i∂ t , we must have: Thus, we introduce the nilpotent operators Q andQ, which anticommute with D andD and such that {Q,Q} = 2i∂/∂t, it is a simple exercise to prove that the classical action (28) is invariant at first order in¯ under the field transformation δΦ i (t, θ,θ) =¯ Q Φ i ; the variation of the kinetic action S kin [Ξ] taking the form of a total derivative. At the quantum level, this invariance translates as exact (i.e. valid to all orders in the perturbative expansion) relations between correlations functions, summarized through the Ward-Takahashi identities. Let us consider the following partition function, including sources for fermionic fields: From an infinitesimal supersymmetric transformation, the partition function becomes: where δ * is defined from linearity of δ as the image of the transformation δ on the effective field M(t). Assuming the integration measure δΞ)] (without anomaly), we must have δZ k [J ] = 0. Therefore, taking into account the explicit expression of the supersource: we get in terms of the first derivative of the classical action the compact relation: Which is the standard expression of Ward-Takahashi identity, from which all the relations between correlation functions coming from supersymmetry can be deduced. One expects that such a relations have to hold along all physical relevant RG candidates. However, the presence of the regulator must break them, if it is not explicitly supersymmetric, adding to the variation (52) the contribution δ∆S k . It is however easy to check that supersymemtry is ensured by time reversal symmetry, and such that δ∆S k = 0. To check this, note that the infinitesimal transformation δΦ i (t, θ,θ) =¯ Q Φ i , reads explicitly as: Thus, computing the variation ∆S k we get: The first and the last terms cancel exactly. Moreover, taking into account the physical condition (45) for the second and the third terms we get: Finally, integrating by parts the first term, we obtain: Because the regulator R k (t) must be even, and the bracket vanish 6 . Therefore, δ∆S k = 0 and supersymmetry is ensured from time 6 This can be specially checked for the regulator (62), which reads explicitly as: reversal symmetry. Moreover, from the way we derive the relation (54), but working with the coarse-grained partition function, we get the identity: ensuring that the relations between correlations functions coming from supersymmetry, and valid for the IR theory, remains true along the RG flow, at least as long as no singularities are encountered.
4. Choice of the regulator. The derivation of the flow equation proceeds in two steps. Following [35] we choose the one-parameter family of regulators: which is compatible with causality and whose Fourier transform can be computed analytically 7 : The Fourier transform (using the same symbol to designate the function and its Fourier transform) being defined as in the rest of this paper as: The parameter α ∈ R can be adjusted numerically to optimize the truncation. Moreover, the condition (45) allows computing R k : where we assumed R (2) k (t) even in time from its definition below (34). Explicitly: Let us examine the boundary conditions. For k → Λ 1, we designed the first regulator R k=Λ (ω) ∼ Λ, ensuring that fluctuations are frozen in the deep UV limit. In the same way, however, R k→Λ (ω) → −Z Λ α; and the initial correlation noise is not exactly recovered. Nevertheless, we expect that such an initial condition, which has only for effect to modify the initial value of D → D(1 − αZ Λ ) does not affect the universal behaviors of the system, at least as long as αZ Λ = 1. Finally, for k → 0, the two regulators both vanish ensuring that all fluctuations are integrated out.

Multi-local truncations for p = 3 models
We denote this section to the construction of approximate solutions of the exact flow equation (42) which cannot be solved exactly excepts for very special cases. We are aiming to construct approximate solutions able to deal both with the nonperturbative regime and non-locality of the interactions. The difficulty to solve the exact RG flow equation can be traced back to the fact that it describes a trajectory into an infinitedimensional functional space, thus we resort to approximations that consist in truncating it to form a finite-dimensional subspace. Here, the non-local structure of the interactions provides a non-conventional difficulty. To deal with it, we introduce a method inspired by the replicated construction considered by [82,86,87,88,89] for the random field Ising model, organizing truncation as a hierarchical multi-local expansion but keeping in mind two differences with their construction: 1. The role of the discrete replica index is now played by the continuous super-coordinate z.
2. The coarse-grained dimension is time rather than space.
We introduce the methodology for the unconstrained model and close this section by considering the spherical model.

Parametrization of the theory space: minimal bi-local truncation
As the perturbative expansion shows straightforwardly, the non-local structure of the interaction term in (38) contaminates the free energy W := lnZ, which becomes a nonlocal function of the external sources J as well. In such a way, we expect an expansion of the form: The origin of this expansion can be easily understood as follows. The potential in the classical action (38) writes schematically as a sum of two terms: , the dotted edge materializing the contraction of internal indices i. Thus, from standard perturbation theory and properties of Gaussian integrations, the loop expansion in power where the thick edges materialize the Wick contraction with propagator defined by the kinetic part of the classical action (38). Note that we voluntarily omitted the combinatorial factors, signs and additional numerical factors as well as classical (i.e. tree) contributions, irrelevant for our purpose. Finally, the upper index p in Now, we are aiming to construct an approximation of the exact equation (42) through a suitable parametrization of the full theory space able to deal with the non-local structure of the interactions. Moreover, we can only consider a finite-dimensional region of this theory space. This can be achieved through an ansatz for the effective averaged action Γ k [M] retaining only a finite set of couplings. This ansatz moreover must be able to deal with some exact constraint coming from the partition function (37). In particular, because we focus on equilibrium, we consider a truncation compatible with the time-reversal symmetry, i.e. with the transformation (44). Moreover, assuming that supersymmetry holds, we have only to retain linear interactions inφ (this can be viewed as a requirement imposed by the fluctuation-dissipation theorem as well, see [35]).
We choose to work with notations that make the supersymmetry explicit. To simplify, we denote the classical fields using the same symbols corresponding to the random variables, and view M = (φ, ψ,ψ,φ) as a superfield. Moreover, we define the operators D = ∂ θ − 2iθ∂ t andD = ∂θ, such that D 2 =D 2 and {D,D} = −2i∂ t . We thus consider the following ansatz for Γ k , keeping only first order derivatives in time and bi-local contributions for the potential: with: To return to the original variables, we expand the functions Y k (M(z)), and use the standard properties of the Grassmann variables. After a tedious calculation, the kinetic part of the truncation (65) (keeping only the quadratic terms involving a first derivative with respect to time) writes as: where (Z k ) 1 2 i must be a function of φ only. Explicitly: This equation means that the own wave function renormalization factors Z k,φ i and Z 1 2 k,ψ i ψ i respectively associated with the fieldsφ i , φ i and the product (ψ i ψ i ) 1 2 (without summation) must be equal: This nonperturbative constraint comes from time reversal symmetry, and drastically reduces the number of independent parameters in the truncation. The remaining part of the truncation, that we call "interaction" writes as: where the symbols (X k ) ijl depend only on φ and are defined as: This is the second constraint coming from time reversal symmetry. In the rest of this paper we focus on the simpler approximation which consist in neglecting the dependence of the wave function renormalization on the classical field φ. In such a way (X k ) ijl (φ) = 0.
Finally, rescaling the fields as M i → (Z k ) 1 2 i ; and using the same symbols to designate the rescaled fields, we get the simplified form: We focus on rank three tensors (p = 3) and quartic deterministic forces (M = 2, see (1) and below). The bare potential W 0 writes as: where the deterministic part V is: In the rest of this paper, we will use the notation h ≡ Z k h 1 . Note that due to the integration over the disorder, the higher powers in fields are always even, and thus have nonvanishing Witten index [81,94], meaning that supersymmetry is expected to be unbroken dynamically. Moreover, we focus on the large-N limit, discarding irrelevant contributions in 1/N . Due to the non-local nature of the interaction, we have to choose a suitable parametrization for the running effective potential U k [M]. As a first approximation we consider sixtic truncation; imposing: As a first step we choose a parametrization build as a sum of monomial interactions reproducing essentially the structure of the first effective loop effects deduced from the potential (72). Explicitly: The initial conditions are given by the potential (72). In particular, for some UV cut-off Λ in frequencies (assuming Λ 1), Physically, the bi-local terms must includeφ-φ products. These terms renormalize the Gaussian measure for the response field, which is nothing but the effective noise correlation. In contrast, the local terms must include only a single fieldφ. They correct the effective driving force. Note moreover that we chose the power on N in front of each interaction accordingly with the extensive nature of the effective action Γ k . Finally, note that we discarded 2-point couplings of the form: With this approximation theφ−φ component of the effective propagator remains diagonal in the internal indices (see (83) below). It can be justified from the observation that the corresponding coupling must be of order (u 4 , a sub-leading order correction with respect to the other contributions. Moreover, even though we provide a formalism including odd vertices (which come from the fact that J 0 = 0), a large part of our numerical analysis focus on the centered distribution J 0 = 0, discarding in that way all the odd contributions as u (1) 3 . This will be especially the case for our investigations about the spherical model, see Section 5. The same observation holds for terms such that: and the truncation (75) can be justified semi-perturbatively in J 0 8 . It reduces the investigations into the interior of a (despite everything vast) subregion of dimension 15 of the full theory space spanned by the different couplings involved in it (including Z k ). In the rest of this section, we will derive and solve numerically the resulting flow equations.

Solving RG flow equation: 0 and 1-point functions
In this section, we illustrate the computations of flow equations for the couplings u 1 and u 0 , before introducing a more systematic procedure in the next subsection.

Computation of the effective propagator
As a first step we derive the explicit expression of the effective propagator at scale k, G k (ω). Denoting as Φ = (φ, φ) the field with 2N components bringing together the entire bosonic sector, the corresponding components of the effective propagator G k must read in Fourier space: each component (G k,ΦαΦ β (ω)) being a N × N matrix. Note that it is suitable to work in the frequency space, and integrate over frequencies ω rather than on time t. To compute each component, we start by computing the inverse matrix from the truncation (75). We have: where A and B are N × N matrices given explicitly by: and: where R (q) The matrix can be inverted as a 2 × 2 block matrix; .
Remark 1 Two-points function for the response field. The fact that Gφφ = 0 in (79) (or equivalently that Γ φφ = 0) is not a simplification but an exact relation [34,4]. It can be traced back to the observation that adding a linear driving force: is equivalent to translating by −ik i (t) the i-th component of the sourcej for the response fieldφ into the partition function (30) before averaging over disorder: Therefore, becauseZ W 0 [j = 0,j = 0] = 1 we must haveZ W 0 [0, −ik] = 1 ∀ k. Thus, for vanishing source j:

Flow equations for u 1 and u 0
We denote asR k the super-matrix in the boson-fermion space: the regulator for fermions R Note that from this point we reserve the dot for derivative with respect to ln(k),Ẋ = kdX/dk and we denote with a " " the time derivatives, X := dX/dt. The flow equation for Γ (1) k α can be deduced from the Wetterich equation (42) taking the first derivative with respect to M α . We get (note that from this point,Ẋ = kdX/dk and X := dX/dt): the notation STr for "supertrace" meaning that we trace over all indices, with a global minus sign in front of fermionic contributions, and Γ k,···α denotes the 3-point function meaning that the third index is fixed to be α. To make the projection on both sides, we notice that: and 9 : In such a way, Γ . Therefore, to deduce the flow equation for u 1 , we have to expand the right-hand side of the equation (89) at the zeroth order in M, and identify the terms of order zero on both sides. Note that α must belong to the auxiliary fieldφ; so that from (91) the two remaining free indices of Γ (3) k,···α have to belong to the bosonic field φ. Thus, because the supermatrixR k is diagonal in the boson-fermion superspace, and that we have to keep only the zeroth order term on fields in the expansion of 1/(Γ (2) k +R k ), the supertrace reduces to a trace over the bosonic sector. Straightforwardly, we get: and denoting as G k,αβ (t) the (symmetric) inverse of the matrix (Γ k +R k ) αβ evaluated for M = 0, we thus obtain: the remaining trace tr running over the internal Latin indices, and V (3) denotes the "vertex matrices" with elements: To go beyond, we have to specify the choice of the regulator. Note that the undefined integral over all times dt in (93) must be formally identified with δ(ω = 0); δ(ω) being the standard Dirac function. Such an undefined factor appears generally on both sides of the flow equations, and are formally canceled as an artifact of the momentum conservation. Thus, accordingly with the definitions (59) and (62), the first integral on the right-hand side of equation (93) reads as follows: where: and following its standard definition we introduced the running anomalous dimension η k as: The derivative of ρ (1) k and ρ (2) k with respect to the flow parameter ln(k) can be easily computed and leads to:ρ and:ρ From dimensional analysis, it is suitable to introduce a dimensionless version of I 2 , sayĪ 2 , such that I 2 =Ī 2 /k, with:Ī where we introduced the notation ρ (i) (x) := ρ (i) k=1 (x). In the same way the explicit expression for G k,φφ (ω) can be straightforwardly computed, and it reads as: with: and g (2) k,φφ (ω) := The second integral on the right-hand side of the equation (93) reads as 4πδ (0) the kernelsJ being defined as follows: 1,φφ (x) + 2πl and:J 1,φφ (x) + 2πl 1,φφ δ(x)) .
In this equation we abusively used of the notation G 1,φφ (x) to denote the diagonal component of the 2-point function, namely: where h =: kh. Moreover, note thatl In terms of these dimensionless couplings, the flow equation forū 1 reads as follows: In the same way, we get for u 0 : The effectiveφ-loop vanishes because G k,φφ = 0. However, the stability of supersymmetry ensured by W-T identities requires u 0 = 0, and we expect that the two remaining contributions cancel exactly. We will check briefly this point. First, the bosonic contribution writes explicitly as: In the same way, because of G k,ψψ (t − t ) = −G k,ψψ (t − t) and taking into account that from (88) R k (−ω), we get for the fermionic loop: .
(110) thus, up to the change of variable ω → −ω into the last integral, we explicitly showed that L B + L F = 0, ensuringu 0 = 0.
Finally, note that we can easily make contact with the choice of the regulator done by some authors who studied stochastic equations from nonperturbative RG, for instance, [81]. The apparent contradiction seems to come from the fermionic regulator, but the relation between our choice and the choice done in the literature can be stressed from the relation (61) coming from time-reversal symmetry:

"φφ " 2-points couplings
Taking the second derivative on both sides of equation (42) with respect to the classical superfield, we get: with: and Because of the constraints coming from relations (69), it is suitable to set α = β =φ.

The relevant component for the 4-point function is
, which can be explicitly computed from the truncation (75) as: From this expression, we deduce that Γ (2) k (4) ϕ kφl must have the following structure: where the solid edges materialize the Kronecker deltas, and the dotted edges weighted with a grey bubble materialize the effective propagator G kṘk G k . Because of the 1/N arising from u 4 ; it is not hard to convince that only the third term, which create a closed trace in the Latin indices, survives in the large N limit. Computing it, we get: In the same way, we can represent the contribution Γ (2) k (3,3) ϕ kφl as: where the white circle materializes the propagator G k , the black dots and square correspond respectively to the fields φ • andφ • and the solid edge materialize the scalar product between fields. In such a way, groups of three isolated dots must correspond to interaction with coupling u 3 and groups with an edge to interaction with coupling u 3 . Because of the way the internal indices are contracted between the remaining free fieldsφ • , the first five diagrams must contribute to the flow of ∆ (3) whereas the last diagram must contribute to the flow of the field strength Z k . Let us consider the first diagram. Explicitly we get: where the traces "tr" run over internal indices and Ω denotes to the external momenta. However, because we truncate around local couplings, we have to keep only the first term in the "local" expansion in power of Ω: where derivative interactions includes derivatives of the Dirac delta. Our truncation scheme (75) focuses on products of local terms; however, we have to be careful with the component G k,φφ of the effective propagator (84). Indeed, omitting the internal index structure, this component reads schematically as: Two local contributions linked by the component denoted as A provide an effective local contribution; when the same local contributions linked together by the component B provides us with a product of local contributions again.
Focusing on the local approximation, from the explicit expressions for the components of the propagator, we get for instance in Fourier variables: 3 + 2(J (for (i, j) ∈ {1, 2}) arising from loop integrals being defined as: where: and: Following the same strategy for the remaining diagrams, we obtain: and finally: Note that we have to be careful with the positions of gray and whites bubbles, which are for instance responsible for a factor 2 in the first contribution to the left-hand side of equation (132) concerning the second contribution.
The flow equations for higher order interactions can be obtained in the same way. The detail being unnecessarily technical, we have placed it in the Appendix (A). We introduce for this purpose a diagrammatic method, which could easily be automated by a numerical routine for further investigations (i.e. for higher truncation).

First numerical investigations for N → ∞
In this section we provide a first look at numerical investigations of the formalism that we introduced in the previous sections, and we plan to devote a companion paper for a deeper numerical analysis. For this first look we focus on the large N limit, and set J 0 = 0. The last condition discard from our analysis all the odd couplings u and so on. The remaining flow equations split in two sets. The first set, which reduces to (ḣ,u (1) 4 ) is such that, even if initially (h, u (1) 4 ) = (0, 0), the corresponding flow equations does not vanish due to u 6 . The second set on the contrary includes the remaining couplings whose flow equations remain zero if we set vanishing initial conditions. For our analysis, we discard this last class of couplings. Moreover we set α = 1, which simplifies even more the original complicated set of flow equations that we previously derived and which reduces now to: Note that the last equations implies that in the large N limit u 6 remains constant. First, let us investigate the fixed point structure of the coupled system (ḣ,u 4 ), by considering u 6 as an external adjustable parameter. This system can be easily triangulated, the conditionḣ = 0 definingū    the fixed point solution for several values ofū 6 . Note that we focused on the region h < 0, aiming to describe the spherical regime; andū (1) 4 > 0 as required by stability. On the left, the Figure 1 shows the behavior of F (h,ū 6 ) for some values ofū 6 , which has to be negative to be interpreted as a disorder effect, see (72).  4 goes to zero. Asymptotically, and exploiting the global invariance reparametrization for fields, it is suitable to fix the asymptotic value of the effective non-zero vacuum as q 2 = N . Figure 1 on the right shows that the Wilson-Fisher fixed point does not exist for all values ofū 6 . Indeed, recalling thatū 6 < 0, we show the existence of a minimumū * 6 ∼ −0.037. Below this value, the fixed point disappears, which is illustrated on the left of Figure 1. Because |ū 6 | increases exponentially with ln(k/Λ), the fixed point is expected to eventually disappear in the deep IR, and the phase transitions will be first-order [26]. The behavior of the RG flow below this limit is illustrated in Figure 2 (on right), where we show that trajectories are not bounded from below. The nature of the transition can be inferred by studying the behavior of such a trajectory, for the nonzero disorder. Figure 3 (on the bottom) shows the typical evolution of mass and quartic coupling starting in the vicinity of the Wilson-Fisher fixed point with very small disorder (ū 6 (Λ) ∼ −10 −4 ). It shows that, while the behavior of the mass h remains quite similar with its behavior without the disorder, the behavior of the coupling, in contrast, diverges at finite RG parameter t = ln(k/Λ). The origin of the divergences can be traced back to the behavior of the integralsĪ 2 ,J 2 ,Ī 3 ,J 3 involved into the flow equations.   4 (t) forū 6 (Λ) = −10 −4 (solid blue curve) andū 6 (Λ) = −10 −5 (dotted yellow curve) respectively, for t := ln(k/Λ). As the disorder decreases, the singularity comes later and ends up disappearing when the disorder is canceled out. Note that this kind of divergence already occurs when the disorder is canceled out, along the attractive direction (in green in the Figure 2) pointing towards the negative masses. The flow diverges in that case because it reaches the singularity of the integrals. As the disorder |ū 6 | remains smaller enough some trajectories have time to escape the divergence ( Figure 5 (on right)). However, as soon as it is large enough, and in particular for |ū 6 | > |ū * 6 |; trajectories fail to converge. Note that this happens at finite scale t * = ln ū 6 (Λ) u * 6 . Moreover, it diverges earlier and earlier as the disorder increases, and the supersymmetric RG flow fails to provide a suitable description (i.e. Ward identities (57) cannot be integrated along with the flow beyond a finite time interval, which collapses to zero). Because of the relation between time reversal symmetry and supersymmetry, we expect that such a divergent behavior signals the end of the equilibrium regime 11 , and a dynamical breaking of ergodicity, which is not controlled by a fixed point and thus have to be of the first order. Note that because the flow for coupling is unbounded from below, the ground state fails to be normalizable (at least for the quartic truncation). Finally, note also that trajectories in the vicinity of the Wilson-Fisher point always converge for the zero disorder case, as Figure 5 (on the left) shows. Moreover, note that, for u 6 = 0, 4 > 0. Such a conclusion obviously has to be confirmed with more sophisticated approaches. For instance, our field expansion around zero classical field does not seem truly suitable to describe the spherical model. Indeed, one could envisage a development around the vacuum of the local part of the effective potential. This, essentially, corresponds to V [Φ] in equation (72). In terms of the field φ, it corresponds to a linear interaction for the response field ∼φV [φ] + (ψψ)V [φ]. Let us focus on the quartic potential (72). It has two minima, for φ 2 = 0 and φ 2 = −2Dh 1 N/h 1 , both making the linear coupling vanish in the response field. The solution φ = 0 moreover extremizes the second derivative V [φ], V [φ = 0] = 0. Because we cannot use O(N ) symmetry to choose the non-zero vacuum along a single axis, we impose φ 2 = r 2 N and N i=1 φ i ∼ 0. With this respect, for h < 0, the vanishing solution becomes a minimum of Thus, the zero vacuum is expected to be well convergent with this respect (see [81,92,93]). Let us briefly consider the non-zero vacuum φ 2 = r 2 N for the following truncation: where Γ k,kin is again given by (67) and with: Taking the second derivative with respect to classical (super-)field of the Wetterich equation (42), we get schematically: From the truncation that we consider, odd functions can be discarded in the large N limit. The computation of G −1 k,φφ follows the strategy explained in the previous section, except that we have to keep into account the contribution of the quartic interaction evaluated for φ 2 = r 2 (k)N . In replacement of (81) and (82), we get for the diagonal parts of A and B: and which, evaluated for r 2 = −h/u 4 reduce to: and The off-diagonal terms behaving like ∼ The resulting flow equations can be derived following the diagrammatic expressions given in A and numerically investigated. Once again, finite time scale divergences appear, as Figure 6 shows explicitly for some initial conditions in the spherical region. Our methods remain highly descriptive, as we mentioned before, and should be improved to make useful physical predictions beyond this finite scale divergence effect, which is however a general feature for disordered systems [40,87,89].

RG for the p = 2 models
In this section we briefly consider the case p = 2 for which equation (1) reduces to: for J ij J kl = λ 2 2N (δ ik δ jl + δ jk δ il ), J ij = 0 and delta-correlated noise (see (2)). We especially focus on the quartic case, which goes toward the spherical model in the large N limit for κ 1 < 0 and κ 2 > 0, κ 1 /κ 2 = −O (1). Note that this model can be solved exactly [26,54,60], and is know to exhibit a weak ergodicity breaking, as long has the laps time in the two point correlation remains smaller than the age of the system. For long time, the system behave likely as a disguised ferromagnet rather than a truly glassy system. Second order phase transitions has been pointed out for this system, investigated for instance through the Ruelle formalism in [90]. We recover it here within our renormalization group formalism. In particular, we show that the divergences exist only for N → ∞, and is governed by a true interacting fixed point.
Integrating out the disorder, we get the following effective potential, in replacement of (72): The truncation that we consider is as soon as given by formula (71), but with the effective potential: The flow equations can be easily derived following the method detailed in the previous section (section 4). Note that because the distribution of the disorder is centered, all the odd functions vanish and η k = 0. For the other functions, we get from the strategy detailed in Appendix A: where J 1,φφ (x) + 2πl Remark 2 It is instructive to compare the equations (201) and (157), corresponding to tensor and matrix models respectively. In the first case, in the large N limit,u 6 → 0, and u 6 can be viewed as a constant parameter along the flow. By contrast,u 42 = O(u 2 42 ) in the large N limit, and the behavior of u 42 remains non-trivial. A moment of reflection helps seeing that this difference can be traced back to the fact that, for n ≥ 5 the flow equation for Γ For numerical investigations, we use the same strategy as for tensors, but in this case we will consider also solutions for finite N . The fixed point solution can be investigated by triangulation; using the first equation, and assuming u 41 = 0, we may extractū 42 ≡ u 42 (ū 2 ,ū 41 ), solving the second solution we get finally:ū 41 ≡ū 41 (ū 2 ),ū 42 ≡ū 42 (ū 2 ); and the fixed point solutions can be found by investigating the zero of the function F (ū 2 ) = u 42 (ū 2 ,ū 41 (ū 2 ),ū 42 (ū 2 )). Figure 7 shows zeros of F for different N . For N large enough, we get two zeros, labeled respectively as FP1 and FP2 on the Figure the remaining coordinate Z along the irrelevant direction v 1 being set to its value at the fixed point. To understand the role played by the fixed point FP1, it is instructive to investigate the behavior of the couplings along a typical trajectory, as we did for p = 3. Figure 9 shows the behavior of the eigencouplings Y and X on the plan of Figure 8 along the red trajectory. As for p = 3, we show that the trajectory exhibit a singularity at a finite time scale. This behavior, typical of disordered systems [40,87,89] is, as for the tensor cases interpreted as the signal that supersymmetry and equilibrium conditions cannot survive for arbitrary large time scale, the lack of analyticity ensuring that Ward identities fail to be integrated along RG trajectories. In order to connect this singularity to the concrete behavior of a vector obeying the Langevin dynamics in Equation 1, we plotted in Figure 11 the two-point correlation in function of the time for N = 1000 (considered as large N) and τ = 50. We observe the existence of two phases, a phase of convergence where the 2-points correlation is O(1), that is followed by a phase where the trajectory becomes divergent. We see in the Figure 11 that the divergent phase appears later for D = 0.01 (left) than for D = 0.1 (right), which is in accordance with the previous insight that a smaller D leads to a later divergence. Such a behavior moreover seems to be a consequence of the large N limit, numerical integration of the flow for small N showing, on the contrary, good convergence properties, as we can see on Figure 9 (on right). We observe a divergence of the two point correlation of q which seems to be related to the singularity identified theoretically.

Conclusions and outlooks
In this work, we provided a nonperturbative functional renormalization group framework based on a coarse-graining in time for a special kind of dynamical systems described by disordered Langevin equations. These results aim to bring new insights and expand our knowledge on the dynamics of gradient descent (in particular Langevin dynamics) within high dimensional non convex landscapes similar to glassy systems. The understanding of these dynamics in these landscapes is central and essential for multiple physical and concrete applications, such as biological [39] and neuroscience systems [11], material physical [51], optimization in data science [68], machine learning algorithms [62], etc. We focused on the equilibrium regime and constructed explicit supersymmetric approximate solutions within a self-consistent framework based on a multi-local expansion. Keeping interactions up to order two in that expansion, we were able to construct flow equations describing the evolution of effective couplings. Planning to provide a deeper numerical investigation of the resulting equation in a companion work, we mainly studied the case of a centered disorder distribution and for a minimal truncation in the large N limit. We moreover essentially studied the symmetric phase expansion on the p = 3 and p = 2 models, where the disorder is respectively incarnated as a rank 3 random tensor and as a random matrix.
In both cases, we showed that effective potential becomes nonanalytic at a finite time scale. This singular behavior implies that Ward-Takahashi identities fail to be continued beyond this point, and supersymmetry breaks down. Because of the close relation between supersymmetry and time reversal symmetry, we interpret this singular behavior as the signal of breaking of ergodicity, where the equilibrium description fails to be suitable for the system. Besides these similarities, the p = 2 and p = 3 distinguish from the nature of the transition. In the p = 2 case, we showed that the transition is controlled by a non-Gaussian fixed point whose disorder coupling is one direction of instability. In contrast, this fixed point disappears at a finite time scale for p = 3, and the transition must be of the first order. These observations moreover seem to be in accordance with literature [12,25,95] that used other methods. However, our conclusions at this stage remain very dependent on the approximations that we have used, and we have already planned several follow-ups to go further. For example, computations that we considered could be quite easily automatized, i.e. performed by a computer program; allowing to investigate larger multi-locals truncations. Another interesting avenue for the case p = 2 would be to work in the space of the eigenvalues of the disorder, the distribution of which at this limit is described by the distribution of Wigner [74,91]. In action, the disorder term will then behave like an effective moment term, for a 1 + 1 dimensional (super-)field theory with finite volume. One could thus consider a division into scale at the same time in frequencies and eigenvalues. All these questions and more will be addressed in our next work. Note that the case of non-equilibrium process and taking into account the case where N is finite will be considered in forthcoming investigation.

A Higher order flow equations: diagrammatic method
In this section we present the diagrammatic method discussed in the paper to obtain higher order flow equations.

A.1 Naive power counting
Let us consider an interaction with coupling g, made as a product of n fields φ, m fieldsφ and l time integrals. To be dimensionless, the dimension of the coupling g must be equals to: Note that if each time integral is the remaining part of an integral over the supercoordinate z, l must be equal to m, In the language of field theory, this means that all the interactions are relevant.

A.2 Graphical representation for higher flow equations
For the derivation of the higher flow equations, it is suitable to use of a graphical representation, like the one used for equations (116) and (118). To this end, we introduce the following rules: • To each field ψ,ψ, φ andφ we associate respectively the symbols , , • and .
• To each sum as i M α i := M α • we associate an isolated vertex of type corresponding to the index α.
• To each scalar product of type M α ·M β we associate a solid edge bounded by vertices corresponding to the indices α and β.
• The vertex corresponding to fields sharing the same superpath coordinate z are surrounded with a closed dotted path materializing the integral dz. We call bubbles these closed paths.
We call local bubble graph such a construction and Figure 12 provides an example of such a bubble graph. As an illustration, let us consider the coupling J = dzdz M 2 • (z) M(z) · M(z ) M • (z ) involved in the truncation (75). We have five vertices, two of them linked with a solid edge and the others being isolated. These vertices moreover leave in the Figure 12: A typical example of bubble graph, corresponding to an interaction involving 11 fields. Note that a bubble cannot contain more than one black square or a pair of triangles.
interior of two distinct bubbles, corresponding to the integrals dt and dt . Using this graphical representation, J decomposes as: where the numerical coefficient in front of each diagram's counts the number of equivalent configurations. In the rest of this paper, we use this graphical notation to represent as well as the effective interactions and the corresponding effective vertices Γ (n) k , the context avoiding ambiguities. Using this graphical notation, flow equations will take the form of an equality between sums of graphs. On the left-hand side the bubble graphs involved in the truncation, weighted with β-functions, and on the right-hand side a series of graphs built as bubble graphs contracted with effective propagators. We will denote as effective edges the effective propagators, and to simplify the notations we materialize them as a solid thick edge. Following these rules, the first term on the right-hand side of equation (118) must we read as: Figure 13: Local structure of a typical effective bubble graph build of four vertices and four effective edges, the local ones being materialized with a solid blue edge whereas the non-local ones are labeled with a dashed blue edge. Each vertex is labeled with the total external momenta Ω i . For effective bubble graphs a and b, all the external momenta satifies a global conservation: Ω 1 + Ω 2 + Ω 3 + Ω 4 = 0. For effective bubble graphs c and d however, external momenta partition in two disjoint conserved sets. Then, Ω 1 = 0 and Ω 2 + Ω 3 + Ω 4 = 0 for graph c, Ω 1 + Ω 2 = 0 and Ω 3 + Ω 4 = 0 for graph d.
By contrast, we call such a diagram an effective bubble graph. Note that we have to keep in mind that into a loop, one of the solid thick edges corresponds to the propagator G kṘk G k when the others correspond only to G k . The next step is to define a rule allowing to identify on both sides the contributions corresponding to the same β-function. To this end, we notice that the components of the effective propagators are of two kinds (see (121)), either local, proportional to δ(ω+ω ) or non-local, proportional to δ(ω)δ(ω ). Now, let us consider an effective bubble graph involving a single loop of length n: L n = { 1 , · · · , n }. Each of these edges i can be local or non-local. If all of them are local, the resulting effective bubble graph is local as well (i.e. all external momenta are constrained by a single global Dirac delta). From connectivity, if one of these edges is non-local, the chain of local edges build a spanning tree, and the resulting diagram remains local. If at least two effective edges become non-local however, the resulting effective bubble graph becomes a product of local effective bubble graphs (i.e. the external momenta splits into two disjoint sets, both of them constrained by a Dirac delta). See Figure 13. We denote as local components each sub-local bubble graph. Finally, each of these local components can be expanded in powers of external momenta, the first term of this expansion corresponding to the ultra-local approximation.
As an example, let us consider the effective loop (168). In this equation, the loop depends on the external momenta Ω (see equation (4.3)), shared by the black squares, and do not correspond with a local bubble graph. From the previous discussion, we must have, formally: where K 1 denotes the local component built with local edges. Expanding the effective loop around vanishing external momenta Ω to obtain the ultra-local approximation, we get: and the ultra-local approximation of the effective bubble graph follows: Note that as soon as we projected into the superfield components, the bubbles materialize time rather than supercoordinate integrations. In general, we introduce the definition: Definition 1 Let us consider a graph G = {G α , a } built as a set of bubble graphs G α and effective edges a . The boundary graph ∂G of G is the bubble graph obtained from G by: 1. Deleting all the effective edges { a } as well as their boundary vertices.
2. We merge all the bubbles linked by at least one effective line.
Within this definition, it is clear that the boundary bubble graph of G is nothing but the local approximation of G. We thus define the following projection rule: From a given flow equation involving a few numbers of effective bubble graphs on the right-hand side, we must have to: 1. Identify on both sides coefficients in front of graphs having the same boundary bubble graph.
2. Expands the loops in the power of the external momenta and keeping only the zeroorder term of the expansion.
As a first application of the graphical method that we propose, we consider the flow equation for the coupling h, which behaves like a mass term in the field theory vocabulary. From equation (114), considering the pairφφ for external fields; we get: Thus, keeping only diagrams having the same boundary graphs on both sides, we arrive to the equation: where on the left-hand side we assume to keep only diagrams having the same boundaries as the bubble graph on the right-hand side. Note moreover that we discarded the combinatorial coefficients in front of each bubble graphs for convenience. The first two diagrams have been computed previously, we get: 2 ) + 16 kJ and finally: where:Ī 3 = iZ 3 k dx 2π (η k ρ (2) (x) +ρ (2) (x))(G 1,φφ (x)) 2 G 1,φφ (−x) , and:J 1,φφ (x) + 2πl (1) 1,φφ δ(x)) .
Remark 3 One can suspect that imaginary contributions occur due to the presence of factors i. However, it is not hard to check that purely imaginary contributions cancel for each loop. This comes from the observation that denominators can be always rewritten as even functions of ω. Because each ω share a factor i, imaginary contributions are even in ω in the numerator and cancels by symmetry.
In the same manner the diagrammatic representation of the flow from u −u (1) 5 3 ) 2 N 4 + .