Well-posedness and numerical schemes for one-dimensional McKean–Vlasov equations and interacting particle systems with discontinuous drift

In this paper, we first establish well-posedness results for one-dimensional McKean–Vlasov stochastic differential equations (SDEs) and related particle systems with a measure-dependent drift coefficient that is discontinuous in the spatial component, and a diffusion coefficient which is a Lipschitz function of the state only. We only require a fairly mild condition on the diffusion coefficient, namely to be non-zero in a point of discontinuity of the drift, while we need to impose certain structural assumptions on the measure-dependence of the drift. Second, we study Euler–Maruyama type schemes for the particle system to approximate the solution of the one-dimensional McKean–Vlasov SDE. Here, we will prove strong convergence results in terms of the number of time-steps and number of particles. Due to the discontinuity of the drift, the convergence analysis is non-standard and the usual strong convergence order 1/2 known for the Lipschitz case cannot be recovered for all presented schemes.


Introduction
In this article, we study the existence and uniqueness of strong solutions for classes of McKean-Vlasov SDEs, where the drift exhibits a discontinuity in the spatial component.We also provide time-stepping schemes of Euler-Maruyama type, for which we show strong convergence of a certain rate.
A McKean-Vlasov equation (introduced in [46,47]) for a d-dimensional process X = (X t ) t∈[0,T ] , with a given finite time-horizon T > 0, is an SDE where the underlying coefficients depend on the current state X t and, additionally, on the law of X t .We consider more specifically the one-dimensional equation of the form dX where L 0 2 (R) denotes the space of real-valued, F 0 -measurable random variables with second finite moments, (W t ) t∈[0,T ] is a one-dimensional standard Brownian motion and L Xt denotes the marginal law of the process X at time t ∈ [0, T ].In particular, we are concerned with the well-posedness of equations of the form (1), where R ∋ x → b(x, µ) is discontinuous in zero, and piecewise Lipschitz on the subintervals (−∞, 0) and (0, ∞).Concerning the measure component of the drift, we will require global Lipschitz continuity with respect to the Wasserstein distance with quadratic cost denoted by W 2 (see below for a precise definition).The diffusion term will be only state dependent and globally Lipschitz continuous.
Our setting contrasts with the standard case with globally Lipschitz continuous coefficients, which is well-studied in the literature, both from an analytic and numerical perspective, e.g., in [48], [60] and [7], respectively.
The study of SDEs and McKean-Vlasov equations with discontinuous drift is motivated by such models in biology (see, e.g., [25]) and financial mathematics (see, e.g., Atlas models in equity markets [3,33] and dividend maximisation problems [61]).Further, in stochastic control a discontinuous control can lead to equations with discontinuous drift (see [59]).In the context of stochastic N -player games, non-smooth cost functions (such as the ℓ 1 -regularisation) or constraints on the size of the control process can result in discontinuous controls (bang-bang type optimal controls) and hence will give controlled state dynamics with discontinuous drift, as in [11].
We start our literature review with some key references on standard SDEs with irregular and discontinuous drift, namely [64,63,40,41], and then proceed to discuss some recent articles on McKean-Vlasov SDEs with non-Lipschitz drift.Zvonkin [64] (for one-dimensional SDEs) and Veretennikov [63] (for the multi-dimensional setting) prove the existence of a unique strong solution for an SDE where the drift is assumed to be measurable and bounded, but the diffusion coefficient σ needs to satisfy rather strong assumptions, namely that it is bounded and uniformly elliptic, i.e., there is a λ > 0 such that for all x ∈ R d and all v ∈ R d , we have v ⊤ σ(x)σ(x) ⊤ v ≥ λv ⊤ v.An interesting addition to the aforementioned results in the case where the diffusion is not uniformly elliptic was established in the one-dimensional case in [40].The authors assume the drift coefficient to be piecewise Lipschitz and σ to be globally Lipschitz with σ(η) = 0 for each of finitely many points of discontinuity η of the drift.This condition guarantees that the process does not spend a positive amount of time in the singularity.Under these assumptions, by explicitly constructing a transformation that removes the singularities, the existence of a unique strong solution can be proven and a numerical procedure for solving this class of SDEs can be constructed.
The main contribution of [41] is the extension of the one-dimensional case to the multi-dimensional setting under the assumption of piecewise Lipschitz continuity of the drift.In [41], the authors introduce a meaningful concept of piecewise Lipschitz continuity in higher dimensions, which is based on the notion of the so-called intrinsic metric.As already indicated by the one-dimensional case, there needs to be an intricate connection between the geometry of the set of discontinuities and the diffusion coefficient.We note that the exceptional set of singularities, denoted by Θ, is assumed to be a C 4 hypersurface and for the diffusion part one requires the following: There exists a constant C > 0 such that |σ(η) ⊤ n(η)| ≥ C for all η ∈ Θ, where n(η) is orthogonal to the tangent space of Θ in η and |n(η)| = 1.Under these assumptions (and some additional technical conditions on the coefficients and on the geometry of Θ) the existence of a unique strong solution for multi-dimensional SDEs with piecewise Lipschitz continuous drift can be proven.
Moving on to McKean-Vlasov equations, the existence and uniqueness theory for strong solutions of such SDEs with coefficients of linear growth and Lipschitz type conditions (with respect to the state and the measure component) is well-established (see, e.g., [60,13]).More general existence/uniqueness results for weak and strong solutions of McKean-Vlasov SDEs can be found in [50,38,6].The article [6] is concerned with the weak and strong existence/uniqueness of one-dimensional equations with additive noise, where the drift is assumed to be measurable, continuous in the measure component with respect to the Monge-Kantorovich metric and further satisfies a linear growth condition.In [50] and [38], a d-dimensional setting is considered, where the drift is assumed to be bounded, measureable (and possibly path-dependent) and Lipschitz continuous in the measure component with respect to the total variation distance.The diffusion is non-degenerate and independent of the measure.Under these assumptions (and some technical conditions) weak existence and uniqueness is proven.
For further recent existence and uniqueness results for strong and weak solutions of McKean-Vlasov SDEs, including results concerning standard Lipschitz assumptions on the coefficients, we refer to [30,44,15,20,21,57,27] and the references given therein.
The numerical analysis of SDEs with discontinuous drifts has received significant attention over the last few years, see, e.g., [40,41,42,53,52,51,29,54,19,39] and the references therein for well-posedness results and for strong and weak convergence rates of numerical schemes.
In particular, in [40] (for the one-dimensional case) and in [41] (for the multi-dimensional setting) the standard strong convergence rate of order 1/2 for a method derived from the Euler-Maruyama scheme was proven.However, the applicability of these schemes is limited as they require the explicit knowledge of a transformation (and its inverse) to map the SDE with discontinuous coefficients into one with Lipschitz continuous coefficients.In [42], an Euler-Maruyama scheme without the aforementioned transformation is introduced (in a multi-dimensional setting).While this scheme is easier to apply, the authors only show a strong convergence rate of order 1/4 − ε for any ε > 0, imposing also the stronger assumption of boundedness for both coefficients of the underlying SDE.The central idea of [42] is to quantify the probability that a multi-dimensional process is in a small neighbourhood of the set of discontinuities, using an occupation time formula.In the one-dimensional case, with coefficients of linear growth, these techniques were refined in [51] and the expected strong convergence rate of order 1/2 was recovered.Other recent works concerned with the numerical approximation of SDEs with discontinuous drifts include [53,52], where a higher order scheme and an adaptive time-stepping scheme were introduced, respectively.In [39] a numerical scheme for classical one-dimensional diffusion processes generated by a differential operator involving discontinuous coefficients is presented.As the generator is non-local for McKean-Vlasov equations it seems a challenging problem to use these techniques in our framework.
The simulation of McKean-Vlasov SDEs typically involves two steps: First, at each time t, the true measure L Xt is approximated by the empirical measure where δ x denotes the Dirac measure at point x and ( ) ⊤ t∈[0,T ] , an interacting particle system, is the solution to the R dN -dimensional SDE with components Here, W i = (W i t ) t∈[0,T ] and ξ i , for i ∈ {1, . . ., N }, are independent Brownian motions (also independent of W ) and independent copies of ξ, respectively.In a next step, one needs to introduce a reasonable timestepping method to discretise the particles (X i,N t ) t∈[0,T ] over some finite time horizon [0, T ].Numerical schemes for interacting particle systems with Hölder continuous coefficients (in the state variable) and with coefficients satisfying certain assumptions on monotonicity (in the state variable) and Lipschitz continuity (in the measure variable), can be found in [4] and [56,5,37] (and the references cited therein), respectively, where a strong convergence analysis is conducted.In [9,1], a quantitative L p -error analysis in terms of density and cumulative distribution function approximation is presented.The survey [7] discusses several examples and numerical schemes for McKean-Vlasov equations involving singular drifts, e.g., a probabilistic interpretation of the Burgers equation, see also [9], of the 2D-incompressible Navier-Stokes equation (see e.g., [22,49]) and turbulent flow models [55].Other examples of McKean-Vlasov equations with singular drifts appear in the Keller-Segel equation [26], the Coulomb gas model [17], the Thomson problem [32], and the Stefan problem [35,36].
Our numerical schemes present an original approximation method which, as of now, is restricted to the specific case of a one-dimensional spatial and one-point discontinuous drift component, but provides, again in this specific framework, a suitable alternative to the methodical mollification/cut-off approximation methods.
In this article, we first focus on the decomposable case, namely that where b 1 is piecewise Lipschitz continuous with a discontinuity in zero, and b 2 satisfies the usual Lipschitz assumptions in both components.This structure allows us to present the main ideas of the analysis to be used later in a more general setting, in particular a transformation of the state variable to remove the discontinuity.In this setting, we prove well-posedness of the McKean-Vlasov equation and the associated particle system.This structure includes the important class of McKean-Vlasov equations of the form where V describes an external potential and β an interaction kernel; see, e.g., [31] and the references cited therein related to mean-field over-damped Langevin equations.These models also embed the class of self-stabilizing diffusions and the McKean-Vlasov model related to the granular media equation.
We then relax the structural assumption on decomposability slightly, but still have to require certain continuity of the measure derivatives at the points of discontinuity, which encompasses the above setting as a special case.The necessity for this condition arises due to the explicit measure or time dependence of the employed transformation.A future research direction concerns a setting where the point of discontinuity is time-dependent, or depends on the distribution of the process (X t ) t∈[0,T ] , which is relevant to study further practically important models, e.g., from [33].
Having established the existence of a unique strong solution with bounded moments, we propose two Euler-Maruyama schemes for the particle systems as numerical approximations to the McKean-Vlasov equations.For an Euler-Maruyama scheme applied to the SDE in the transformed state, strong convergence of order 1/2 follows immediately, while for a direct time-discretisation of the particle system without transformation, we are only able to show order 1/9.Numerical tests indicate that this order is in general not sharp.We will discuss the reasons for this gap and possible improvements later.
The main contributions of the present article are as follows.First, we establish the well-posedness of McKean-Vlasov SDEs (with a certain discontinuity) and of their associated particle systems.Techniques from variational calculus on the measure space P(R) equipped with the Wasserstein distance W 2 will be essential in the proofs, due to the possible measure dependence of the transformation applied to the processes as described above.The second central contribution of the present paper is the development of numerical schemes for approximating such McKean-Vlasov SDEs and their associated particle systems.
Here, a non-standard strong convergence analysis based on occupation time estimates of the discretised processes in a neighbourhood of the discontinuity will be presented.
The remainder of the paper is organised as follows: In Section 2, we collect all preliminary tools and notions needed throughout the paper.The precise problem description and the main results are presented in Section 3.Then, Section 4 discusses numerical schemes for McKean-Vlasov SDEs with discontinuous drift.We show strong convergence of certain orders with respect to the number of particles and time-steps, respectively.In Section 5, we apply our numerical scheme to a model problem arising in neuroscience [25] and to a slight modification of a mean-field game in systemic risk [28,16].

Preliminaries
In the sequel, we will introduce several concepts and notions, which will be needed throughout this article.In addition, we will give a brief introduction to the so-called Lions derivative (abbreviated by Lderivative), which allows us to define a derivative with respect to measures of the space P 2 (R) (see below for a precise definition).Also, we recall the transformation used to cope with drifts having discontinuities in a given finite number of points and first developed in [40].We give a summary of important properties of this mapping.Note that generic constants used in this article are denoted by C > 0. They are independent of the number of particles and number of time-steps, and might change their values from line to line.

Notions and notation
We start with introducing some notions and fixing the notation.
• Throughout this article, (Ω, F , (F t ) t∈[0,T ] , P) will denote a filtered probability space, where is the natural filtration of W augmented with an independent σ-algebra F 0 and (Ω, F , P) is assumed to be atomless.
Euclidean space.As a matrix-norm, we will use A := sup v∈R d |Av|, for any A ∈ R d×d .
• We use P(R) to denote the family of all probability measures on (R, B(R)), where B(R) denotes the Borel σ-field over R and define the subset of probability measures with finite second moment by • We recall the definition of the standard Wasserstein distance with quadratic cost: For any µ, ν ∈ P 2 (R), we define , where Π(µ, ν) denotes the set of all couplings between µ and ν.
• For a given p ≥ 2, L 0 p (R) refers to the space of real-valued, F 0 -measurable random variables X satisfying E[|X| p ] < ∞ and for a terminal time T > 0, S p ([0, T ]) refers to the space of realvalued continuous, F-adapted processes, defined on the interval [0, T ], with finite p-th moments, i.e., processes We briefly introduce the L-derivative of a functional f : P 2 (R) → R, as it will appear in the proofs presented in the main section.For further information on this concept, we refer to [12] or [10].Here, we follow the exposition of [14].We will associate to the function f a lifted function f , defined by f (X) = f (L X ), where L X is the law of X, for X ∈ L 2 (Ω, F , P; R).
This will allow us to introduce L-differentiability as Fréchet derivative on the lifted space.In particular, a function f on P 2 (R) is said to be L-differentiable at µ 0 ∈ P 2 (R) if there exists a random variable X 0 ∈ L 2 (Ω, F , P; R) with law µ 0 , such that the lifted function f is Fréchet differentiable at X 0 .Now, the Riesz representation theorem implies that there is a (P-a.s.) unique Φ ∈ L 2 (Ω, F , P; R) with with the standard inner product and norm on L 2 (Ω, F , P; R).If f is L-differentiable for all µ 0 ∈ P 2 (R), then we say that f is L-differentiable.
It is known (see, e.g., [14,Proposition 5.25]) that there exists a Borel measurable function χ : R → R, such that Φ = χ(X 0 ) almost surely, and hence Note that χ only depends on the law of X 0 , but not on X 0 itself.We define ∂ µ f (L X0 )(y) := χ(y), y ∈ R, as the L-derivative of f at µ 0 .If, in addition, for a fixed y ∈ R, there is a version of the mapping We require a definition describing regularity properties of a function f : P 2 (R) → R in terms of the measure derivative (see [14,18]).Definition 2.1.Let f : P 2 (R) → R be a given functional.
• We say that f is an element of the class C (1,1) b , if f is continuously L-differentiable, for any µ, there is a continuous version of the mapping R ∋ y → ∂ µ f (µ)(y) and the derivatives exist, are bounded and jointly continuous in the variables (µ, y) such that y ∈ Supp(µ).
• We say that f is an element of the class and in addition the second order Lions derivative ∂ 2 µ f (µ)(y, y ′ ) exists, is bounded and is again jointly continuous in the corresponding variables.Also, the joint continuity of all derivatives is here required globally, i.e., for all (µ, y, y ′ ).
We give the following additional remark, which links the L-derivative of functions of empirical measures to the standard partial derivatives of its empirical projections.For a functional f : P 2 (R) → R, we associate with it the finite dimensional projection f N : R N → R defined as , then f N is twice differentiable (in a classical sense) and where δ i,k is the Kronecker delta, see, e.g., [14,Proposition 5.35].

Properties of the transformation G
In [40], the authors consider one-dimensional SDEs of the form with a piecewise Lipschitz continuous drift coefficient b that is discontinuous in K ∈ N points η 1 , . . ., η K , and a Lipschitz diffusion coefficient σ that does not vanish in any η k .A mapping G : R → R is defined to transform the SDE into one for Z = G(X) with globally Lipschitz continuous coefficients.For simplicity, we restrict the discussion to K = 1 with η 1 = 0. We define the mapping G by where and c is a constant satisfying 0 < c < 1/|α|.The choice of α yields a Lipschitz continuous drift coefficient for the SDE of Z = G(X), in particular, it removes the discontinuity in 0 from the drift.The restriction on c guarantees that G possesses a global inverse.
It is known from [41] that G satisfies the following properties: Therefore, G is Lipschitz continuous and has an inverse G −1 : R → R that is Lipschitz continuous as well.

Existence and uniqueness results
The following subsections are devoted to proving well-posedness results for certain classes of onedimensional McKean-Vlasov SDEs with a drift having a discontinuity in zero.In a first step, we study a simple class where the resulting transformation will not depend on the measure.Here, the transformation techniques developed in [40] will allow us to prove existence and uniqueness of a strong solution.The second class of McKean-Vlasov SDEs investigated below has the intrinsic difficulty that the required transformation will depend on the measure (i.e., will be time-dependent).Hence, a fixed-point iteration in the measure component will be required and we need to use techniques from variational calculus on the measure space P 2 (R), in particular an Itô formula for functionals acting on this space.
For each of these classes of McKean-Vlasov SDEs, we will additionally study the well-posedness of their associated interacting particle system.Although they can be considered as N -dimensional classical SDEs, with N denoting the number of particles, the resulting set of discontinuities of the N -dimensional drift cannot be handled by the main results of [41].
Future work is needed to extend the methods developed in this article to a multi-dimensional framework.In particular, it seems that the decomposable case can be generalised when discontinuities of the form discussed in [41] are considered.

McKean-Vlasov SDEs and interacting particle systems with decomposable drift
For a given terminal time T > 0 and given p ≥ 2, we consider a one-dimensional McKean-Vlasov SDE of the form dX where b : R × P 2 (R) → R and σ : R → R are measurable functions.
In the following, we state the model assumptions which will specify the set-up for this subsection: We have σ(0) = 0 and there exists a constant L > 0 such that ( The drift is decomposable in the following sense: where b 1 : R → R is Lipschitz continuous on the subintervals (−∞, 0) and (0, ∞) and there exists a constant L 1 > 0 such that We now state the main results of this section:  The interacting particles (X i,N t
Proof.In contrast to the set-up in [41], the set of discontinuities, denoted by Θ, is not a differentiable manifold, but has the form

However, we may define G
where G is as in Proposition 3.1 and x N = (x 1 , . . ., x N ) ⊤ , which allows us to transform the particle system into a new particle system with globally Lipschitz continuous coefficients.Now, G N has a global inverse, as the mapping G has a global inverse, due to the choice of c (see Section 2.2).Therefore, applying Itô's formula to the inverse allows to deduce the claim.

McKean-Vlasov SDE with non-decomposable drift
Here, we consider again a one-dimensional McKean-Vlasov SDE of the form (3), However, in contrast to the above setting, we will not assume that b can be decomposed in two parts as in Assumption (H.1(2)) from the previous section, and therefore the transformation will also depend on the measure.To be precise, for any (x, µ) ∈ R × P 2 (R), we define where and c > 0 is a constant small enough to guarantee the invertibility of G.When we speak of an 'inverse' of G(x, µ), we mean 'inverse with respect to x', i.e., the inverse is a function , G may also be viewed as a mapping G : R × [0, T ] → R.
In the following, we state the model assumptions which will specify the set-up for this subsection: ) is satisfied and we require: (1) There exist constants L, L 1 > 0 such that sup Additionally, for all µ ∈ P 2 (R), R ∋ x → b(x, µ) is piecewise Lipschitz continuous on the subintervals (−∞, 0) and (0, ∞), uniformly with respect to µ.
(2) α ∈ C (1,1) b , and the mapping (3) For any µ ∈ P 2 (R), the mapping R ∋ y → ∂ µ α(µ)(y) vanishes in zero, i.e., is needed to apply an Itô formula for α (see [14,Proposition 5.102]).Remark 3.2.Note that (H.2(4)) could also be replaced by the following alternative set of assumptions: On each of the two subintervals (−∞, 0) and (0, ∞), b(•, µ) is a C 1 (R, R) function with bounded derivative and, additionally, for any µ ∈ P 2 (R), the mapping R ∋ y → ∂ y ∂ µ α(µ)(y) vanishes in zero, i.e., The following proposition shows the Lipschitz continuity of the mapping R× P 2 (R) ∋ (x, µ) → G −1 (x, µ).Proposition 3.3.Let the function G be defined as in (6) with c < 1/ sup µ∈P2(R) |α(µ)| and (7) and let Assumption (H.2(1)) be satisfied.Then, there exists a constant L(c) satisfying L(c) → 0 as c → 0, such that for any x, y ∈ R and µ, ν ∈ P 2 (R) Proof.First, we note that by differentiating (6), we get for It is easy to verify that sup which implies that for all x ∈ R and µ ∈ P 2 (R) In particular, if c < 1/ sup µ∈P2(R) |α(µ)|, then for all x, y ∈ R and µ ∈ P 2 (R), we have ).It is easy to show that the mapping x → ∂ x G(x, µ) is also Lipschitz continuous.Denote the Lipschitz constant of ∂ x G with respect to the first and second argument by L x and L µ , respectively.Writing we obtain For |x| < c, we have and hence Gronwall's inequality implies ) by the definition of G.We finally obtain, for all x, y ∈ R and all µ, ν ∈ P 2 (R), with L(c) := 4L µ ce 4Lxc , In what follows, we will assume that c < 1/ sup µ∈P2(R) |α(µ)| and is small enough such that the Lipschitz constant of the mapping P 2 (R) ∋ µ → G −1 (x, µ) (i.e., the constant L(c) a few lines above) is less than a half.The reason for this requirement will become clearer in the proof of Theorem 3.1.
Similar to the previous section, we aim to recover a unique strong solution of ( 5) by setting X t = G −1 (Z µ t , µ t ), where µ t = L Xt for t ∈ [0, T ], and (Z µ t ) t∈[0,T ] is the process obtained by applying the transformation G to X.Even though G is not twice continuously differentiable in the state variable, Itô's formula is still applicable due to the special form of the discontinuity (see [24,Theorem 2.1] and the comments after the proof of this theorem).Now, observe that Itô's formula along a flow of measures (µ t ) t∈[0,T ] ∈ C([0, T ], P 2 (R)) (see, e.g., [14,Proposition 5.102]) implies where we recall that ∂ y ∂ µ G(x, µ t )(y) denotes the derivative of the mapping R ∋ y → ∂ µ G(x, µ t )(y) and Hence, we define In the following, we will show that the decoupled SDE (9) where the flow is fixed has Lipschitz continuous coefficients.Note that for a such a flow of measures the process X µ in (5) (interpreted as classical SDE) has bounded moments uniformly in (µ t ) t∈[0,T ] , which is a consequence of (H.1(1)) and (H.2(1)).In particular, for ξ ∈ L 0 2 (R), we have the a-priori estimate ]) =: C, where C > 0 only depends on T, p and the constants appearing in the model assumptions.We will introduce the following subspace of C([0, T ], P 2 (R)): We define where C is defined as above and complete this space with the metric sup t∈ Proof.Using (H.1(1)), the Lipschitz continuity of z → G −1 (z, ν t ) and the uniform boundedness of α, see (H.2(1)), in combination with [41, Lemma 2.5], gives the Lipschitz continuity of z → σ(z, ν t ), with a Lipschitz constant independent of ν t .Similarly, we can deduce that ν t → σ(z, ν t ) is Lipschitz continuous, due to Proposition 3.3 and the Lipschitz continuity of ν t → α(ν t ).
We are now ready to present the main result and its proof of this section: Theorem 3.1.Let Assumption (H.2) be satisfied, let ξ ∈ L 0 p (R) for a given p ≥ 2 and assume that the constant c is sufficiently small (as in Remark 3.3).Then, the McKean-Vlasov SDE defined in (5) has a unique strong solution in S p ([0, T ]).
Proof.First, we remark that for any given flow of measures (µ t ) t∈[0,T ] ∈ C([0, T ], P 2 (R)) the SDE defined in (9) has a unique strong solution by Lemma 3.1.Now, for (µ t ) t∈[0,T ] , (ν t ) t∈[0,T ] ∈ P b , we obtain, for any t ∈ [0, T ], using Lemma 3.1, BDG's inequality and Hölder's inequality along with Gronwall's inequality For k ≥ 0 and t ∈ [0, T ], we define the Picard iteration with Z µ 0 t = G(ξ, δ ξ ) and µ 0 t = δ ξ .Note that by Itô's formula (applied to G −1 ), the process defined by Recall that (X k+1 t ) t∈[0,T ] has uniformly bounded moments (uniformly in k), due to (H.1(1)) and (H.2(1)), i.e., we have The applicability of Itô's formula for G −1 is a consequence of the fact that the inverse inherits the regularity of G, in particular the mapping P 2 (R) ∋ µ → G −1 (y, µ) is still an element of the class C where L := 2L(c) < 1 due to the choice of c and C > 0 is a constant depending on the constants appearing in Proposition 3.3 and Lemma 3.1.Let now 0 < T 0 < T such that CT 0 + L < 1.With this choice the uniqueness of ( 5) on [0, T 0 ] follows from the estimate in (15) by assuming there exist two solutions (X, µ) and (Y, ν) to (5), with µ t and ν t the marginal laws of X t and Y t , respectively, for t ∈ [0, T 0 ].In addition, we observe that the sequence of flows (µ k ) k , for µ k = (µ k t ) t∈[0,T0] , is a Cauchy sequence in the complete metric space C([0, T ], P 2 (R)) equipped with the Wasserstein distance sup t∈[0,T0] W 2 (µ t , ν t ).Hence, (13) has a fixed point, in particular we have X t = G −1 (Z µ t , µ t ), where µ t = L Xt .Itô's formula applied to G −1 yields the claim for the time interval [0, T 0 ].Repeating the above procedure starting at T 0 , we can extend the solution to the interval [T 0 , T 1 ], for some T 0 < T 1 < T .This is possible as the choice of T 1 depends on X T0 only through the second moment of X T0 , for which we have a uniform bound for the entire interval [0, T ], see (14).Proceeding in such a manner, we can obtain well-posedness of ( 5) on [0, T ].

Interacting particle system with non-decomposable drift
In the following, we state the model assumptions which will specify the set-up for this subsection: H.3. Assumptions (H.2(3)) and (H.2(4)) are satisfied and we require: (1) Assumption (H.1(1)) holds and there exists a constant L > 0 such that |σ(x)| ≤ L for all x ∈ R.
is a bounded function and the mappings are bounded and Lipschitz continuous.
Remark 3.4.Note that, compared to (H.2(1)), we do not require the drift to be uniformly bounded in the measure component.
In contrast to the case of particle systems with decomposable drift, we set α : (which could also be interpreted as a mapping α N : R N → R) and apply to each particle the following transformation G : R × P 2 (R) → R, We set R N ∋ x N → G i (x N ) := G x i , µ x N , and use these mappings to define G N : R N → R N by To obtain the transformed process ( we proceed as follows: For any t ∈ [0, T ] and i ∈ {1, . . ., N }, we have, using [14, Proposition 5.35] (see also Section 2) The applicability of Itô's formula for the function G is guaranteed by [41,Theorem 3.19].Note that Assumption 3.4 therein is imposed to guarantee Lipschitz continuity of second order derivatives of G outside the set of discontinuities.Assumption (H.3(3)) and the definition of G substitute this condition.
Assuming for now the global invertibility of G N , we may introduce where and In the following lemma, we will prove the invertibility of G N .
Lemma 3.2.Let Assumption (H.3(3)) be satisfied and assume that the constant c in (17) satisfies c < min 1, sup Then, G N has a global inverse.
Proof.We will employ Hadamard's global inverse function theorem (see, e.g., [58, Theorem 2.2]) to prove that G N has a global inverse.To do so, we need to verify the following properties of The first two mentioned conditions are obvious, due to the definition of G N and the uniform boundedness of α and φ(x) := x|x|φ x/c .Hence, we need to prove that G ′ N (x N ) is invertible.First, note that where N is a row vector.Now, we define and remark that G ′ N (x N ) can be identified with the linear operator I N ×N + A(x N ) : R N → R N .Therefore, succeeding in showing that c can be chosen (uniformly in x N ) in a way such that the operator norm of A(x N ) is smaller than one would yield the claim, as in this case G ′ N (x N ) is close to the identity (see, [41,Lemma 3.17]).We compute 3)) guarantees that α and its derivatives are uniformly bounded, i.e., c can be chosen uniformly in x N .Therefore, Hadamard's global inverse function theorem proves that, for each given N ≥ 1, G N : R N → R N is a diffeomorphism.
We proceed by showing that the transformed SDE has (locally) Lipschitz continuous coefficients.Lemma 3.3.Let Assumption (H.3) be satisfied.Then, the coefficients B N and Σ N introduced in ( 18) and ( 19), respectively, are locally Lipschitz continuous with linear growth.
Proof.For x N , y N ∈ R N , we obtain using ( 19) where we used (H.Noting that we further obtain for the drift That the terms Π 3 , Π 4 and Π 5 allow a Lipschitz bound is a consequence of (H.3 (1)) and (H.3(3)).
The linear growth of B N and Σ N , i.e., that there exists a constant C > 0 such that is a direct consequence of the growth conditions on b and σ along with the bounds for the derivatives of G and α.Theorem 3.2.Let Assumption (H.3) be satisfied, let ξ ∈ L 0 p (R) for a given p ≥ 2 and assume that the constant c in (17) satisfies (20).Then, the interacting particle system defined in ( 16) has a unique strong solution in S p ([0, T ]).
Proof.From Lemma 3.3 and the linear growth of B N and Σ N , we can deduce that the SDE for Z N has a unique strong solution (see, [45, Chapter 5, Theorem 2.5]).Applying now Itô's formula to G −1 N (Z N t ) proves the strong uniqueness of the particle system defined in (16).Note that G −1 N exists due to Lemma 3.2 and Itô's formula is applicable for G −1 N as it inherits the regularity of G N (see, Appendix A.2 for details).

Euler-Maruyama scheme with and without transformation
In this section, we restrict most of our discussion to the case of a McKean-Vlasov SDE with decomposable drift, due to the simpler structure of the underlying transformation and the particle systems.
In the following subsections, we will present two Euler-Maruyama schemes to discretise the particle system defined in (4) in time.
For the first scheme (Scheme 1), we will discretise the transformed (continuous) particle system in time and then exploit the global inverse G −1 to obtain approximations of the original (discontinuous) particle system.A slight modification of this scheme will also be applied in the non-decomposable case.An approximation result with respect to the number of particles will also be presented.
The second scheme (Scheme 2) will be defined by directly discretising the discontinuous particle system, without making use of the transformation G.We give strong convergence rates in terms of the number of time-steps and pathwise strong propagation of chaos results in order to obtain quantitative L 2 -approximations for the underlying McKean-Vlasov SDE.

Scheme 1: Euler-Maruyama after transformation (decomposable case)
We define the following explicit Euler-Maruyama scheme to discretise the particle system (4) in time.In a first step, we partition a given time interval [0, T ] into subintervals of equal length h = T /M , for some integer M > 0, and define t n := nh.Then, we simulate the transformed particle system by for n ∈ {0, . . ., M − 1}, where We introduce the notation η(t) := sup{s ∈ {0, h, . . ., M h} : s ≤ t}, for t ∈ [0, T ], which allows us to define the continuous time version of ( 21) Then, we propose an Euler-Maruyama approximation to X i,N t , for i ∈ {1, . . ., N } and t ∈ [0, T ], by The convergence of this algorithm is proven in the following theorem: Theorem 4.1.Let Assumption (H.1) be satisfied, let ξ ∈ L 0 p (R) for some p > 4 and assume c < 1/|α|.For i ∈ {1, . . ., N }, let (X i t ) t∈[0,T ] be the unique strong solution of (3) driven by the Brownian motion (W i t ) t∈[0,T ] with initial data ξ i , and (X i,N,M t ) t∈[0,T ] be given by ( 22) and (23).Then, there exists a constant C > 0 (independent of N and M ) such that and recall that the dynamics of where μZ . Therefore, we obtain for some constant C > 0 where the last inequality can be derived similarly to the propagation of chaos results for equations with Lipschitz continuous coefficients as in, e.g., [13] for d = 1 (note that the rate N −1/2 can be improved into N −1 in case where b 2 (x, µ) = R β(x, y) µ(dy), with β Lipschitz continuous, see [48]).To apply the aforementioned propagation of chaos result, we need the requirement that ξ ∈ L 0 p (R) for p > 4 (see, also [14,Theorem 5.8]).Furthermore, we have since the SDEs for (Z i,N t ) t∈[0,T ] and (Z i,N,M t ) t∈[0,T ] have globally Lipschitz continuous coefficients.From these two estimates the claim follows.

Scheme 2: Euler-Maruyama without transformation (decomposable case)
As G and G −1 may be difficult to construct in multi-dimensional settings, and since the evaluation for the inverse at each time point can be computationally expensive, it would be preferable to discretise the particle system (X i,N ) i∈{1,...,N } in time directly, without the use of the transformation G.In addition, a drawback of Scheme 1 is that an SDE with additive diffusion term will be transformed into one with multiplicative noise and therefore the Euler-Maruyama scheme no longer coincides with the Milstein scheme.We employ an Euler-Maruyama scheme to discretise the particle system (4) and compute an approximate solution by The following results are concerned with moment stability of the discretised particle system and estimates for the occupation time of the particle system in the neighbourhood of the set of discontinuities.

Moment stability:
We first remark that due to the linear growth of the coefficients in the state component, the Lipschitz continuity in the measure variable and the fact that all particles are identically distributed, we have the following result (see, e.g., [45] for details): Proposition 4.1.Let Assumption (H.1) be satisfied, and let ξ ∈ L 0 p (R) for p ≥ 2.Then, there exist constants and for all i ∈ {1, . . ., N } and for all t ∈ [0, T ], Occupation time formula for (24): Below, we will show an estimate of the expected occupation time of a fixed particle of the system defined by ( 24) in a neighbourhood of zero.
Proposition 4.2.Let Assumption (H.1) be satisfied and let i be an arbitrary but fixed particle index.Further, let ξ ∈ L 0 p (R) for p ≥ 2 and let (X i,N,M t ) t∈[0,T ] be given by (24).Then, there exists a constant C > 0 such that for all N, M ∈ N and all sufficiently small ε > 0 Proof.We aim to apply [42, Theorem 2.7], which states the following: Let (X t ) t∈[0,T ] be an R d -valued Itô process with progressively measurable processes A = (A t ) t∈[0,T ] and B = (B t ) t∈[0,T ] , where A is R d -valued and B is R d×d -valued.The set of discontinuities, Θ, is assumed to be a C 3 hypersurface of positive reach.Namely, there exists ε > 0 such that p(x) = arg min y∈Θ |x − y| is a single valued function of class C 3 on the tubular neighbourhood Θ ε := {x ∈ R d : inf y∈Θ |x − y| < ε} (see Definition 2.4 in [42] for details).Then, there exists a constant C > 0, such that for all sufficiently small ε > 0 assuming, additionally, that 1. the processes A and B are almost surely bounded by a constant C if X t (ω) is in a small neighbourhood of Θ, and 2. there exists a constant C > 0 such that for almost all ω ∈ Ω, we have: where n(x) has length one and is orthogonal to the tangent space of Θ in x.
We now return to our particular model problem.First, we remark that Θ i satisfies all regularity conditions of [42, Theorem 2.7], i.e., it is a C 3 hypersurface of positive reach.We then observe that the N -dimensional particle system can be rewritten as where W N t = (W 1 t , . . ., W N t ) ⊤ and B N : R N → R N and Σ N : R N → R N ×N are defined by Further, we observe that there is a constant C > 0 such that: If, for any t ∈ [0, T ] and ω ∈ Ω, as σ is continuous and σ(0) = 0. Also, note that a normal vector of the tangent space of Θ i is e i , i.e., the i-th unit vector.Further, a close inspection of the proof of [42,Theorem 2.7], shows that the boundedness assumption on the coefficients in a neighbourhood of Θ i is not needed in our case, due to the moment bound of an individual particle established in Proposition 4.1.Hence, where the constant C > 0 is independent of N , due to the fact that the normal vector is the i-th unit vector.

Auxiliary proposition:
Based on the occupation time estimate from Proposition 4.2, we will prove the following result, which is needed in the proof of Theorem 4.2 given below.
Proposition 4.3.Let Assumption (H.1) be satisfied, and let ξ ∈ L 0 p (R) for p ≥ 8. Furthermore, let (X i,N,M t ) t∈[0,T ] be given by (24).Then, there exists a constant C > 0 (independent of N and M ) such that for any t ∈ [0, T ], we have Proof.First, observe that the linear growth of σ and the piecewise Lipschitz continuity of G ′′ imply that there exists a constant C > 0 such that where ε > 0 will be specified later.With this at hand, we derive ds , where we used Hölder's inequality, Proposition 4.1 and Proposition 4.2 in the last display.Markov's inequality along with Proposition 4.1 imply that there exists a constant C > 0 such that Choosing ε = h 4/9 gives the result.
We are now ready to present our main convergence result.In this case, we only obtain the strong convergence rate of order 1/9: Theorem 4.2.Let Assumption (H.1) be satisfied, let ξ ∈ L 0 p (R) for some p ≥ 8 and assume c < 1/|α|.Furthermore, let (X i t ) t∈[0,T ] be the unique strong solution of (3) driven by the Brownian motion (W i t ) t∈[0,T ] with initial data ξ i and (X i,N,M t ) t∈[0,T ] given by (24).Then, there exists a constant C > 0 (independent of N and M ) such that Proof.Note that where in the last display, we used the pathwise propagation of chaos result as in the previous subsection.Further, we have where in the last estimate, we employed standard strong convergence results for the Euler-Maruyama scheme applied to SDEs with globally Lipschitz continuous coefficients.Following similar arguments to [42] or [51], one further obtains, by applying Itô's formula to G(X i,N,M t ), where the second summand on the right side is of order h 2/9 due to Proposition 4.3.Hence, Gronwall's inequality yields and the claim follows combining (25) and (26).
Remark 4.1.The convergence rate in terms of the number of particles in the above theorem can again be improved to 1/2 if the drift has the form (2), see [48].The convergence rate in terms of number of time-steps established in Theorem 4.2 could be improved by employing exponential tail estimate techniques, as in [42].The resulting strong convergence rate would be 1/4 − ε, for an arbitrarily small ε > 0. However, to achieve this, one would need to assume boundedness of the coefficients in equation (3).Another possibility to recover a better convergence rate in our setting would be to require that the initial data ξ ∈ L 0 p (R) for any p ≥ 2 and that σ is uniformly bounded.This would enable us to obtain sharper estimates, when employing Markov's inequality in the proof of Proposition 4.3.If we assume moment boundedness of the initial data of any order, but allow σ to grow-linearly, we would obtain a rate of order 1/8 − ε.
Moreover, although we expect that the optimal convergence rate of the Euler-Maruyama scheme applied to the interacting particle system is 1/2 (as for equations with Lipschitz coefficients), we only achieve the order 1/9 (or 1/4 − ε under stronger assumptions on the initial data or the coefficients of the underlying McKean-Vlasov SDE), due to the estimate of the probability that X i,N,M η(s) and X i,N,M s have a different sign, i.e., that the term |G ′′ (X i,N,M s ) − G ′′ (X i,N,M η(s) )| in the proof of Proposition 4.3 does not allow a Lipschitz type estimate.Refined estimates of the aforementioned expected sign change, as derived in [51] for one-dimensional SDEs, are not easy to prove for an interacting particle system.The proof of the main result in [51] is not applicable to our setting as an individual particle (seen as a one-dimensional equation) does not satisfy Markov properties (due to the dependence of interaction terms), which are key in [51].
Remark 4.2.The implicit function theorem applied to the function implies that we can express y N in terms of x N .The applicability of the implicit function theorem follows from similar arguments to the ones presented in Lemma 3.2 along with Proposition A.1.].Similar arguments as for Scheme 2 could possibly be used here, but our current analysis does not allow us to derive an optimal convergence rate in h.Theorem 4.3.Let Assumption (H.3) hold, let ξ ∈ L 0 p (R) for p > 4 and assume the mapping R \ {0} × P 2 (R) ∋ (x, µ) → b(x, µ) to be uniformly bounded.Let (X i t ) t∈[0,T ] be the unique strong solution of (5) driven by the Brownian motion (W i t ) t∈[0,T ] with initial data ξ i , and X i,N,M tn for n ∈ {0, . . ., M } be defined by the above algorithm.Then, there exists a constant C > 0 (independent of N and M ) such that Proof.We start with the observation that for any n ∈ {0, . . ., M } Using the arguments from the previous lemma, the first term can be shown to be of order O(N −1 ).For the second term, we derive the estimate where L(c) → 0 as c → 0, as in Proposition 3.3.Therefore, we get Using now the definition of Xi,N,M tn and Lemma 3.3, one shows that there exists a constant C > 0 such that E |G(X i,N tn , µ X N tn ) − Xi,N,M tn | 2 ≤ Ch.This together with Lemma 4.1 yields the claim.

Numerical illustration
In the following, we will present two examples of McKean-Vlasov SDEs and interacting particle systems exhibiting discontinuous drifts in order to motivate the theoretical study of such equations and to numerically illustrate the strong convergence behaviour of an Euler-Maruyama scheme.These models find applications in biology and mathematical finance, in particular systemic risk.
As we do not know the exact solution of the considered equations, the convergence rate (in terms of number of time-steps) was determined by comparing two solutions computed on a fine and coarse grid, respectively, where the same Brownian motions were used.In order to illustrate the strong convergence behaviour in the uniform time-step h, we compute the root-mean-square error (RMSE) by comparing the numerical solution at level l of the time discretisation with the solution at level l − 1 at the final time T = 1.To be precise, as error measure we use the quantity where M l = 2 l T and by X i,N,M l T we denote the approximation of X at time T computed based on N particles and 2 l T time-steps.The number of particles used in the tests will be specified below.where a ≥ 0 is the mean-reversion rate, κ 1 < 0, κ 2 > 0 and σ > 0. The strong well-posedness of (30) follows from Proposition 3.1.This equation can be linked to a model of systemic risk in [16], where a mean-field game of N banks borrowing from, and lending to, a central bank is proposed.The banks control the rate of their borrowing depending on their (log-)monetary reserves, which are modelled by a system of SDEs with interaction through their average.In this setting flocking, and thus systemic default events, may occur.Here, we give a slight reformulation of this problem following and μt = L X μ, β t for all t ∈ [0, T ].
For simplicity, we did not add the common noise term as in [16].In addition, we modified the diffusion to allow it to be degenerate.In [28], a constraint β t ∈ [κ 1 , κ 2 ] on the borrowing/lending rate is imposed for all t ∈ [0, T ].The minimiser of the objective function for this mean-field game (with constant diffusion term) is shown to be a control of bang-zero-bang type (see [28, equation (24)] for an analytic expression).Written in feedback form, this optimal control strategy has discontinuities that are time-dependent, i.e., the zero-control region changes over time, a setting not covered by our current analysis.
Using instead the special bang-bang type control (β * t ) t∈[0,T ] of the form and plugging β * t back into (31) results into an equation of the form (30) with a one-point discontinuity.In our numerical experiments, we set κ 1 = −0.5, κ 2 = 0.5, x = 0 and σ = 0.7.Further, we consider three different choices for the mean-reversion rate, i.e., we set a = 1, 5, 10.The expected value is approximated by the empirical mean of N = 10 4 particles.For a larger value of a, we expect sample paths of (30) to be more concentrated around the mean of the samples; see Fig. 2a and 2b (with T = 1 and M = 2 7 ).We can also observe that for a stronger concentration effect the strong approximation behaviour becomes better; see Fig.

A Auxiliary results
A.1 L-differentiability of G −1 Proposition A.1.Let G : R × P 2 (R) → R be defined by (7) with c > 0 as in Remark 3.3 and let the assumptions of Proposition 3.3 be satisfied.Assume that G is L-differentiable, and R ∋ y → ∂ µ G(x, µ)(y) is in C 1 (R, R) for all (x, µ) ∈ R × P 2 (R).Then, the inverse G −1 : R × P 2 (R) → R is continuously Ldifferentiable with for all (x, µ, y) ∈ R × P 2 (R) × R.
0) , in order to eliminate the discontinuity in zero, yields a McKean-Vlasov SDE with globally Lipschitz continuous coefficients.This can be shown in a similar manner to [40, Theorem 2.5]).Moreover, G has a global inverse due to the choice c < 1/|α| (see [41, Lemma 2.2]), and Itô's formula can be applied to G −1 , which allows to deduce the claim.

Remark 4 . 3 .= G − 1 (
We could also define an explicit scheme by setting X i,N,M tn Xj,N,M tn (dx).However, to derive a strong convergence rate for the resulting scheme, one has to analyse the quantity E[|X i,N,M tn − Xi,N,M tn | 2

1 Figure 1 :
Figure 1: Strong convergence of the Euler-Maruyama scheme applied to the particle system obtained by approximating the equation for the action potential of the neurons.

Figure 2 :
Figure 2: Sample trajectories of the particle system associated with (30) for N = 10 and a = 1 (left) and a = 10 (right).

Figure 3 :
Figure 3: Strong convergence of the Euler-Maruyama scheme applied to the particle system obtained by approximating the equation (30) with a = 1, 5, 10.