Mean-field limits for interacting diffusions with colored noise: phase transitions and spectral numerical methods

TL;DR

This paper investigates how colored noise influences phase transitions in interacting particle systems, employing spectral methods to solve McKean-Vlasov equations and validating results with simulations.

Contribution

It introduces a spectral method for solving mean-field Fokker-Planck equations with colored noise, analyzing bifurcations and phase transitions.

Findings

Colored noise affects bifurcation diagrams of particle systems.

Spectral method accurately solves complex Fokker-Planck equations.

Small correlation time regime analyzed via perturbation theory.

Abstract

In this paper we consider systems of weakly interacting particles driven by colored noise in a bistable potential, and we study the effect of the correlation time of the noise on the bifurcation diagram for the equilibrium states. We accomplish this by solving the corresponding McKean-Vlasov equation using a Hermite spectral method, and we verify our findings using Monte Carlo simulations of the particle system. We consider both Gaussian and non-Gaussian noise processes, and for each model of the noise we also study the behavior of the system in the small correlation time regime using perturbation theory. The spectral method that we develop in this paper can be used for solving linear and nonlinear, local and nonlocal (mean-field) Fokker-Planck equations, without requiring that they have a gradient structure.

Tables

Table 1. Table 1 : Value of ζ 𝜁 \zeta

Model	OU	H	B	NS
$\zeta$	$1/\sqrt{2}$	$1/\sqrt{2}$	$0.624$	$0.944$

Equations

$$\frac{\mathrm{d}{X}_t^i}{\mathrm{d}t}=-\left(V^{\prime }(X_t^i)+\theta \left(X_t^i-\frac{1}{N}\sum _{j=1}^N{X}_t^j\right)\right)+\sqrt{2{\beta }^{-1}}{\xi }_t^i,i=1,\dots ,N,$$

$$\frac{\partial \rho }{\partial t}=\frac{\partial }{\partial x}\left(V^{\prime }(x)\rho +\theta \left(x-{\int }_{\mathbf{R}}x\rho (x,t)\mathrm{d}x\right)\rho +{\beta }^{-1}\frac{\partial \rho }{\partial x}\right).$$

$$\mathcal{F}[\rho ]={\beta }^{-1}{\int }_{\mathbf{R}}\rho (x)\ln \rho (x)\mathrm{d}x+{\int }_{\mathbf{R}}V(x)\rho (x)\mathrm{d}x+\frac{\theta }{2}{\int }_{\mathbf{R}}{\int }_{\mathbf{R}}F(x-y)\rho (x)\rho (y)\mathrm{d}x\mathrm{d}y,$$

$$\frac{\partial }{\partial x}\left(V^{\prime }(x){\rho }_{\infty }(x)+\theta \left(x-m\right){\rho }_{\infty }(x)+{\beta }^{-1}\frac{\partial {\rho }_{\infty }}{\partial x}(x)\right)=0,$$

$$m={\int }_{\mathbf{R}}x{\rho }_{\infty }(x)\mathrm{d}x.$$

$$V_{\text{eff}}(x;m,\theta )=V(x)+\frac{\theta }{2}(x-m{)}^2,$$

$${\rho }_{\infty }(x;m,\beta ,\theta ):=\frac{1}{\mathcal{Z}(m,\beta ,\theta )}{e}^{-\beta V_{\text{eff}}(x;m,\theta )},$$

$$m={\int }_{\mathbf{R}}x{\rho }_{\infty }(x;m,\beta ,\theta )\mathrm{d}x=:R(m,\beta ,\theta ).$$

$$\displaystyle\mathcal{F}[\rho_{\infty}(\,\cdot\,;m,\beta,\theta)]=-\beta^{-1}\,\ln\mathcal{Z}(m,\beta,\theta)-\frac{\theta}{2}\big{(}R(m,\beta,\theta)-m\big{)}^{2},$$

$$\displaystyle\mathchoice{\frac{\partial{}}{\partial m}}{\partial_{m}}{}{}\,\mathcal{F}[\rho_{\infty}(\,\cdot\,;m,\beta,\theta)]=-\beta\theta^{2}\big{(}R(m,\beta,\theta)-m\big{)}\,\mathrm{Var}(\rho_{\infty}(\cdot\,;\,m,\beta,\theta)),$$

$$\mathrm{Var}(\psi ):={\int }_{\mathbf{R}}{\left(x-{\int }_{\mathbf{R}}\psi (x)\mathrm{d}x\right)}^2\psi (x)\mathrm{d}x.$$

$$\mathrm{d}{X}_t^i=-\left(V^{\prime }(X_t^i)+\theta \left(X_t^i-\frac{1}{N}\sum _{j=1}^N{X}_t^j\right)\right)\mathrm{d}t+\sqrt{2{\beta }^{-1}}{\eta }_t^i\mathrm{d}t,1\le i\le N,$$

$$\mathrm{d}{\mathbf{Y}}_t^i=A{\mathbf{Y}}_t^i\mathrm{d}t+\sqrt{2}D\mathrm{d}{\mathbf{W}}_t^i,i=1,\dots ,N,$$

$$\mathrm{d}{\eta }_t^i=-{\eta }_t^i\mathrm{d}t+\sqrt{2}\mathrm{d}{W}_t^i.$$

$$K_{OU}(t)=e^{-|t|}.$$

$$A=\left(\begin{array}{cc}0& 1\\ -1& -\gamma \end{array}\right),D=\left(\begin{array}{cc}0& 0\\ 0& \sqrt{\gamma }\end{array}\right),{\mathbf{y}}_{\eta }=\left(\begin{array}{c}1\\ 0\end{array}\right).$$

$$K_H(t)=e^{-\frac{|t|}{2}}\left(\cos \left(\frac{\sqrt{3}}{2}t\right)+\frac{\sqrt{3}}{3}\sin \left(\frac{\sqrt{3}}{2}t\right)\right).$$

$$\mathrm{d}{\eta }_t^i=-V_{\eta }^{\prime }({\eta }_t^i)\mathrm{d}t+\sqrt{2}\mathrm{d}{W}_t^i,$$

$${\int }_{\mathbf{R}}\eta e^{-V_{\eta }(\eta )}\mathrm{d}\eta =0.$$

$$V_{\eta }(\eta )=\frac{{(\eta -\alpha )}^4}{4}-\frac{{(\eta -\alpha )}^2}{2}+(\eta -\alpha ),$$

$$\frac{\partial \rho }{\partial t}=\frac{\partial }{\partial x}\left(V^{\prime }\rho +\theta (x-m(t))\rho -\sqrt{2{\beta }^{-1}}⟨{\mathbf{y}}_i,{\mathbf{y}}_{\eta }⟩\rho \right)+{\mathcal{L}}_{\mathbf{y}}^{\ast }\rho ,$$

$$m(t)={\int }_{\mathbf{R}}{\int }_{{\mathbf{R}}^n}x\rho (x,\mathbf{y},t)\mathrm{d}\mathbf{y}\mathrm{d}x.$$

$${\mathcal{L}}_{\mathbf{y}}^{\ast }\rho =\{\begin{array}{ll}{\partial }_{\eta }\left(\eta \rho +{\partial }_{\eta }\rho \right),& \text{for scalar OU noise,}\\ {\partial }_{\lambda }\left(\lambda \rho +{\partial }_{\lambda }\rho \right)+\left(\eta {\partial }_{\lambda }\rho -\lambda {\partial }_{\eta }\rho \right),& \text{for harmonic noise,}\\ {\partial }_{\eta }\left(V_{\eta }^{\prime }\rho +{\partial }_{\eta }\rho \right),& \text{for non-Gaussian noise}.\end{array}$$

$$\frac{\partial }{\partial x}\left(V^{\prime }(x)\rho +\theta (x-m)\rho -\sqrt{2{\beta }^{-1}}⟨{\mathbf{y}}_i,{\mathbf{y}}_{\eta }⟩\rho \right)+{\mathcal{L}}_{\mathbf{y}}^{\ast }\rho =0,$$

$$m={\int }_{\mathbf{R}}{\int }_{{\mathbf{R}}^n}x\rho (x,\mathbf{y})\mathrm{d}\mathbf{y}\mathrm{d}x.$$

$${\eta }_t^i\to \zeta {\eta }_{t/{\epsilon }^2}^i/\epsilon ,$$

$${\int }_0^{\infty }{\zeta }^{2}K(t/{\epsilon }^2)/{\epsilon }^2\mathrm{d}t={\int }_0^{\infty }{\zeta }^{2}K(t)\mathrm{d}t=\frac{1}{2}.$$

$$\mathrm{d}{X}_t^i=\left(-V^{\prime }(X_t^i)-\theta \left(X_t^i-\frac{1}{N}\sum _{j=0}^N{X}_t^j\right)\right)\mathrm{d}t+\sqrt{2{\beta }^{-1}}\mathrm{d}{W}_t^i,i=1,\dots N,$$

$${\rho }_{\infty }^{\epsilon }(x,\mathbf{y};m,\beta ,\theta )=p_0(x,\mathbf{y};m,\beta ,\theta )+\epsilon p_1(x,\mathbf{y};m,\beta ,\theta )+{\epsilon }^2{p}_2(x,\mathbf{y};m,\beta ,\theta )+\cdots ,$$

$${\rho }_{\infty }^{\epsilon }(x;m,\beta ,\theta )$$

$$=:{\rho }_{\infty }(x;m,\beta ,\theta )+\epsilon p_1(x;m,\beta ,\theta )+{\epsilon }^2{p}_2(x;m,\beta ,\theta )+\cdots ,$$

$${\rho }_{\infty }^{\epsilon }(x;m,\beta ,\theta )={\rho }_{\infty }(x;m,\beta ,\theta )\left[1+{\epsilon }^2\left(C_{OU}-\frac{\beta }{2}{\left(V_{\text{eff}}^{\prime }\left(x\right)\right)}^2+V_{\text{eff}}^{\prime \prime }\left(x\right)\right)\right]+\mathcal{O}({\epsilon }^4)$$

Full text

Keywords: McKean–Vlasov PDEs, Nonlocal Fokker–Planck equations, Interacting particles, Desai–Zwanzig model, Colored noise, Hermite spectral methods, Phase transitions.

AMS: 35Q70, 35Q83, 35Q84, 65N35, 65M70, 82B26,

1 Introduction

Systems of interacting particles appear in a wide variety of applications, ranging from plasma physics and galactic dynamics [5] to mathematical biology [13, 34], the social sciences [17, 37], active media [4], dynamical density functional theory (DDFT) [19, 18] and machine learning [33, 42, 45]. They can also be used in models for cooperative behavior [9], opinion formation [17], and risk management [16], and also in algorithms for global optimization [41].

In most of the existing works on the topic, the particles are assumed to be subject to thermal additive noise that is modeled as a white noise process, i.e. a mean-zero Gaussian stationary process that is delta-correlated in time. There is extensive literature studying the behavior of these systems; we mention for example works on the rigorous passage to the mean field limit [38], the long-time behavior of solutions (see [9, 44] for a case of a ferromagnetic (quartic) potential, and [20] for more general potentials), multiscale analysis [21], and phase transitions [47].

In a more realistic scenario, the system has memory and the hypothesis of Markovianity does not hold [26, 27, 28]. This memory can be modeled by using colored noise, i.e. noise with a nonzero correlation time (or, more precisely, a nonsingular autocorrelation function), which is the approach we take in this paper. For simplicity, we will assume that the noise is additive and that it can be represented by a finite-dimensional Markov process, as in the recent study [11] on mean field limits for non-Markovian interacting particles.

In this paper we will study the dynamics of a system of interacting particles of the Desai–Zwanzig type, interacting via a quadratic Curie–Weiss potential. The system of interacting particles is modeled by a system of stochastic differential equations (SDEs):

[TABLE]

where $N$ is the number of particles, $V(\cdot)$ is a confining potential, $\theta$ is the interaction strength, $\beta$ is the inverse temperature of the system, and $\xi_{t}^{i}$ are independent, identically distributed (i.i.d.) noise processes.

Before discussing the Desai–Zwanzig model with colored noise, we present a brief overview of known results [9, 44] for the white noise problem. When $\xi_{t}^{i}$ are white noise processes, we can pass to the mean field limit $N\to\infty$ in eq. 1.1 and obtain a nonlinear and nonlocal Fokker–Planck equation, known in the literature as a McKean–Vlasov equation, for the one-particle distribution function $\rho(x,t)$:

[TABLE]

The McKean–Vlasov equation (1.2) is a gradient flow with respect to the quadratic Wasserstein metric for the free energy functional

[TABLE]

where $F(x):=x^{2}/2$ is the interaction potential. The long-time behavior of solutions depends on the number of local minima of the confining potential $V$ [47]. It follows directly from eq. 1.2 that any steady-state solution $\rho_{\infty}(x)$ solves, together with its first moment, the following system of equations:

[TABLE]

Since eq. 1.4a is, for $m$ fixed, the stationary Fokker–Planck equation associated with the overdamped Langevin dynamics in the confining potential

[TABLE]

solutions can be expressed explicitly as

[TABLE]

where $\mathcal{Z}(m,\beta,\theta)$ is the normalization constant (partition function); see [9, 20, 21] for more details. By substitution in eq. 1.4b, a scalar fixed-point problem is obtained for $m$, the self-consistency equation:

[TABLE]

The stability of solutions to eq. 1.4 depends on whether they correspond to a local minimum (stable) or to a local maximum/saddle point (unstable) of the free energy functional. The free energy along the one-parameter family (1.6), with parameter $m$, can be calculated explicitly [21],

[TABLE]

from which we calculate that

[TABLE]

where, for a probability density $\psi$,

[TABLE]

Though incomplete, this informal argument suggests that the stability of a steady-state solution can also be inferred from the slope of $R(m,\beta,\theta)-m$ at the corresponding value of $m$: if this slope is positive, the equilibrium is unstable, and conversely. The self-consistency map and the free energy of $\rho_{\infty}(x;m,\beta,\theta)$, for a range of values of $m$, are illustrated in fig. 1 for the bistable potential $V(x)=\frac{x^{4}}{4}-\frac{x^{2}}{2}$.

It is well-known that, when $V(\cdot)$ is an even potential, eq. 1.2 possesses a unique, mean-zero steady-state solution for sufficiently large temperatures (i.e., small $\beta$). As the temperature decreases, this solution loses its stability and two new solutions of the self-consistency equation emerge, corresponding to a pitchfork bifurcation; see [9, 21] for details.

As mentioned above, in this paper we focus on the case where the noise processes $\xi_{t}^{i}$ in eq. 1.1 have a nonzero correlation time, and in particular we assume that each noise process can be represented using a (possibly multi-dimensional) SDE, in which case eq. 1.1 leads to a Markovian system of SDEs in an extended phase space. The colored noise will be modeled by either an Ornstein–Uhlenbeck process, harmonic noise [39, Example 8.2], or a non-Gaussian reversible diffusion process.

Though more realistic, the use of colored noise presents us with some difficulties. First, the introduction of an extra SDE for the noise breaks the gradient structure of the problem; while we can still pass formally to the limit $N\to\infty$ in eq. 1.1 and obtain a McKean–Vlasov equation for the associated one-particle distribution function, it is no longer possible to write a free energy functional, such as eq. 1.3, that is dissipated by this equation. Second, the McKean–Vlasov equation is now posed in an extended phase space, which increases the computational cost of its numerical solution via PDE methods. And third, it is no longer possible to obtain an explicit expression for the one-parameter family of (possible) stationary solutions to the mean field equation, as was possible in eq. 1.6, which renders the calculation of steady states considerably more difficult.

When the correlation time of the noise is small, the latter difficulty can be somewhat circumvented by constructing an approximate one-parameter family of solutions through appropriate asymptotic expansions in terms of the correlation time, from which steady-state solutions of the McKean–Vlasov dynamics can be extracted by solving a self-consistency equation similar to (1.7), see eq. 2.10. Outside of the small correlation time regime, however, finding the steady-states of the McKean–Vlasov equation requires a numerical method for PDEs in all but the simplest cases.

In this work, we propose a novel Hermite spectral method for the time-dependent and steady-state equations, applicable to the cases of both white and colored noise. Discretized in a basis of Hermite functions, the McKean–Vlasov equation becomes a system of ordinary differential equations with a quadratic nonlinearity originating from the interaction term. In contrast with other discretization methods for PDEs, the use of (possibly rescaled) Hermite functions for the problem under consideration leads to an efficient numerical method, first because Hermite functions have very good approximation properties in $L^{2}$, but also because all the differential operators appearing in the McKean Vlasov equation lead to sparse matrices in Hermite space, with a small bandwidth related to the polynomial degree of $V$ (provided that a suitable ordering of the multi-indices is employed). To solve the finite-dimensional system of equations obtained after discretization of the time-dependent equation, we employ either the Runge–Kutta 45 method (RK45) or a linear, semi-implicit time-stepping scheme.

To assess the performance of our numerical method, we compare its efficiency in the white noise case with that of the finite volume scheme developed in [8], the applicability of which depends on the existence a gradient structure of eq. 1.2. We also verify that our results agree with known analytical solutions in simple settings, and with explicit asymptotic expansions in the small correlation time regime. We then use our spectral method, together with asymptotic expansions and Monte Carlo (MC) simulations of the particle system, to construct the bifurcation diagram of the first moment of the steady-state solutions as a function of the inverse temperature.

For the reader’s convenience, we summarize here the main results of this paper:

The systematic study of the effect of colored noise, both Gaussian and non-Gaussian, on the long-time behavior of the McKean–Vlasov mean field equation, including the effect of colored noise on the structure and properties of phase transitions. 2. 2.

The development and analysis of a spectral numerical method for the solution of linear or nonlinear, local or nonlocal Fokker–Planck-type equations. In particular, our method does not depend on an underlying gradient structure for the PDE.

The rest of the paper is organized as follows. In Section 2, we present the models for the colored noise and we derive formally the mean field McKean–Vlasov equation associated with the interacting particle system. In Section 3, we present the numerical methods used to (a) solve the time-dependent and steady-state Fokker–Planck (or McKean–Vlasov) equations and (b) solve the finite-dimensional system of interacting diffusions (1.1). In Section 4, we study the performance of our numerical method in the small correlation time regime, and we verify numerically the convergence rates to the white noise solution in the limit where the correlation time tends to 0. In Section 5, we describe our methodology for constructing the bifurcation diagrams and we present the associated results. Section 6 is reserved for conclusions and perspectives for future work.

2 The model

We consider the following system of weakly interacting diffusions,

[TABLE]

where the noise processes $\eta^{i}_{t}$ are independent, mean-zero, second-order stationary processes with almost surely continuous paths and autocorrelation function $K(t)$. In the rest of this paper, we will assume that the interaction strength $\theta$ is fixed and equal to 1 and we will use the inverse temperature $\beta^{-1}$ as the bifurcation parameter. We will consider two classes of models for the noise: Gaussian stationary noise processes with an exponential correlation function, and non-Gaussian noise processes that we construct by using the overdamped Langevin dynamics in a non-quadratic potential.

Gaussian noise

Stationary Gaussian processes in $\mathbf{R}^{n}$ with continuous paths and an exponential autocorrelation function are solutions to an SDE of Ornstein–Uhlenbeck type:

[TABLE]

where $\mathit{A},\mathit{D}$ are $n\times n$ matrices satisfying Kalman’s rank condition [32, Chapter 9], and $\boldsymbol{\mathbf{W}}_{t}^{i}$, $1\leq i\leq n$, are independent white noise processes in $\mathbf{R}^{n}$. We assume here that the noise is obtained by projection as $\eta_{t}^{i}=\left\langle{\boldsymbol{\mathbf{Y}}_{t}^{i},\boldsymbol{\mathbf{y}}_{\eta}}\right\rangle$, where $\left\langle{\cdot,\cdot}\right\rangle$ denotes the Euclidean inner product, for some vector $\boldsymbol{\mathbf{y}}_{\eta}\in\mathbf{R}^{n}$. Throughout this paper we will consider two particular examples, namely the scalar OU process and the harmonic noise [39, Chapter 8].

(OU)

Scalar Ornstein–Uhlenbeck process:

[TABLE]

The associated autocorrelation function is

[TABLE]

(H)

Harmonic noise:

[TABLE]

In this case the noise is the solution to the Langevin equation, with the first and second components of $\boldsymbol{\mathbf{Y}}$ corresponding to the position and velocity, respectively. Throughout this paper we will assume $\gamma=1$ for simplicity. The associated autocorrelation function of $\eta^{i}$ is given by

[TABLE]

Non-Gaussian noise

In this case, instead of eq. 2.2 we consider

[TABLE]

where now $V_{\eta}$ is a smooth non-quadratic confining potential satisfying the mean-zero condition:

[TABLE]

We consider the following choices for $V_{\eta}$:

(B)

The bistable potential $V_{\eta}(\eta)=\eta^{4}/4-\eta^{2}/2$.

(NS)

The shifted tilted bistable potential

[TABLE]

with the constant $\alpha\approx 0.885$ such that eq. 2.3 is satisfied.

Remark 2.1**.**

For the two Gaussian noise processes we consider, it would have been equivalent (by a change of variables) to include the inverse temperature $\beta$ in the noise equation eq. 2.2 rather than in eq. 2.1. This is not the case for non-Gaussian noise processes, for which including the temperature in the noise equation leads to an effective diffusion coefficient, in the limit as the correlation time tends to 0, with a nonlinear dependence on $\beta$.

2.1 Mean field limit

For weakly interacting diffusions, the derivation of the mean field McKean–Vlasov PDE is a standard, well-known result [9, 44, 38]. When $\xi_{t}^{i}$ in (1.1) are colored noise processes, it is also possible to pass to the mean field limit $N\to\infty$ in eq. 1.1 and to obtain a McKean–Vlasov equation for the one-particle distribution function $\rho(x,\boldsymbol{\mathbf{y}},t)$:

[TABLE]

Here $\boldsymbol{\mathbf{y}}$ are the noise variables, denoted by $(\eta,\lambda)$ in the case of harmonic noise and just $\eta$ otherwise, and

[TABLE]

A formal derivation of the mean field limit is presented in [49], and this derivation can be justified rigorously using the results in [12, 36].

The main goal of this paper is the study of the effect of colored noise on the structure of the bifurcation diagram for the McKean–Vlasov equation with colored noise, eqs. 2.5a and 2.5b. In other words, we want to gain insight into the number of solutions to the following stationary PDE and associated constraint (self-consistency equation):

[TABLE]

Although there still exists, for fixed $\beta$ and fixed $\theta$, a one-parameter family of solutions to (2.6a) (with parameter $m$), which we will denote by $\{\rho_{\infty}(x,\boldsymbol{\mathbf{y}};m,\beta,\theta)\}_{m\in\mathbf{R}}$, no closed form is available for these solutions. This is because the detailed balance condition no longer holds in the presence of colored noise, i.e. the probability flux at equilibrium does not vanish. Here, by probability flux, we mean the argument of the divergence in the Fokker–Planck operator;

see [39, Section 4.6].

2.2 The white noise limit

To study the limit of small correlation time, it will be convenient to rescale the noise as

[TABLE]

where $\varepsilon$ is a time scaling parameter, and $\zeta$ is a model-dependent parameter ensuring that the autocorrelation function of the rescaled noise, given by $\zeta^{2}\,K(t/\varepsilon^{2})/\varepsilon^{2}$, satisfies

[TABLE]

Then the autocorrelation of the noise converges to a Dirac delta when $\varepsilon\to 0$, and it can be shown that, in this limit, the solution of eq. 2.1 converges to that of

[TABLE]

where $W^{i}$, $i=1,\dotsc N$, are independent Wiener processes; see [6] and [40, Chapter 11]. While not strictly necessary, including the parameter $\zeta$ is convenient to obtain simpler formulas. The value of $\zeta$ for each of the noise models considered in this paper is presented in table 1. For the models B and NS, $\zeta$ was calculated numerically and rounded to three significant figures in this table.

In view of the convergence of the solution of the finite-dimensional particle system when $\varepsilon\to 0$, we expect that also the $x$-marginals of the steady-state solutions to the McKean–Vlasov equation with colored noise, obtained by solving eqs. 2.6a and 2.6b, should converge to their white-noise counterparts as $\varepsilon\to 0$. It turns out that this is the case and, using asymptotic techniques from [26], it is possible to approximate the solutions $\rho_{\infty}(x,\boldsymbol{\mathbf{y}};m,\beta,\theta)$ to eq. 2.6a by a power series expansion in $\varepsilon$; using a superscript to emphasize the dependence on $\varepsilon$,

[TABLE]

From eq. 2.7, we obtain a power series expansion for the $x$-marginal by integrating out the noise variable:

[TABLE]

The methodology to obtain expressions for the terms works by substituting eq. 2.7 in eq. 2.6a and grouping the terms in powers of $\varepsilon$ in the resulting equation. This leads to a sequence of equations that can be studied using standard techniques. Details of the analysis leading to an explicit expression of the first nonzero correction in (2.8) can be found in [26, Section 8] for the particular case of the OU noise, and in [49] for the other noise models we consider.

The order of the first nonzero correction in this expansion depends on the model: it is equal to 1 for model NS, to 2 for models OU and B, and to 4 for model H. In all cases, the first nontrivial term in the series expansion (2.8) can be calculated explicitly (possibly up to constant coefficients that have to be calculated numerically). For scalar Ornstein–Uhlenbeck noise, for example, we have, omitting the dependence of $V_{\textrm{eff}}$ (the effective potential defined in eq. 1.5) on $m$ and $\theta$ for notational convenience,

[TABLE]

Here $C_{OU}$ is a constant such that the correction integrates to [math]. Similar expressions can be obtained for the other models; see [49].

Taking into account only the first nontrivial correction, the order of which we denote by $\delta$, the steady-state solutions to the McKean–Vlasov equation with colored noise can be approximated by solving the approximate self-consistency equation

[TABLE]

We show in fig. 2 that the equation $R_{0}(m,\beta,\theta)+\varepsilon^{2}R_{2}(m,\beta,\theta)=m$, for fixed $\beta=10$, $\theta=1$ and $\varepsilon=0.1$, admits three solutions in the case of OU noise, similarly to the case of white noise. This figure was generated using the asymptotic expansion (2.9).

3 The numerical method

In this section, we describe the spectral numerical method that we will use in order to solve the time-dependent McKean–Vlasov equation, eqs. 2.5a and 2.5b, as well as the steady-state equation, eqs. 2.6a and 2.6b. Before looking at colored noise, we consider the case of white noise, for which our method can be tested against the results in [21], which were obtained using the finite volume scheme developed in [8].

3.1 Linear Fokker–Planck equation with white noise

We start by presenting the methodology used in the absence of an interaction term, in which case eq. 2.5a reduces to a linear Fokker–Planck equation:

[TABLE]

We assume that $V(\cdot)$ is a smooth confining potential and, consequently, the unique invariant distribution is given by $\rho_{s}=\frac{1}{\mathcal{Z}}\operatorname{e}^{-\beta V}$, where $\mathcal{Z}$ is the normalization constant [39, Proposition 4.2]. The Fokker–Planck operator in eq. 3.1 is unitarily equivalent to a Schrödinger operator; see [1] and [39, Section 4.9]. Defining $u=\rho/\sqrt{\rho_{s}}$, the function $u$ satisfies

[TABLE]

with the initial condition $u(x,t=0)=\rho_{0}/\sqrt{\rho_{s}}=:u_{0}$. Several works made use of Hermite spectral methods to study equations of this type, e.g.  [1, 14, 15]. The Schrödinger operator on the right-hand side of eq. 3.2 is selfadjoint in $L^{2}(\mathbf{R})$ and it has nonpositive eigenvalues. Under appropriate growth assumptions on the potential $V(x)$ as $x\to\infty$, it can be shown that its eigenfunctions decrease more rapidly than any exponential function in the $L^{2}(\mathbf{R})$ sense, in that they satisfy $\operatorname{e}^{\mu\left|x\right|}\,\varphi(x)\in L^{2}(\mathbf{R})$ for all $\mu\in\mathbf{R}$; see [15] and also [2] for a detailed study. Under appropriate decay assumptions at infinity on the initial condition, we expect the solution to eq. 3.2 to also decrease rapidly as $\left|x\right|\to\infty$.

We denote by $\mathbf{P}(d)$ the space of polynomials of degree less than or equal to $d$, and by $\left\langle{\cdot,\cdot}\right\rangle$ the usual $L^{2}(\mathbf{R})$ inner product. For a quadratic potential $V_{q}=\frac{1}{2}\left(\frac{x}{\sigma}\right)^{2}$, with $\sigma$ a scaling parameter, the Galerkin method we employ consists in finding $u_{d}(t)\in\operatorname{e}^{-V_{q}/2}\,\mathbf{P}(d)$ such that

[TABLE]

Here the subscript $\hat{d}\geq d$ on the right-hand side of eqs. 3.3a and 3.3b indicates that the inner product is performed using a numerical quadrature with $\hat{d}+1$ points. With appropriately rescaled Gauss–Hermite points, inner products calculated using the quadrature are exact for functions in $\operatorname{e}^{-V_{q}/2}\,\mathbf{P}(\hat{d})$,

[TABLE]

which is why we did not append the subscript $\hat{d}$ to the inner products in the left-hand side of eqs. 3.3a and 3.3b. When $V$ is a polynomial, it is possible to show using the recursion relations eqs. A.2 and A.1 in appendix A, that the inner product $\left\langle{\mathcal{H}_{x}u_{d},w_{d}}\right\rangle_{\hat{d}}$ on the right-hand side of eq. 3.3a is exactly $\left\langle{\mathcal{H}_{x}u_{d},w_{d}}\right\rangle$ when $\hat{d}\geq d+\deg(|V^{\prime}|^{2})$. This is the approach we take in all the numerical experiments presented in this paper, and we will therefore omit the subscript $\hat{d}$ in eq. 3.3a from now on.

The natural basis of $\mathbf{P}(d)$ (from which a basis of $e^{-V_{q}/2}\,\mathbf{P}(d)$ follows) to obtain a finite-dimensional system of differential equations from the variational formulation (3.3a) is composed of rescaled Hermite polynomials $H_{i}^{\sigma}(x):=H_{i}(x/\sigma)$, $0\leq i\leq d$, where $H_{i}(x)$ are the Hermite polynomials orthonormal for the Gaussian weight $\mathcal{N}(0,1)$; the corresponding basis functions of $\operatorname{e}^{-V_{q}/2}\,\mathbf{P}(d)$ are then rescaled Hermite functions. The fundamental results on Hermite polynomials, Hermite functions and the related approximation results that are used in this paper are summarized in appendix A. In general, in addition to the rescaling parametrized by $\sigma$, a translation could be applied in order to generate more suitable basis functions (e.g. when most of the mass of $\operatorname{e}^{-\beta V}$ is itself localized away from $x=0$), but for simplicity we confine ourselves in this work to the case where $V_{q}$ is symmetric around $x=0$.

Remark 3.1**.**

Although eq. 3.2 and the associated variational formulation (3.3a) are convenient for analysis purposes, for numerical purposes it is useful to perform a second unitary transformation; defining $v=\operatorname{e}^{V_{q}/2}\,u=:u/\sqrt{\rho_{q}}$, the function $v$ satisfies

[TABLE]

with the initial condition $v(x,t=0)=\rho_{0}/\sqrt{\rho_{s}\,\rho_{q}}$, and to approximate the solution to this equation in $\mathbf{P}(d)$. While clearly equivalent, this approach enables us to work directly with Hermite polynomials, for which a range of free and open-source software tools are available, e.g. in the NumPy package for scientific computing in Python [29]. Building upon the tools provided by NumPy, we have developed a thin Python library, available online [48], that offers the possibility of fully automating spectral methods based on Hermite polynomials. We note that, when $V_{q}=\beta V$, eq. 3.4 is merely the backward Kolmogorov equation corresponding to eq. 3.1.

It is possible to prove the convergence of the method presented above in the limit as $d\to\infty$ given appropriate additional assumptions on the confining potential $V(\cdot)$. For simplicity we will make the following assumption, which is satisfied for the bistable potential that we consider in this work, but we note that less restrictive conditions would be sufficient.

Assumption 3.1**.**

The confining potential $V(\cdot)$ is a polynomial of (even) degree greater than or equal to 2. Consequently, it satisfies

[TABLE]

for constants $C_{1},C_{2},C_{3}>0$ and a natural number $k\geq 1$.

We will denote by $H^{m}(\mathbf{R};\mathcal{H}_{x})$ the Hilbert space obtained by completion of $C^{\infty}_{c}(\mathbf{R})$, the space of smooth compactly supported functions, with the inner product

[TABLE]

The norm associated with this Sobolev-like space will be denoted by $\|\cdot\|_{m,\mathcal{H}_{x}}$.

Theorem 3.1**.**

Suppose that 3.1 holds and that the initial condition $u_{0}$ is smooth and belongs to $H^{m}(\mathbf{R};\mathcal{H}_{x})$ for some natural number $m\geq 2k$, where $k$ is as in 3.1. Then for any $d\geq m-1$, any final time $T$ and for all $\alpha>0$, it holds that

[TABLE]

for a constant $C_{\alpha}$ not depending on $d$, $u_{0}$, or $T$, and where $\|\cdot\|$ denotes the $L^{2}(\mathbf{R})$ norm.

Proof.

See appendix B. ∎

Remark 3.2**.**

Theorem 3.1* is not optimal. One one hand, it overestimates the error for large times: both the numerical and exact solutions converge to stationary solutions as $t\to\infty$, so we expect the error $\|u(t)-u_{d}(t)\|^{2}$ to tend to finite limit when $t\to\infty$. Although the error between the stationary solutions can be bounded similarly to the transient error, see remark B.1, we have not obtained a result that combines both errors; we plan to return to this interesting question in future work. On the other hand, the bound on the transient error of theorem 3.1 is probably not sharp. Indeed, when the initial condition $u_{0}$ is smooth and, together with all its derivatives, decreases exponentially as $x\to\infty$, theorem 3.1 implies only that the error decreases faster than any negative power of $d$. In most practical examples, however, we observed numerically that the convergence is in fact exponential.*

Remark 3.3**.**

The condition that $u_{0}\in H^{m}(\mathbf{R};m,\mathcal{H}_{x})$ is quite restrictive. It requires in particular that $u_{0}\in L^{2}(\mathbf{R})$, which is equivalent to requiring that $\rho_{0}\in L^{2}(\mathbf{R};\rho_{s}^{-1})$, because $u_{0}=\rho_{0}/\sqrt{\rho_{s}}$ by definition. Though natural from an $L^{2}$-theory perspective, see [39, Sec. 4.5] and [35], this condition excludes a large class of initial conditions. If $V(x)$ behaves as $|x|^{4}$ as $|x|\to\infty$, then it excludes Gaussian initial conditions, for example.

Remark 3.4**.**

An alternative manner of solving eq. 3.1 numerically is to approximate the eigenvalues and eigenfunctions of the operator $\mathcal{H}_{x}$ in terms of Hermite functions, after which integration in time becomes trivial. This approach requires finding (or approximating) the $d+1$ solutions $(\varphi_{d}^{i},\lambda_{d}^{i})$, $\varphi_{d}^{i}\in\operatorname{e}^{-V_{q}/2}\,\mathbf{P}(d)$, $0\leq i\leq d$, of the eigenvalue problem

[TABLE]

Once these have been calculated, an approximation of the solution to eq. 3.2 is obtained simply as

[TABLE]

It is readily seen that $u_{d}$, defined by this equation, is also the unique solution of (3.3a). This is the approach taken in [7] for the calculation of drift and diffusion coefficients in the diffusion approximation of kinetic equations. While equivalent, this methodology is more computationally expensive because it requires the full solution of the eigenvalue problem (3.5).

3.2 McKean–Vlasov equation with white noise

In the presence of an interaction term, the Fokker–Planck equation becomes nonlinear:

[TABLE]

For this equation the weighted $L^{2}(\mathbf{R};\operatorname{e}^{V})$ energy estimate of the linear case (B.7) does not hold, and there is therefore no longer a natural space for the Galerkin approximation.

Because of this, and since we would like to employ the spectral numerical method with Gaussian initial conditions, which is not possible with a variational formulation of the type (3.3) in view of remark 3.3, we will use Hermite functions to approximate the solution to eq. 3.7 directly, i.e. we will look for an approximate solution in the space $\operatorname{e}^{-V_{q}/2}\,\mathbf{P}(d)$.

The variational formulation corresponding to the Galerkin approximation is then to find $\rho\in\operatorname{e}^{-V_{q}/2}\,\mathbf{P}(d)$ such that

[TABLE]

Dividing by $\left\langle{1,\rho_{d}}\right\rangle_{\hat{d}}$ in eq. 3.8b is useful to account for changes in the total mass of $\rho_{d}$, which can compromise the accuracy of the method when $d$ is low, but doing so becomes unnecessary for large enough $d$. In contrast with the operator $\mathcal{H}_{x}$ in eq. 3.3a, the operator $(\mathcal{L}_{x}^{m_{d}})^{*}$ is not selfadjoint in $L^{2}(\mathbf{R})$, and therefore the associated stiffness matrix is not symmetric. In addition, the quadratic form $\left\langle{(\mathcal{L}_{x}^{m})^{*}\cdot,\cdot}\right\rangle$ is not necessarily negative for the usual $L^{2}(\mathbf{R})$ inner product, and indeed we observe numerically that the eigenvalue with smallest real part of the discrete operator is often negative, although small when $d$ is large enough. This is illustrated in fig. 3 for the same parameters as in the subsequent convergence study.

For the integration in time, we used either the RK45 method (using the solve_ivp method from the SciPy integrate module), or a linear semi-implicit method obtained by treating $m_{d}$ explicitly and the other terms implicitly at each time step. The former is most useful when an accurate time-dependent solution is required, while the latter enables the use of larger time steps and is therefore more convenient when only the steady-state solution is sought, as will be the case for the construction of bifurcation diagrams. Denoting the time step by $\Delta t$ and the Galerkin approximation of $\rho_{d}(n\,\Delta t)$ by $\rho_{d}^{n}$, the semi-implicit method is based on obtaining $\rho_{d}^{n+1}$ by solving:

[TABLE]

Convergence study

The analysis of the Hermite spectral method for general types of McKean–Vlasov equations will not be presented here. For the purposes of this work, it will be sufficient to present a detailed numerical study of the convergence of the method. To study empirically the validity of the Galerkin method (3.8) and of the associated time-stepping scheme (3.9), we compare our method with the positivity preserving, entropy decreasing finite volume method proposed in [8] for nonlinear, nonlocal gradient PDEs111We reiterate the fact that one of the main advantages of our numerical method is that it does not require that the PDE has a gradient flow structure. . The parameters used here are $\beta=3$, $\theta=1$, and the initial condition was the Gaussian $\mathcal{N}(10^{-1},1)$. The same time points were used for the finite volume method and semi-implicit Galerkin method (with a mean time step of approximately 0.002), and for RK45 the absolute and relative tolerances were both set to $10^{-11}$. For the finite volume method, 600 equidistant mesh points were used between $x=-6$ and $x=6$.

Figure 3 presents the $L^{\infty}(0,T;L^{1}(\mathbf{R}))$ norm of the errors associated with the solutions obtained, for values of $d$, the degree of Hermite polynomials used, ranging from $10$ to $80$. A very accurate solution, obtained by using our spectral method with $d=120$, was employed for the calculation of the errors. We observe that, as $d$ increases initially, the solutions obtained using the semi-implicit (3.9) and the RK45 methods are indistinguishable and converge exponentially fast. From $d\approx 40$, the accuracy of the semi-implicit method no longer improves, indicating that the error introduced by the time-stepping scheme dominates from that point on. From $d\approx 50$, the Galerkin/RK45 approximation becomes more precise than the finite volume method. We therefore conclude that an accuracy as good as that obtained using the finite volume scheme can be reached with roughly ten times fewer unknowns using the spectral discretization (3.8). Our spectral method also enjoys a low computational cost: it ran in only about a minute with a Intel i7-3770 processor, even for a value of $d$ as high as 80, whereas the finite volume simulation took over an hour.

Figure 4 presents snapshots of the solutions at different times. We observe that, although the number of Hermite functions employed in the expansion is relatively low (=25), the solutions are in extremely good agreement.

In the simulations presented in this section, the scaling factor was set to $\sigma^{2}=\frac{1}{10}$. As discussed in appendix A, choosing this factor appropriately can significantly improve the accuracy of the method. In particular, given that the solution to eq. 3.7 decreases rapidly as $\left|x\right|\to\infty$, $\sigma$ should decrease with $d$, with the optimal scaling being $\sigma\propto\sqrt{d}$, as demonstrated in [46]. For convergence studies, however, it is convenient to use a fixed $\sigma$, first because this is assumed by most convergence results (such as theorem 3.1) and, second, because this simplifies the calculation of the matrices involved in the Galerkin formulation (only the last row and the last column have to be calculated upon incrementing $d$).

Remark 3.5** (Computational considerations).**

Discretizing the operators appearing in the Galerkin approximations (3.3a) and (3.8) requires the calculation of multiple matrices corresponding to operators of the type $\Pi_{d}\,\left(f\,\partial_{x}\right)\Pi_{d}$, where $f$ is a polynomial and $\Pi_{d}$ is the $L^{2}(\mathbf{R};\operatorname{e}^{-V_{q}})$ projection operator onto $\mathbf{P}(d)$. These calculations can be carried out by noticing that

[TABLE]

The matrix representation of the first operator on the right-hand side, in a basis of Hermite polynomials, can be obtained from the Hermite transform of $f$. The matrix representation of the second operator, on the other hand, is a matrix with zero entries everywhere except on the $m$-th superdiagonal, in view of the recursion relation (A.2).

3.3 Linear Fokker–Planck equation with colored noise

In this section, we turn our attention to the case of Gaussian or non-Gaussian colored noise given in terms of overdamped Langevin dynamics. The case of harmonic noise can be treated in a similar fashion, and for conciseness we do not present the associated Galerkin formulation explicitly here. We start by considering the linear (without the interaction term) Fokker–Planck equation with colored noise:

[TABLE]

We recall that $\varepsilon^{2}$ controls the correlation time of the colored noise and $\zeta$ is a parameter such that the white noise limit is recovered (with inverse temperature $\beta$) when $\varepsilon\to 0$. We include $\varepsilon$ in eq. 3.10 because, although we do not consider the white noise limit in this section, large values of $\varepsilon$ are in general more difficult to tackle numerically, and it will be therefore convenient to use smaller correlation times in the numerical experiments below. The problem is now two-dimensional and the operator on the right-hand side of eq. 3.10 is no longer elliptic. In contrast with the white noise case, there does not exist an explicit formula for the steady-state solution for eq. 3.10.

The procedure for obtaining a Galerkin formulation is the same as in Subsection 3.1, except that we now use tensorized Hermite polynomials/functions. To retain some generality, we will consider that the Galerkin approximation space is of the form $S_{d}=\operatorname{e}^{-U(x,\eta)/2}\,\operatorname{e}^{-V_{q}(x,\eta)/2}\,\mathbf{P}(\mathcal{I}_{d})$ for some function $U:\mathbf{R}^{2}\mapsto\mathbf{R}$, a nondegenerate quadratic potential $V_{q}$ to be determined, and where $\mathbf{P}(\mathcal{I}_{d}):=\operatorname{span}\left\{x^{\alpha_{x}}\,\eta^{\alpha_{\eta}}:(\alpha_{x},\alpha_{\eta})\in\mathcal{I}_{d}\right\}$ for some index set $\mathcal{I}_{d}\subset\mathbf{N}^{2}$ that grows with $d\in\mathbf{N}$. Compared to the one-dimensional case, there are now two scaling parameters, $V_{q}:=x^{2}/2\sigma_{x}^{2}+\eta^{2}/2\sigma_{\eta}^{2}$. The Galerkin approximation we propose consists in finding $\rho_{d}\in S_{d}$ such that

[TABLE]

with appropriate initial conditions. The choice of the weight $\operatorname{e}^{U}$ in the inner products of eq. 3.11 is motivated by the fact that differential operators admit sparse representations in the Hermite-type basis naturally associated with $S_{d}$, and we note that $\operatorname{e}^{-U(x,\eta)/2}\,\operatorname{e}^{-V_{q}(x,\eta)/2}\,\mathbf{P}(\mathbf{N}^{2})$, where $\mathbf{P}(\mathbf{N}^{2})$ is the space of polynomials in two dimensions, is dense in $L^{2}(\mathbf{R}^{2};\operatorname{e}^{U})$. In practice, we obtain $\rho_{d}$ as $\operatorname{e}^{-U(x,\eta)/2}\,\operatorname{e}^{-V_{q}(x,\eta)/2}\,v_{d}$, where $v_{d}$ is obtained by solving

[TABLE]

where, for a test function $\varphi$, $\mathcal{H}_{\varepsilon}\varphi:=(\operatorname{e}^{U/2}\,\operatorname{e}^{V_{q}/2})\mathcal{L}_{\varepsilon}^{*}\,(\operatorname{e}^{-U/2}\,\operatorname{e}^{-V_{q}/2}\,\varphi)$, and the basis functions used for eq. 3.12 are Hermite polynomials orthonormal with respect to the Gaussian weight $\operatorname{e}^{-V_{q}}$. Regarding the index set, several choices are possible, with the simplest ones being the triangle $\{\alpha\in\mathbf{N}^{2}:\left|\alpha\right|_{1}\leq d\}$ and the square $\{\alpha\in\mathbf{N}^{2}:\left|\alpha\right|_{\infty}\leq d\}$, see figs. 15 and 16 below. We demonstrate in Section 4 that, in order to study the limit $\varepsilon\to 0$, a rectangle-shaped index set is usually the only suitable choice. When studying the behavior as $d$ increases, however, we observed spectral convergence irrespectively of the index set utilized.

Clearly, it is necessary that $\rho\in L^{2}(\mathbf{R}^{2};\operatorname{e}^{U})$ for the Galerkin discretization (3.11) to produce good results. Since the $1/\varepsilon^{2}$ part of the operator on the right-hand side of eq. 3.10, $\mathcal{L}_{0}^{*}\cdot=\mathchoice{\frac{\partial{}}{\partial\eta}}{\partial_{\eta}}{}{}(V_{\eta}^{\prime}(\eta)\cdot+\,\mathchoice{\frac{\partial{}}{\partial\eta}}{\partial_{\eta}}{}{}\cdot)$, is selfadjoint in $L^{2}(\mathbf{R};\operatorname{e}^{-U_{x}(x)/2-V_{\eta}(\eta)/2})$ for any choice of $U_{x}$, it is natural to choose $\operatorname{e}^{-U(x,\eta)/2}=\operatorname{e}^{-U_{x}(x)/2-V_{\eta}(\eta)/2}$ for some one-dimensional potential $U_{x}$. This guarantees that the matrix representation of $\mathcal{L}_{0}^{*}$ is symmetric and negative semi-definite, but this is not a requirement.

The performance of the Galerkin approximation (3.11) is investigated through numerical experiments in appendix C. An asymptotic analysis of the numerical method in the limit as $\varepsilon\to 0$ is presented in [22], which is a longer version of this paper.

3.4 McKean–Vlasov equation with colored noise

We consider now the nonlinear McKean–Vlasov initial value problem with OU noise: recalling that $\zeta=1/\sqrt{2}$ in this case,

[TABLE]

for some initial distribution $\rho_{0}(x,\eta)$ such that the noise is not necessarily started at stationarity. The method that we use in this case, which applies mutatis mutandis to the other noise models, is the same as in eq. 3.11, with the addition of the interaction term, and we use the same time-stepping schemes as in Subsection 3.1.

Numerical experiments, testing the convergence of the method for the two time-stepping schemes, are presented in appendix C.

3.5 Monte Carlo simulations

We will compare the bifurcation diagrams obtained using the spectral method described above to those obtained by direct MC simulations of the system of interacting particles (2.1). We use the Euler–Maruyama method:

[TABLE]

where $\eta_{k}^{i}$ is the appropriate projection of the stochastic process $\boldsymbol{\mathbf{Y}}_{t}$. In the case of Gaussian noise, this is discretized as follows

[TABLE]

where $\xi\sim N(0,1)$, and $X_{k}$, $\boldsymbol{\mathbf{Y}}_{k}$ and $\eta_{k}$ are the approximations to $X(k\Delta t),\,\boldsymbol{\mathbf{Y}}(k\Delta t)$ and $\eta(k\Delta t)$, respectively. The time step used was always $\mathcal{O}(\varepsilon^{2})$, to ensure the accurate solution of the equation. This scheme has weak order of convergence one, see [24, 25], and we find that we capture the correct behavior as long as the time step is sufficiently small.

4 Asymptotic analysis for the Galerkin formulation

In Section 5, we will construct bifurcation diagrams of $m$ as a function of $\beta$ for different values of $\varepsilon$, and we will verify that the bifurcation diagram of the white noise case is recovered when $\varepsilon\to 0$. Since the spectral method presented in Section 3 will be used to that purpose, it is useful to study the behavior of the solution to the Galerkin formulation (3.11) in the limit $\varepsilon\to 0$, which is the purpose of this section. We will then confirm numerically the rates of convergence to the white noise limit presented in Subsection 2.2, i.e. $\mathcal{O}(\varepsilon^{2})$ for Ornstein–Uhlenbeck noise and $\mathcal{O}(\varepsilon^{4})$ for harmonic noise.

For simplicity, we confine ourselves for the analysis to the case where the noise process is one-dimensional and the weight function $\operatorname{e}^{-U(x,\eta)/2}$ can be decomposed as $\operatorname{e}^{-U(x,\eta)/2}=\operatorname{e}^{-U_{x}(x)/2}\,\operatorname{e}^{-V_{\eta}(\eta)/2}$. As before, $\hat{\Pi}_{d}$ denotes the $L^{2}(\mathbf{R}^{2};\operatorname{e}^{U})$ projection operator on the space of Hermite functions (with appropriate scalings). Decomposing the operator $(\hat{\Pi}_{d}\,\mathcal{L}_{\varepsilon}^{*}\,\hat{\Pi}_{d})$ in eq. 3.11 in powers of $\varepsilon$, we obtain the equation

[TABLE]

As a consequence of the choice of $\operatorname{e}^{U}$, the largest (sign included) eigenvalue of $\hat{\mathcal{L}}_{0}$ is a nonpositive, nondecreasing function of $d$. Since we cannot expect the leading order term of the discrete generator to have an eigenvalue exactly equal to 0, we look for a solution of the form $\rho_{d}=\operatorname{e}^{-\left|\lambda_{0,d}\right|\,t/\varepsilon^{2}}(\varrho_{0}+\varepsilon\,\varrho_{1}+\varepsilon^{2}\,\varrho_{2}+\dotsb)$, where $\lambda_{0,d}$ is the largest eigenvalue of $\hat{\mathcal{L}}_{0}$. Denoting by $\operatorname{id}$ the identity operator and gathering equal powers of $\varepsilon$, we obtain the equations:

[TABLE]

4.1 Suitable index sets

Let $H^{x}_{i}$ and $H^{\eta}_{j}$ denote the (possibly rescaled) Hermite functions in the $x$ and $\eta$ directions, respectively. Let also $\mathcal{I}_{i,\eta}$ be a slice of the index set, $\mathcal{I}_{i,\eta}:=\{j:(i,j)\in\mathcal{I}_{d}\}$, and $\Pi_{x}\mathcal{I}_{d}$ be the projected index set given by $\{i:(\exists j\in\mathbf{N})[(i,j)\in\mathcal{I}_{d}]\}$. Expanding $\varrho_{0}$ in terms of the basis functions used for the Galerkin discretization in the first equation, we obtain

[TABLE]

where $L_{jk}:=\left\langle{\mathcal{L}_{0}^{*}(\operatorname{e}^{-V_{\eta}/2}\,H^{\eta}_{j}),\operatorname{e}^{-V_{\eta}/2}\,H^{\eta}_{k}}\right\rangle_{\operatorname{e}^{V_{\eta}}}$. Since $\{\operatorname{e}^{-U_{x}/2}\,H^{x}_{i}\}_{i\in\Pi_{x}\mathcal{I}_{d}}$ are linearly independent, this implies that

[TABLE]

implying that, for all $\forall i\in\Pi_{x}\mathcal{I}_{d}$, the vector of coefficients ${(c_{ij})}_{j\in\mathcal{I}_{i,\eta}}$ is in the kernel of the matrix ${(L_{jk}+\left|\lambda_{0,d}\right|)}_{j,k\in\mathcal{I}_{i,\eta}}$. Therefore, if for some $i$ this matrix has full rank, then the corresponding Hermite coefficients must be [math]. This is of particular relevance when the eigenfunction in the kernel of $\mathcal{L}_{0}^{*}$ cannot be exactly represented in terms of a finite number of the approximating basis functions, as is the case with the noise processes B and NS considered in Section 2.

Consider for example the case of a triangular index set. In this case the matrix $M_{i}:={(L_{jk}+\left|\lambda_{0,d}\right|)}_{j,k\in\mathcal{I}_{i,\eta}}$ is of shape $d+1-i\times d+1-i$: incrementing $i$ by one corresponds to removing the last line last and last column of the matrix. Consequently, the maximal eigenvalue of $M_{i}$ (sign included), which we denote by $\lambda_{\max}(i)$, is a nonincreasing function of $i$ such that $\lambda_{\max}(0)=0$. In the case of noise model B, it holds that $\lambda_{\max}(i)<0$ for $i>0$ or $i>1$, depending on whether $d$ is even or odd. Consequently, in this case eq. 4.3 implies that $c_{ij}=0$ for these values of $i$ and all $j\in\mathcal{I}_{i,\eta}$.

From this we conclude that, in order to capture the correct solution as $\varepsilon\to 0$, it is necessary to choose a rectangularly shaped index set, which is consistent with the fact that, in the limit $\varepsilon\to 0$, the solution can be expressed as a tensor product $\rho_{\infty}(x,\eta)=\rho^{x}(x)\,\rho^{\eta}(\eta)$.

To illustrate the point made in the previous paragraph, we present side by side in fig. 5 the results of numerical experiments performed using either a triangular index set or a square index set, for the parameters $\varepsilon=0.01$, $d=20$, $\beta=15$, $\theta=0$, $\sigma_{x}^{2}=\sigma_{\eta}^{2}=1/15$, $\operatorname{e}^{-U_{x}(x)}=1$. While the probability density obtained using a square index set is close to the exact solution and clearly exhibits four local maxima, the solution obtained using a triangle index set is concentrated around $x=0$, and all the associated Hermite coefficients $c_{ij}$ with $i>0$ are very close to zero.

4.2 Effective drift and diffusion coefficients

We assume from now on that the index-set has a rectangular shape, $\mathcal{I}_{d}=\{0,1,\dotsc,d_{x}\}\times\{0,1,\dotsc,d_{\eta}\}$. At the continuous, infinite dimensional level, the absence of an effective drift term for the models of the noise B and NS as $\varepsilon\to 0$ is ensured by the centering condition (2.3). Any deviation from zero would lead to an effective drift term, scaling as $1/\varepsilon$ and proportional to

[TABLE]

At the finite-dimensional, numerical level, a parasitic effective drift can arise even when $V_{\eta}$ satisfies (2.3), as we demonstrate below. This can occur when $V_{\eta}$ is not an even function, and it is especially critical when the number of basis functions used to approximate the solution in the $\eta$ direction is relatively low, leading to a nonzero first moment of the approximate equilibrium probability density of the noise process. In these cases, it is useful to introduce an artificial drift term $\mu_{d}/\varepsilon$ in Galerkin formulation, for some constant $\mu_{d}$ to be determined. To formulate result 4.1 below, let $\varphi_{0,d}$ denote the (assumed unique up to a sign) normalized one-dimensional eigenfunction associated with $\lambda_{0,d}$,

[TABLE]

Result 4.1**.**

Let $\mu_{d}$ be defined by

[TABLE]

Then, when $\varepsilon\ll 1$, the solution $\rho_{d}$ of

[TABLE]

with $\mathcal{L}_{\varepsilon}^{*}$ as in eq. 3.10, can be approximated by $\rho_{d}^{x}(x,t)\,\varphi_{0,d}(\eta)$, where $\rho_{d}^{x}$ satisfies

[TABLE]

Here the effective diffusion $A_{d}$ is equal to

[TABLE]

with $b_{\eta}:=\left(\mu_{d}+\sqrt{2\,\beta^{-1}}\,\zeta\,\eta\right)$ and where $\mathcal{L}_{0}$ is the formal $L^{2}$ adjoint of $\mathcal{L}_{0}^{*}$.

Proof.

The argument below is formal, but it can be turned into a rigorous proof using standard methods in multiscale analysis; see e.g. [39]. Expanding the solution in powers of $\varepsilon$ and gathering terms multiplying equal powers of $\varepsilon$, a system of equations similar to eqs. 4.2a, 4.2b and 4.2c can be obtained, differing only by the presence of the corrective drift term next to $\hat{\mathcal{L}}_{1}$. The solvability condition for the first equation implies that $\varrho_{0}=\rho^{x}_{d}(x,t)\,\varphi_{0,d}(\eta)$. For the second equation, we see from the definition of $\mu_{d}$ and using the symmetry of $\hat{\mathcal{L}}_{0}$ in $L^{2}(\mathbf{R}^{2};\operatorname{e}^{V_{\eta}(\eta)+U_{x}(x)})$, that the Fredholm solvability condition is automatically satisfied, which enables solving for $\varrho_{1}$:

[TABLE]

Writing out the solvability condition for the third equation, we obtain the effective equation for $\rho^{x}_{d}$:

[TABLE]

which after expansion of the terms is eq. 4.5. ∎

In the examples we considered in Section 3, both $\varphi_{0,d}$ and $V_{\eta}$ were even functions, so the corrective drift term defined by eq. 4.4 was zero. This is not the case for the noise model NS, where the noise is confined by the uneven potential (2.4). To illustrate the importance of including the corrective drift term (4.4), we present simulation results with and without it in the case of the noise model NS. We consider the following parameters for the equation: $\varepsilon=2^{-5}$, $\theta=0$, $\beta=1$, $V(x)=x^{4}/4-x^{2}/2$, and for the numerics we choose $d=20$, $\sigma_{x}^{2}=\sigma_{\eta}^{2}=0.1$, $\operatorname{e}^{-U_{x}(x)}=1$. The solutions obtained are presented in figs. 6 and 6.

We note however that, for fixed $\varepsilon$, $\mu_{d}\to 0$ as $d\to\infty$, by convergence of $\varphi_{0,d}$ to $C\,\operatorname{e}^{-V_{\eta}}$. This is illustrated in fig. 7, where the same parameters as those used in fig. 6 were used. While the value of the bias decreases exponentially with the degree of approximation, we have found through numerical experiments that, when $\varepsilon$ is of the order of $0.01$, failure to account for the parasitic drift is very detrimental to the accuracy of the solution for $d$ as high as 30.

4.3 Numerical verification of the rates of convergence

In this section, we verify numerically the rate of convergence of the stationary solution to the Galerkin formulation eq. 3.11 in the limit $\varepsilon\to 0$. We recall that this Galerkin formulation is associated to the linear Fokker–Planck equation (3.11). In addition to verifying the rates of convergence, examining the limit $\varepsilon\to 0$ numerically will enable us to gain insight into the accuracy of the asymptotic expansions for moderate values of $\varepsilon$.

One-dimensional Ornstein–Uhlenbeck noise

For this test we use the same parameters as in the convergence study for the bistable potential in Subsection 3.3. We verify the accuracy of the asymptotic expansion up to order $\varepsilon^{2}$ of the full solution in the $x$–$\eta$ plane, given in [49, Chapter 3], by comparing it to numerical results obtained using the spectral method introduced in Section 3. The convergence is presented in fig. 8. We notice that, even for the smallest value of $\varepsilon$ considered ($2^{-6}$), the norm of the difference between the asymptotic and spectral solutions appears to be roughly constant for $d\geq 20$. This is because, beyond this point, the spectral method is more accurate than the asymptotic expansion.

Harmonic noise

The second case we consider is that of harmonic noise, for which the order of the first nontrival correction in the expansion of the solution is $\mathcal{O}(\varepsilon^{4})$; see [49]. We confirmed this numerically for $V(x)=x^{4}/4-x^{2}/2$ and $\beta=5$ using 50 basis functions in each direction, with scaling factors $\sigma_{x}^{2}=1/30,\sigma_{p}^{2}=\sigma_{q}^{2}=1$. The results are illustrated in fig. 9.

5 Results: effect of colored noise on bifurcations

In this section we present the bifurcation diagrams corresponding to the four models of the noise introduced in Section 2. We begin with the case of Gaussian noise, and later move to the case of non-Gaussian noise.

5.1 Construction of the bifurcation diagrams for the mean field equation

We constructed the bifurcation diagrams using three different approaches:

Monte Carlo simulations

We solved the system of interacting particles (2.1) with a sufficiently large number of particles, and we approximated the first moment by ergodic average over an interval $(T,T+\Delta T)$, where $T$ is sufficiently large to guarantee that the system has reached its stationary state and $\Delta T$ is sufficiently large to ensure that the ergodic averages are accurate. By applying this procedure for a range of inverse temperatures, $\beta=0.1,0.15,0.2,\dotsc,10$, we obtained the desired bifurcation diagram.

Perturbation expansions

This approach, which we already outlined in Section 2, relies on the fact that the self-consistency map can be approximated as $R(m,\beta)\approx R_{0}(m,\beta)+\varepsilon^{\delta}R_{\delta}(m,\beta)$, with good accuracy when $\varepsilon\ll 1$. Here we used the same notation as in Section 2, and in particular $\delta$ denotes the order of the first nontrivial correction in eq. 2.8. Using arclength continuation222We do this using the Moore–Penrose quasi-arclength continuation algorithm. The rigorous mathematical construction of the arclength continuation methodology can be found, e.g., in [30] and [3]. Some useful practical aspects of implementing arclength continuation are also given in [10]. See also [20].

for the resulting approximate self-consistency equation, $m=R_{0}(m,\beta)+\varepsilon^{\delta}R_{\delta}(m,\beta)$, we can plot the first moment $m$ as a function of $\beta$ for a fixed value of $\varepsilon$. We note that, in view of the typical shape of the self-consistency map, depicted in a particular case in fig. 2, a standard root finding algorithm can be employed to initiate the arclength continuation at some initial inverse temperature $\beta_{0}$.

The spectral method

Finally, we employed the Galerkin method presented in Subsection 3.3. We considered two different methodologies: on the one hand, by calculating numerically an approximation $\rho_{d,\infty}(x,\eta;\beta,m)$ of the steady-state solution of the linear Fokker–Planck equation (2.6a) with fixed $m$ and $\beta$, we approximated the self-consistency map as $R(m,\beta)\approx\int_{\mathbf{R}}\int_{\mathbf{R}^{n}}x\,\rho_{d,\infty}(x,\eta;\beta,m)\,\mathrm{d}x\,\mathrm{d}\boldsymbol{\mathbf{y}}$, after which a bifurcation diagram can be constructed by using the same method as in the previous paragraph. Each evaluation of the self-consistency map requires the computation of the eigenvector associated with the eigenvalue of smallest magnitude of the discretized operator, which can be performed efficiently for sufficiently small systems using the SciPy toolbox. On the other hand, the time-dependent (nonlinear) McKean–Vlasov equation can be integrated directly using our spectral method. Since only the final solution is of interest to us, the semi-implicit time-stepping scheme (3.9) can be used with a large time step, which enables a quick and accurate approximation of the steady-state solutions. While both methodologies work well in the two-dimensional case, in three dimensions (harmonic noise) solving the McKean–Vlasov equation directly proved more efficient, so this is the approach we employed for all the tests presented in this section.

5.2 Gaussian case

The one-dimensional Ornstein–Uhlenbeck noise provides an ideal testbed for the three methods we use to construct bifurcation diagrams. Figure 10 below plots the bifurcation diagram of the first moment $m$ as a function of $\beta$ for $\varepsilon=0.1,0.2,\dotsc,0.5$. Three different initial conditions ($\boldsymbol{\mathbf{X}}_{0}\sim N(0,0.1),\,\boldsymbol{\mathbf{X}}_{0}\sim N(0.1,0.1)$, and $\boldsymbol{\mathbf{X}}_{0}\sim N(-0.1,0.1)$) were used for the MC simulations. Although we observe that the results of MC simulations tend to be less precise around the bifurcation point, the agreement between the three methods overall is excellent for $\varepsilon=0.1,0.2$. For the other values of $\varepsilon$, while the results of MC simulations and of our spectral method continue to agree, those obtained from the asymptotic expansion are significantly less accurate, which is consistent with the observations presented in fig. 8.

The case of harmonic noise, corresponding to a three-dimensional McKean–Vlasov equation, is more challenging to tackle using our spectral method. When using 40 basis functions in each direction, the CPU time required to construct the full bifurcation diagram was of the order of a week. As a consequence of the lower number of basis functions used in this case, we observe a small discrepancy between the results of the spectral method and those of MC simulations for large $\beta$ in the case $\varepsilon=0.4$. Nevertheless, as can be seen in fig. 11, for small $\varepsilon$ the overall agreement between the three methods is excellent. We note in particular that, as suggested by the asymptotic expansions, the use of harmonic noise produces results much closer to the white noise limit than scalar OU noise.

5.3 Non-Gaussian noise

For the non-Gaussian noise processes we consider, the $x^{4}$ asymptotic growth of the confining potentials in both directions causes the McKean–Vlasov equation to be stiffer than in the cases of OU and harmonic noise, especially for large values of $\varepsilon$. Consequently, we were not able to consider as wide a range of $\varepsilon$ as in the previous subsection using the spectral method. Since, on the other hand, MC simulations become overly computationally expensive for small $\varepsilon$, the comparisons in this section comprise only results obtained using our spectral method and asymptotic expansions. Results of simulations for the bistable noise (model B) are presented in fig. 12, in which a very good agreement can be observed.

For nonsymmetric noise (model NS), the two branches in the bifurcation diagram are separate, as illustrated in fig. 13. Here too, the agreement between the spectral method and the asymptotic expansion is excellent. In contrast with the other models considered, the first nonzero term in the asymptotic expansion is of order $\varepsilon$, which is reflected by the manifestly higher sensitivity to the correlation time of the noise. In the right panel of fig. 13, we present the graph of $R_{0}(m;\beta)+\varepsilon R_{1}(m;\beta)$ for a value of $\beta$ close to the point at which new branches (one stable and one unstable) emerge.

5.4 Dependence of the critical temperature on $\varepsilon$

For the noise models OU, H and B, the effect of colored noise on the dynamics is a shift of the critical temperature: the pitchfork bifurcation occurs for smaller values of $\beta$ (i.e., larger temperatures) as the correlation time increases. In order to further investigate the effect of the correlation time on the long time behavior of the system of interacting particles, we will compute the critical temperature as a function of $\varepsilon$ based on the asymptotic expansions and compare with the results of spectral and MC simulations, see fig. 14. Rather than finding the critical inverse temperature $\beta_{C}$ for a range of values of $\varepsilon$ (and for a fixed $\theta$), it is convenient to fix $\beta_{C}$ and find the corresponding $\varepsilon$, satisfying

[TABLE]

which is merely a polynomial equation in $\varepsilon$, the coefficient of which can be calculated by numerical differentiation. With this procedure, the dependence of the critical $\beta$ upon $\varepsilon$ can be calculated on a fine mesh. In the case of OU noise, for example, both coefficients on the left-hand side of eq. 5.1 are positive, implying that the equation has a solution (in fact, two, but one of them negative) only if $\beta_{C}$ is lower than the inverse critical temperature in the white noise case.

Of the three methods employed in fig. 14, the approach based on the asymptotic expansions has the lowest computational cost: calculating all the solid curves took only about a couple of minutes on a personal computer with an Intel i7-3770 processor. The data points associated with the spectral method and the MC simulations were obtained from the bifurcation diagrams presented above.

6 Conclusions

In this paper, we introduced a robust spectral method for the numerical solution of linear and nonlinear, local and nonlocal Fokker–Planck-type PDEs that does not require that the PDE is a gradient flow. We then used our method to construct the bifurcation diagram for the stationary solutions of the mean field limit of a system of weakly interacting particles driven by colored noise.

To verify our results, we also constructed the bifurcation diagrams by using two other independent approaches, namely by MC simulation of the $N$-particle system and by using explicit asymptotic expansions with respect to correlation time of the noise. In the small correlation time regime, we observed a very good agreement between all three methods. For larger values of the correlation time, the asymptotic expansions become inaccurate, but the results obtained via the spectral method and MC simulations continue to be in good agreement.

It appeared from our study that, unless the potential in which the noise process is confined is asymmetric, the correlation structure of the noise does not influence the topology of the bifurcation diagram: the mean-zero steady-state solution, which is stable for sufficiently large temperatures, becomes unstable as the temperature decreases below a critical value, at which point two new stable branches emerge, in the same manner as reported in [9, 44]. The correlation structure does, however, influence the temperature at which bifurcation occurs, and in general this temperature increases as the correlation time of the noise increases. In the presence of an asymmetry in the confining potential of the noise, on the other hand, the two stable branches in the bifurcation diagram are separate, indicating that the system always reaches the same equilibrium upon slowly decreasing the temperature. This behavior is similar to what has been observed previously in the white noise case when a tilt is introduced in the confining potential $V(\cdot)$, see [20, 21].

Several problems remain open for future work. On the theoretical front, we believe that the analysis we presented in Subsections 3.1 and B for the linear Fokker–Planck equation can be extended to both the linear Fokker–Planck equation with colored noise and the nonlinear McKean–Vlasov equation. Another direction for future research could be the rigorous study of bifurcations and, more specifically, of fluctuations and critical slowing down near the bifurcation point. On the modeling front, it would be interesting to consider more general evolution equations for the interacting particles, such as the generalized Langevin equation, and also to study systems of interacting particles subject to colored noise that is multiplicative.

Appendix A Hermite polynomials and Hermite functions

In one dimension, the orthonormal Hermite polynomials can be defined by:

[TABLE]

They form a complete orthonormal basis of the weighted space $L^{2}(\mathbf{R};g)$, where $g:=\operatorname{e}^{-x^{2}/2}/\sqrt{2\,\pi}$, and they satisfy the following recursion relations (eq. A.3 can be obtained by combining eqs. A.1 and A.2):

[TABLE]

From eqs. A.1 and A.2, we see that Hermite polynomials are the eigenfunctions of a second-order operator:

[TABLE]

which is essential to proving approximation results. We note that the eigenvalues grow linearly, which explains the square root in the rate of convergence in the results presented below. Since Hermite polynomials constitute an orthonormal basis of $L^{2}(\mathbf{R};g)$, any function $u$ in that space can be expanded as a series of Hermite polynomials, and

[TABLE]

We will call Hermite transform the operator:

[TABLE]

Denoting by $\Pi_{d}$ the $L^{2}(\mathbf{R};g)$ projection operator on $\operatorname{span}\{H_{0},H_{1},\dotsc,H_{d}\}$, the following theorem follows from (A.2); see [43] for details.

Theorem A.1** (Approximation by Hermite polynomials).**

For any $u\in H^{m}(\mathbf{R};g)$ with $0\leq m\leq d+1$, the following inequality holds:

[TABLE]

The Hermite polynomials, as defined above, are suitable for the approximation of functions with respect to the norm of $L^{2}(\mathbf{R};g)$, which assigns a significant weight only to the region around $x=0$. For the approximation with respect to the flat $L^{2}(\mathbf{R})$ norm, or with respect to other norms that penalize growth as $x\to\infty$, such as the weighted $L^{2}(\mathbf{R};e^{x^{2}/2})$ norm, one can use the basis functions $(\operatorname{e}^{-U/2}\,H_{i})_{i=0}^{\infty}$ for some function $U$, which constitute an orthonormal basis of $L^{2}(\mathbf{R}^{n};\operatorname{e}^{U}g)$. For $u\in L^{2}(\mathbf{R};\operatorname{e}^{U}g)$, we define the generalized Hermite transform associated with the factor $\operatorname{e}^{-U/2}$, which we denote by $\mathcal{T}_{U}$:

[TABLE]

Note that $\mathcal{T}_{U}(u)=\mathcal{T}(\operatorname{e}^{U/2}u)$, i.e. $T_{U}(u)$ is the usual Hermite transform of $\operatorname{e}^{U/2}u$. This formalism enables us to treat in a unified manner the case of Hermite polynomials ($U=0$), of Hermite functions ($\operatorname{e}^{-U/2}=\sqrt{g}$), as well as other useful cases. As an example of why such generality can be useful, it has been shown in [14] that the choice $\operatorname{e}^{-U/2}=g$ leads to basis functions, referred to as generalized Hermite functions in that paper, that can be used to design an efficient numerical method for the solution of the Kramers Fokker–Planck equation. By the property (A.4), we see that $(\operatorname{e}^{-U/2}H_{i})_{i=0}^{\infty}$ are the eigenfunctions of the operator:

[TABLE]

which is of Schrödinger type when $\operatorname{e}^{-U}=g$, see [39]. Introducing the notations $\Pi_{d}^{U}:=\operatorname{e}^{-U/2}\Pi_{d}\operatorname{e}^{U/2}$ and $\partial_{x}^{U}:=\partial_{x}+U^{\prime}/2=\operatorname{e}^{-U/2}\partial_{x}(\operatorname{e}^{U/2}\cdot)$, we have the following immediate corollary of theorem A.1.

Corollary A.2** (Approximation by generalized Hermite functions).**

For any $u$ such that $(\partial_{x}^{U})^{m}u\in L^{2}(\mathbf{R};\operatorname{e}^{U}g)$ with $0\leq m\leq d+1$,

[TABLE]

In the case $\operatorname{e}^{-U}=g$ (orthonormal Hermite functions in $L^{2}(\mathbf{R})$), one can prove a similar statement with the usual derivative instead of $\partial_{x}^{U}$ in the left-hand side, see [43, Theorem 7.14]. We note that theorems A.1 and A.2 can be extended to the multi-dimensional case, see e.g. [1].

In addition to the function $\operatorname{e}^{-U/2}$ multiplying the Hermite polynomials in the definition of basis functions, it is usual in numerical simulations to introduce a scaling factor, which can be chosen appropriately depending on how localized the function to be approximated is. We define $H^{\sigma}_{i}(x):=H_{i}(x/\sigma)$, and note that these polynomials form an orthonormal basis in $L^{2}(\mathbf{R};g_{\sigma})$, where $g_{\sigma}$ is the normal distribution with mean 0 and variance $\sigma^{2}$. Although we do not present them explicitly, approximation results similar to theorems A.1 and A.2 can be proved in the presence of this scaling factor, with the only difference being the presence of additional constant factors on the right-hand side; see, for example, [1].

The importance of choosing this scaling factor in a suitable manner was demonstrated rigorously in [46], where the author shows that, for Gaussian-type functions, the optimal value of $\sigma$ should depend on the number of basis function used with a dependence of the form $1/\sqrt{d}$. This can be justified intuitively by taking into account that, on the one hand, the behavior of Hermite functions as $x\to\infty$ is well-understood, with the final inflection point occurring at $x\propto\sqrt{d}\,\sigma$ and a rapid decrease to 0 beyond that point [31], and that, on the other hand, scaled Hermite functions are eigenfunctions of the Fourier transform operator, up to dilations/contractions. Since a contraction in real space results in a dilation in Fourier space, by choosing a scaling factor that decreases with $d$ one effectively favors exploration in Fourier space. In particular, choosing $\sigma\propto 1/\sqrt{d}$ leads to a situation where the position of the final inflection point of the Hermite function of highest order remains approximately constant (while the position of the final inflection point of its Fourier transform grows linearly), so it is crucial in that case to ensure that $\sqrt{d}\,\sigma$ is large enough to ensure that the basis functions cover the support of the target function.

In practice, calculating the Hermite transform numerically requires the introduction of a quadrature. To bound the associated error, results similar to theorems A.1 and A.2, with the projection operators replaced by interpolation operators, can be proved; see [43, Theorems 7.17, 7.18].

Appendix B Proof of theorem 3.1

Using the same notation as in appendix A, we let $\Pi_{d}$ be the $L^{2}(\mathbf{R};e^{-V_{q}})$ projection operator on $\mathbf{P}(d)$ and $\hat{\Pi}_{d}:=e^{-V_{q}/2}\,\Pi_{d}\,e^{V_{q}/2}$. The solution $u_{d}$ of eq. 3.3a satisfies $\partial_{t}u_{d}=\hat{\Pi}_{d}\,\mathcal{H}_{x}\,\hat{\Pi}_{d}\,u_{d}=:\mathcal{H}_{d}\,u_{d}$. Clearly, the operator $\mathcal{H}_{d}$ is selfadjoint on $\operatorname{e}^{-V_{q}/2}\mathbf{P}(d)$ with the $L^{2}(\mathbf{R})$ inner product, and it is also negative, by negativity of $\mathcal{H}_{x}$:

[TABLE]

To prove the convergence of $u_{d}$ when $d\to\infty$, we will rely on the following lemma.

Lemma B.1**.**

Let $\hat{\partial}_{x}:=\partial_{x}+x/2$, and assume that $\hat{\partial}_{x}^{n}u\in L^{2}(\mathbf{R})$ for $n=0,\dots,m$. Then for all natural numbers $m_{1},m_{2}$ such that $m_{1}+m_{2}\leq m$, it holds that $x^{m_{1}}\,u^{(m_{2})}\in L^{2}(\mathbf{R})$ and

[TABLE]

where $K_{1}(m),K_{2}(m)$ are positive constants depending only on $m$ and $\|\cdot\|$ is the usual $L^{2}(\mathbf{R})$ norm.

Proof.

We denote by $H^{m}(\mathbf{R};\operatorname{e}^{-x^{2}/2})$ the Sobolev space weighted by $\operatorname{e}^{-x^{2}/2}$,

[TABLE]

and by $\|\cdot\|_{m,\operatorname{e}^{-x^{2}/2}}$ the associated norm: $\|v\|_{m,\operatorname{e}^{-x^{2}/2}}^{2}=\sum_{i=0}^{m}\|v^{(i)}\|_{\operatorname{e}^{-x^{2}/2}}^{2}$. For the first inequality, we know from [43, Lemma B.6] that

[TABLE]

Applying this inequality repeatedly, we obtain

[TABLE]

for a constant $C(m)$ depending only on $m$. By definition, $\hat{\partial}_{x}u=\operatorname{e}^{-x^{2}/4}\partial_{x}(\operatorname{e}^{x^{2}/4}u)$, so the assumption implies that $\operatorname{e}^{x^{2}/4}u\in H^{m}(\mathbf{R};\operatorname{e}^{-x^{2}/2})$, from which we obtain using eq. B.3 that, for $0\leq m_{1}\leq m$,

[TABLE]

This proves the first inequality of eq. B.2 in the case $m_{2}=0$. We assume now that the statement is proved up to $m_{2}-1$, and we show that it is valid for $m_{2}$. Using the triangle inequality we obtain

[TABLE]

The derivatives in the first term are of order strictly lower than $m_{2}$, and therefore this term can be bounded by the induction assumption. The second term is bounded by applying the base case to $\hat{\partial}_{x}^{m_{2}}u$: introducing $v:=\hat{\partial}_{x}^{m_{2}}u$, we notice that $\hat{\partial}_{x}^{m-m_{2}}v=\hat{\partial}_{x}^{m}u\in L^{2}(\mathbf{R})$ by assumption, so we can apply the first inequality in eq. B.2, without any derivative of $v$ in the left-hand side, to deduce

[TABLE]

The second inequality in (B.2) then holds trivially by expanding $\hat{\partial}_{x}$ and applying a triangle inequality. ∎

With 3.1, we can show that the two norms in lemma B.1 can be bounded from above by the norm $\sqrt{\left\langle{{(-\mathcal{H}_{x}+1)}^{m}u,u}\right\rangle}$ for appropriate $m$.

Lemma B.2** (Bound by alternative norm).**

If 3.1 holds, then

[TABLE]

for any smooth $u$ for which the right-hand side is well-defined. Here $C$ is a positive constant that depends on $\beta$, $m$, and on the particular expression of the potential $W$ defined in 3.1.

Proof.

Below $C_{1}$ and $C_{2}$ denote the same constants as in 3.1. First we notice that, for any constant $K>1$,

[TABLE]

because $\mathcal{H}_{x}$ is a negative operator. Since $W$ is a polynomial, its derivatives grow asymptotically more slowly that $W$ itself, and so it is possible for any $\varepsilon>0$ to find $K\geq C_{2}$ large enough that

[TABLE]

For this proof to go through, it is in fact sufficient that this inequality be satisfied for $i=1,\dotsc,m$.

We decompose $-\mathcal{H}_{x}+K$ as $(-\beta^{-1}\mathchoice{\frac{\partial^{2}{}}{\partial x^{2}}}{\partial_{x}^{2}}{}{})+(W(x)+K)$. The two operators in this sum are positive because $K\geq C_{2}$ and by assumption $W(x)+C_{2}\geq C_{1}(1+\left|x\right|^{2})$. Expanding the inner product in the left-hand side of eq. B.4 and using integration by parts,

[TABLE]

where the remainder terms originate from the fact that the operators $\mathchoice{\frac{\partial{}}{\partial x}}{\partial_{x}}{}{}$ and $(W(x)+K)$ do not commute. By eq. B.5, these terms can be bounded for sufficiently large $K$ by half the leading term in eq. B.6. To conclude, we further expand this leading term:

[TABLE]

from which lemma B.1 allows us to conclude. ∎

Proof of theorem 3.1.

We assume for simplicity that $\sigma=1$, and we begin by splitting the error as $u_{d}-u=(u_{d}-\hat{\Pi}_{d}u)+(\hat{\Pi}_{d}u-u)=:e_{d}+\delta_{d}$. The first term is related to the so-called consistency error, and the second to the approximation error. We obtain from eqs. 3.2 and 3.3a

[TABLE]

Taking the inner product with $e_{d}$ and using (B.1), this implies

[TABLE]

where we used Young’s inequality. We see from this equation that $e_{d}$ can be controlled if one can bound the second inner product on the right-hand side. For this we use arguments similar to the ones employed in [1, 15]. Since $\mathcal{H}_{x}$ is negative and selfadjoint, we notice

[TABLE]

for $i=1,2,\dots$, which implies that the inner products $\left\langle{(-\mathcal{H}_{x})^{i}u,u}\right\rangle$ remain bounded for all positive times. We can now apply corollary A.2 to obtain, using lemmas B.1 and B.2 and 3.1,

[TABLE]

We note that when $V$ is quadratic, $k=1$ is a valid choice in 3.1, and the bound above can be obtained by simply expanding $u$ in terms of the eigenfunctions of $\mathcal{H}_{x}$, which in that case are just rescaled Hermite functions. Using Grönwall’s inequality, we finally obtain

[TABLE]

The first term, proportional to $\|e_{d}(0)\|^{2}$, depends only on the interpolation error of the initial condition, which is nonzero when using a Gauss–Hermite quadrature. It was proved that this error term also decreases spectrally, see e.g. [43, Theorems 7.17, 7.18], and in our case faster than the second error term. For the approximation error $\delta_{d}$, similar inequalities to the ones used above can be used to obtain a bound of the type (B.8), which leads to the conclusion. ∎

Remark B.1**.**

As mentioned in remark 3.2, theorem 3.1 is not optimal. It leaves open, in particular, the question of precisely how the error behaves as $t\to\infty$. In this remark, we give a partial answer to the question: we show how a bound on the stationary error can be obtained, under the assumption that the solution to the discretized equation is rescaled in time in such a way that the integral of $\rho_{d}:=\sqrt{\rho_{s}}\,u_{d}$ remains equal to $1$. With this rescaling and with the notation $S_{d}:=\hat{\Pi}_{d}(L^{2}(\mathbf{R}))$, the stationary solution of the discretized (in space) equation is

[TABLE]

Taking

[TABLE]

we deduce

[TABLE]

Since $\mathcal{H}_{x}$ is self adjoint in $L^{2}(\mathbf{R})$ and $\mathcal{H}_{x}\sqrt{\rho_{s}}=0$, it follows that

[TABLE]

An approximation argument similar to the one above can be employed to show that the right-hand side, and therefore also the left-hand side, decrease to zero as $d\to\infty$ faster than $d^{-n}$ for any $n>0$. To conclude, a Poincaré-type inequality can be invoked, which is justified because

[TABLE]

to obtain a bound of the type

[TABLE]

Appendix C Benchmark tests for the spectral numerical method

In this section, we investigate the performance of the spectral method through numerical experiments.

C.1 Linear Fokker–Planck equation with colored noise

We focus first on the Galerkin approximation (3.11). Here we consider only the cases where $V(\cdot)$ is a quadratic or a bistable potential and where the noise is described by an OU process, but results of additional numerical experiments, corresponding to harmonic noise and non-Gaussian noise, are presented in [49].

We start with the case $V(x)=x^{2}/2$, for which the exact solution to the Fokker–Planck equation (3.10) can be calculated explicitly by substitution of a Gaussian ansatz, see [39, Section 3.7]. We study the convergence of the steady-state solution, obtained by calculating the eigenfunction associated with the eigenvalue of lowest magnitude of $\hat{\Pi}_{d}\,\mathcal{L}_{\varepsilon}^{*}\,\hat{\Pi}_{d}$, where $\hat{\Pi}_{d}$ is the $L^{2}(\mathbf{R}^{2};\operatorname{e}^{U})$ projection operator on $S_{d}$, directly using the method eigs from the SciPy toolbox. The parameters used for this simulation are the following: $\beta=\varepsilon=1$, $\sigma_{x}^{2}=\frac{1}{10}$, $\sigma_{\eta}^{2}=1$, $\operatorname{e}^{-U(x,\eta)/2}=\operatorname{e}^{-V_{\eta}(\eta)/2}=\operatorname{e}^{-\eta^{2}/4}$. With these parameters, the steady-state solution to eq. 3.10 is equal to $\rho_{\infty}(x,\eta)=\operatorname{e}^{-2x^{2}+2x\eta-\eta^{2}}/\pi$, and clearly $\rho_{\infty}\in L^{2}(\mathbf{R}^{2};\operatorname{e}^{U})$. Figure 15 presents the steady-state solution, obtained using the spectral method with Hermite polynomials up to degree 100 ($d=100$) and a triangular index set, and fig. 15 presents the convergence of the method. Since the solution satisfies $\rho_{\infty}(x,\eta)=\rho_{\infty}(-x,-\eta)$, the Hermite coefficients corresponding to even values of $i+j$ are zero, where $i$ and $j$ are the indices in the $x$ and $\eta$ directions, respectively.

Now we consider that $V$ is the bistable potential $x^{4}/4-x^{2}/2$, which was solved numerically in [23] using generalized Hermite functions and a variation of the matrix continued fraction technique. For this case an explicit analytical solution is not available. The parameters we use are the following: $\beta=1$, $\varepsilon=\frac{1}{2}$, $\sigma_{x}^{2}=\frac{1}{20}$, $\sigma_{\eta}^{2}=1$. Through numerical exploration, we noticed that a good convergence could be obtained by using the multiplier function $\operatorname{e}^{-U(x,\eta)/2}=\operatorname{e}^{-\beta V(x)/2-\eta^{2}/4}$, rather than just $\operatorname{e}^{-\eta^{2}/4}$ in the previous paragraph. We note that this would have been the natural choice if the noise in the $x$ direction had been white noise. The solution obtained using a square-shaped index set and $d=100$, as well as the corresponding Hermite coefficients up to degree 10, is illustrated in fig. 16. We observe that the Hermite coefficients corresponding to the degree 0 in the $\eta$ direction (i.e. to the basis function $\operatorname{e}^{-\eta^{2}/2}$) are significantly larger than the other coefficients, which is consistent with the fact that, as $\varepsilon\to 0$, the steady-state solution approaches $\operatorname{e}^{-\beta V(x)}\,\operatorname{e}^{-\eta^{2}/2}$ (up to a constant factor). The associated convergence curves are presented in fig. 16.

We observe that the convergence is exponential, which is better than the rate of convergence predicted in theorem 3.1.

C.2 McKean–Vlasov equation with colored noise

We focus now on the Galerkin approximation to (3.13). When the potential $V(\cdot)$ is quadratic and the initial condition is Gaussian, it is well-known that the McKean–Vlasov equation has an explicit solution and that this solution is Gaussian. We assume that $V(x)=x^{2}/2$ and we rewrite eq. 3.13 in the formalism of [11], as

[TABLE]

where $\boldsymbol{\mathbf{x}}=(x,\eta)^{T}$ and

[TABLE]

Adapting [11, Proposition 2.3] to our case, we deduce that the solution is of the type

[TABLE]

where $\mu(t)$ and $\Sigma(t)$ are given by

[TABLE]

This solution can be obtained by introducing $g=-\ln\rho$, rewriting eq. 3.13 as an equation for $g$, and using a quadratic ansatz for $g$. The eigenvalue decomposition of $B+K$ is

[TABLE]

which enables the explicit calculation of the integral in the expression of $\Sigma(t)$. From eq. C.1 and the structure of $B$ and $K$, we notice that, as $t\to\infty$, $\mu\to 0$ and

[TABLE]

which coincides with the solution of the steady state linear Fokker–Planck equation corresponding to the McKean–Vlasov equation when $m$ is a parameter equal to 0. For this test case, we use the following parameters: $\beta=\theta=1$, $\varepsilon=1/2$, $\sigma_{x}^{2}=\sigma_{\eta}^{2}=1/5$, $\operatorname{e}^{-U(x,\eta)/2}=\operatorname{e}^{-V(x)/2}\,\operatorname{e}^{-V_{\eta}(\eta)/2}$. The initial condition is taken to be the Gaussian density $\mathcal{N}\left((1,1)^{T},I_{2\times 2}\right)$. The evolution of the probability density is illustrated in fig. 17, and the convergence of the method, in the $L^{\infty}(0,T;L^{1}(\mathbf{R}^{2}))$ norm, is illustrated in fig. 18.

Acknowledgments

The authors are grateful to J.A. Carrillo for making available the finite volume code employed in Subsection 3.1 for the comparison with our spectral method for the McKean–Vlasov dynamics. This work was supported by EPSRC through grants number EP/P031587/1, EP/L024926/1, EP/L020564/1 and EP/K034154/1. The work of SG is supported by the Leverhulme Trust through Early Career Fellowship ECF-2018-536.