Slepian models for Gaussian Random Landscapes

Phenomenologically interesting scalar potentials are highly atypical in generic random landscapes. We develop the mathematical techniques to generate constrained random potentials, i.e. Slepian models, which can globally represent low-probability realizations of the landscape. We give analytical as well as numerical methods to construct these Slepian models for constrained realizations of a full Gaussian random field around critical as well as inflection points. We use these techniques to numerically generate in an efficient way a large number of minima at arbitrary heights of the potential and calculate their non-perturbative decay rate. Furthermore, we also illustrate how to use these methods by obtaining statistical information about the distribution of observables in an inflationary inflection point constructed within these models.


I. INTRODUCTION
The low energy description of many higher dimensional theories involve a large number of fields (moduli fields) that need to be stabilized. This is normally achieved by the existence of a potential that fixes the values of these fields to a local minimum of that potential function. A typical example of this procedure can be found in String Theory compactification scenarios. In particular, models of flux compactification have been shown to lead to an enormous set of possible 4d potentials that can have many local minima. The typical number of moduli fields in these cases is quite large, reaching often the order of a few hundred. This makes prohibitively difficult to study these potentials in detail and one is forced to look for simple models where the field space has been truncated to a small subset of fields. Alternatively, one can try to study these models by taking a more statistical approach, where the scalar potential is regarded as a random field whose sample space is the set of 4d low-energy effective potentials. These ideas have been pursued in relation to the study of the stability of critical points in these potentials in [1][2][3], as well as the description of cosmological models for the early universe in [4][5][6].
In many of these studies one is interested in particular points of the landscape such as, for example, a minimum with some value of its cosmological constant, or an inflection point with a particular set of conditions in its derivatives necessary for it to sustain inflation. However, depending on the restrictions imposed, it may be very difficult to obtain an example of the potential with these characteristics by producing random realizations of the scalar potential.
Indeed, metastable de Sitter vacua and inflationary points compatible with observations are very rare in generic landscapes, with probabilities scaling as P ∼ exp(−N p f ), where N f is the number of scalar fields in the theory, and p > 0 is a number of order one [7][8][9][10][11][12]. To obtain realizations with the desired properties, one can of course use a Taylor expansion around the point in question and take into account the probability distribution for its coefficients [13,14]. However this becomes quite complicated as one increases the number of fields and the field range that one is interested in 1 . Moreover, with this type of procedures it is not possible to capture correctly the global properties of the scalar potential, which are essential to study quantum decay processes in the landscape. Here we present a different strategy to generate these potentials that locally will be constrained to have a particular form, but that globally will still represent a faithful realization of the random landscape, the so-called Slepian models [16].
Several different methods have been suggested as a way to represent these random potentials in the landscape. In this paper we will concentrate on potentials described by Gaussian 1 For another method of generating a specific class of constrained Gaussian random fields, see [15].
Random Fields (GRFs). This is based on the assumption that the 4d potential can be thought of as a sum of many different terms, of classical and quantum origin, coming from the compactification mechanism rendering the final result a Gaussian random field. This type of models have also been studied in connection to the distribution of vacua and its stability [9,17,18] as well as inflation [13-15, 19, 20] in the landscape. As an illustration of the mathematical techniques presented here for the construction of constrained GRFs we develop Slepian models that are locally described by critical points (maxima, minima and saddle points) as well as inflection points and use these realizations to extract important statistical information about them.
In particular, we will first study the quantum mechanical stability of local minima in these landscapes. In order to do so, we will compute numerically the decay rate of these minima using the quantum tunneling techniques first described in a series of papers in [21,22]. The result of this quantum instability is the creation of a bubble instanton that interpolates between the false vacuum and the true vacuum states. Using these Euclidean methods one can evaluate the probability of this decay channel and therefore estimate the lifetime of any specific vacuum. The calculation of these tunneling events in a multidimensional potential is however notoriously difficult. Recently some work on this direction has been done in relation to the stability of vacua in models with large number of dimensions in field space.
It has been argued that the probability of the decay depends exponentially on the number of fields although the particular scaling is still uncertain [23][24][25][26].
In this paper we will study these tunneling events in models of Gaussian random potentials. In particular we are interested in studying the dependence of the tunnelling rate with the height of the potential at the false vacuum. For large values of the cosmological constant this calculation would be impossible without constraining methods, since the number of these minima is negligible compared to the minima at lower values of the field.
Our techniques allowed us to efficiently generate the same number of minima for different heights and have a good sample of cases from where we can extract statistical information.
The obtained distribution for the instanton actions S E (which determines the decay rate Γ ∼ e −S E ) is displayed in Figs. 5 and 6, where we found that the average dependence of the decay rate on the false vacuum height V fv is given by where U 0 and Λ are the characteristic energy and length scale (in field space) of our potential respectively. The distribution for the Euclidean action becomes increasingly peaked around its mean, and thus more predictive, for larger values of V fv . As we show in the main text, this enhancement of the predictability can be explained using Slepian models for very atypical extrema of the potential, such as high minima.
Our second application involves the generation of inflection points. These are some of the most likely points in the landscape where cosmological inflation can happen. However this does not mean that an arbitrary inflection point would lead to inflation. Obtaining a successful inflationary period consistent with the current cosmological observations still requires some amount of fine tuning of the potential around the inflection point. Therefore, to characterise the distribution of observables for these inflationary models in the landscape one should again use some sort of constraining method, and look at a particular set of non-generic inflection points. In the present paper we will explore the dependence of the observable parameters of inflation to its initial conditions in the landscape. In particular we will take the initial conditions for the fields to be the ones determined by the exit point of an instanton describing the transition from a nearby parent false vacuum. Note that in order to perform this analysis, one requires not only the knowledge of the potential around the inflection point but also its relation to nearby minima. Hence our method, which accurately captures the global statistical properties of the potential, is particularly suitable to carry out this investigation. It is worth noting that, to the best of our knowledge, this is the first time that an Slepian model for inflection points is presented in the literature. The effect of the tunneling in the initial stages of inflation has also been discussed in [14,[27][28][29].
The remaining of the paper is organized as follows. In section II we introduce the notation that we will be using for describing our random potential function as a GRF. In section III we will outline the method for generating constrained random potentials as Slepian models.
In Section IV, we implement these ideas for a 2d field space landscape and generate a large set of random potentials with a minimum at a specific point in field space. This allows us to compute the tunneling paths from these minima and determine the statistics of the decay rate. In section V, we condition the random potential to have an inflection point suitable for inflation, and study the effect of the initial conditions set by the tunneling process from a nearby minimum. We conclude in Section VI with some comments on the results and some further ideas that can be implemented with these numerical techniques. Some of the mathematical details and numerical proofs have been left for the Appendices. In the present work, unless otherwise stated, we will use reduced Planck units M −2 pl = 8πG/( c) = 1.

II. PRELIMINARIES FOR GAUSSIAN RANDOM FIELDS
In this paper we will take our random potential, V (φ), to be a Gaussian random field defined over a N -dimensional field space, which we will parametrize with the vector φ = {φ i }, with i = 1, . . . , N . Furthermore, we will consider the probability distribution for the random potential to be homogeneous and isotropic, so its covariance function will only depend on the distance between the points at which it is evaluated, in other words it is of the form We will additionally require the potential to have a null mean: In the rest of the paper we will evaluate our expressions using the following simple covariance function: for the case of N = 2 field space dimensions. The parameter U 0 sets the energy scale of the potential while Λ represents the correlation length in field space. Generalizing this construction to other covariance functions, or a different number of dimensions in field space is straightforward 2 .
In the following we will be interested in the value of the field and its derivatives at a particular point in field space, which we can take to be φ = 0 without loss of generality, and we will refer to it as the center of field space. In order to simplify the notation we introduce the following definitions for the value of the potential and its derivatives: Furthermore, we will denote the eigenvalues of the Hessian matrix by λ i with i = 1, 2 which will single out the directions 1, 2 in our field space. Note that the derivatives of the scalar potential are also Gaussian random variables, and therefore any collection of the previous quantities forms a Gaussian random vector. In Appendix A 4 we will give the expressions for the correlators between these different derivatives of the potential as a function of the derivatives of the covariance function C(φ). These correlations will play an important role in some parts of our discussions.

III. SLEPIAN MODELS FOR CONSTRAINED GAUSSIAN RANDOM FIELDS
A key point in our construction of the GRF rests on the fact that a conditioned GRF maintains its Gaussian nature. More specifically, homogeneous and isotropic processes (such as the GRFs we are dealing with) can be conditioned using the Kac-Rice formula [30] in order to obtain new mean and covariance functions which generate GRFs with the required 2 Note that this covariance function leads to a somewhat special form of the Hessian matrix for the minima in this GRF (See for example the discussion of this point in [18].) It would be interesting to check whether this could have any quantitative effect on the conclusions of our paper. constraints 3 . The models for stochastic processes dealing with conditional events and crossings where pioneered by David Slepian [16], and have thus been coined in the mathematical literature as Slepian models.
We can describe these constrained processes in a generic form in the following way. For simplicity, let us consider first a Gaussian random p-dimensional vector, composed of jointly Gaussian variables, x T = (x 1 , . . . , x p ), whose probability distribution function (PDF) is given by, where µ = x is the mean vector and Σ is the covariance matrix, whose elements are given by with a, b = 1, . . . , p.
Let us now consider the following decomposition of the random vector where x 2 are p c components of the vector x that will be constrained by a condition x 2 =x, and x 1 are the remaining p − p c unconstrained elements. Then one can show [30,31] that the distribution probability for x 1 holding x 2 fixed to the desired values is given by, which shows that the distribution for the variables x 1 is indeed a Gaussian distribution but now with a mean and covariance functions given in terms of the original ones as where µ 1 and µ 2 are the means of the vectors x 1 and x 2 respectively, and This is possible because one can always find a new Gaussian random vector x = (x 1 , x 2 ), connected to the original one with a non-singular linear transformation x = A · x, such that We show in Appendix A 2 a proof of this statement. In the rest of the paper we will use this fact in several different ways, applying this technique for Gaussian random vectors made of different quantities of our potential.
In this section we will use the methods described earlier to generate a Gaussian random field with a critical point with a specific height at the center, φ = 0. In other words, we will find a description of the new GRF conditioned so that the point at its center satisfies the following properties: V (0) = u and V i (0) = η i = 0 for i = 1, 2. In order to do this we will follow the prescription used in the mathematical literature for maxima in GRF [32] and adapt it to our case. Let us start by introducing the following Gaussian random vector: where we denote by φ a , with a = 1, . . . , q, the position in field space of a discrete set of q points. One can show that the Gaussian random vector x has zero mean, and a probability distribution that can be readily computed using the form of the covariance function and its derivatives. This is a somewhat lengthy calculation and we have given the general expression in Appendix A 6. According to the description for constrained Gaussian random vectors given above this is all we need to obtain the new mean and covariance function for the new conditioned vector (and thus, also for the constrained GRF).
Using the results in Appendix A 6, one can show that the new mean function for the GRF with the constrained conditions is given by, This result corresponds to the particular choice of covariance function in Eq. (3), and is written in terms of the the value of the field V (0) = u and the eigenvalues of the Hessian matrix at the center, λ i , which are to be drawn from the distribution in Eq. (12) below.
The new covariance function is which is no longer homogeneous, but it is still isotropic.
It is important to note that the eigenvalues of the Hessian are not statistically independent of the height of the potential. This is intuitively clear since, for example, one would expect the typical minimum at a large height to be quite shallow compared to the minima situated well bellow the mean value of the potential. This expectation can be translated to the existence of important correlations between the field and its second derivatives at a point, and in particular at critical points. In order to take this effect into account one can calculate the joint probability distributions for the Hessian eigenvalues (λ i ) and heights (u) at critical points to obtain 4 where N is a normalizing constant. This distribution includes all types of critical points, namely maxima, minima and saddle points. Depending on the kind we are interested in, we simply need to impose positivity or negativity conditions on the values of each λ i .
Using these results we can generate a Gaussian random field with a critical point with the desired properties by the following procedure. Let us consider for example a minimum with fixed height u. Our first step will be to generate a set of eigenvalues drawn from the distribution (12) taking into account the value of u, imposing the non-negativity condition λ i ≥ 0, and fixing the normalization factor accordingly.
Using these values for λ i we can then generate realizations of the potential using the expression where we have denoted by ∆(φ) an inhomogeneous, zero-mean Gaussian random field whose covariance function is given byC(φ 1 , φ 2 ) in Eq. (11). We show in Fig. (1) an example of the different ingredients that make up a Slepian model for a local minimum in a 1d GRF.
We can use a similar procedure to generate other critical points, such as saddle points with different number of negative eigenvalues, by generating the appropriate samples of λ i 's.
An important conclusion that can be derived from the Slepian model (13), first noticed in [32], is that for highly non-generic extrema |u| U 0 (such as very low maxima or high minima), the shape of this GRF becomes very deterministic around the critical point, and it is described very accurately by the first two terms in Eq. (13). Indeed, one can see from Eq.
(11) that the standard deviation of the random component ∆(φ) is always smaller than U 0 , and that it approaches zero near the extremum located at φ = 0 (see also Fig. 1). Therefore the last contribution in (13) can be neglected in a neighbourhood of the extremum where On the other hand, in the limit |u| U 0 the eigenvalue distribution of the Hessian (12) is approximately given by 5 which indicates that in this limit the magnitude of the eigenvalues is very suppressed |λ i | U 0 /Λ 2 . Then, as we mentioned above, for highly non-generic extrema the decomposition (13) 4 See the calculation in Appendix A 6. 5 Note that for very high minima u > 0 and λ i > 0, while for very low maxima all signs are reversed. the study of Large Scale Structure formation in the universe (see e.g. [33][34][35][36][37]).

B. Slepian models for inflection points
As we discussed in the Introduction, we are also interested in inflection points in the landscape. The reason is that in a cosmological context these points could be one of the regions of the potential that give rise to a cosmological inflationary period. However, in order to be compatible with the latest cosmological observations, one needs to restrict the form of these inflection points. This leads us to consider an inflection point at φ = 0 as a realization of the GRF with a small gradient of the potential in the φ 1 direction, denoted by η 1 , and the rest of the coefficients of the Taylor expansion of the field around that point of the form The intuitive picture of these choices is clear, we are looking for a one dimensional inflection point that allows the slow-roll conditions to be satisfied along the direction φ 1 while the perpendicular directions have positive curvature. In other words, we are looking for a potential where inflation is effectively one dimensional locally. This also explains the last condition, which is imposed in order to allow for enough slow-roll inflation in the vicinity of this inflection point. This is admittedly a very particular form of the potential around the inflection point and, even though it could be interesting to identify this type of points in a GRF in other contexts, we have not seen any studies of this class of constrained points on GRFs in the mathematical literature. However, it is not difficult to follow a similar procedure to the one for critical points in order to obtain Slepian models in this case. The first thing we should do is to enlarge the form of our initial Gaussian random vector (9), since we now want to constrain not only first derivatives but second derivatives as well. This suggests that we should take the vector of the form, which, similarly to the critical point case, can now be conditioned to have the desired properties given in Eq. (15) Following the computations given in the Appendix A 7 one arrives to the result that a GRF with an inflection point at φ = 0 is described by the expression where Γ(φ) is an inhomogeneous zero-mean GRF with covariance functioñ In these expressions u, λ i and ρ ijk should be drawn from the joint probability distribution for heights, first, second and third derivatives of the potential at inflection points 6 where In the last distribution, the condition η 1 · ρ 111 > 0 should also be imposed if one is interested in 'inflationary' inflection points.
We have checked the accuracy of these distributions by numerically computing them from a large set of generic (unconstrained) GRF examples. We have identified all the inflection points of our sample, and used this information to compute the distributions of the parameters of the inflection points we are interested in. See Appendix B for the details of these numerical checks, which are summarised in figure 14.

C. 2D numerical implementation
All GRFs generated for this work were constructed following the Karhunen-Love expansion (see e.g. [31]), which is briefly described in Appendix B. This algorithm generates values for a GRF discretized over a lattice which is to be interpolated afterwards.
Based in the criteria developed in [14], we used 5 lattice points per correlation length (25 per length squared). The resulting grid was then interpolated with fourth-order splines in order to analyse up to third-order derivatives of the field as faithfully as possible. The generated GRFs were found to follow successfully the initial mean and covariance function, as well as other properties such as the distribution of critical points and eigenvalues thereof.
Two examples of (rather extreme) GRFs generated following the steps in this section have been plotted in figure 2.

IV. TUNNELING IN A GAUSSIAN RANDOM LANDSCAPE
A Gaussian random landscape possesses a large number of perturbatively stable minima.
However, we know that quantum mechanically these vacua are not completely stable and can decay by the nucleation of a bubble of the new state. This means that a typical vacuum in our landscape will have many channels to decay into, each of them with a different probability. Here we would like to study the statistics of these decay channels in a controlled way by generating a large number of GRF realizations, and analyse their dependence on the parameters of the central minimum.
In order to do that we will use the instanton techniques first described by Coleman and collaborators [21] where it was shown that for a given minimum of the potential the decay probability per unit time and per unit volume is given by where S E is the Euclidean action for the bounce solution that interpolates between the new state and the original one 7 .
In the absence of gravity, one can show that the most likely decay channel is given by the O(4)-symmetric instanton solution in a 4-dimensional Euclidean spacetime; therefore we will be interested in solving the following set of Euclidean equations of motion where the prime denotes a derivative with respect to the radial coordinate in 4-dimensional Euclidean spacetime, r, and we have assumed that the fields φ(r) = {φ 1 (r), . . . , φ N (r)} are canonically normalized. Finally the boundary conditions are where φ F V is the location of the false vacuum in field space, the minimum of the potential from which the decay happens. Once the field equations have been solved, the action in the exponent of (21) reads Computing the coupled system of the instanton equations (22) is no easy task; particularly, as the dimensionality of the field space grows, the solutions tend to be increasingly unstable. There are, however, several publicly available algorithms in the literature to tackle the problem (see, e.g., [38,39]); additionally, some alternative methods have been recently proposed to find the action and escape point for the instanton, as in [40,41].
In this work, we use AnyBubble [42] to compute the instanton actions for our realizations.
AnyBubble is a Mathematica Package based on efficient numerical methods for the solution and optimization of the tunneling equations, see [42] for details.
In order to obtain statistics of the tunneling action in terms of the properties of the central minimum, we sampled false vacua with heights between -2 and 5 (in units of U 0 , see Eq. (3)) in uniform intervals. As explained in [14], we can write the Euclidean action as so thatS corresponds to the Euclidean action of a potential with covariance function (3) with U 0 = Λ = 1. Unless otherwise specified, all histograms corresponding to the action are given in terms ofS due to numerical simplicity.
Following the procedure of the Slepian models described the previous sections, for each value of the false vacuum height, we generated 2 · 10 4 Gaussian random field realizations 7 Here we will not be concerned with the pre factor A. See [22] for a detailed description of its computation.  We can readily see the power of the machinery described in the previous section when constraining the field to have a minimum with a vacuum energy higher than 1.5U 0 . If we tried to find a minimum higher than that drawing samples from an unconstrained GRF, we would need to generate tens (if not hundreds) of random fields before finding a single minimum satisfying that condition, see figure 13(a) in Appendix B. For example, from equations (B5), we can easily check that the probability of any minimum being higher than 5U 0 is O(10 −16 ), so finding one by chance happens to be quite remarkable. With the aid of conditioning methods, we are able to construct very efficiently large samples of random fields subject to a condition as difficult to meet as this one.
In order to study tunneling processes on each generated example, we identified all the minima near the center of field space and computed the tunneling rate between the central minimum (which always acts as a false vacuum, in our analysis) to all lower minima. An example of this procedure is plotted in figures 3 and 4, where we show the paths followed in field space by the different instanton decay channels 8 . We have only considered tunneling to minima around the center to avoid problematic issues with minima close to the boundaries of our realizations.
A. Statistics of the instanton action 1. Dependence with the height Figure 5 shows the resulting distributions 9 for the tunneling action, for different values of the false vacuum height. There is an interesting correlation between the mean and width of this distribution and the height of the false vacuum. Namely, we find that the higher the false vacuum is the lower the action and thus, the higher the probability of tunneling is.
This behaviour is quite intuitive; as we can see from the examples in Fig. 2, tunneling from a minimum high up in field space requires crossing a lower barrier to the true vacuum, which in turn results in a lower action for those transitions. Figure 6 (blue dots) shows the median of each distribution along with the range of actions between the first and third quartiles. and third quartiles of each distribution, for the optimal path, the linear (straight-path) and the thin wall approximations, along with a fitting curve (see (26)).
We see, once again, that higher false vacua lead to accumulation over lower values of the action.
The obtained data for each potential height was found to be easily fitted to a log-normal distribution. More specifically, the logarithm of the median of each distributionS med (which, in this case, is very similar to the mean of log 10S ) can be fitted by the following expression where V fv stands for the height of the false vacuum. As we see from figure 6, increasing V fv reduces the width of the distribution significantly, thus increasing the predictive power of (26) for the expected value of the action. This enhancement of the predictability of the Slepian model for large values of V fv corresponds precisely to what we anticipated in the previous section. Indeed, there we showed that near high minima the random potential becomes dominated by the first term in the decomposition (13), and therefore the landscape is very deterministic in a neighbourhood of false vacua with large V fv . Consistent with this result, when studying the non-perturbative stability from these vacua we observe a reduction of the variance of tunneling actions for large heights of the false vacuum. This agreement also suggests that in the case of minima with a large V fv the value of the instanton action is dominated by the local structure of the minimum. We will provide further evidence for this claim below.

B. Approximations for the calculation of the action
Due to the inherent instability of the equations to be solved to compute tunneling profiles, it is clear that as we increase the domain and dimensionality of the potential under study, the required computational time to solve the system will grow accordingly. Evidently, this makes the study of higher-dimensional GRFs and their tunneling properties almost prohibitive in this sense. Motivated by these limitations, we turn to computing several different approximations of tunneling actions suggested in the literature, and compare them with our exact results.

Thin wall approximation
The thin-wall prescription was already discussed in the original papers by Coleman in [21]. In this approximation the instanton action is given in terms of the difference between potential at the false vacuum (V fv ) and true vacuum (V tv ) and σ, the tension of the wall interpolating between them, namely, This approximation is accurate as long as the difference between V fv and V tv is small.
We evaluated (27) for each bounce we previously found with AnyBubble in order to check this expression and its predictive power for GRFs. In the computation we restricted the field to a straight line in field space connecting the true and false vacua. Figure 6 shows the evolution of the median ofS tw as a function of the false-vacuum height. While the width and median of the distribution in this case follow the same pattern as the optimal action, the values diverge rapidly from the optimal ones as the false vacuum height increases. This is not too surprising since, as one increases the height of the false vacuum minimum, the field can tunnel to a minimum with quite different values of the potential, what violates one of the premises of the thin wall approximation.

Straight-path approximation
While the thin-wall prescription provides a solid upper bound on the bounce action [43], it does not provide any useful estimation on the actual value on the bounce in our case.
This fact calls for an alternative way to estimate the action, mostly for higher-dimensional landscapes.
A straightforward simplification to this problem was introduced in [44], which we will denote by straight-path approximation. This prescription is based on reducing the field space to a single straight line connecting the false and true vacua, thus making the problem of tunneling effectively one-dimensional. As can be seen from figure 3, this approximation may not be too unreasonable. Even though there are some paths which do curve over the field space, many (if not most) of them follow a straight trajectory in field space. Note, however, that this restriction in field space may yield effective potentials where the bounce does not exist or might even correspond to a different bounce in the full theory. For more details on the properties of this approximation, see [45].
For each optimal path, we considered a straight line in the two-dimensional GRF connecting the true and false vacua, and computed the corresponding estimate of the action, S sp , in each case. In principle,S sp represents an upper bound on the optimal actionS, as the former only considers variations of the action in the direction of the straight path [44].
It is thus expected (and explicitly shown in [45]) that this approximation will diverge from the full solution as the dimensionality of the potential is increased.
We found that in this case the distribution of actions in terms of false vacuum height is identical to the optimal one shown in Fig. 5, though slightly shifted to higher values. As we can see from Fig. 6, the change in the median is minimal when the straight-path approximation is considered. Although, as we just mentioned, the straight-path approximation is not expected to give precise results for potentials in a higher field space dimension, this result suggests that it would be interesting to explore the validity of this method with GRFs in higher dimensions. Indeed, due to the computational complexity of such an analysis, a rough statistical estimate of the decay rate obtained with this approximation would still be very valuable.

C. The lowest action
In many circumstances one will be interested in the lowest action for a particular kind of minima. This will of course correspond to the path that would dominate the decay for those minima. In this subsection we will investigate the characteristics of such trajectories in field space.

Exit angle
An intuitive way to think about the most likely decay process would be to imagine that the tunneling occurs along the trajectory with the lowest barrier. One can check this idea in our case by first identifying the angle (in our 2d field space), θ, that the instanton trajectory makes with respect to the direction of the lowest eigenvalue of the Hessian at the minimum.
A distribution of such angles obtained for different values of the height is plotted in figure 7.
We see that there is a clear tendency of the tunnelings to occur around θ ≈ 0 but the correlation is not very strong.

Estimating the lowest action
The correlation of the instanton path with the lowest eigenvalue direction at the false vacuum suggests that one can try to estimate the lowest action by analyzing the potential along the lowest eigenvalue direction alone. This has been recently proposed in the context of the landscape in [26]. In the following we will use our large sample of realizations to test this idea in detail in our 2d GRF model of the landscape.
In order to evaluate the instanton action along the lowest eigenvalue direction we first take a slice of the potential along that direction and fit it to be of the form, Note that this procedure does not guarantee that the resulting one-dimensional potential is suitable for a tunneling process. In fact, in many cases the potential constructed this way does not have a lower minimum along this direction and therefore it cannot be used to estimate the decay rate. In the following we will only compute the instanton action in the successful cases where this 1d truncation gives an acceptable form, what in particular requires ρ 111 < 0.
Considering this simple form of the potential as the most likely exit path for the decay transition we can estimate the instanton action. In order to do that we will use the parametrization of the Euclidean action for the bounce that was obtained by Sarid in [46].
In our notation this becomes, We show in figure 8   all high minima should look very similar in all realizations, with its shape dominated by the first term in (13). This explains why the distribution of instanton actions becomes more deterministic (Fig. 5) for larger values of V fv , and therefore also the agreement between the

V. INFLATION IN A SLEPIAN RANDOM LANDSCAPE
Up to now we have been using all the software and mathematical tools described above for the computation of bounce profiles and actions with Gaussian random fields conditioned to have a minimum at φ = 0. In this section, we turn to studying constrained GRFs with inflection points at the origin of field space focusing on their application to cosmological inflation.
More specifically, inflation around inflection points has received special attention for being capable of sustaining enough e-folds to make contact with observations, while taking place in a small region of field space with an effectively one-dimensional potential.
While most of the obtained results and distributions seem promising, they have only been tested within Taylor expansions around these points, instead of using full GRFs. As we mentioned before, such methods do not capture correctly the global features of the potential, what is essential for characterising the non-perturbative stability of vacua. Therefore, this procedure is unsuitable for studying models of inflation where the initial conditions are determined by the decay of a parent false vacuum.
In this section we will apply Slepian models to constrain Gaussian random fields to have inflection point with the desired properties to sustain inflation, and then we will study the dependence of its cosmological observables on the initial conditions, set by different realizations of the parent vacuum.

A. 1D Inflection point inflation
Let us briefly review the main results for one-dimensional inflection-point inflation (see [20,48] and references therein for more details). Let us consider a potential of the form, where, in order to satisfy the slow-roll conditions around the inflection point, we will assume that η u. Note that we do not need to assume that the third derivative is too small. In fact, following typical conditions for a GRF we will consider the case where u ρ. Taking this into account one can show that slow-roll inflation conditions will be satisfied in the which together with the condition u ρ implies that we are describing small field inflation.
Using the slow-roll conditions, it is easy to check that the expected number of e-folds, N exp , that can be achieved within that region is where = (V (φ)/ √ 2V (φ)) 2 and N max is the maximal number of e-folds achievable in the whole potential. Moreover, defining where N CMB is the e-fold number at which the CMB scales leave the horizon, the spectral index of scalar perturbations can be shown to be given by Finally, the amplitude of scalar perturbations can be expressed as where f (x, y) = cos 2 (x)(y tan(x) + 1) 2 x 2 (y 2 + 1) .
With these expressions at hand, we can easily obtain a set of parameters for the inflection point (u, η and ρ) that are in agreement with the current cosmological observations, namely, N exp > N CMB ≈ 50, n s ≈ 0.965 and ∆ 2 R ≈ 2 × 10 −9 (see Eq. (38) below).

B. Numerical inflection points in a 2D Landscape
We now want to embed 1d inflection-point inflation in our 2d GRF landscape. In order to do that we can follow the procedure explained in Section IIIB for Slepian models in the case of inflection points. In the notation introduced earlier, the 1d parameters η = η 1 and ρ = ρ 111 , correspond to the derivatives along the flat direction of the multidimensional inflection point. Note that, in principle, u and ρ 111 (when evaluated at the same point) are uncorrelated, but the same is not true for u and the second derivative along the inflaton direction λ 1 ; similarly η 1 and ρ 111 are also correlated, see Eq. (20). Here we are interested in studying the global properties of the landscape on the cosmological observables so we will focus on a particular type of inflection point where we have fixed its 1d parameters 11 .
Following the steps from the previous section, we built two-dimensional GRFs with an inflection point whose inflating direction has fixed features. In the forthcoming sections we where U 0 = 6.0 · 10 −16 M 4 Pl and Λ = 0.5 M Pl define the energy scale and correlation length respectively, with the Planck masses written explicitly for clarity.
Once u, η and ρ have been fixed, using the probability distributions listed in (19) and (20), we can obtain the remaining parameters of the two-dimensional inflection point set at 11 It is also interesting to study the effects of varying these parameters together with the global properties of the GRF. We leave the details of this calculation for a future publication. the origin of field space φ = 0, and generate in a efficient way a large sample of GRFs with the listed properties. 12 As an example, we show in figure 9 a field constructed with the above constraints. We then used AnyBubble to tunnel from a higher false vacuum to the central inflection point.
We note that even though in every realization the inflection point has the same properties along the φ 1 direction up to third order, the potentials are different away from that point.
This means that the false vacuum, which decays to the region around the inflection point, is located in a different place and it also has different features in each realization, e.g. vacuum energy and barrier height. Using AnyBubble we computed the exit points of a large set of realizations. After that we used these exit points of the instanton decay as the starting points of a Lorentzian evolution of a FRW universe with this potential.
In order to study the inflationary trajectory we used mTransport [49], a Mathematica code developed to compute inflationary observables using the transport method. The cosmological evolution inside of a bubble universe created from tunneling is described by an open FRW universe [50]. Here, for simplicity, we used the flat-space approximation for the evolution of the cosmological interior of the bubble 13 .
In the example from figure 9, the dashed line represents the tunneling trajectory, while the solid one marks the inflationary one. We found this path to sustain a total of 124.1 e-folds and a spectral index of n s = 0.964 at the observable scale.

C. Statistics of inflationary parameters
In order to test the method described above to generate inflationary random fields, we generated 5000 GRFs constrained to have an inflection point with the same properties as the one in the example of figure 9 (see Eq. (38)). Next, in each of these realizations, we found all minima lying above the central inflection point and used anyBubble to compute the tunneling trajectory from the former to the latter in each case. Considering the exit point as the starting point of an inflationary phase, we used mTransport to find the number of e-folds, power spectrum, tensor-to-scalar ratio, spectral index and its running. The distributions of the e-fold number and the spectral index are shown in figure 10, for a pivot scale of 50 13 Note that in reality the initial cosmological evolution is dominated by the spatial curvature of the open FRW slices that describe the bubble interior. This will have some effect on the initial stages of the evolution of the scalar field in a multidimensional potential. See [19,28] for a discussion of these effects. e-folds, whereas the action associated to the tunneling to the inflection point is shown in Fig. 11. This is a different distribution than the ones we found earlier, since the common factor in these decays is the final point and we do not impose anything about the initial (false vacuum) state. It is interesting to see that this distribution is quite peaked around an action of the order of 10 3 .
We have also obtained the distributions for the amplitude of scalar perturbations, tensorto-scalar ratio and running of spectral index which turned out the be centered around the values ∆ 2 R = (2.02 ± 0.04) · 10 −9 , r = (8.0 ± 0.1) · 10 −9 and α = (−2.49 ± 0.02) · 10 −3 , (39) respectively 14 . Our results in this section are fully compatible with the 1d studies in [14]. 14 The cosmological evolution of these Lorentzian trajectories continue after inflation until they reach a lower minimum. We have not fine-tuned this minimum to be in Minkowski space, so in general the evolution leads to eternal de Sitter or to an Anti-deSitter crunch. We are only interested in the statistics of the inflationary period so we have stopped this evolution after the field leaves the slow-roll regime. It is important to remember that all these realizations have the same 1d inflection point parameters. In order to extract the complete statistical information about the predictions of a particular GRF we should combine these results with the ones obtained from inflection points with other parameters with their correct statistical weight. This is a much more numerically intensive problem and we leave it for a future publication.

VI. SUMMARY AND CONCLUSIONS
Slepian models are a powerful mathematical technique for modelling the statistics of random landscapes conditioned to satisfy a certain set of constraints. For this reason they are particularly useful to characterise phenomenologically interesting corners of the landscape, e.g. de Sitter vacua or inflationary regions consistent with the cosmological data, which are known to have highly suppressed probability to occur in generic random potentials. On the one hand, Slepian models provide a way to generate numerically large samples of a random landscape containing the region of phenomenological interest to be studied, regardless of the low probability of the realizations. On the other hand, this technique can also be used as an analytical description of conditioned random potentials, and thus to obtain valuable insight about properties of the landscape around these regions of interest. A particularly attractive feature of Slepian models, as opposed for example to the use of Taylor expansions, is that they can capture the global features of the random potential, and therefore they are specially useful for studying quantum mechanical instabilities in the landscape. In this paper we have presented the mathematical techniques for studying conditioned Gaussian random landscapes. We have applied these method to condition a 2d random potential to have a de Sitter minimum with a specific vacuum energy and also to study 2d landscapes containing an inflection point capable of sustaining a period of inflation compatible with the data.
More specifically, regarding our discussion of de Sitter minima, we have considered the non-perturbative decay of these vacua to lower minima, and characterised the statistical distribution of their decay rate as a function of the height of the false vacuum. For this purpose we have used our Slepian model to generate numerically large samples of vacua with varying values of the vacuum energy, and then computed the corresponding decay rates both solving the full instanton equations, and using various approximate methods present in the literature: the thin-wall approximation [21], the straight-path approximation [44], and the estimate proposed by Sarid [46] for the lowest instanton action (see Eq. (29)).
Our analysis shows that the thin-wall approximation is in good qualitative agreement with the numerical results, but only provides an accurate estimate of the instanton action action for minima with a relatively small vacuum energy. Indeed, consistently with the thin-wall prediction of the instaton action, we observe that the decay rate increases (on average) for increasing values false vacuum height. This can be understood noticing that, in a Gaussian random landscape, the barrier height that needs to be crossed to escape from the vacuum decreases when the vacuum energy of the minimum increases. However, for minima with a large vacuum energy the tunneling typically occurs to much lower vacua, what violates the assumptions of the thin-wall approximation, and thus it cannot provide a good quantitative estimate of the decay rate.
In the straight-path approximation one assumes the decay is effectively one-dimensional, so that it occurs along the line connecting the false and true vacua. We have shown this simplification agrees remarkably well with the results of our full numerical analysis in all cases we studied in a 2d Gaussian landscape. It is interesting to check if this simplification still provides a rough estimate (see [45] for a discussion) for the instanton action in higher dimensional landscapes where the numerical resolution of the full instanton equations becomes prohibitively difficult. In the particular case of a Gaussian random landscape this approach is specially attractive, since the statistics of the random field along the straight line connecting the false and true vacua can be fully described by simply restricting its covariance function to that line. Therefore, if this method would prove useful to estimate the non-perturbative stability of vacua in large-dimensional Gaussian landscapes, it would not be necessary to produce a sample the full higher dimensional GRF, it would suffice to generate one dimensional realizations of the random field with the same covariance.
Regarding the estimate of Sarid [46] for the lowest action (the most likely decay chan- exhibit a very deterministic shape in a large region around it, which is dominated by the first term in equation (13). As we argued in the main text, combining the estimate of [46] for the lowest action, with the Slepian analysis one concludes that the distribution for the instanton actions should become increasingly peaked and deterministic for higher minima. Finally, it is worth mentioning that the present work has potentially very interesting applications to characterise the landscape of 4d effective field theories in String Theory flux compactifications at tree-level. Actually, as was discussed in [1], the superpotential defining the effective supergravity description of flux compactifications can be modelled as a (complex) Gaussian random field with a specific covariance function determined by the geometry of the compact dimensions. The superpotential encodes a large amount of information about the low energy theory: the critical points of the superpotential represent minima of the tree-level moduli potential; the supersymmetry breaking scale is given by its absolute value; and the eigenvalues of its Hessian encode the mass spectrum of the moduli fields and their fermionic superpartners. Thus, the conditioning methods presented in this paper can be immediately translated into this context, allowing to study the statistical properties of the 4d effective theory when constrained to satisfy one or various conditions (see [2,3]), e.g. the existence of a vacuum with a particular supersymmetry breaking scale, or to have a mass spectrum containing a certain number of light modes. 15 A realistic study of the observable parameters of inflation in this model should also include a prescription to calculate their probability distribution in the multiverse. This will require the introduction of a measure.
Here we have not discussed this issue any further. Throughout this Appendix, we will give a detailed description of the tools and derivations needed in order to generate conditioned Gaussian random fields, such as the ones we have been using throughout the main text. We will be mainly following [31,32].

Introductory remarks and some properties of Gaussian random variables
A random variable x is said to follow a normal or Gaussian distribution if its probability distribution function (PDF) is given by where µ = x and σ = x 2 are the mean and variance of the distribution, respectively.
Likewise, a p-dimensional vector x T = (x 1 , . . . , x p ) is defined as a Gaussian random vector (composed of jointly Gaussian variables) if every linear combination satisfies that is, it follows a normal distribution. The PDF of the whole vector is where µ = x is the mean vector and Σ is the (non-degenarate) covariance matrix, whose elements are given by

Conditioned Gaussian random vectors
Let A be a p × p matrix and x T = (x 1 , . . . , x p ) a Gaussian random vector. Then, by definition, is also a Gaussian random vector with mean µ and covariance matrix Σ . Since (A5) is a linear transformation, the new mean is given by whereas the new covariance matrix is or, more compactly, In order to introduce conditional probability notions to jointly Gaussian random variables, let us discuss some interesting properties of grouped random variables. If we split some Gaussian vector x into two parts, namely, then the mean vector and covariance matrix will also split accordingly: each block in Σ having the proper dimensions to accommodate the covariances among the vectors x 1 and x 2 .
With these remarks at hand, let us perform a linear transformation on x, choosing After some straightforward algebra, one can show that the new Gaussian vector y is whose associated mean vector and covariance matrix are meaning that the new y 1 and x 2 are uncorrelated and, therefore, independent.
Given a bivariate joint probability distribution function f (x 1 , x 2 ), the conditional probability f (x 1 |x 2 =x) is defined by [52] f ( Let x be a Gaussian random vector, a subset of which has been set to x 2 =x. We could, in principle, substitute the value of the variables x 1 into (A3) and proceed with the remaining (and normalized) expression. However, more interesting conclusions can be drawn if the above results are applied. Instead of working with x = (x 1 , x 2 ), let us use the PDF associated to y = Ax, where A is given by (A12): Fixing x 2 =x and applying (A16) to the resulting probability distribution function, we From the expression above, we can conclude that conditioned Gaussian random vectors retain their Gaussian nature with mean and covariance matrix given byμ andΣ respectively.

Gaussian random fields
The idea of Gaussian random vectors can be generalized to random variables dependent on a certain set of parameters. Instead of having p Gaussian variables, we will have an infinite amount of them; the mean vector and covariance matrix will thus transform into a mean and covariance functions.
A Gaussian random field (GRF) {V (t), t ∈ R n } is defined as a function satisfying at every point of its domain. The mean function will be given by µ(t) = V (t) whereas the covariance function must satisfy C(t, s) = V (t)V (s) . If C(t, s) = f (t − s) the GRF is said to be homogeneous; if, on the other hand, C(t, s) = g(t · s, |t|, |s|) the field is isotropic.
GRFs which are both homogeneous and isotropic are referred to as stationary, and satisfy In the main text, we will we working with this last type of covariance function.
Finally, note that any GRF V (t) with mean µ(t) can always be decomposed as where W (t) is a mean-zero GRF sharing the same covariance function V (t). This construction will be useful to construct GRFs numerically (see Appendix B).

Useful correlations
Since linear combinations of Gaussian variables are Gaussian as well, it is straightforward to see that the derivatives of Gaussian random fields at any point of their domain are Gaussian too. Some of the most important covariance functions relating different Gaussian variables are the following [31, sect. 5.5]: Let us change the notation to ∂ φ j V (0) = V j (0) and evaluate the previous expression for some useful cases: In the above expressions, α i , α ij and α ijk are numerical constants which depend only on the covariance function of the (unconstrained) Gaussian random field. Note that in the two-dimensional case α 222 will be absent from all derivations, since the indices appearing in the correlation function between the third derivatives can only take two different values.
Note also that odd derivatives of the GRF are uncorrelated with even ones when they are evaluated at the same point in field space. This is due to the isotropy of the covariance function: if it is written as a power series, only even powers such as φ 2 i , φ 2 i φ 2 j will be involved. Therefore, only those correlations which end up involving even derivatives of the covariance function are non-zero. This however, does not mean the fields V (φ) and, say, V i (φ) are completely uncorrelated.
If we evaluate them at different points in field space, it can be shown [30, theorem 2.3 therefore, a GRF and any of its derivatives are correlated as processes.
The multidimensional 16 Kac-Rice formula for this field gives us the expected number of 16 Note that this formula is only valid for fields mapping R n → R n .
times a certain event, say, V(φ) = u, happens in an interval φ ∈ I of volume V: where det V (φ) stands for the Jacobian determinant of the vector field 17 , that is, If the field is stationary, that is, homogeneous and isotropic, we can simplify the expression above. Denoting V 0 = V(0) and V 0 = V (0), we find, assuming ergodicity, where the integral is performed over the whole domain of V 0 and V 0 and P (V 0 , V 0 ) is the joint PDF of V 0 and its derivatives.
More than one simultaneous event can be considered in the expressions above by enlarging the vector V and introducing more Dirac deltas representing each event 18 .
While the above expression can certainly be used to obtain the number of times a certain event happens in a given interval, it can also be used to obtain distribution functions. More specifically, applying ergodicity theorems, it can be shown [30] that the probability of an event A happening, given that B has happened, that is, P (A|B), can be obtained by If A depends on continuous parameters (such as the position in field space of the GRF), then the expression above represents a probability distribution function.

Conditioned Gaussian random field for a critical point
With the tools presented in the sections above, we are now ready to begin conditioning GRFs. We can begin applying (A38) and specializing it for critical points. We denote by A the event describing the field V (φ) taking a particular configuration, while B imposes V (0) ≡ V 0 = u and V i (0) ≡ η i = 0, that is, a critical point lying in the center of field space at height u. In order to proceed more easily, we shall discretize V (φ) as In this case, the conditioning event involves the Gaussian random vector field V = ∇V , whose Jacobian is the Hessian of the original field V evaluated at φ = 0. Therefore, its determinant is simply the product of the eigenvalues of the Hessian evaluated at the origin, Applying the Kac-Rice formula (A37) into (A38) yields where the integration domain will depend on the kind of critical point we are working with.
∆(λ) ∝ i<j |λ i −λ j | is the Jacobian of the variable change from components of the Hessian matrix to its eigenvalues, the proportionality constant depending on the dimensionality of the field space. For simplicity, the denominator in (A39) has been considered as a normalization factor for the distribution in the numerator.
We can rewrite (A39) in a more useful way: We can now see the power of this method. Assuming we have discretized our field space, we can readily compute the conditional probability distributions in (A41) and (A42) using the results from section A 2. This leads, together with (A42), to a distribution from which we can draw eigenvalues for a minimum of height u. These can be plugged in (A41) to generate iterations of GRFs with a minimum (or any other critical point) at their origin.
In order to apply all this machinery, let us introduce the following Gaussian random vector: where we denote by φ q the position in field space of a discrete set of points whose center is located at 0, V i (0) describes the first derivative along φ i and V ij (0) is the (i, j)-th element of the Hessian matrix. In order to unclutter the notation, we will compactify the previous vector as which has dimension q +1+n+n+ 1 2 n(n−1). The mean of (A43) is zero, and the covariance matrix of these quantities can be computed from the results in section A 4: In order to simplify the notation, since the jointly Gaussian probability distribution in the end depends on two-point functions, we can actually write 19 (A45) in the following way: With these arrangements, the Gaussian random vector corresponding to (A53) is 19 We basically have evaluated the first row for a given φ 1 and the first column for a given φ 2 , just as in [32]. Doing so allows us to treat the independent variable as a continuous one, rather than a discrete one.
obtain the mean function and covariance matrix of the conditioned process 20 . Using the results given above, one gets that the expectation value for the GRF around a critical point where V 0 = u and V 0 = 0, is given by, where h = h 11 , . . . , h nn , h 12 , . . . , h (n−1)n represents a certain configuration of the Hessian components of the field around the origin.
Furthermore, the covariance function for the conditioned GRF is now We can also obtain (A42), the distribution of eigenvalues at a critical point of a given height u, following the same steps as above, using as initial covariance matrix the bottomright block of (A53).

a. Analysis of a conditioned 2D Gaussian field
Let us apply these expressions to a two-dimensional isotropic and homogeneous GRF with covariance function and zero mean. For this case, we obtain the conditioned mean from (A57), which gives 20 Strictly speaking, we should be getting the mean and covariance of the random vector {V (φ 1 ), V (φ 2 )}.
Due to the isotropy of the GRF, φ 1 and φ 2 can be any points in field space. Thus, in order to unclutter the notation, we will only keep track of a single component of the resulting mean vector. Likewise, we will only keep the V (φ 1 )V (φ 2 ) component of the covariance matrix.
where h 21 = h 12 , by definition. Since we are free to choose the basis of φ, in order to simplify the expression we will employ the eigenvector basis of the Hessian matrix, therefore transforming (A60) toμ where λ i denote the two eigenvectors, drawn from (A42) specialized to this case (see below).
As for the conditioned covariance, from (A58) we obtaiñ Note that the covariance function of the conditioned process is not homogeneous anymore! This, however, makes complete sense. We have actually made the center of every realization special, meaning that homogeneity is broken in this sense. In fact, the new covariance is isotropic with respect to φ = 0, further stating that the center of the GRF is somehow different from the rest of the points.
All the presented machinery works not only for minima, but also for maxima and saddle points as well; the only difference among these being the sign of each λ i .

b. Distribution of heights and eigenvalues of the Hessian at a critical point
In order to calculate the probability distribution of the eigenvalues of the Hessian at a certain height of the potential at critical points we should pay attention to two ingredients.
The first one is the fact that the height and the second derivatives are correlated, so we need to calculate the multivariate covariance function for these quantities together. Furthermore, we also want to calculate this at critical points which can be done with the use of the generalized Kac-Rice formula.
Assuming a critical point located at φ = 0, the probability distribution to be computed is We can easily compute the PDF by conditioning the following random vector: of mean zero and covariance matrix    and plugging them into (A42), we get where N is a normalization factor and, in this two-dimensional example, ∆(λ) = |λ 1 − λ 2 | · π/2.
Setting V 0 to a constant value, say V 0 = u, in (A67) yields the distribution q u (λ 1 , λ 2 ), defined in (A42). On the other hand, integrating out either V 0 or the eigenvalues, gives the marginal distribution for the remaining variables in critical points (see Appendix B for more detail).
Another interesting application of (A66) is that it can be used to count the expected number of critical points in a certain region of field space. For example, to compute the expected number of minima per correlation volume Λ 2 in the example above, a direct application of (A37) yields In this case, the eigenvalues have been assumed to be positive. Setting other integration limits can give the expected number of maxima and saddle points, for example.

Conditioned Gaussian random field for an inflection point
We shall define an inflection point on our GRF as a point where the gradient of the field points in the direction of a Hessian eigenvector whose corresponding eigenvalue is zero.
Furthermore, we will also demand that the non-zero eigenvalue of the Hessian to be positive at this point.
In order to do this we can expand the discussion of the previous section by taking into account the third derivatives of the GRFs along with the lower ones. In order to simplify this description we will give a detail account of this construction for a 2d GRF only. Extending this to higher dimensions is straightforward. In particular we will be interested in the Gaussian random vector whose components have zero mean. As for the covariance matrix, it can be expressed as where (for the 2D case) Following the same steps as in the critical point case, we can obtain (for the covariance function (A59)) the expression for a GRF once we conditioned everything up to the third derivative. In order to do this we can first compute the mean value of the GRF in the vicinity of our inflection point, which is given bỹ where the basis of φ has been chosen to be the eigenbasis of the Hessian matrix (whose components are described by h and its eigenvalues by λ i ) and we have denoted by η and ρ the components of the first and third derivatives at the origin along the eigenbasis.
The conditioned covariance, on the other hand, reads which, once again, is isotropic around the origin of the field.

a. Probability distribution for the inflection point parameters
We can extend the treatment for the eigenvalues of the hessian that we did for the critical points to inflection points. The difference is that we will now impose that one of the eigenvalues vanishes while the other one is positive. Furthermore we will also impose that the gradient in the second eigenvalue direction also vanishes. These conditions have to be included in the calculation of the PDF of the parameters of the inflection points (V 0 , η 1 , λ 2 , ρ). Using a generalized version of the Kac-Rice procedure we arrive to, P inf dV 0 dλ 2 dη 1 dρ = N |λ 2 | 2 |ρ 111 | P V 0 , λ 2 | λ 1 = 0 P (η 1 , ρ ijk | η 2 = 0) (A77) where In (A77), one of the |λ 2 | factors comes from the Jacobian of the variable change to the eigenbasis of the Hessian (though with λ 1 = 0); the remaining |λ 2 ||ρ 111 | factor is just the determinant appearing in Kac-Rice's expression.
These last expressions can be used as in (A68) to compute the expected number of inflection point per correlation volume Λ 2 , which yields, for our choice of covariance function, .
In order to generate realizations of two-dimensional Gaussian random fields, we resorted to the so-called spectral or Karhunen-Love decomposition, due to its mathematical and computational simplicity.
Given a certain mean function µ(t), covariance function C(t, s) and a discretized space {t a } (where a runs over all n points in the lattice space) of a GRF, we can build the matrix C ab = C(t a , t b ), which by construction is symmetric and positive definite; therefore, we can always decompose C ab as where Λ = diag(λ 1 , . . . , λ n ) is the diagonal eigenvalue matrix, consisting of non-negative entries, and U is constructed by inserting all eigenvectors along its rows. Since Λ > 0, we can further decompose C as This procedure is tantamount to performing a Cholesky decomposition [53] on C; which is by far the most expensive step in this algorithm, in terms of computational cost.
Once we have computed L, constructing the GRF on the discretized space is straightforward. We only need to construct a random vector ξ of length n whose entries are independently distributed as Gaussian variables of zero mean and unit variance, and introduce the following variables: where µ a = µ(t a ). It can be easily shown that this gives the correct correlations among the values of the GRF evaluated at different points t a , (V a − µ a )(V b − µ b ) = L ac ξ c L bd ξ d = L ac L bd ξ c ξ d = L ac L bd δ cd = L ac L bc = L ac L T cb = (LL T ) ab = C ab = C(t a , t b ). (B4) The main advantage of using this procedure to generate GRFs is that the main computationally costly step, constructing the L matrix, needs to be performed only once. The rest of the algorithm is highly trivial from this perspective and allows for further simplification, as we have seen.

Numerical evaluations of Critical points
Using the expressions above we can compute the normalized distribution of heights of minima, maxima and saddle points for a 2d GRF, Furthermore, we can also compute the marginal distribution for the Hessian eigenvalues at critical points regardless of their height. This distribution is given by, We have checked the distributions above with numerical realizations of unconstrained Gaussian random fields in Mathematica. Regarding the heights of critical points, the numerical results fit the analytical prediction perfectly, as shown in figure 13(a).
As for the eigenvalue distribution, Fig. 13(b)-(d) shows that the histograms fit the analytical predictions perfectly once again. An important feature of these distributions is the fact that critical points with one of the eigenvalues close to zero or both eigenvalues close to each other are very rare; this effect (referred to as eigenvalue repulsion) is a direct consequence of the presence of the Vandermonde determinant in the distributions, as well as the Jacobian of the gradient field in the Kac-Rice formula.

Numerical evaluations of Inflection points
Using the results given above, we can obtain the following distributions for the parameters of the inflection points in a typical GRF.
where the complementary error function is defined as Once again, we found these expressions to be fully consistent with the numerical results, as shown in figure 14.
In order to find inflection points in our numerically generated potentials, we looked for

roots of the system
where H is the Hessian matrix, η T = (η 1 , η 2 ) represents the gradient at any point of the field and η T ⊥ = (−η 2 , η 1 ). It can be easily shown that simultaneous roots of Eq. (B11) are either critical or inflection points.
Finding inflection points numerically is quite tricky and the algorithm sometimes incorporates spurious points that, upon further study, are proven to be fictitious inflection points.
In order to make a proper comparison to the general expressions we have found analytically and avoid the inclusion of those spurious inflection points, we only considered those points which satisfied λ 2 /η 1 > 4 . This cut removes around 30% of the potential inflection points. Note that even though we might be removing a portion of real inflection points, the