Mean Field Limits for Interacting Diffusions in a Two-Scale Potential

In this paper, we study the combined mean field and homogenization limits for a system of weakly interacting diffusions moving in a two-scale, locally periodic confining potential, of the form considered in Duncan et al. (Brownian motion in an N-scale periodic potential, arXiv:1605.05854, 2016b). We show that, although the mean field and homogenization limits commute for finite times, they do not, in general, commute in the long time limit. In particular, the bifurcation diagrams for the stationary states can be different depending on the order with which we take the two limits. Furthermore, we construct the bifurcation diagram for the stationary McKean–Vlasov equation in a two-scale potential, before passing to the homogenization limit, and we analyze the effect of the multiple local minima in the confining potential on the number and the stability of stationary solutions.

(1.1) denotes the position of the interacting agents, V (·) a confining potential, θ the strength of the interaction between the agents, {B i t } N i=1 standard independent one-dimensional Brownian motions and β the inverse temperature. The total energy (Hamiltonian) of the system of interacting diffusions (1.1) is (1.2) Passing rigorously to the mean field limit as N → ∞ using, for example, martingale techniques (Dawson 1983;Gärtner 1988;Oelschläger 1984), and under appropriate assumptions on the confining potential and on the initial conditions (propagation of chaos), is a well-studied problem. Formally, using the law of large numbers we deduce that where the expectation is taken with respect to the "1-particle" distribution function p(x, t). 1 Passing, formally, to the limit as N → ∞ in the stochastic differential equation (1.1), we obtain the McKean SDE (1. 3) The Fokker-Planck equation corresponding to this SDE is the McKean-Vlasov equation (Frank 2005;McKean 1966McKean , 1967) (1.4) The McKean-Vlasov equation is a nonlinear, nonlocal Fokker-Planck type equation that we will sometimes refer to as the McKean-Vlasov-Fokker-Planck equation. It is a gradient flow, with respect to the Wasserstein metric, for the free energy functional where we write the interaction potential as F(x) = 1 2 x 2 . Background material on the McKean-Vlasov equation can be found in, e.g., Carrillo et al. (2006), Frank (2005) and Villani (2003).
The finite-dimensional dynamics (1.1) has a unique invariant measure. Indeed, the process x t defined in (1.1) with V being a confining potential is always ergodic, and in fact reversible, with respect to the Gibbs measure (Pavliotis 2014, Ch. 4), On the other hand, the McKean dynamics (1.3) and the corresponding McKean-Vlasov-Fokker-Planck equation (1.4) can have more than one invariant measures, for nonconvex confining potentials and at sufficiently low temperatures (Dawson 1983;Tamura 1984). This is not surprising, since the McKean-Vlasov equation is a nonlinear, nonlocal PDE and the standard uniqueness of solutions for the linear (stationary) Fokker-Planck equation does not apply (Bogachev et al. 2015).
The density of the invariant measure(s) for the McKean dynamics (1.3) satisfies the stationary nonlinear Fokker-Planck equation Based on earlier work (Dawson 1983;Tamura 1984), it is by now well understood that the number of invariant measures, i.e., the number of solutions to (1.7), is related to the number of metastable states (local minima) of the confining potential-see Tugaut (2014)  This one-parameter family of probability densities is subject, of course, to the constraint that it provides us with the correct formula for the first moment: (1.9) We will refer to this as the self-consistency equation and it will be the main object of study of this paper. Once a solution to (1.9) has been obtained, substitution back into (1.8) yields a formula for the invariant density p ∞ (x; θ, β, m). Clearly, the number of invariant measures of the McKean-Vlasov dynamics is determined by the number of solutions to the self-consistency equation (1.9). It is well known and not difficult to prove that for symmetric nonconvex confining potentials a unique invariant measure exists at sufficiently high temperatures, whereas more than one invariant measure exists below a critical temperature β −1 c (Dawson 1983, Thm. 3.3.2;Tamura 1984, Thm. 4.1, Thm. 4.2); see also Shiino (1987). In particular, for symmetric potentials, m = 0 is always a solution to the self-consistency equation (1.9). Above β c , i.e., at sufficiently low temperatures, the zero solution loses stability and a new branch bifurcates from the m = 0 solution (Shiino 1987). This second-order phase transition is similar to the one familiar from the theory of magnetization and the study of the Ising model. In Fig. 1, we present the solution to the self-consistency equation and the bifurcation diagram for stationary solutions of the McKean-Vlasov equation for the standard bistable-Landau-potential V (x) = x 4 4 − x 2 2 . To compute the critical temperature, we need to solve the equation obtained by differentiating the self-consistency equation with respect to the order parameter m at m = 0 (see Shiino 1987;Frank 2005, Sec 5.1.3 for more details): (1.10) The number of times that m and R(m; θ, β) cross, i.e., the number of stationary measures, depends on the slope of R(m; θ, β) at the origin. This is given precisely by Eq. (1.10).
The main purpose of this paper is to study the dynamics and, in particular, bifurcations and phase transitions for a system of interacting diffusions moving in a rugged energy landscape, coupled through the Curie-Weiss interaction. We are particularly interested in understanding the combined effect of the presence of several local minima (metastable states) in the confining potential and of the passage to the mean field limit. We will study the problem for a system of interacting diffusions of the form (1.1) moving in a two-scale, locally periodic confining potential 2 . This class of potentials provides us with a natural testbed for testing several techniques and methodologies for the study of multiscale diffusions such as maximum likelihood estimation (Papavasiliou et al. 2009;Pavliotis and Stuart 2007), particle filters and filtering (Imkeller et al. 2013;Papavasiliou 2007), importance sampling and large deviations (Spiliopoulos 2013) and optimal control (Hartmann et al. 2014).
Of particular relevance to us is the multiscale analysis presented in Duncan et al. (2016a), Duncan et al. (2016b). 2 In these works, the homogenized SDE for a Brownian particle moving in a two-scale potential in R d valid in the limit of infinite scale separation → 0 was obtained and the effect of the multiscale structure on noise-induced transitions was investigated. It was shown, in particular, that the homogenized SDE is characterized by multiplicative noise. For a single Brownian particle in R d moving in a two-scale potential (1.11) (or, equivalently, for a system of d noninteracting Brownian particles in a two-scale potential), the homogenized equation reads where M(·) denotes the diffusion tensor and (·) the free energy-see Sect. 2. It is important to note that, in addition to the presence of multiplicative noise, the potential energy driving the dynamics is not simply the average of the two-scale potential over its period, but, rather, the free energy = −β −1 ln e −βV (x,y) dy . Since the dynamics (1.13) is finite-dimensional, no phase transitions can occur. In fact, the 2 In fact, in these papers a potential with N microscales and one macroscopic scale of the form where V is periodic in all the microvariables is studied. For the purposes of this work, it is sufficient to consider a potential with two characteristic, widely separated, length scales. homogenized dynamics is reversible with respect to the thermodynamically consistent Gibbs measure; see the discussion in Sect. 2. It is well known, however, that multiplicative noise can lead to noise-induced transitions, i.e., to changes in the topological structure of the invariant measure (Horsthemke and Lefever 1984;Pavliotis 2014, Sec. 5.4). Such phenomena, including multiscale-induced hysteresis effects, for a one-dimensional Brownian particle moving in a multiscale potential, were studied in detail in Duncan et al. (2016a). Our goal is to study mean field limits for multiscale interacting diffusions of the form (1.14) where the two-scale potential is given by (1.11). The interaction potential F(·) is assumed to be a smooth even function, with F(0) = 0 and F (0) = 0. All of the numerical experiments that we will present will be for the Curie-Weiss quadratic interaction potential F(x) = 1 2 x 2 . The main issues that we address in this work are: 1. What is the effect of the presence of (infinitely) many local minima in the locally periodic confining potential on the bifurcation diagram? In other words, how do the bifurcation diagrams for 1 but finite and → 0 differ? 2. Do the homogenization and mean field limits commute, in particular when also passing to the long time limit T → +∞? In other words: are the bifurcations diagrams corresponding to the N → ∞, T → ∞, → 0 and → 0, N → ∞, T → ∞ limits the same?
Two typical examples of the type of locally periodic potentials that we will study in this paper are shown in Fig. 2: It should be clear from these two figures that the homogenization and mean field limits, when also combined with the long time limit, do not necessarily commute. First, the homogenization process tends to smooth out local minima and to even "convexify" the confining potential-think of a quadratic potential perturbed by fast periodic fluctuations. This implies, in particular, that even though many additional stationary solutions, i.e., branches in the bifurcation diagram may appear for all finite values of , most, if not all, of them may not be present in the bifurcation diagram for the homogenized dynamics. Furthermore, multiplicative/nonseparable fluctuations of the type presented in Fig. 2b tend to flatten the potential around x = 0. As we will see in Sect. 4, this phenomenon is very much related to the lack of commutativity of the limits N → ∞, T → ∞, → 0 and → 0, N → ∞, T → ∞. We will study these problems using a combination of formal multiscale calculations, (some) rigorous analysis and extensive numerical simulations. There are many technical issues that we do not address, such as the rigorous homogenization study of the McKean-Vlasov equation and the rigorous study of bifurcations in the presence of infinitely many local minima. We will address these in future work.
The rest of the paper is organized as follows. In Sect. 2, we study the mean field limit for a system of homogenized interacting diffusions, i.e., the first → 0, then N → ∞ limit. In Sect. 3, we study the homogenization problem for the McKean-Vlasov equation in a two-scale potential. In Sect. 4, we present extensive numerical simulations. Section 5 is reserved for conclusions.

Mean Field Limit of the Homogenized Interacting Diffusions: First
→ 0, then N → ∞ In this section, we consider the one-dimensional version of the system of SDEs (1.14). We first take the homogenization limit ( → 0) and then the mean field limit (N → ∞). The homogenization theorem for a system of finite-dimensional interacting diffusions moving in a two-scale confining potential is presented in Duncan et al. (2016b). The mean field limit of the homogenized SDE system can be obtained by using the results of Gärtner (1988), Oelschläger (1984).

Homogenization for Finite System of Interacting Diffusions in a Two-Scale Potential
We consider the system of interacting diffusions where F is a smooth even function with F(0) = 0 and F (0) = 0 and V is a smooth locally periodic potential of the form (2.7). We introduce the notation x t = (X 1 t , . . . , X N t ), so that we have and B t is a standard Brownian motion in R N . This equation is of the same form as Duncan et al. 2016b, Eq. (1) and Duncan et al. (2016a), Since F does not depend on the fast scale, the results of Duncan et al. (2016b) apply directly to (2.2) and we deduce that the sequence x t converges, as → 0, to the solution of the homogenized equation where W N (x, y) is defined as in Eq. (1.2), The convergence is in the sense of weak convergence of probability measures, i.e., the law of the process x t converges weakly to the law of the limiting process x t . The proof of this result, which is quite standard, is based on the application of Itô's formula to the solution of an appropriate Poisson equation, Eq. (2.12), the decomposition of the rescaled process into a martingale part and a remainder part, and the use of the martingale central limit theorem. The details can be found in Duncan et al. (2016b). It will be useful to decompose the two-scale potential into its large-scale confining part and the modulated, mean-zero, fluctuations: Notice that this decomposition is not unique, since we can define the average of the two-scale potential over the unit cell with respect to a different, e.g., Gibbs, weight.
However, the choice of the weight does not affect our results. See, e.g., the proof of Proposition 3.1. We note that the free energy N is of the form and, for fixed x ∈ R d , is the unique weak solution in H 1 (Y) to ∇ y · e −βV 1 (x,y) (I + ∇ y (x, y)) = 0, y ∈ Y, (2.12) such that the centering condition Y (x, y)e −βV 1 (x,y) dy = 0, for all x ∈ R d is satisfied. The proof of uniqueness of centered solutions to this equation is based on the Lax-Milgram lemma and can be found in Duncan and Pavliotis (2016b, Thm. 2.3).
To compute the diffusion tensor (see Duncan et al. 2016a, Appendix A for a similar computation), we observe that (2.13) Since there is no coupling between the different y i components od Eq. (2.12), it follows that (x, y) can be written in the form ( Substituting in (2.13), we obtain (2.15) and the diffusion tensor is diagonal, with As it is well known (Pavliotis and Stuart 2008, Sec 13.6.1), the one-dimensional Poisson equation (2.14) can be solved explicitly, up to quadratures. We can then obtain formulas for the diagonal elements M ii of the diffusion tensor M(x): (2.16) We can write the system of stochastic differential equations for the homogenized system of interacting particles: (2.16), prime denotes differentiation with respect to x and is given by Equations (2.8)-(2.9).
We note that the homogenized system of SDEs (2.17) is characterized by multiplicative noise. 3 Furthermore, the diffusion coefficient of the ith particle depends only on the position of the particle itself, and not of the other particles. The dynamics (2.17) is reversible with respect to the Gibbs measure (2.18)

Mean Field Limit for the Homogenized SDE
We can now pass to the mean field limit N → ∞. The system of SDEs (2.17) is of the form which is in the same form to the one considered in Gärtner (1988), Oelschläger (1984), with slightly different drift and diffusion coefficients. 4 It is straightforward to check that the homogenized equation satisfies the conditions in the aforementioned papers. 5 Taking the mean field limit of (2.17), we obtain the following nonlinear Fokker-Planck equation: where denotes the convolution operator in x, and M(x) is defined in (2.10). We note that the solution of Eq. (2.19) represents the density of the empirical measure of the process in the limit N → ∞. The McKean stochastic differential equation corresponding to (2.19) is We reiterate that the correction to the drift β −1 M (X t ) dt is not the Stratonovich correction, but rather the Klimontovich (kinetic) one. This interpretation of the stochastic integral ensures that the homogenized dynamics is reversible with respect to the (thermodynamically consistent) Gibbs measure(s) that we can calculate by solving the stationary Fokker-Planck equation.
The (one or more) stationary distributions p ∞ (x; θ, β, m) are solutions to the stationary Fokker-Planck equation The detailed balance condition implies that and since M(x) is strictly positive, a simple variant of Tamura (1984, Lemma 4.1) enables us to obtain an integral equation for the invariant distribution: where ψ(x) is given by Eq. (2.20). In particular, p ∞ is independent of the diffusion tensor M(x). For the particular case of a quadratic interaction potential F(x) = x 2 2 , which is the case that we will study here, all stationary solutions are given by the one-parameter family of Gibbs states of the form (1.8) and the integral equation (2.23) reduces to a nonlinear equation, the self-consistency equation (Shiino 1987) By solving this equation, we can construct the full bifurcation diagram of the stationary Fokker-Planck equation. This will be done in Sect. 4. We are also interested in the equation for the critical temperature (1.10), which in this case is given by (2.25) Assuming that the large-scale part of the potential is symmetric, we have that x p ∞ (x; β, θ, m = 0) dx) = 0 and the equation above simplifies to From the definition of ψ(x) in Eq. (2.20)), we can conclude that for separable potentials, i.e., when V 1 (x, y) is independent of x, then ψ(x) becomes a constant. This, in turn, means that the stationary solutions to the homogenized McKean-Vlasov equation are the same to the ones for the system without fluctuations (V 1 (x, y) = 0)see Corollary 3.2 in Sect. 3. For example, when the large-scale part of the potential V 0 (x) is convex, there are no phase transitions for the homogenized dynamics. We will show in Sections 3 and 4 that this is not the case if we take the limits in different order.

Multiscale Analysis for the McKean-Vlasov Equation in a Two-scale Potential
In this section, we consider the homogenization problem for the McKean-Vlasov equation in a locally periodic potential for the case of a quadratic (Curie-Weiss) interaction. In particular, we first pass to the mean field limit (i.e., send N → ∞) in Eq. (1.14) with F(x) = x 2 2 and study the effects of finite (but small) on the bifurcation diagram, before sending → 0.

Mean Field Limit for Interacting Diffusions in a Two-Scale Potential: N → ∞, > 0 Finite
We start with the system of interacting diffusions (3.1) The notation is the same as in Sect. 2, i.e., V (x) := V x, x is a smooth confining potential that is L-periodic in its second argument, θ > 0 is the interaction strength, β the inverse temperature and {B i t , i = 1, . . . , N } are standard independent onedimensional Brownian motions.
Taking the limit as N → ∞, we obtain the McKean-Vlasov-Fokker-Planck equation: The equilibrium solutions, i.e., stationary states, of this equation are given by a oneparameter family of two-scale Gibbs distributions-see Eq. (1.8): Our goal now is to study the → 0 limit of the self-consistency equation-see and also the equation for the critical temperature, Proof The proof of this result follows from properties of periodic functions (Pavliotis and Stuart 2008, Thm. 2.28). Consider u ∈ L 2 (R; C per (Y)), > 0 and define (3.6) We will use this fact to identify the limits as → 0 of p ∞ (x; θ, β, m ) and Z (m ; θ, β), in order to obtain the limits of the first and second moments. First, we note that both the invariant density p and the first moment m depend on . For a fixed > 0, it is straightforward to check that the two-scale potentials verify the conditions presented in Arnold et al. (1996, Eqs. (3.1), (3.2)), as long as the nonseparable fluctuations are truncated outside the interval [−a, a]-this is the case for us; see Table 1 in Sect. 4. This, by estimates (Arnold et al. 1996, Eqs. (3.7), (3.8)), implies uniform boundedness of the first moment, m , as well as existence of a unique global weak solution for the McKean-Vlasov equation. We can therefore extract a converging subsequence that converges to some m ∈ R. We use the notation (2.7). We note that V e f f depends smoothly on m. We use the convergence of m to m and (3.6) to deduce: Combining (3.7) and (3.8), we obtain (3.9) Arguing in a similar way for the variance, we conclude that We conclude that Eqs. (3.9) and (3.10) are different, from (2.24) and (2.26).
The two limits → 0, N → ∞, T → ∞ and N → ∞, T → ∞, → 0 commute in the case where the fluctuations in the potential are independent of the macroscale x, V 1 = V 1 (y) in (2.7). An immediate corollary of the above proposition is the following. Proof When the fluctuations are separable (i.e., V 1 (x, y) does not depend on x), ψ(x, β) in (2.24), (2.26) becomes a constant that we can ignore since it also appears in the partition function and they cancel out. Similarly, the terms of the form R e −βV 1 (x,y) dy in Eqs. (3.9) and (3.10) become constants independent of x and cancel with the corresponding terms in the partition function (3.7).
To illustrate the fact that the two limits do commute when the fluctuations are independent of the macroscale, we present in Figs. 3 and 4 the plots of R(m ; θ, β) for various values of and fixed β and θ , which we compare with the solution of the homogenized self-consistency equation R(m; θ, β) = m. We present results both for a convex and nonconvex confining potential, with periodic fluctuations. More details about the two-scale potentials that we use for the numerical simulations will be given in Sect. 4.
As is evident from Fig. 2a, the oscillatory part of the potential introduces (infinitely many) additional local minima. Consequently, Tugaut (2014), the self-consistency equation R(m ; θ, β) = m has multiple solutions. Furthermore, as shown in Fig. 3b, in the limit → 0, the curves R(m ; θ, β) (various dashed lines) approach those given by R(m; θ, β) computed from Eq. (2.24) (full black line), in accordance with Corollary 3.2, showing the commutativity of the two limits.
Let us consider now the case of nonseparable fluctuations. As we have already discussed, see Fig. 2b and also the inside panels of Figs. 7a and 9a, the resulting two-scale potential does not only contain many additional local minima, but it is also

Multiscale Analysis for the McKean-Vlasov Equation in a Two-Scale Confining Potential
In this section, we study the problem of periodic homogenization for the McKean-Vlasov equation in a locally periodic confining potential, for the Curie-Weiss quadratic interaction and in one dimension. We only present formal arguments. The rigorous analysis of this problem will be presented elsewhere. We consider the nonlinear Fokker-Planck equation (3.2) with F(x) = x 2 2 : with initial conditions p (x, 0) = p in (x), independent of and where the prime denotes differentiation with respect to x. The PDE (3.11) is coupled to the selfconsistency equation This homogenization problem is (slightly) different from the standard one for the Fokker-Planck equation in a two-scale potential that was studied in Duncan et al. (2016a), Duncan et al. (2016b) due to the self-consistency equation (3.12). In particular, in addition to the standard two-scale expansion for the solution of the Fokker-Planck equation (3.11), we also need to expand the solution of (3.12) into a power series in : where, as usual (Pavliotis and Stuart 2008), we take { p j = p j (x, ·, t) , j = 0, 1, . . . } to be L-periodic in their second argument. Substituting (3.13) into (3.11) and (3.12) and using the standard tools from the theory of periodic homogenization, e.g., Fredholm's alternative, we obtain the homogenized equation (2.19), satisfied by the marginal of the first term in the two-scale expansion p(x, t) = L 0 p(x, y, t) dy and with the partial free energy ψ(x) given by (2.20) and with The convergence of m (t) to m(t) can be justified using the a priori estimates on moments of the solution to the McKean-Vlasov equation that were derived in Arnold et al. (1996), in particular (Arnold et al. 1996, Eqs. (3.1), (3.2)).
Alternatively, we can work with the backward Kolmogorov equation: We recall that Eq. (3.11) corresponds to the McKean SDE (3.14) with V (x) = V x, x . We introduce the auxiliary variable y t = x t , see, e.g., Pavliotis and Stuart (2007), and using the chain rule, we can write (3.14) as a system of interacting diffusions across scales, driven by the same Brownian motion, (3.16) We start by expanding the first moment m in powers of as in (3.13b). The backward Kolmogorov equation for the observable u (x, y, t) = E( f (x t , y t )|x 0 = x, y 0 = y) reads (neglecting terms of O( ) that are due to the expansion of m ) We can now proceed with the analysis of (3.17a), first for the choice f (x) = x, i.e., the evolution of the first moment, and then for arbitrary observables. We obtain, thus, the homogenized backward Kolmogorov equation, from which we can read off the homogenized McKean SDE and the corresponding Fokker-Planck equation: where ψ(x) = −β −1 ln L 0 e −βV 1 (x,y) dy and M(x) is defined in (2.10). For the sake of brevity, we will omit the details.

Numerical Simulations
In this section, we construct the bifurcation diagram for the stationary McKean-Vlasov equation (both for finite values of and in the homogenization limit), present the results of Monte Carlo (MC) simulations based on the numerical solution of the particle/SDE approximation and also solve the time-dependent McKean-Vlasov PDE. Our goal is to investigate numerically the issue of (lack of) commutativity of the mean field and homogenization limits. We consider interacting diffusions (and the corresponding McKean-Vlasov) in one dimension and we study two types of large-scale and fluctuating parts of the potential. We consider both convex and nonconvex potentials, and both additive (separable) and multiplicative (nonseprarable) fluctuations. The four potentials that we use for our simulations are tabulated in Table 1. We remark that the nonseparable fluctuations V × 1 (x) are truncated outside the interval [−a, a] in order to prevent the oscillations from growing as |x| → +∞. 6 We note that this is necessary for the proof of the homogenization theorem in Duncan et al. (2016b) and that, furthermore, it ensures that the a priori estimates on the moments from Arnold et al. (1996) hold. 7 Throughout this section, we consider fluctuations which have period L = 2π . In all cases, we will consider the Curie-Weiss interaction potential F(x) = x 2 2 , and throughout Sections 4.1 and 4.2, we will fix the interaction strength to be θ = 5. We choose this value because larger values of θ allow for bifurcations to occur at higher temperatures, i.e., lower β, which is easier to handle numerically. In fact, the relevant bifurcation parameter for our problem is given by the combination βθ; see Eq. (1.10). Fixing θ allows us to construct the bifurcation diagram by varying only the temperature. It is also clear from Eq. (1.10) that this equation has no solutions for negative values of θ , i.e., that no (pitchfork) bifurcations can occur for θ < 0.
Using Eq. (2.16), we note that the diffusion coefficient for separable fluctuations in the potential is independent of x and is given by (4.1) where I 0 (·) is the modified Bessel function of the first kind (Duncan et al. 2016a). On the other hand, for nonseparable fluctuations (cases 2 and 4 in Table 1) we obtain (4.2) Furthermore, we obtain the following formulas for the partition functions (4.3) We can now solve the self-consistency equation (1.9) and the equation for the critical temperature (1.10) for the various potentials given in Table 1. We will track each branch of the bifurcation diagram using arclength continuation, which will enable us to plot the first moment m as a function of the inverse temperature β for a fixed value of the interaction strength θ . We do this using the Moore-Penrose quasi-arclength continuation algorithm. 8 The stability of each branch was determined in two different ways: First, we checked whether it corresponded to a local minimum or maximum of the confining potential. Second, we solved the time-dependent McKean-Vlasov equation-see details in Sect. 4.5-using a perturbation of the steady state belonging to each branch (for a particular value of β and θ ) as initial condition. Finally, we have confirmed the stability of each branch by computing the free energy (1.5) of a steady state from that branch at a particular value of β, chosen so that all the branches plotted were present. Stable branches, plotted in blue in all the figures presented in this section, correspond to local minimizers of the free energy functional; unstable branches, plotted in red, correspond to local maxima of the free energy.

Mean Field Limit of the Homogenized System of SDEs: The → 0, N → ∞ Limit
As discussed before (see discussion of Corollary 3.2), when the fluctuations are separable the partial free energy ψ(x) defined in Eq. (2.20) drops out from the homogenized stationary Fokker-Planck equation. This implies, in particular, that the invariant measure(s) of the homogenized dynamics is(are) independent of the fluctuating part of the potential. In particular, there are still no phase transitions when the large-scale part of the potential is convex and still only one pitchfork bifurcation for the bistable potential case-see Fig. 5-where two new, stable, branches emerge from the zero mean solution. We note that in this case the homogenized confining potential in the homogenized equation depends on the inverse temperature β; see the inside panels in Fig. 5. In particular, the values of the local minima of the effective potential are shifted, although their location remains the same, and there are no changes in the topology of the bifurcation diagrams. For nonseparable fluctuations, the mean field and homogenization limits do not commute (see Prop. 3.1). In fact, the homogenization procedure can convexify the effective potential, and we still see no bifurcations when the large-scale part of the potential is convex, while for the bistable potential there is still only one phase transition. The effect of fluctuations on the bifurcation diagram is visible by a shift of the critical temperature at which the phase transition occurs.
where p ∞ (x; θ, β, m) is a stationary solution of the McKean-Vlasov equation. We start the algorithm at a sufficiently large β 0 , i.e., at a sufficiently low temperature for which we have a good initial guess for the value of the order parameter. Since there are no phase transitions for the convex potential (cases 1 and 2 in Table 1), we do not present numerical results for this case. We present in Fig. 5 the plots of R(m; θ, β) and the bifurcation diagrams for the bistable potential with separable and nonseparable fluctuations (cases 3 and 4, respectively). We observe that, for nonseparable fluctuations, the function R(m; θ, β) is flat around m = 0; see Fig. 5b. As we have already mentioned, the topology of the bifurcation diagram does not change, in comparison with that of the bistable potential V b 0 (x) with no fluctuations (see Fig. 1b for this case); thus, the effect of fluctuations is only observed by a shift in the critical temperature.

Mean Field Limit of the Multiscale System of SDEs: Effects of Finite
In this section, we present numerical results on the bifurcation diagram when we first pass to the mean field limit, while keeping small but finite. We are particularly interested in the finite effects on the bifurcation diagrams for the two-scale potentials presented in Table 1.

Convex Confining Potential with Separable and Nonseparable Fluctuations
We first consider Case 1 in Table 1: a convex large-scale potential with separable fluctuations. We present in Fig. 6 the solution to the self-consistency equation R(m; θ, β) = m, the two-scale potential, and the bifurcation diagram for this case. For all finite values of , the resulting potential is nonconvex. This results in the selfconsistency equation having multiple solutions (in fact, as → 0, there are infinitely many solutions). In addition to the emerging pitchfork bifurcation (second-order, or continuous, phase transition), we observe the emergence of discontinuous branches that correspond to metastable states, since they are not (global) minimizers of the free energy; see the results presented in Table 2.
Next, we consider the second case in Table 1: a convex large-scale potential V c 0 (x) with nonseparable fluctuations. Similarly, we present in Fig. 7 the solution to the selfconsistency equation R(m ; θ, β) = m , the two-scale potential and the bifurcation diagram. We note that, as we mentioned before, we restrict the nonseparable fluctuations to a finite interval. In our computations, we use a = 5, in the characteristic function in Table 1.
We observe in Fig. 7b that no pitchfork bifurcations appear; all new branches that appear do not emerge continuously from the mean-zero solution. This is due to the flatness observed in the potential around m = 0 (see Fig. 7a). Furthermore, the mean-   Table 7c, where they are listed in the same way as in Table 2, i.e., in decreasing order of nonnegative m. The free energies of the different branches are presented in Fig. 7d. These new branches correspond to metastable states.
We have checked the stability of each branch by computing the free energy (1.5) of a steady state from that branch at a particular value of β, chosen so that all the branches plotted were present. We summarize the results in Table 2. Since we only consider symmetric potentials, it is sufficient to calculate the free energy for the branches with, say, nonnegative values of m. In each column of Table 2, the values of the free energy are presented from the branch with largest value of m to the lowest; the last value presented in each column corresponds to the branch with m = 0. We summarize the results in Table 2.
We observe that the branch corresponding to a pitchfork bifurcation (i.e., secondorder phase transition), when present, has the lowest value of the free energy, i.e., it is the globally stable one. Furthermore, when a pitchfork bifurcation does not occursee Fig. 7-the branch corresponding to m = 0 is the one with the lowest value of the free energy. Finally, we observe that the stability of the branches in Fig. 9b does not alternate in the same manner as in the previous figures. This is due to the flatteness of the potential around x = 0 for nonseparable oscillations.
The results on the stability of the different branches that are reported in this section are preliminary. A more thorough study of the local (linear) and global stability of the stationary states of the McKean-Vlasov dynamics in multiwell potentials will be presented elsewhere. We mention in passing the early rigorous work on the global stability of the steady states for the McKean-Vlasov equation in Tamura (1987) as well as the careful study of the connection between the loss of linear stability of the uniform state and phase transitions for the McKean-Vlasov equation on the torus (without a confining potential) and with finite-range interactions in Chayes and Panferov (2010).

Bistable Confining Potential with Separable and Nonseparable Fluctuations
Here we consider cases 3 and 4 in Table 1, the bistable potential V b 0 (x). In this case, the large-scale potential exhibits a second-order phase transition even in the absence of small-scale fluctuations (see the pitchfork bifurcation in Fig. 1b) due to the existence of two local minima for V b 0 (x). We are interested in analyzing the topological changes that rapid oscillations in the potential induce to the bifurcation diagram.
We start with separable potentials-see Fig. 8. We observe that the self-consistency equation R(m ; θ, β) = m exhibits a larger number of solutions for finite , which, as for the convex case, result in the emergence of metastable states that are not continuously connected with the mean-zero Gibbs state.
Similarly, for the last case in Table 1, case 4 (bistable potential V b 0 (x) and nonseparable fluctuations), there are more solutions to the self-consistency equation. However, the flatness of the potential (and therefore of the curves R(m; θ, β) near m = 0) reduces the number of additional branches. Moreover, the topological structure of the bifurcation diagram changes, and we now observe a nonparabolic curve for the main branch, which bifurcates from the mean-zero solution via a pitchfork bifurcation.

Numerical Study of the Critical Temperature as a Function of
Here we study the influence of finite on the critical temperature β C , the solution of (3.10) for two-scale potentials, after which continuous phase transitions (pitchfork bifurcations) occur. We do this by solving the equation (we only consider symmetric potentials) for the various potentials in Table 1.  Table 1 We present in Fig. 10 plots of the critical temperature, β C as a function of for a fixed θ = 5. The results are presented for cases 1 (Fig. 10a), 3 (Fig. 10b) and 4 ( Fig. 10c) from Table 1. We do not present the remaining case because, as shown in Fig. 7b, there is no pitchfork bifurcation from the mean-zero solution for case 2. The dependence of the critical temperature on is different for separable and nonseparable potentials. It appears that the critical temperature can change considerable by varying , which implies that a different number of branches might be present in the bifurcation diagram at a fixed temperature, for different values of . This issue will be studied in detail in future work.

Simulations of the Interacting Particles System
In this section, we present the results of Monte Carlo (MC) simulations for the system of interacting diffusions, both for the full, i.e., -dependent, (2.1) and for the homogenized dynamics (2.17). Our focus is on the study of the convergence of the interacting particles system to their equilibrium state. It should be emphasized that no phase transitions occur for the finite-dimensional particles system. However, the numerical simulation of the two interacting particles systems, (2.1) and the homogenized particle system (2.17) clearly exhibit the lack of commutativity between the mean field and homogenization limits.
For the full dynamics (2.1), we used δ = 1 and = 0.1. We solved the SDEs using the Euler-Maruyama scheme. For the homogenized dynamics (2.17), since the noise is multiplicative (for nonseparable potentials), we used the Milstein scheme. In both cases, the time step used was dt = 0.01, which is of O( 2 ). Finally, in both cases we initialized the N particles as being normally distributed, with mean zero and variance 4, which was large enough so that all the local minima were contained within two standard deviations of the Gaussian distribution.
In Figs. 11, 12 and 13, we present the results of our simulations for case 1 in Table 1 the results for = 0.1, while the right panels show the results for the homogenized system. In Fig. 12, we present snapshots of the histogram for the N = 1000 particles for the same time and parameter values, which are δ = 1, = 0.1, θ = 2 and β = 8. On the t = 5000 snapshot, we superpose the corresponding invariant measure, rescaled for comparison, and we observe that the empirical density of the system of interacting diffusions converges to the steady-state solution computed by solving the stationary McKean-Vlasov equation. We also calculate the empirical average of the interacting particle system as a function of time t. We observe that in both cases, the average converges to 0 as expected, but that the convergence for the homogenized SDE (2.17) is slower. The position of the N particles follows approximately the same qualitative behavior (with the particles clustering close to 0), but as we can see from the corresponding histogram there exist additional wells (nonconvexity) for the finite case.
We performed similar experiments for case 4 in Table 1 Here we used N = 500 particles, and smaller values of θ and β. The parameters used were θ = 0.5, β ≈ 5.6, δ = 1 and = 0.1, and the results are plotted in Figs. 14, 15 and 16. In Fig. 14, we present snapshots of the position of each of the N = 500 particles for t = 0 (top panels), t = 100 (middle panels) and t = 5000 (bottom panels). The left panels show the results for = 0.1, while the right panels show the results for the homogenized SDE (2.17). Here we can observe the noncommutativity of the limits: The particles evolve toward different steady states, which shows the effect of the fluctuations on the critical temperature β C at which phase transitions occur. This will be confirmed below when we present the mean value of the solution.
We present in Fig. 15 snapshots of the histogram of the N = 500 particles at t = 0, t = 100 and t = 5000. Again, we observe that the particles converge to different equilibria, the homogenized system converging to a mean-zero distribution with peaks at 1 and − 1, while for positive values of the parameter the system converges to a distribution withX t = − 1. Similarly to the previous case, we superpose the corresponding invariant measure, rescaled for comparison, for this parameter regime on the t = 5000 snapshot, and again we observe that the empirical density of the system of interacting diffusions converges to the steady-state solution computed by solving the stationary McKean-Vlasov equation, which is also obtained by time evolution of the Fokker-Planck equation (see Figs. 17,18 in Sect. 4.5).
Finally, we plot in Fig. 16 the averageX N t of the N = 500 particles for the case of a bistable large-scale potential with nonseparable fluctuations. We observe here that the critical temperature for the homogenized dynamics is different than that for the full dynamics. In particular, the phase transition occurs for β ≈ 10.4 > 5.6 for the homogenized problem, while for finite values of there already exist several branches at this value of β.

Time-Dependent McKean-Vlasov Evolution
We performed time-dependent simulations of the evolution of the nonlinear McKean-Vlasov equation both for the full and for the homogenized dynamics. We present below the results corresponding to the cases presented for the Monte Carlo simulations. We recall that, for the case when we take N → ∞ first while keeping > 0 fixed, the McKean-Vlasov-Fokker-Planck equation that we need to solve is whereas for the case when we first homogenize the dynamics and then pass to the mean field limit the McKean-Vlasov equation becomes with ψ(x) and M(x) given by (2.20) and (2.16), respectively. To solve the McKean-Vlasov evolution PDE, we used the positivity preserving, entropy decreasing finite volume scheme from Carrillo et al. (2015). We point out that this scheme solves the equations using no-flux boundary conditions. We use these boundary conditions and a sufficiently large domain. We used the same initial conditions for the time-dependent Fokker-Planck simulations as the ones used for the Monte Carlo simulations, i.e., the initial condition was the PDF for a normal distribution with mean zero and variance 4. However, for the bistable large-scale potential with nonseparable fluctuations in the finite but positive case-see left panel in Fig. 18-we needed to use a different initial condition: Here we used a normal distribution with mean − 0.1 and variance 4. This is likely because the value of β we chose here was close to the bifurcation point and the mean-zero solution was still being picked up on the time evolution.
We present below the results for the case of a convex large-scale potential V c 0 with separable fluctuations-the same case presented in Figs. 11, 12 and 13. The parameters used were θ = 2, β = 8, δ = 1, = 0. its steady state faster for the full dynamics than for the homogenized equation. This observation can be quantified by comparing the convergence rates in the weighted L 2 or relative entropy exponential estimates, in particular by comparing the constants in the Poincaré and logarithmic Sobolev inequalities for the full and for the homogenized dynamics. A preliminary study of this-for the Fokker-Planck operator of the finitedimensional dynamics-was presented in Duncan et al. (2016b).
Finally, we present numerical results for the case of a bistable large-scale potential V b 0 with nonseparable fluctuations-the same case presented in Figs. 14, 15 and 16. The parameters used were θ = 0.5, β ≈ 5.6, δ = 1, = 0.1.
As expected, the solutions converge to those computed by solving the stationary McKean-Vlasov equation and are qualitatively similar to those obtained from the particle system simulations; see Fig. 15. In this case, the solution of the time-dependent McKean-Vlasov PDE converges to a steady state slower for the full dynamics, in comparison with the homogenized dynamics. We believe that this is related to the phenomenon of critical slowing down (Shiino 1987) when the dynamics is close to a bifurcation, since the inverse temperature β −1 that we use for the simulations is close to the critical temperature β −1 C for the full dynamics.

Conclusions and Further Work
The combined mean field and homogenization limit for a system of interacting diffusions in a two-scale confining potential was studied in this paper. In particular, the homogenized McKean-Vlasov equation was obtained and studied and the bifurcation diagram for the stationary states was considered. It was shown, by means of analysis and extensive numerical simulations, that the homogenization and mean field limits, at the level of the bifurcation diagram (i.e., when combined with the long time limit), do not commute for nonseparable two-scale potentials. Furthermore, it was shown that the bifurcation diagrams can be completely different for small but finite and for the homogenized McKean-Vlasov equation.
It should be emphasized, as is clearly explained in Chayes and Panferov (2010), see in particular the remarks at the end of Sec. 2 of this paper, that the connection between bifurcations and phase transitions for the McKean-Vlasov dynamics is not entirely straightforward. In particular, in order for a bifurcation point to correspond to a genuine phase transition, it is not sufficient to have the emergence of a new branch of solutions, but these emergent solutions should have a lower free energy. More precisely, it was shown in Chayes and Panferov (2010) for the McKean-Vlasov dynamics on the torus and with a finite-range interaction potential that the loss of linear stability of the uniform state-which corresponds to the mean-zero Gibbs state in our setting-does not imply a second-order phase transition. Furthermore, the critical temperature (or, equivalently, critical interaction strength) at which first-order phase transitions occur is lower than the temperature at which the pitchfork bifurcation happens. For the problem that we studied, supercritical pitchfork bifurcations occur which correspond to second-order (continuous) phase transitions. On the other hand, when only saddle node bifurcations are present, e.g., in Fig. 7b, then the mean-zero solution is still the global minimizer of the free energy; see Fig. 7d. In particular, no first-order phase transitions seem to appear in the McKean-Vlasov model that we studied in this work.
There are many open questions that are not addressed in this work. First, the rigorous multiscale analysis for the McKean-Vlasov equation in locally periodic potentials needs to be carried out. Perhaps more importantly, the rigorous construction of the bifurcation diagram in the presence of infinitely many local minima in the confining potential, thus extending the results presented in, e.g., Dawson (1983), Tamura (1984), Tugaut (2014), appears to be completely open. Furthermore, the study of the stability of stationary solutions to the McKean-Vlasov equation in the presence of a multiscale structure, as well as the analysis of the problem of convergence to equilibrium in this setting is an intriguing question. Finally, the extension of the work presented in this paper to higher dimensions presents additional challenges. We mention, for example, that the corresponding nonlinear diffusion process does not have to be reversible (Lelievre et al. 2013;. We believe that the results reported in this work open up a new exciting avenue of research in the study of mean field limits for interacting diffusions in the presence of many local minima, with potentially interesting applications to the study of McKean-Vlasov-based mathematical models in the social sciences.