SDE with oblique reflection on domains defined by multiple constraints

We present simple assumptions on the constraints defining a hard core dynamics for the associated reflected stochastic differential equation to have a unique strong solution. Time-reversibility is proven for gradient systems with normal reflection, or oblique reflection with a fixed oblicity matrix. An application is given concerning the clustering at equilibrium of particles around a large attractive sphere.


Introduction
Since the first works of Skorokhod [11] on existence and uniqueness for strong solutions of reflected stochastic differential equations, many authors have investigated this type of equation and extended his results on half-spaces to more general domains: convex sets (Tanaka [12]), admissible sets (Lions-Sznitman [6]), domains satisfying only the Uniform Exterior Sphere and the Uniform Normal Cone conditions (Saisho [7]) or some weaker version of these conditions (Dupuis and Ishii [2]). The question of equilibrium states of the reflected process (construction of time-reversible initial measures) has also been investigated (see e.g. [10]).
All these existence, uniqueness and reversibility results were obtained under some smoothness assumptions on the boundary of the domain. Typically, the existence of at least one normal inward vector at each point of the boundary is a necessary condition to define the normal reflection direction.
In most cases, the domain in which the process has to live is defined by constraints which are physically natural rather than by its geometrical properties as a subset of some Euclidean space. For example, consider a system of n identical hard spheres with radius r in R d . The domain in which they evolve is the set of configurations (x i ) 1≤i≤n satisfying the constraints |x i − x j | > 2r (i.e. the distance between the centers of any two spheres is larger than twice their radius); the geometrical description is much more complicated: the complementary set in R nd of some starconvex subset whose boundary can be locally approximated by a tangent sphere and a cone.
Unfortunately, for reflected processes in dimension larger than three, the geometrical properties of the domain are not that obvious from the physical constraints. In the already mentioned nddimensional example of a finite system of hard spheres, the main part of the paper [8] is devoted to proving that the set of allowed ball configurations satisfies the Uniform Exterior Sphere and Uniform Interior Cone property. In [9] and [3] too, a meticulous and extensive geometrical study has to be performed before the stochastic analysis of the dynamics.
We present in this note a constraint-based assumption for existence, uniqueness and reversibility results for Skorokhod equations.
Our aim is to deal with assumptions as simple and physically natural as possible, even if they are not the weakest ones. Moreover, we treat the delicate case of oblique reflection. It will be the basis of a general approach for analysing stochastic dynamics under constraints. We also prove that for well-chosen diffusion coefficients the oblique reflected gradient diffusion admits a reversible measure.
This note is divided in two parts. The first part (section 2) exhibits a new compatibility criterion for constraints. If it is satisfied, then the reflected stochastic differential equation admits a unique strong solution, and this solution is time-reversible in the case of a gradient system whose reflection direction is consistent with its diffusion coefficient.
In section 3 we present an application. We consider the behaviour of many particles around a sphere called planet. The particles are spherical with random radii oscillating between a minimum and a maximum value. Their motion is perturbated by collisions into other particles and into the planet. The planet also generates a gravitational field which has a smooth attractive influence on the particles. We are interested in the equilibrium situation: Are the particles scattered in a wide area as a typical configuration, or do they tend to cluster around the planet ? We prove that at equilibrium and for low temperature all particles are as close as possible to the planet, all located beneath some altitude, with high probability.

Reflected stochastic differential equation under multiple constraints
We are interested in a process living in the closure of a domain D. This domain is defined by a finite set F of smooth R-valued constraint functions on R d : D is an intersection of smooth sets (arbitrary many of them, provided they are in finite number), so its boundary is a finite union of smooth boundaries: Since we want the process to be reflected on the boundary of D, we have to assume some regularity on the functions in F. The reflection at any point x ∈ ∂D occurs either in the inward normal direction ∇f (x) or in an oblique direction depending only on the normal direction. So we have to suppose the existence of a direction which is normal to the boundary: ∇f (x) = 0 for each x ∈ D such that f (x) = 0. We actually assume something more: the first derivatives of the functions of F should admit some positive uniform lower bound, their second derivatives should be uniformly bounded and, most important, we suppose that the boundary of each single-constraint set {x ; f (x) > 0} crosses the boundaries of the other single-constraint sets at "not too sharp an angle". To be more precise, we have to exclude infinitely sharp "thorns" which contain a point admiting inward normal vectors with opposite directions. This is what we call compatibility between the constraints: Definition 2.1 Let F be a finite set of R-valued C 2 -functions on R d . These functions are called compatible constraints if > 0 for each f ∈ F is a non-empty connected set ; • for each f ∈ F, inf{|∇f (x)|; x ∈ D, f (x) = 0} > 0 and where Conv(x) is the convex hull of the unit normal vectors to the boundaries at point x: Here and in the sequel, δ denotes the Euclidean distance in R d , |y| denotes the Euclidean norm of vector y, and |M | = sup{|M y|/|y| ; y ∈ R d } denotes the norm of the matrix M . Lebesgue measure is denoted by dx. The next main theorem states that our compatibility definition provides a convenient assumption to ensure the existence of a reflected process within a set defined by constraints. In most models, for the sake of simplicity, the reflection direction is the inward normal direction on the boundary. Here we consider a general oblique reflection, as in the case treated in [4] or in Section 3. We state the result with a fixed deviation θ t θ from the normal direction. t θ denotes the transposed matrix. Normal reflection corresponds to the special case θ = I d .
has for each starting point x ∈ D a unique strong solution in D, where the local times L f satisfy In this theorem "Strong uniqueness of the solution" stands for strong uniqueness, in the sense of [5] chap.IV def.1.6, of the process X, not of the local times L f . If σ is constant and the drift b is a gradient, then the equation also admits a time-reversible measure µ (i.e. the distribution of the solution with initial measure µ is invariant under the transformation X(·), (L f (·)) f ∈F −→ X(T − ·), (L f (T − ·) − L f (T )) f ∈F for each T > 0): Theorem 2.3 (Reversibility of the reflected gradient process) Let θ be a fixed d × d invertible matrix and F a set of compatible constraints. If Φ is a C 2 -function on R d with bounded derivatives, then the solution of where the dot denotes the Euclidean scalar product. Though this statement is longer and apparently more difficult to obtain than an uniform lower bound on the norms of the convex combinations, it is in some sense more intuitive. It states the existence of cones (with vertex x, axis v and aperture 2 arccos β 0 ) which contain all the inward normal vectors given by the constraints at point x. The positivity condition ensures that these cones do not degenerate into half-spaces. This condition is easier to check in some concrete situations (see section 3).

Remark 2.5: (Stability of the compatibility property)
Let F be a set of compatible constraints on R d .
• If all constraints disregard one of the d coordinates, then F induces a set of compatible The remaining of this section is devoted to the proofs of the above results. We first prove remarks 2.4 and 2.5 which will be useful in the other proofs, and then proceed to theorems 2.2 and 2.3.

Proof of remark 2.4
The third compatibility condition is ∃β 0 > 0 ∀x ∈ ∂D δ(0, Conv(x)) ≥ β 0 . The condition in remark 2.4 can be rewritten as Thus it suffices to prove that for each x ∈ ∂D |∇f (x)| .v, which holds for every unit vector v and every y ∈ Conv(x) because families (c f ) of non-negative numbers summing up to 1 satisfy: Since the convex hull Conv(x) is a closed set, it contains an element z with minimal norm: |z| = δ(0, Conv(x)). For each f satisfying f (x) = 0 and for each positive ε, the convex combination |∇f (x)| belongs to the convex hull, hence its norm can not be smaller than |z|: ) and provides the upper bound on δ(0, Conv(x)).

Proof of remark 2.5
Let us prove the compatibility of the set ∀f ∈ F f (θy) > 0 is a non-empty connected set as continuous image of the non-empty connected set D. θ also transforms the bounds on the f 's into bounds on the g's. Remark 2.4 with v replaced by θv provides the existence of some positive β 0 such that: Replacing x by θy we obtain: Thanks to remark 2.4 with β 0 = β 0 |θ −1 | | t θ| , this proves that F θ is a set of compatible constraints. In order to prove the second part of remark 2.5, we now assume that f (x, is a non empty connected set as a projection of a non-empty connected set. The lower bound on ∇f and the upper bound on D 2 f transfer to f because ∇f = (∇f , 0) and From the compatibility of F, we also get the existence of a positive β 0 such that for each Proof of theorems 2.2 and 2.3 The case of normal reflection: we assume here that θ = I d . According to corollary 3.6 of [3], equation (1) has a unique strong solution as soon as D satisfies the four assumptions of the inheritance criterion for Uniform Exterior Sphere and Uniform Normal Cone conditions (proposition 3.4 in [3]). We will check these four assumptions in the unusual order (i) (ii) (iv) (iii) because some parameter appearing in (iii) depends on a parameter defined in (iv). We use the notations: ∇f For simplicity sake, we assume that ∇ d f (x) > 0 (the idea easily adapts to k = d and to negative partial derivatives). Applying the implicit function theorem to the C 2 -function f , we obtain the existence of a neighborhood V of (x 1 , . . . , x d−1 ), a neighborhood U of x d and an increasing C 2 -function h such that the Assumption (ii): Let us prove that x ∈ R d ; f (x) ≥ 0 satisfies the Uniform Exterior Sphere property restricted to D: according to definition 3.1 in [3], we have to prove that there exists some positive α f such that, for each x ∈ D satisfying f (x) = 0, one has Let us fix x ∈ D on which f vanish. Taylor formula gives: for each y ∈ R d , with some c * ∈ [0; 1] depending on y and x. In particular, for y such that Assumption (iv): We have to prove the existence of some β 0 > 0 such that for each x ∈ ∂D there exists a unit vector l 0 x satisfying l 0 x .∇f (x) ≥ β 0 |∇f | for each constraint such that f (x) = 0. But this has already been done in remark 2.4 with β 0 = inf x∈∂D d(0, Conv(x)) and l 0 x = z |z| for some z with minimal norm in Conv(x).

Assumption (iii):
We have to prove that each set x ∈ R d ; f (x) ≥ 0 satisfies the Uniform Normal Cone property restricted to D with constant β f smaller than β 2 0 /2. Using Taylor formula for the derivative of f yields to ∇f (y) = ∇f (x)+D 2 f (x+c * (y−x))(y−x) for some c * ∈ [0; 1] depending on y and x. We obtain for x and y on which f vanish: )(y − x)| the right hand side is not smaller than: As a consequence, for any β f ∈]0, 1[ one can choose a δ f > 0 small enough such that for each x ∈ D satisfying f (x) = 0 and each y ∈ D satisfying f (y) = 0 and |y − x| ≤ δ f one has ∇f (y) This proves that x ∈ R d ; f (x) ≥ 0 satisfies the Uniform Normal Cone property restricted to D with any constant β f ∈]0; 1[, in particular with β f < β 2 0 /2 as requested. To complete the proof of theorems 2.2 and 2.3 for θ = I d , we proceed as in the proof of theorem 3.3 in [3], replacing the probability measure dµ(x) = 1 Z 1 I D (x)e −Φ(x) dx in that proof by the (σfinite but maybe unbounded) measure µ defined by dµ(x) = 1 I D (x)e −Φ(x) dx. Girsanov theorem yields the density of the distribution of the process with initial measure µ with respect to the distribution of reflected Brownian motion with Lebesgue measure as initial measure. Since both this density and the distribution of reflected Brownian motion starting from Lebesgue measure are time-reversal invariant, we obtain the reversibility of the solution with initial measure µ.
The case of oblique reflection: Let us check that the results obtained in the normal reflection case θ = I d transfer to the case of any invertible matrix θ. Using the notation X θ = θ −1 X, existence and uniqueness for equation (1) is equivalent to existence and uniqueness for (4) in the closure of the set θ −1 D = y ∈ R d ; ∀f ∈ F f (θy) > 0 with local times satisfying the condition L f (·) = · 0 1 I f (θX θ (s))=0 dL f (s).

An application : cluster of particles near the surface of an attractive planet
As an application of the previous results, let us study the configuration of a large number of particles around a fixed sphere called planet. As in [3], we consider spherical hard particles with random radii oscillating between a minimum and a maximum value. Each particle is driven by a Brownian motion, and undergoes the influence of the gravitational attraction generated by the planet. The motion is perturbated as the particles bump into each other and into the planet. What we are interested in is the existence and uniqueness of such a dynamics, and the equilibrium distribution of the particles: Do they tend to be scattered in a large area, or are they typically close together, all located beneath some altitude, with high probability ? That is, do they form some sort of "ocean" around the planet ? Using the results of section 2, we prove in proposition 3.3 that the particles eventually tend to cluster at the surface of the planet when the temperature (represented by the diffusion coefficient) tends to zero. More precisely, the planet is the d-dimensional closed ball B(0, R) centered at the origin with radius R. The space around this planet contains a large number n of particles, which are spheres with random positions and random radii. Each one will be represented by the position x i of its center in R d and the valuex i of its radius. Thus configurations are vectors x = (x 1 ,x 1 , . . . , x n ,x n ) in R n(d+1) .
To prevent negative radii, we enforcex i ∈ [r − , r + ] for some fixed values 0 < r − < r + . The random oscillations of the positions of the particles are not on the same scale as the random oscillations of their radii. We take this into account via the elasticity coefficientσ > 0 of their surface.
We assume that the gravity field ϕ generated by the planet is isotropic: it only depends on the norm |x| (more general cases would be mathematically tractable, but without physical meaning). As usual (see e.g. [1]) the gravitational attraction appears through a drift in the dynamics: Function ϕ is an increasing function which is C 2 on ]0; +∞[. The drift decreases with the distance, but not too fast (or equivalently the potential energy −ϕ increases not too fast) in the sense that ϕ ≤ 0 and lim inf ρ→+∞ ρϕ (ρ) > 0.
An important example in dimension d = 3 is ϕ(ρ) = C st ln(ρ) which gives the drift −ϕ (ρ) = − C st ρ corresponding to the usual gravitational acceleration −ϕ (ρ) = C st ρ 2 . At temperature θ > 0, the random motion of particles is modelized by the stochastic differential system: In this equation, vector (X i (·),X i (·)) 1≤i≤n represents the positions and radii of the n particles, the W i 's are independent R d -valued Brownian motions and theW i 's are independent one-dimensional Brownian motions, also independent from the W i 's. The amplitude of the Brownian oscillation of the positions depends on temperature θ, while the amplitude of the radii oscillation depends on both the temperature θ and the surface elasticityσ. Since ϕ is positive, the drift of X i is always directed toward the origin, i.e. toward the planet as expected. The local time L R i represents the repulsion received by the i th particle when it collides with the planet (impulsion away from the origin, in the direction of the unit vector X i R+X i ), and the local times L ij represent the collisions between particles, which tend to move the involved particles away from each other: unit direction Collisions between particles are symmetric (L ij ≡ L ji ).These local times also appear in the dynamics of the radii, because particles are like bubbles, and collisions have a deflating effect on them (the radii decrease). The local times L + i and L − i are here to comply with the conditionx i ∈ [r − , r + ] and give a positive (resp. negative) impulsion to the radius if it becomes too small (resp. too large). The impulsions are only given on the boundary of the set of allowed configurations, therefore the L R i 's, L + i 's, L − i 's and L ij 's should satisfy: The set of constraints is therefore: (the particles do not intersect the planet); • f + i (x) = r + −x i > 0 for 1 ≤ i ≤ n (radii are smaller than the maximum value); (radii are larger than the minimum value); . . , n} (particles do not overlap).
. We want to study the effect of the attraction of the planet on a "typical" particle configurations, i.e. at equilibrium.
Our aim, once the existence and uniqueness of the dynamics is proved, is to verify that, at low temperature, all particles cluster around the planet with high probability. In other words, there exists with high probability an interface between two regions around the planet: no particle over a certain altitude, and beneath this altitude a particle density so high that one cannot add any particle more (see figure 1). Figure 1: A configuration with an interface between high particle density and empty space.

Proposition 3.3
For each positive ε, let A ε be the set of configurations which do not pack into a minimal volume: A ε = {x ∈ D; ∃y ∈ D ∃k ≤ n s.t. ∀i = k y i = x i and |y k | < |x k | − ε} The probability that A ε occurs at equilibrium tends to zero as the temperature tends to zero: The end of the paper is devoted to the proofs of the three above propositions.

Proof of proposition 3.1
The are C ∞ and the corresponding set D of possible configurations is obviously a non-empty connected set. The first derivative of each constraint function is uniformly positive on its vanishing set because: In order to check the condition inf x∈∂D d(0, Conv(x)) > 0, in the form given in remark 2.4, we have to find some positive β 0 , and some non-vanishing vector v depending on x ∈ ∂D, such that From an intuitive point of view, v is the "shortest way to go back" into D from the point x on the boundary of D, i.e. the quickest way for colliding particles to go apart, for particles with maximum (resp. minimum) radius to become smaller (resp. larger) and for particles touching the planet to go away. C R will denote the indices of these globules: Intuitively, the best way to separate colliding particles is to move them away from the center of gravity of the cluster, that is to give each center x i an impulsion in the direction . . , n} is the cluster of colliding particles around x i (i.e. C i is the set containing i and all indices connected to i in the graph constructed on the vertices {1, . . . , n} by the edges j ∼ j ⇐⇒ |x j − x j | =x j +x j ). Similarly, the best way for particles touching the planet to go away is for each center x i with i ∈ C R to receive a small impulsion proportional to x i (this impulsion will also separate clusters of colliding particles). So a convenient v should be: Let us prove that the above vector v satisfies the desired inequalities: |∇f (x)| is bounded from below, uniformly in x ∈ ∂D and f ∈ F such that f (x) = 0. To complete the proof of proposition 3.1, it only remains to find a uniform upper bound for |v|: If C i is any cluster of colliding globules: Similarly, if C i is a cluster with at least one globule at distance R +x i of the origin: and the same upper bound holds for a sum over a union of such clusters. Consequently: Since the sum over {i; C i ∩ C R = ∅} is smaller than n(n − 1)(2n − 1), the norm of v is uniformly bounded from above as a function of x. This completes the proof.

Proof of proposition 3.2
We use theorem 2.2 with the n(d + 1) × n(d + 1) diagonal matrix θ which has n times the sequence (θ, . . . , θ, θσ) as its main diagonal entries. Since the constraints are compatible on R n(d+1) , for any C 2 -function Φ on R n(d+1) with bounded derivatives, the equation: has a unique strong solution in the closure of the set D defined by the constraints. Choosing The property L f (·) = · 0 1 I f (X(s))=0 dL f (s) implies that condition (E θ ) is satisfied for these new local times. Then the solution of equation (5) is the solution of (E θ ). Thanks to theorem 2.3, the measure 1 dx is a timereversible measure for the solution. To complete the proof, let us check that this measure can be renormalized as a probability measure for θ small enough.
-If, among the x k 's such that there exists y ∈ D satisfying |y k | < |x k | − ε and y i = x i for i = k, at least one has a norm smaller than K, then n i=1 ϕ(|x i |) > ϕ D + εϕ (K) > ϕ D ; -If no "ε-pushable" x k has a norm smaller than K, then all particles in x are at distance at least K from the origin (because K ≥ R + 2nr + + nε and it is impossible to fill so large a volume with only n particles without one of them to be "ε-pushable"). When the x k which has the largest norm is shifted at distance R + r + from the origin, becoming y k and giving the configuration y, we have So we obtain n i=1 ϕ(|x i |) > ϕ D + ε for all x ∈ A ε , with a positive ε equal to ε lim +∞ ϕ if this limit does not vanish and to min(εϕ (K), 3 ) otherwise.
The normalization constant of the probability measure µ θ is larger than: The denominator does not depend on θ. Dominated convergence theorem ensures that the numerator converges to zero when θ tends to 0. So we obtain lim θ→0 µ θ (A ε ) = 0.