Introduction

Within the field of Partial Differential Equations (PDEs), a large body of literature is concerned with investigating well-posedness of PDEs. That is, with studying the existence and uniqueness of solutions, in appropriate spaces. It is however the case that many nonlinear PDEs of pivotal importance in applications admit more than one solution. This paper is devoted to the problem of determining the number of solutions of stationary nonlinear PDEs (and then possibly finding properties of such solutions), i.e. solving problems of the form

where ℒ is a suitable (nonlinear) operator acting on a given function space. Studying the set of solutions of

royalsocietypublishing.org/journal/rsta Phil. Trans. R. Soc. A 383: 20240239 ......................................................................................................................... (1.1) is central in at least three broad contexts. First, solutions of (1.1) are candidate equilibria of the associated time-dependent PDE, 𝜕 t u = ℒu . (1.2)

This is relevant for both deterministic and stochastic problems. Indeed, when ℒ is the Fokker-Planck operator associated with a given Stochastic (ordinary) Differential Equation (SDE), the solutions of (1.1) are the invariant measures of the SDE. Second, at least when the solutions of (1.1) are stable-in the sense that they are stable equilibria of (1.2)-one expects them to correspond to metastable states of the associated Stochastic Partial Differential Equation (SPDE), i.e. of the evolution

where 𝜖 > 0 is a small parameter and 𝜂 is spatio-temporal noise. Finally, if the problem has a variational structure, solutions of (1.1) are related to extrema of associated functionals, often referred to as energy functionals or potentials. Hence the connection with the many applications of interest in calculus of variations and with energy landscape exploration problems. Many problems in both natural and applied sciences can be recast as one of (1.1), (1.2) and (1.3). Let us give some examples. Stochastic processes with multiple invariant measures naturally arise in the study of interacting particle systems and associated mean field limits. Such systems are paradigmatic models in statistical mechanics and kinetic theory, in connection with the study of phase transitions, as well as in the study of collective navigation and consensus formation [1,2]. The importance of understanding processes with multiple invariant measures cannot be understated. Indeed, while statistical sampling has been one of the main motivations for the study of processes with one invariant measure-ergodic processes-many processes in natural and engineered systems are likely to exhibit multiple invariant measures. This includes the tendency of fireflies to synchronize their flashing or not, the microstructure of nematic crystals aligning to several distinct equilibrium configurations (giving rise to different material properties), the opinion formation process in social media, which can converge towards various possible outcomes, etc. [3][4][5][6][7][8].

Related to problems of the form (1.3), metastability theory and landscape exploration are central to molecular dynamics, computational chemistry and condensed matter physics, particularly for the calculation of reaction rates. Indeed, a large number of methods have been proposed to understand transitions between metastable states, including the string method in finite [9] and infinite dimensions [10], the nudged elastic band method [11], surface walking approaches such as the dimer method, the gentlest ascent dynamics [12], the activation-relaxation technique, and more. All such methods require some prior knowledge of the metastable states (so that one can e.g. 'initialize the string'), which are a subset of the solutions of stationary problems like (1.1).

Despite the pervasiveness of PDEs with multiple solutions, the methods available to find all such solutions are very limited. To the best of our knowledge, aside from the brute force approach of starting the Newton algorithm from a large number of different initial points, there are few systematic strategies that work in full generality. The available tools include deflation [13] (which only requires pre-knowledge of one solution and can find any stationary solution, irrespective of its stability), continuation and homotopy-continuation [14] (which have similar requirements). These approaches are all numerical or computational as, except for rare cases, one does not expect to be able to find solutions analytically 1 . We mention also so-called 'landscape exploration methods' such as the eigenvector-following method [16] and minimax approaches [17] to find saddle points, which are also relevant to the study of reaction rates. However, these approaches mostly apply to finite dimensional landscapes.

We emphasize that the complexity of the problems that arise in this context is high; examples are illustrated in the works [18,19] concerning the study of nematic crystals (different solutions of the associated stationary PDE problem here correspond to distinct equilibrium configurations of the crystal), in the setting of de Gennes theory. In [19], 28 solutions were found for the problem considered there, with no guarantee that all had been found. For semi-linear Schrödinger equations, there are (at least) countably many soliton type solutions, even after accounting for symmetries [20,21].

There are two challenges to determining the set of solutions to (1.1): determining the number of solutions and then finding the solutions themselves (or at least some properties of such solutions). These challenges should be regarded as distinct. Deflation is a very powerful method, proven useful in many settings. But even by deflation one cannot guarantee all solutions have been found. This limitation is true of all the other methods we have mentioned, and it will certainly be true, also, of the method we propose in this paper. Thus, it is pragmatic to have at one's disposal an array of methods and to exploit synergies between them, to make progress on the specific problem at hand. For example, if by deflation one has found n solutions but some 'solution counting method' reveals that we should expect at least say n + 2 then one can redouble the search. Knowledge of such solutions can, in turn, provide the basis to use paradigms such as e.g. the string method, which requires pre-knowledge of the metastable states (to initialize the string).

This paper attempts to put forward an approach to address the problem of counting the number of (stable) solutions of a given stationary PDE of the form (1.1). The purpose of this short paper is to present the main idea of this approach though it is undoubtedly the case that further investigation and (non-trivial) analysis will be needed in order to understand its strengths and limitations. We comment on this in §1a below, §2 and at various points throughout.

The paper is organized as follows: in §1a we present our method and then, in subsequent sections, we apply this method to find the number of stable stationary solutions of the McKean-Vlasov (or porous media) equation. Such stationary solutions correspond to invariant measures of stochastic processes that are nonlinear in the sense of McKean. We do this both for parameter regimes where the number of solutions is, a priori, known and also in settings where the issue of finding the number of solutions is yet to be settled. In §2, we provide background material on the McKean-Vlasov PDE and on the associated McKean-Vlasov SPDE; the latter will be the central tool we employ to find the number of stable solutions of the stationary McKean-Vlasov PDE. Section 2 also serves the purpose of illustrating how determining the number of stationary solutions of a given PDE can be rather non-trivial, even in simple situations. Sections 3 and 4 are devoted to numerical and simulation aspects, respectively (more detail below).

(a) Our approach

We propose that in order to count the number of stable solutions of the stationary problem (1.1) one can sample from the invariant measure (assuming it exists and is unique) of the SPDE obtained by perturbing the time-dependent problem (1.2) by additive noise; that is, one can sample from the invariant measure of the SPDE

where 𝜂 is a suitable space-time noise. At least in the case in which the solutions of (1.1) are isolated, we expect the number of modes of the invariant measure of (1.4) to coincide with the number of stable solutions of (1.1) and the number of stable equilibria of (1.2). Further, we expect a correspondence between the metastable states of (1.4), the modes, and the stable stationary solutions of (

1.1) and (1.2). The noise 𝜂 must be appropriately chosen not only to make sure that (1.4) is well posed and has a unique invariant measure, but also to, potentially, satisfy other constraints royalsocietypublishing.org/journal/rsta Phil. Trans. R. Soc. A 383: 20240239 ......................................................................................................................... (a) (b) Figure 1. (a and b) A simple scalar multi-well potential, U, and the associated trajectory of (1.5). As expected, the trajectory tends to persist near one of the two minima. specific to the problem under consideration, see §2 for more detail on this pointfoot_4 . Intuitively, if the noise 𝜂 is strong enough we would expect (1.4) to have at most one invariant measure.

Let us explain the rationale of the approach in a simpler, finite dimensional setting. Consider one-dimensional Langevin dynamics in multi-well potential:

where Y t ∈ ℝ, U ∶ ℝ → ℝ is a confining potential, 𝛼 > 0 and 𝛽 t is one dimensional standard Brownian motion. Under very general assumptions on U, as soon as 𝛼 > 0, the SDE (1.5) has only one invariant measure, the Maxwellian e -2U∕𝛼 (after choosing U to include the correct normalization constant, which can always be done). By contrast, when 𝛼 = 0 this equation is a simple ODE with as many stable stationary states as the number of minima (wells) of the potential U (assuming such minima are isolated). Dirac measures centred around such stationary states can be viewed as the (multiple) stable invariant measures of the ODE. For a double well potential, one has three stationary states, two stable and one unstable. When U has many wells, a naive way to count them is to simulate (1.5) for long enough and count the number of modes of the measure e -2U∕𝛼 , as this measure will concentrate around the minima of U. This is shown in figure 1. This approach comes with limitations. First and foremost, one never knows if the whole state space has been adequately sampled; hence there are no guarantees that all the minima of U have been found. This problem is unavoidable, whether in finite or infinite dimensions, and irrespective of the method used to sample. The potential landscape might also be more complex with non-isolated minima, elusive saddle points, etc. All these problems persist in the infinite dimensional setting. This is why we propose this method to be used in conjunction with other methods, as an aid for exploration. We also emphasize, again, that the method we propose will only help count the number of stable solutions of (1.1), in the sense that they are stable in (1.2). It will be far less informative about unstable solutions.

Analogously with the finite dimensional case, one might see the idea of using the SPDE (1.4) as a landscape exploration method, especially if the deterministic part of the equation has some gradient structure. Some comments on this-and words of caution-can be found in Note 2.1.

While infinite dimensional Markov Chain Monte Carlo (MCMC) sampling algorithms are now well-developed for many problems of interest, one might not know the invariant measure of (1.4) explicitly, making the use of MCMC algorithms impractical. The option-which we take in the example of this paper-is to simulate the SPDE for a sufficiently long time. Despite all the caveats that we have highlighted, this approach is simple to understand and general, in the sense that, since it relies on a concept that does not hinge on the specific structure of the equation at hand, it can be also applied to equations different from our test case. More comments on this in §2. In the next section, we provide an illustration of its effectiveness in the context of the search of stationary solutions of porous media equations.

Background on the McKean-Vlasov PDE and associated SPDE

As a test case, we consider the problem of finding the number of stationary solutions of the McKean-Vlasov PDE,

for the unknown 𝜌 = 𝜌(t, x) ∶ ℝ + × 𝕋 → ℝ, where 𝕋 denotes the one-dimensional torus given by the interval [0, 2𝜋] with periodic boundary conditions 3 . That is, we want to find the stable solutions of the problem (1.1) when the operator ℒ is the nonlinear differential operator given by the righthand side of (2.1), which we henceforth refer to as the McKean-Vlasov differential V, F ∶ 𝕋 → ℝ are given functions, commonly referred to as the environmental and inter-particle potential, respectively, ′ denotes derivative with respect to the argument of the function and * denotes convolution. Let us briefly recap what is known about the evolution (2.1).

The PDE (2.1) is a popular and flexible model in a range of applications and has been studied for decades. It models the evolution of the density of a system of interacting particles. More precisely it is associated to the following interacting particle system

where, for every i = 1 … N, X i,N t ∈ 𝕋 is the position of the i-th particle at time t ≥ 0. Using mean field techniques [2] one can show that the empirical measure

of the system converges, as N tends to infinity, to a deterministic probability measure. Such a measure has a density, which is precisely the solution 𝜌 of the PDE (2.1). The PDE (2.1) can also be viewed as the Fokker-Planck equation of the following nonlinear SDE

In the above 𝜌 t is the density of the law of X t -so that the SDE is nonlinear in the sense of McKean-and, by Itô's formula, one can see that 𝜌 t solves (2.1). Hence the stationary solutions of (2.1) are the invariant measures of (2.3). A key fact to observe is that the evolution equation (2.1) preserves mass and positivity; that is, ∫ 𝕋 𝜌 t (x)dx is a constant, independent of time, and 𝜌 0 ≥ 0 implies 𝜌 t ≥ 0. This is coherent with the probabilistic interpretation of the PDE. The problem of determining the set of stationary solutions of (2.1) is far from trivial, even when the spatial variable is one-dimensional. Following the seminal paper [1], this problem has been tackled in a variety of settings [3,[22][23][24][25]. In particular, the number of stationary solutions of (2.1) depends on the choice of potentials V and F and on the strength of the noise. Here, to fix ideas, we consider the case where V and F are given by

(2.4)

The relevant properties of this choice are that V is non-convex (more specifically, double well and symmetric) and F introduces an attractive force between particles. One can show that with this choice of V and F, there exists a critical value of the noise, 𝜎 c > 0, such that if 𝜎 > 𝜎 c equation (

2.1) admits exactly one stationary solution, 𝜌 0 , which is bimodal and equally concentrated around the bottom of two wells. The two wells have minima at x 1 = 𝜋∕2 and x 2 = 3𝜋∕2. For 1∕2 < 𝜎 < 𝜎 c and royalsocietypublishing.org/journal/rsta Phil. Trans. R. Soc. A 383: 20240239 .........................................................................................................................

for 𝜎 small enough there are exactly three stationary solutions of (2.1)-the conjecture being that for any 0 < 𝜎 ≤ 𝜎 c there should be three stationary solutions. One is 𝜌 0 itself; the other two, 𝜌 ± , are still bimodal with modes around x 1 and x 2 but 𝜌 + concentrates most of the mass in x 1 while 𝜌 -concentrates most of its mass around x 2 . The reader may find intuition for this property in, among other references, [26]. Here we briefly recall (as it will be relevant to §4 and for later comments) that this result is obtained by first observing that, irrespective of the choice of V and F, any stationary solution 𝜌 ∞ of (2.1) solves the following fixed point problem

.

(2.5)

Once V and F are chosen as in (2.4), the above, a priori, infinite dimensional fixed point problem, can be reduced to a finite-dimensional, two-parameter fixed point problem. One can indeed observe (see [26, §2]) that in this case 𝜌 ∞ solves (2.5) if and only if it has the following structure

where m 1 , m 2 ∈ ℝ are fixed points of the following equations

respectively. To recap, when ℒ is the McKean-Vlasov differential operator, solving (1.1) corresponds to solving the fixed point problem (2.5). The results we have summarized above on the number of stationary states of the McKean-Vlasov PDE are obtained by using this strategy. More comments on this in Note 2.1. We note in passing that the fact that the PDE (2.1) can admit several stationary solutions is in contrast to the behaviour of the particle system (2.2), which, irrespective of the strength of the noise 𝜎 > 0, always admits just one invariant measure (although clearly this invariant measure is made of multiple metastable states). So the particle system is not necessarily helpful to our goals (see Note 2.1).

While (2.1) has many stationary states, it is natural to imagine that, after adding a strongenough noise to (2.1), we will obtain an SPDE with a unique invariant measure (this invariant measure is a measure on the space of measures). Technical details aside, this turns out to be the case. Indeed, consider the SPDE

where W t is cylindrical Wiener noise and Q is a trace class operator, so that

and the w k t 's are independent one-dimensional standard Brownian motions, {e k } k is a basis of L 2 (𝕋; ℝ) (in which Q is assumed to diagonalize) and {𝜆 2 k } k∈ℕ is the sequence of non-negative eigenvalues of Q. We will assume that the eigenvalues of Q are of the form

for some 𝛾 > 0 and s > 1∕2. Well-posedness and ergodicity of the SPDE (2.8) have been studied in [26]. Using the results of [26] one can easily see that with the above choice of Q, the SPDE (2.8)

, while it is strong Feller and irreducible, respectively, for

, respectively. This implies that for 1∕2 < s < 1 the SPDE admits at most one invariant measure. Whether such an invariant measure exists or not is still an open problem, but we believe it to be the case. We proceed by assuming (2.8) admits a unique invariant measure, at least under (2.10). Our simulations in §4 will, as a by-product, support this conjecture.

To our purposes, a more pressing issue is the fact that the SPDE (2.8) does not preserve mass in general and it does not preserve positivity either. This means that sampling from the invariant measure of (2.8) may also result in exploring parts of state space in which we are not interested, as we only want stationary solutions of (2.1) that have mass one and are non-negative. The choice 𝜆 0 = 0 in (2.10) has been made to remedy this issue and guarantee mass conservation-by formally integrating the equation in space and using the fact that ∫ e k (x) dx = 0 for every k ≥ 1 one can see that, if 𝜆 0 = 0, then the quantity ∫ u t (x) dx is constant in time. As for the positivity, a straightforward calculation shows that stationary solutions of (2.1) with mass one are always non-negative 4 . Therefore, once condition (2.10) is enforced, we expect the solution of (2.8) to concentrate only around the desired (stable) stationary solutions of (2.1), i.e. those that have mass one and are hence positive.

To the best of our knowledge, when 0 < 𝜎 < 𝜎 c , there are no results on the global stability of the three stationary states of (2.1) (with V and F as in (2.4))-though some partial results in this direction have been obtained in [27] and later in [28]. When 0 < 𝜎 < 𝜎 c one expects that two of the equilibria (𝜌 ± ) will be stable and one unstable (𝜌 0 ). In appendix A, we numerically show that this is indeed the case. Hence, when sampling from the invariant measure of (2.8), we expect to see bimodal behaviour in our measure. This is confirmed by the simulations of §4; more detail is provided in that section. In general the linearly unstable solutions may be akin to saddle points; they will be challenging to see directly in the simulation of (2.8), but may be critical to understanding transition rates among metastable solutions.

An explicit form for the invariant measure of (2.8) is not known. This precludes the use of many MCMC sampling strategies as they often depend on evaluating the (unnormalized) invariant density in a Metropolis ratio. Moreover, since we are looking to sample from a measure supported on an infinite dimensional space, it would be advisable to use an MCMC method that is well posed in infinite dimension [29][30][31]. But, again, we do not know of any MCMC method that is well posed in infinite dimension and that can be implemented in absence of information on the analytical form of the target measure [29,32,33].

We therefore take the approach of simulating the SPDE (2.8) for sufficiently long. As the simulation of this SPDE needs care, we present in §3 a numerical scheme that is used for its study. We reserve for future work the proof of the convergence of this numerical scheme. Note 2.1. Some further comments on the SPDE (2.8) and its simulation.

-The SPDE (2.8) can be obtained as a limit of an appropriate (weighted) interacting particle system. This has been proved in [34]. Such a particle system could be used as a way of simulating the SPDE. However, similarly to what happens for the particle system (2.2) and its limiting PDE (2.1), it may be the case that the ergodic properties of the SPDE (2.8) are different from those of the associated particle system. Indeed, the work [34] proves that the particle system converges to the SPDE over finite time horizons, but there is a priori no reason to believe that convergence will hold uniformly in time, i.e. over infinite time horizons. This is why we opted to simulate the SPDE with a numerical discretization scheme. -As we mentioned at the end of §2a, the deterministic part of the evolution (2.8) is in gradient form (see [3]). So one would expect, in analogy to the finite dimensional setting discussed around equation (1.5), that the approach we take here could be interpreted as a landscape exploration method. Heuristically this is perhaps the case. However, there is a subtlety to this point. As we mentioned, there are no rigorous results on the existence

royalsocietypublishing.org/journal/rsta Phil. Trans. R. Soc. A 383: 20240239 .........................................................................................................................

of the invariant measure of (2.8). In particular, we do not know if such a measure is of Maxwellian type, and this may be challenging, in view of recent observations [35]. Interpreting the method we consider here as an infinite dimensional landscape exploration method hinges on making progress on this aspect of the theory. Nonetheless, whether the invariant measure is of Maxwellian type or not, one still expects that such a measure will concentrate around the stable stationary solutions of the McKean-Vlasov PDE (as confirmed by our simulations), and this is all we truly need. -Still, on the matter of potentially using particle systems to solve the problem of interest, it should be noted that the particle system (2.2) has, for each fixed N, a unique invariant measure. Such an invariant measure is metastable and the number of metastable states of this measure coincides with the number of stable stationary states of (2.1) [36]-and it is expected to be the same as the number of metastable states of the invariant measure of the SPDE (2.8). However the exit time from a potential well for the particle system (2.2) grows exponentially with N [36,37], making this approach computationally impractical (if not infeasible), especially in the presence of many wells. By contrast, using the metastability of the SPDE (2.8) is computationally much faster. -Reducing the infinite dimensional fixed point problem (2.5) to a finite dimensional one is only practical when V and F are sufficiently simple. If V and F allow for this drastic reduction in dimensionality (with our choices of potentials this is due to the fact that such potentials are each given by a single Fourier mode), then it is usually possible to solve, numerically, the resulting finite dimensional fixed point problem. For us, this is problem (2.7). However, even in these simple cases, the analytical study of the finite dimensional fixed point problem is challenging, and it tends to offer only partial results (as for the case we have detailed above). See [26], references therein and the list of references at the start of this section.

Even slightly complicating potentials V and F makes the reduction of the infinite dimensional fixed point problem to a finite dimensional one substantially more difficult. Should it succeed, the numerical and analytical study, respectively, of the resulting finite dimensional problem becomes tricky and prohibitive to pursue. In §4, we choose V and F as in (2.4); this provides a good test case, as the number of stable stationary states is known in this case (modulo the comments we made above). However, applying our method is not needed in this relatively simple case. We expect our method will be more useful in more complex situations-including complex choices of V and F or, more so, when the PDE part of the evolution is not of McKean-Vlasov type. Indeed, if the PDE part of the evolution is not of McKean-Vlasov type then it is not necessarily the case that the problem (1.1) can be reduced to some form of fixed point equation. In such instances, the method we have presented allows for greater flexibility. An interesting avenue could be the one of using the method we present here for different nonlinear PDEs (using the associated SPDEs), such as Allen-Cahn type equations, Navier-stokes and equations appearing in nonlinear elasticity, such as those in [19]. This will however require new theoretical efforts in the ergodic theory for SPDEs.

Numerical methods

Here, we review our computational methods for exploring (2.8).

(a) Discretization of the SPDE

To generate sample paths of (2.8), we employ a spectral Galerkin spatial discretization together with Euler-Maruyama time stepping. We refer the reader to, for instance, [38]. Denoting the total number of resolved Fourier modes J, and using complex trigonometric functions, we will achieve spectral accuracy with O(J ln J) scaling through the Fast Fourier Transform (FFT). Euler-Maruyama time stepping for this problem (having only additive noise) will yield a strong error of O(∆t) in the sample paths. This strategy is entirely standard, and though we do not provide a full numerical analysis of it in the present work, our numerical tests are consistent with what would be expected from this method. We defer such numerical tests to §4b. We expect that a complete analysis of the numerical scheme we present here should be accessible and we hope to pursue this in future work.

While our example computations will be for the case that F and V are given by (2.4), these methods are expected to work for all sufficiently regular V and F.

(i) Space discretization: spectral Galerkin method

The fundamental assumption is that the covariance, Q, of the noise process in (2.8) diagonalizes against the complex trigonometric functions on the 𝕋 = [0, 2𝜋)

This allows us to express the noise as (2.9). We now represent the solution, u, of (2.8) by

Substituting in and projecting onto e k , (2.8) becomes

Symbols ℱ and ℱ -1 denote the Fourier transform and its inverse; in software, these are carried out by the FFT algorithm. The above notation is to emphasize that V ′ ⋅ u is computed in real space, while F ′ * u is first computed in Fourier space, and then (F ′ * u) ⋅ u is computed in real space. Using an integrating factor on the linear term, we can write the mild solution (in Fourier coordinates) as ûk (t) = e -𝜎k 2 t ûk (0)

Thus far, we have only formally manipulated (2.8); no approximation has been made. We now truncate to J (assumed to be even) modes, so that

and let P (J) be the projection onto the J lowest modes. Substituting (3.5) into (3.4) and applying this projection we obtain our finite dimensional system:

As is sufficient for problems with a quadratic nonlinearity, the projection operator can be implemented in code by:

-Simulating with 2J modes, if we wish to resolve J of them; -Zeroing out the J highest frequency modes.

This corresponds to anti-aliasing, [39].

Downloaded from https://royalsocietypublishing.org/ on 18 October 2025 royalsocietypublishing.org/journal/rsta Phil. Trans. R. Soc. A 383: 20240239 ......................................................................................................................... x. On the vertical axis of each graph there is time t while on the horizontal axis lies the spatial variable x. The colour bar is common across the three plots. These plots have been obtained by employing CairoMakie [40].

(ii) Time discretization: Euler-Maruyama

The fully discretized problem is obtained via an Euler-Maruyama temporal discretization of (3.6). Hence, on top of the spatial approximation, we make the temporal approximation

with t n = n∆t. After a bit of stochastic calculus on the noise term, the update is

2 ∆t 2𝜎k 2 𝜉 k,n+1 , 𝜉 k,n+1 i.i.d. ∼ N(0, 1). (3.8) 4. Simulations (a) Computed stationary solutions We now use the numerical discretization method presented in §3 to simulate the dynamics (2.8) for long times. When the stationary solutions of (2.1) are isolated, we expect the process {u t } t≥0 solution of (2.8) to exhibit metastable behaviour, hopping between different stable stationary solution of (2.1). To observe this fact numerically, we show a heat plot of the solution u t , see figure 2; moreover, inspired by (2.7), we consider the quantities royalsocietypublishing.org/journal/rsta Phil. Trans. R. Soc. A 383: 20240239 ......................................................................................................................... u t (x) sin x dx for the solution u t to (2.8) with 𝜎 = 0.2, 0.6 and 1. These simulations are run with the same parameters as in figure 2. When V, F are as in (2.4) and 𝜎 is either small enough or 1∕2 < 𝜎 < 𝜎 c , we know that (2.1) has two stable stationary states, 𝜌 ± . Such stable states correspond to the values of (m 1 , m 2 ) given by the points (±1, 0) (for detail see [26]). We expect I 1 (t) to be metastable around the values of m 1 = ±1 and m 2 to fluctuate around the origin. This is consistent with figures 3 and 4, respectively. When 𝜎 > 𝜎 c , which is the case for the right-most plot in figure 2, we know that (2.1) has a unique stationary solution, which we expect to be bimodal and equally concentrated around the two minima. This is confirmed by the simulation. The value of 𝛾 = 10 -2 used in these experiments was found empirically. As we had some foreknowledge of the stationary solutions in the system, we found this value to be a sufficiently small value such that these solutions were metastable on the timescale of the simulation. It is possible that we may miss solutions should this noise level be too large.

Let us now make a further remark about figure 2. The left (𝜎 = 0.2) and right (𝜎 = 1) panels in that figure are consistent with our expectations. Since the value 𝜎 = 0.6 is also below the critical threshold, we would have expected to see metastability also in the panel in the middle. That is, we would have expected the central panel to look a bit more like the one on the left. Upon running the simulation for longer, the expected metastability appears also for the value 𝜎 = 0.6 (not shown in figures). Nonetheless it is at first sight counterintuitive that a higher value of 𝜎 should be slowing down the hopping between metastable states. We do not have a complete explanation of this phenomenon yet, and we intend to carry out further study of this. However, we point out that this result is not wholly unreasonable: in this infinite dimensional scenario the frequency of hopping between metastable states is governed by the parameter 𝛾, which controls the intensity of the infinite dimensional noise, not by 𝜎. Increasing 𝛾 does indeed result in faster hopping between metastable states (not shown in figures). We do believe that understanding why an increase of 𝜎 slows down metastability is related to the role that 𝜎 plays in the (weighted) interacting particle system associated to (2.8)-see [34].

Finally, in figure 5 we run the simulation for different potentials V and F. In particular the potential V we choose in figure 5 has four wells. For these potentials no rigorous proofs are available, though the result we find is in line with heuristics. u t (x) cos x dx for the solution u t to (2.8) with 𝜎 = 0.2, 0.6 and 1. These simulations are run with the same parameters as in figure 2. (b) Verification of algorithm 𝔼[‖u (J,∆t) (t max ) -u(t max )‖ 2 L 2 ] ≈ 1 n trials n trials ∑ i=1 ‖u (J,∆t,i) (t max ) -u ref (t max )‖ 2 L 2 (4.2)

as we vary J and ∆t independently, with n trials = 10 4 independent trials in each case. The L 2 norms are computed in Fourier space. For reference solutions, we used ∆t ref = 10 -2 and J ref = 2 12 to check J convergence; ∆t ref = 10 -6 and J ref = 2 7 to check ∆t convergence. The convergence rates of both the number of Fourier modes J the and time-step ∆t is consistent with analogous discretizations of semi-linear SPDE (i.e. stochastic Allen-Cahn) [38]. The spectral differentiation matrix, D x is dense Toeplitz [41]. F ′ corresponds to the discretization of F ′ * , and 𝝆 is the discrete approximation of 𝜌 ∞ . For the values of J we consider, a standard eigenvsolver suffices. As we see in figure 7-whose numerical simulations have been carried out for different values of 𝜎 with the choice (m 1 , m 2 ) = (0, 0), i.e. 𝜌 ∞ (x) = 1 Z 𝜎 e -1 𝜎 cos 2x -there is an eigenvalue of positive

real part for 𝜎 in the subcritical regime. This corresponds to a linear instability about this solution, and explains why we do not observe it in our time-dependent simulations. By contrast, for 𝜎 in the supercritical regime, the eigenvalue nearest the imaginary axis has negative real part, and it is linearly stable; we do see this solution in the time-dependent simulations. Though the eigenvalues appear stable to numerical refinement, it remains to establish these properties for (A 1). We note that we filtered out zero eigenvalues of L as these corresponded to eigenstates that violated the mean zero condition of (A 1).

For the non-trivial solution (i.e. (m 1 , m 2 ) ≠ (0, 0)), in the subcritical regime, figure 8 implies these are stable solutions. As the images show, the eigenvalue nearest the imaginary axis is negative. Thus, these will become metastable when stochastic forcing is added.