Jingyi Fu

Benoit Perthame
email: benoit.perthame@sorbonne-universite.fr.

Min Tang

Fokker-Plank system for movement of micro-organism population in confined environment

Keywords: 2010 Mathematics Subject Classification. 35Q84, 35Q92, 35K65, 65C05, 65M06, 92C17 Piecewise deterministic Markov process, Fokker-Planck equation, Relative entropy, Monte-Carlo simulations, Semi-discrete numerical schemes

We consider self-propelled particles confined between two parallel plates, moving with a constant velocity while their moving direction changes by rotational diffusion. The probability distribution of such micro-organisms in confined environment is singular because particles accumulate at the boundaries. This leads us to distinguish between the probability distribution densities in the bulk and in the boundaries. They satisfy a degenerate Fokker-Planck system and we propose boundary conditions that take into account the switching between free-moving and boundary-contacting particles. Relative entropy property, a priori estimates and the convergence to an unique steady state are established. The steady states of both the PDE and individual based stochastic models are compared numerically.

Introduction

Self-propelled micro-organisms, such as C. Crescentus [START_REF] Li | Accumulation of swimming bacteibria near a solid surface[END_REF][START_REF] Li | Accumulation of microswimmers near a surface mediated by collision and rotational brownian motion[END_REF] and different types of Escherichia coli [START_REF] Berke | Hydrodynamic attraction of swimming microorganisms by surfaces[END_REF][START_REF] Bianchi | Holographic imaging reveals the mechanism of wall entrapment in swimming bacteria[END_REF][START_REF] Molaei | Failed escape: solid surfaces prevent tumbling of escherichia coli[END_REF] accumulate near boundaries. Micro-swimmers bump into various obstacles and boundaries, alter the motility and are trapped in the near surface region. This phenomena strongly relates to the first stages of biofilm formation and is important when micro-organisms are moving inside micro-fluidic devices [START_REF] Diluzio | Escherichia coli swim on the right-hand side[END_REF][START_REF] Figueroa-Morales | Living on the edge: transfer and traffic of E. coli in a confined flow[END_REF][START_REF] Galajda | A wall of funnels concentrates swimming bacteria[END_REF] or micro-structures [START_REF] Di Leonardo | Bacterial ratchet motors[END_REF][START_REF] Sokolov | Swimming bacteria power microscopic gears[END_REF]. Two possible mechanisms are attributed to the near surface behavior, one is long range hydrodynamic interactions between swimming cell and the surface [START_REF] Berke | Hydrodynamic attraction of swimming microorganisms by surfaces[END_REF][START_REF] Spagnolie | Hydrodynamics of self-propulsion near a boundary: predictions and accuracy of far-field approximations[END_REF], the other is steric interactions [START_REF] Li | Accumulation of swimming bacteibria near a solid surface[END_REF][START_REF] Li | Accumulation of microswimmers near a surface mediated by collision and rotational brownian motion[END_REF][START_REF] Ezhilan | Transport of a dilute active suspension in pressure-driven channel flow[END_REF]. Models based on either of them can lead to accumulation at the boundaries.

Cell-wall interactions have attracted a lot of interest from micro-biologists [START_REF] Berke | Hydrodynamic attraction of swimming microorganisms by surfaces[END_REF][START_REF] Bianchi | Holographic imaging reveals the mechanism of wall entrapment in swimming bacteria[END_REF][START_REF] Drescher | Fluid dynamics and noise in bacterial cell-cell and cell-surface scattering[END_REF][START_REF] Frymier | Three-dimensional tracking of motile bacteria near a solid planar surface[END_REF][START_REF] Li | Accumulation of swimming bacteibria near a solid surface[END_REF][START_REF] Li | Accumulation of microswimmers near a surface mediated by collision and rotational brownian motion[END_REF][START_REF] Molaei | Failed escape: solid surfaces prevent tumbling of escherichia coli[END_REF]. In the experiment in [START_REF] Berke | Hydrodynamic attraction of swimming microorganisms by surfaces[END_REF], the authors measured the distribution of non-tumbling E. coli cells swimming between two parallel plates. They measured the cell density as a function of distance away from the plate. However, which of the two mechanisms, hydrodynamic or steric interactions, dominates the reorientation is not clear.

Motivated by the experiment in [START_REF] Berke | Hydrodynamic attraction of swimming microorganisms by surfaces[END_REF] and the theoretical models in [START_REF] Li | Accumulation of microswimmers near a surface mediated by collision and rotational brownian motion[END_REF], we consider here self-propelled micro-organisms confined between two parallel plates. We assume they move with a constant velocity V , while their moving direction changes by rotational diffusion and we ignore their interaction with the external fluid. More precisely, we consider a 2-dimensional stochastic model for particles confined between two horizontal boundaries {y = ±L}. When the microorganism is away from boundaries, it swims at a constant speed V and changes its orientation with white noise; and when it touches boundary, it keeps changing orientation with the same rotational Brownian motion but stays at the boundary until it is oriented away from the boundary. The evolution of its vertical position Y t and the angle between its orientation and horizontal plane Θ t , is described by a stochastic differential equation (SDE) entering the class of Piecewise-Deterministic Markov Processes [START_REF] Davis | Piecewise-deterministic Markov processes: a general class of nondiffusion stochastic models[END_REF], Here Ω = {(y, θ)| -L < y < L, -π < θ < π} in which cells can move freely with constant velocity V , i.e. cells are at free-swimming phase (FSP). When cells move towards the boundary and touch it, they keep propelling themselves but are not able to move forward due to the steric forces of the boundary. This state is called boundary contacting phase (BCP) for which (Y t , Θ t ) ∈ {(y, θ)| y = ±L, 0 < ±θ < π} =: Ω ± .

       dY t = V sin Θ t dt
Then, cells return to FSP when their orientations are parallel to the boundary due to rotational diffusion. The mechanism of how cells change their movement directions near the boundary is more complex than rotational diffusion, for example the flow field close to the boundary plays an important role and some bacteria may tumble [START_REF] Bianchi | Holographic imaging reveals the mechanism of wall entrapment in swimming bacteria[END_REF][START_REF] Drescher | Fluid dynamics and noise in bacterial cell-cell and cell-surface scattering[END_REF][START_REF] Elgeti | Self-propelled rods near surfaces[END_REF][START_REF] Li | Accumulation of swimming bacteibria near a solid surface[END_REF]. No-flux boundary condition for Fokker-Planck model is proposed in [START_REF] Ezhilan | Transport of a dilute active suspension in pressure-driven channel flow[END_REF], and independently derived in [START_REF] Berlyand | A kinetic approach to active rods dynamics in confined domains[END_REF] as the vanishing inertia limit of a elastic collision boundary condition. However, the no-flux condition represent the balance of translational diffusion and self-propulsion in the wall-normal direction. While in our model, only rotational diffusion is considered for simplicity, and the collision is assumed to be inelastic. The SDE model (0.1) exhibits accumulations at the boundaries due to the direction persistence. The probability density function of the particles can be described by a Fokker-Plank equation inside |Y t | < L, but particles behave differently at the boundaries and stop moving until the Brownian motion generates an angle forcing cells to re-enter the domain. The purpose of this work is to propose appropriate boundary conditions for the Fokker-Plank equation in confined environment, so that it can describe the correct probability density distribution of the SDE model (0.1). The boundary conditions are deduced from the switching between FSP and BCP which gives a coupled system of three equations, one for free swimming cells, the other two for boundary-contacting cells at Y t = L and Y t = -L. Similar models have been proposed and studied before in [START_REF] Lee | Active particles under confinement: aggregation at the wall and gradient formation inside a channel[END_REF][START_REF] Wagner | Steady-state distributions of ideal active brownian particles under confinement and forcing[END_REF], in which semi-analytical solutions are derived for the steady state distribution. Here in this paper, we study the time evolutionary system and establish the relative entropy estimate and give its long term convergence to steady state solution.

The rest of this paper is organized as follows. In section 1, we introduce the Fokker-Plank model with appropriate boundary conditions when the micro-organisms are confined in two parallel plates. We give the relative entropy inequality and a priori estimate in section 2 and the weak convergence of the model system to steady state under certain assumptions in section 3. Finally, in section 4, we employ Von Mises iterations based on an entropy decreasing semi-discrete scheme to get the numerical solution of the steady state problem. We compare the numerical results with Monte-Carlo simulations.

The Fokker-Planck system

We consider the following degenerate Fokker-Planck system:

                 ∂p ∂t + V sin θ ∂p ∂y -D θ ∂ 2 p ∂θ 2 = 0, (y, θ) ∈ Ω, ∂p + ∂t -D θ ∂ 2 p + ∂θ 2 = V sin θp(t, L, θ), (L, θ) ∈ Ω + , ∂p - ∂t -D θ ∂ 2 p - ∂θ 2 = -V sin θp(t, -L, θ), (-L, θ) ∈ Ω -. (1.1a) (1.1b) (1.1c)
Here p(t, y, θ) represents the probability distribution of free-swimming cells at time t, position y and with moving direction θ; p ± (t, θ) are the probability density distributions of boundary-contacting cells at time t, position ±L and with moving direction ±θ ∈ (0, π). The diffusion in θ terms in (1.1a)-(1.1c) describe the rotational diffusion of cells while the y-derivative term models the free transport with velocity V sin θ in y direction when cells are at FSP. Notice the degeneracy of the system (the missing diffusion in y) which is standard in kinetic theory [START_REF] Villani | Hypocoercivity[END_REF][START_REF] Dolbeault | Hypocoercivity for linear kinetic equations conserving mass[END_REF][START_REF] Bessemoulin-Chatard | Hypocoercivity and diffusion limit of a finite volume scheme for linear kinetic equations[END_REF][START_REF] Hérau | Introduction to hypocoercive methods and applications for simple linear inhomogeneous kinetic models[END_REF]. The terms on the right hand side of (1.1b)-(1.1c) describe the transitions from FSP to BCP. Since when cells at FSP reach and move towards the boundary, they will switch to BCP, the changes p ± (t + dt, θ) dθ -p ± (t, θ) dθ, with dt being a small time interval, are composed by two parts, one is due to the rotational diffusion and the other is from cells at FSP near ±L. More precisely, let L 0 = L -V sin θ dt, the number of cells locating in the interval [L 0 , L) (or in (-L, -L 0 ] for p -) and moving forward in direction θ with velocity V is L L0 p(t, s, θ) ds. When dt → 0, we find that the source terms for p ± are as in (1.1b) and (1.1c). On the other hand, cells at BCP can switch to FSP when they are attached to the boundary at y = L (y = -L) and change their swimming direction from θ ∈ (0, π) (θ ∈ (-π, 0)) to 0 or π by rotational diffusion. Therefore the source term for p inside Ω is 0 but there is flux coming from the boundary point y = ±L, θ = 0 and y = L, θ = π, y = -L, θ = -π.

Boundary conditions. The Fokker-Planck system (1.1) comes with boundary conditions which we describe now. We begin with the equation for p(t, y, θ). On the one hand for θ, cells at FSP moving with direction θ = π and θ = -π are the same and thus periodic boundary conditions are used

p(t, y, -π) = p(t, y, π), ∂p ∂θ (t, y, -π) = ∂p ∂θ (t, y, π), -L < y < L. (1.2)
On the other hand for y we have a transport term and we need boundary conditions only when cells enter the domain (-L, L), and since there is no BCP cell in the entering boundary we impose Dirichlet condition

p(t, L, θ) = 0, θ ∈ (-π, 0), p(t, -L, θ) = 0, θ ∈ (0, π), (1.3)
Next, we turn to p ± and their boundary condition related to diffusion. In order to express that BCP cells p ± can only leave the boundary at θ = 0 or ±π, we specify Dirichlet boundary conditions

p + (t, 0) = p + (t, π) = 0, p -(t, 0) = p -(t, -π) = 0. (1.4)
It remains a difficulty related to the degeneracy of the Fokker-Planck equation. We need to express that BCP cells turning to the angles 0 or ±π will re-enter the domain, i.e., will change their status to FSP.

To determine p at ±L, we need four conditions for p(t, y, θ) at the points (±L, 0), (±L, ±π). Since cells switch from BCP to FSP immediately when θ changes from positive to negative at y = L, from (1.1b), (1.1a), we have that, at θ = 0,

∂ ∂t 0+ 0- p(t, L, θ) dθ -D θ ∂p ∂θ θ=0+ θ=0- = - ∂ ∂t 0+ 0- p + (t, θ) dθ -D θ ∂p + ∂θ θ=0+ θ=0- δ y=L p + (t, θ
) is defined for θ ∈ (0, π), but it can be extended to (-π, π) with p + (t, θ) = 0 for θ ∈ (-π, 0), which indicates that ∂p+ ∂θ (t, 0 -) = 0. Noting (1.4), we find

D θ ∂p + ∂θ (t, 0 + ) + D θ L L- ∂p ∂θ dy θ=0+ θ=0- = 0
Thus ∂p ∂θ jumps at the point (y, θ) = (L, 0) while it is continuous at (y, 0) with y ∈ (-L, L). Similar calculations give the jumps of ∂p ∂θ at (y, θ) = (-L, 0). Thus, we find ∂p ∂θ (t, y, θ)

θ=0+ θ=0- = - ∂p + ∂θ (t, 0 + )δ y=L + ∂p - ∂θ (t, 0 -)δ y=-L . (1.5)
On the other hand, due to the periodic boundary conditions for θ in (1.2), by letting p(t, y, θ + 2π) = p(t, y, θ) for θ ∈ [-π, 0), we can extend p(t, y, θ) from θ ∈ (-π, π) to θ ∈ (-π, 2π). Then ∂p ∂θ (t, y, θ) is continuous at θ = π and y ∈ (-L, L). Similar calculations as for θ = 0 yields ∂p ∂θ (t, y, θ)

θ=π+ θ=π- = ∂p + ∂θ (t, π -)δ y=L - ∂p - ∂θ (t, -π + )δ y=-L . (1.6)
Remark 1.1. Another way of taking into account the degeneracy of the Fokker-Plank equation is to add a source term on the right hand side of the density equation for cells at FSP such that

∂p ∂t + V sin θ ∂p ∂y -D θ ∂ 2 p ∂θ 2 = D θ S p , (1.7)
with the source term S p being given by

S p (t, y, θ) = δ y=L δ θ=0 ∂p + ∂θ (t, 0) -δ θ=π ∂p + ∂θ (t, π) + δ y=-L -δ θ=0 ∂p - ∂θ (t, 0) + δ θ=-π ∂p - ∂θ (t, -π) . (1.8)
The source term can be understood as follows. Cells at BCP can switch to FSP when they are attached to the boundary at y = L (y = -L) and change their swimming direction from θ ∈ (0, π) (θ ∈ (-π, 0)) to 0 or π by rotational diffusion. The flux of a diffusion process is driven by the negative gradient. Therefore, ∂p+ ∂θ (t, 0) is the number of cells that go outside of the interval (0, θ) by rotational diffusion at y = L, θ = 0. Cells switch from BCP to FSP immediately when θ becomes negative at y = L, which yields the source term δ y=L δ θ=0 ∂p+ ∂θ (t, 0) in (1.8). Other terms in (1.8) can be understood similarly.

Mass conservation, weak form. Finally, we equip this system (1.1) with the following nonnegative bounded initial data p(0, y, θ) = p 0 (y, θ), (y, θ) ∈ Ω;

p ± (0, θ) = p 0,± (θ), (±L, θ) ∈ Ω ± , (1.9)
The transfer terms and the boundary conditions are compatible with mass conservation, that means

M (t) = L -L π -π p dθ dy + π 0 p + dθ + 0 -π p -dθ := M (0). (1.10)
We impose that the initial data satisfy M (0) = 1. More generally, mass conservation is a special case of the weak form of the equation. We define admissible test functions as smooth functions ϕ(t, y, θ) defined in [0, T )×Ω which satisfy the periodicity condition ϕ(t, y, -π) = ϕ(t, y, π),

∂ϕ ∂θ (t, y, -π) = ∂ϕ ∂θ (t, y, π), ϕ + (t, θ) defined on [0, T )×[0, π] and ϕ -(t, θ) defined on [0, T )×[-π, 0] which satisfy ϕ(t, L, θ) = ϕ + (t, θ), for θ ∈ [0, π] and ϕ(t, -L, θ) = ϕ -(t, θ) for θ ∈ [-π, 0]
. Then, the weak form is written, for all admissible test functions, as

T 0 Ω p - ∂ϕ ∂t -V sin θ ∂ϕ ∂y -D θ ∂ 2 ϕ ∂θ 2 + T 0 π -π p + - ∂ϕ + ∂t -D θ ∂ 2 ϕ + ∂θ 2 + T 0 π -π p -- ∂ϕ - ∂t -D θ ∂ 2 ϕ - ∂θ 2 = = Ω p 0 ϕ(t = 0) + π -π [p 0,+ ϕ + (t = 0) + p 0,-ϕ -(t = 0)] - Ω p(T )ϕ(T ) + π -π [p + (T )ϕ + (T ) + p -(T )ϕ -(T )].
We notice that the formal integration by parts uses all the boundary conditions. This formulation and the a priori bounds in next section, also tells us that the probability distribution of the process (0.1) is composed of p(t, y, θ) which is absolutely continuous with respect to Lebesgue measure dy dθ, and singular parts, concentrated on the boundaries, p ± which are absolutely continuous with respect to dθ.

Relative entropy and a priori estimates

Since the model in (1.1) is linear, we take for granted existence and non-negativity of weak solutions even if its degeneracy and boundary conditions yield interesting mathematical questions. We derive here formally the relative entropy estimate and we deduce some a priori bounds (see [START_REF] Michel | General relative entropy inequality: an illustration on growth models[END_REF] for the general form of entropy structure).

The relative entropy uses the stationary state distributions q(y, θ) and q ± (θ) which satisfy:

                 -D θ ∂ 2 q ∂θ 2 + V sin θ ∂q ∂y = 0, (y, θ) ∈ Ω, -D θ d 2 q + dθ 2 = V sin θq(L, θ), (L, θ) ∈ Ω + , -D θ d 2 q - dθ 2 = -V sin θq(-L, θ), (-L, θ) ∈ Ω -, (2.1a)
(2.1b)

(2.1c)
with the boundary conditions

                         q(y, θ) = q(y, θ + 2π), ∂q ∂θ (y, θ) = ∂q ∂θ (y, θ + 2π), (y, θ) ∈ Ω, q(L, θ) = 0, θ ∈ (-π, 0), q(-L, θ) = 0, θ ∈ (0, π), ∂q ∂θ (y, θ) θ=0+ θ=0- = - ∂q + ∂θ (0 + )δ y=L + ∂q - ∂θ (0 -)δ y=-L , ∂q ∂θ (y, θ) θ=π+ θ=π- = ∂q + ∂θ (π -)δ y=L - ∂q - ∂θ (-π + )δ y=-L , q + (0) = q + (π) = 0, q -(0) = q -(-π) = 0. (2.2)
Since the model is linear, we can normalize (q, q ± ) and let it satisfy

L -L π -π q dθ dy + π 0 q + dθ + 0 -π q -dθ = 1. (2.3)
Thanks to the linearity of the model, and despite its degeneracy, we admit well-posedness and non-negativity of the solution to (2.1)-(2.3). Moreover, we also admit the strong maximum principle and get that q is strictly positive in the domain (-L, L) × (-π, π). Therefore, Ω 0 = {(y, θ) ∈ Ω| q(y, θ) = 0} is restricted to the incoming part of the boundary. Moreover, concavity of q ± guarantees their strictly positivity in Ω ± . Before establishing the relative entropy estimate, we define formally the relative gap ω and ω ± as

                 ω(t, y, θ) = p(t, y, θ) q(y, θ) -1, (y, θ) ∈ Ω \ Ω 0 , ω + (t, θ) = p + (t, θ) q + (θ) -1, (L, θ) ∈ Ω + , ω -(t, θ) = p -(t, θ) q -(θ) -1, (-L, θ) ∈ Ω -.
(2.4)

For the sake of simplicity, hereafter we denote ω(t, ±L, θ) by ω ±L (t, θ), p(t, ±L, θ) by p ±L (t, θ), and q(t, ±L, θ) by q ±L (θ).

We have the following relative entropy inequality:

Theorem 2.1. Consider a solution p to equation (1.1) with the boundary conditions in (1.2) and q is a solution to (2.1) with the boundary conditions in (2.2). For any convex function H ∈ C 2 (R), the relative gaps ω and ω ± as in (2.4) satisfy

d dt L -L π -π qH(ω) dθ dy + π 0 q + H(ω + ) dθ + 0 -π q -H(ω -) dθ = -D (1)
H (t) -D (2)
H (t) -D (3)
H (t) ≤ 0, (2.5)
where the dissipation from diffusion is, D

H (t) = D θ π 0 q + H (ω + ) ∂ω + ∂θ 2 dθ + D θ 0 -π q -H (ω -) ∂ω - ∂θ 2 dθ + D θ q L -L π -π qH (ω) ∂ω ∂θ 2 dθ dy, (2.6) (1)
the dissipation of switching from FSP to BCP is

D (2)
H (t) = π 0 V sin θq +L [H(ω +L ) -H(ω + ) -H (ω + )(ω +L -ω + )] dθ + 0 -π -V sin θq -L [H(ω -L ) -H(ω -) -H (ω -)(ω -L -ω -)] dθ, (2.7)
and the dissipation of switching from BCP to FSP is

D (3)
H (t) = -D θ [H(ω + ) -H(ω +L ) -H (ω +L )(ω + -ω +L )] ∂q + ∂θ π 0 -D θ [H(ω -) -H(ω -L ) -H (ω -L )(ω --ω -L )] ∂q - ∂θ 0 -π .
(2.8)

Proof. Substituting p + = q + ω + + q + into (1.1b) and using (2.1b), we get

q + ∂ω + ∂t =D θ (1 + ω + ) ∂ 2 q + ∂θ 2 + 2D θ ∂q + ∂θ ∂ω + ∂θ + D θ q + ∂ 2 ω + ∂θ 2 + V sin θq +L (1 + ω +L ) =D θ (ω + -ω +L ) ∂ 2 q + ∂θ 2 + 2 ∂q + ∂θ ∂ω + ∂θ + q + ∂ 2 ω + ∂θ 2 . . (2.9)
Multiplying both sides of (2.9) by H (ω + ) yields

q + D θ ∂H(ω + ) ∂t =H (ω + )(ω + -ω +L ) ∂ 2 q + ∂θ 2 + 2 ∂q + ∂θ ∂H(ω + ) ∂θ + q + H (ω + ) ∂ 2 ω + ∂θ 2 =H (ω + )(ω + -ω +L ) ∂ 2 q + ∂θ 2 + 2 ∂q + ∂θ ∂H(ω + ) ∂θ + q + ∂ 2 H(ω + ) ∂θ 2 -H (ω + ) ∂ω + ∂θ 2 .
(2.10)

Combining the above equation with

∂ 2 q + H(ω + ) ∂θ 2 = H(ω + ) ∂ 2 q + ∂θ 2 + 2 ∂q + ∂θ ∂H(ω + ) ∂θ + q + ∂ 2 H(ω + ) ∂θ 2 , we find ∂ ∂t -D θ ∂ 2 ∂θ 2 q + H ω + =D θ [H (ω + )(ω + -ω +L ) -H(ω + )] ∂ 2 q + ∂θ 2 -D θ q + H (ω + ) ∂ω + ∂θ 2 =V sin θq +L H(ω +L ) -V sin θq +L [H(ω +L ) -H(ω + ) -H (ω + )(ω +L -ω + )] -D θ q + H (ω + ) ∂ω + ∂θ 2 .
(2.11)

Similarly, the equation for q -H(ω -) is

∂ ∂t -D θ ∂ 2 ∂θ 2 q -H ω -= -V sin θq -L H(ω -L ) + V sin θq -L [H(ω -L ) -H(ω -) -H (ω -)(ω -L -ω -)] -D θ q -H (ω -) ∂ω - ∂θ 2 .
(2.12)

Substituting p = qω + q into (1.1a) and using (2.1a), we have

q ∂ω ∂t =(1 + ω) D θ ∂ 2 q ∂θ 2 -V sin θ ∂q ∂y + 2D θ ∂q ∂θ ∂ω ∂θ + D θ q ∂ 2 ω ∂θ 2 -V sin θq ∂ω ∂y =2D θ ∂q ∂θ ∂ω ∂θ + D θ q ∂ 2 ω ∂θ 2 -V sin θq ∂ω ∂y .
(2.13)

Multiplying both sides of the above equation by H (ω) gives

∂ qH(ω) ∂t =2D θ ∂q ∂θ ∂H(ω) ∂θ + D θ qH (ω) ∂ 2 ω ∂θ 2 -V sin θq ∂H(ω) ∂y =2D θ ∂q ∂θ ∂H(ω) ∂θ + D θ q ∂ 2 H(ω) ∂θ 2 -H (ω) ∂ω ∂θ 2 -V sin θq ∂H(ω) ∂y =D θ ∂ 2 qH(ω) ∂θ 2 -D θ H(ω) ∂ 2 q ∂θ 2 -D θ qH (ω) ∂ω ∂θ 2 -V sin θ ∂ qH(ω) ∂y -H(ω) ∂q ∂y =D θ ∂ 2 qH(ω) ∂θ 2 -V sin θ ∂ qH(ω) ∂y -D θ qH (ω) ∂ω ∂θ 2 .
(2.14)

The equality in (2.5) can be obtained by rearranging π-

0+ (2.11) dθ+ 0- -π+ (2.12) dθ+ L -L 0- -π+ (2.14) dθ dy+ L -L π- 0+
(2.14) dθ dy. The two terms D

H and D

H are obvious and D

H is from the diffusion terms. After integrating and using (1.3), the second order derivative in θ terms becomes

D θ ∂ q + H(ω + ) ∂θ θ=π- θ=0+ + ∂ q -H(ω -) ∂θ θ=0- θ=-π+ + L -L ∂ qH(ω) ∂θ θ=π- θ=0+ dy + L -L ∂ qH(ω) ∂θ θ=0- θ=-π+
dy .

(2.15)

From the definitions of ω, ω ± in (2.4), we have

p = ωq + q, p ± = ω ± q ± + q ± ,
Then using the continuity of q and ω for (y, θ) ∈ (-L, L) × {0}, (1.3), (1.5) and second and third boundary conditions in (2.2), we find

L -L ∂ qH(ω) ∂θ θ=0+ θ=0- dy = L -L H (ω) ∂(ωq) ∂θ θ=0+ θ=0- dy + L -L H(ω) -ωH (ω) ∂q ∂θ θ=0+ θ=0- dy = L -L H (ω) ∂p ∂θ θ=0+ θ=0- dy + L -L H(ω) -(ω + 1)H (ω) ∂q ∂θ θ=0+ θ=0- dy = - L -L H ω(t, y, 0) ∂p + ∂θ (t, 0 + )δ y=L dy + L -L H ω(t, y, 0) ∂p - ∂θ (t, 0 -)δ y=-L dy - L -L H ω(t, y, 0) -ω(t, y, 0) + 1 H ω(t, y, 0) ∂q + ∂θ (0 + )δ y=L dy + L -L
H ω(t, y, 0) -ω(t, y, 0) + 1 H ω(t, y, 0) ∂q - ∂θ (0 -)δ y=-L dy which we also write

= -H ω +L (t, 0) q + (0 + ) ∂ω + ∂θ (t, 0 + ) + ω + (t, 0 + ) + 1 ∂q + ∂θ (0 + ) + H ω -L (t, 0) q -(0 -) ∂ω - ∂θ (t, 0 -) + ω -(t, 0 -) + 1 ∂q - ∂θ (0 -) -H ω +L (t, 0) -ω +L (t, 0) + 1 H ω +L (t, 0) ∂q + ∂θ (0 + ) + H ω -L (t, 0) -ω -L (t, 0) + 1 H ω -L (t, 0) ∂q - ∂θ (0 -) =-H ω +L (t, 0) q + (0 + ) ∂ω + ∂θ (t, 0 + ) -ω + (t, 0 + ) -ω +L (t, 0 + ) H ω +L (t, 0) ∂q + ∂θ (0 + ) -H ω +L (t, 0) ∂q + ∂θ (0 + ) +H ω -L (t, 0) q -(0 -) ∂ω - ∂θ (t, 0 -) + ω -(t, 0 -) -ω -L (t, 0 -) H ω -L (t, 0) ∂q - ∂θ (0 -) + H ω -L (t, 0) ∂q - ∂θ (0 -).
Then using (1.4) and (2.2), we obtain

- ∂ q + H(ω + ) ∂θ (t, 0 + ) + ∂ q -H(ω -) ∂θ (t, 0 -) - L -L ∂ qH(ω) ∂θ θ=0+ θ=0- dy = H ω +L (t, 0 + ) -H ω + (t, 0 + ) -H (ω +L (t, 0 + )) ω +L (t, 0 + ) -ω + (t, 0 + ) ∂q + ∂θ (0 + ) -H ω -L (t, 0 -) -H ω -(t, 0 -) -H (ω -L (t, 0 -)) ω -L (t, 0 -) -ω -(t, 0 -) ∂q - ∂θ (0 -).
The other terms in (2.15) can be reformulated similarly and we obtain D

H . Due to the positivity of q ± and the boundary conditions that q ± (0) = q ± (±π) = 0, we have ∂q+ ∂θ (0), ∂q- ∂θ (-π) > 0 and ∂q+ ∂θ (π), ∂q- ∂θ (0) < 0. Then convexity of H and the positivity q ± , q give the positivity of D

(1)

H , D (2)
H , D (3)
H . This concludes the derivation of the relative entropy inequality.

The immediate consequence of relative entropy inequality (2.5) is the following a priori estimate which indicates that if initially p, p ± are bounded, they will not blow up.

Corollary 2.2. Assume there exists a constant Λ > 0 such that

ω(0, •, •) L ∞ (Ω) ≤ Λ, ω ± (0, •) L ∞ (Ω±) ≤ Λ, (2.16
)

then for all t ≥ 0, ω(t, •, •) L ∞ (Ω) ≤ Λ, ω ± (t, •) L ∞ (Ω±) ≤ Λ. (2.17)
Proof. By choosing the convex function

H(x) = [(x -Λ) + ] 4 ∈ C 2 (R) in (2.5), we find d dt L -L π -π q(ω -Λ) 4 + dθ dy + π 0 q + (ω + -Λ) 4 + dθ + 0 -π q -(ω --Λ) 4 + dθ ≤ 0.
Then, we conclude that

L -L π -π q(ω(T, y, θ) -Λ) 4 + dθ dy + π 0 q + (ω + (T, θ) -Λ) 4 + dθ + 0 -π q -(ω -(T, θ) -Λ) 4 + dθ ≤ L -L π -π q(ω(0, y, θ) -Λ) 4 + dθ dy + π 0 q + (ω + (0, θ) -Λ) 4 + dθ + 0 -π q -(ω -(0, θ) -Λ) 4 + dθ = 0.
Since q is strictly positive for all (y, θ) in Ω except a zero measure subset, and q ± are strictly positive for θ ∈ Ω ± , we conclude the proof.

Long time convergence

Another consequence of the relative entropy is the weak convergence to steady state for solutions of (1.1), under certain assumptions. The approach is as follows. For k ∈ N, we define the time shifts

ω k (t, y, θ) = ω(t + k, y, θ), (y, θ) ∈ Ω; ω k,± (t, θ) = ω ± (t + k, θ), (±L, θ) ∈ Ω ± . (3.1)
From Corollary 2.2, {ω k }, {ω k,± } are uniformly bounded in L ∞ . By weak compactness of L ∞ space, there exists a subsequence {n k } k∈N and limit functions ω

∞ ∈ L ∞ R × Ω , ω ∞,± ∈ L ∞ R × Ω ± , such that for all functions φ(t, y, θ) in L 1 (R × Ω) and φ ± (t, θ) in L 1 (R × Ω ± ), we have for all T ∈ R, lim n k →∞ ∞ T L -L π -π φω n k dθ dy dt = ∞ T L -L π -π φω ∞ dθ dy dt, lim n k →∞ ∞ T ±π 0 φ ± ω n k ,± dθ dt = ∞ T ±π 0 φ ± ω ∞,± dθ dt. (3.2)
We prove in the subsequent part that ω ∞ , ω ∞,± are zero functions.

Theorem 3.1. Assume that the initial data ω(0, •, •) and ω ± (0, •) are uniformly bounded in L ∞ , that is (2.16) are satisfied. Then

ω ∞ = 0, ω ∞,± = 0, a.e. (3.3)
and thus for any test function φ(y, θ)

∈ L 1 (Ω), φ + (θ) ∈ L 1 (0, π), φ -(θ) ∈ L 1 (-π, 0), lim t→+∞ π -π L -L ω(t, y, θ)φ(y, θ) dy dθ = lim t→+∞ 0 -π ω -(t, θ)φ -(θ) dθ = lim t→+∞ π 0 ω + (t, θ)φ + (θ) dθ = 0. (3.4)
Proof. The proof is divided into five steps, we first show that ω ∞ , ω ∞,± are independent of θ by the relative entropy estimate, then prove their independence of time and y by the equation and finally we use the conservation of mass to show that ω ∞ , ω ∞,± have to be zero almost everywhere.

First step: strong convergence from relative entropy inequality (2.5). Using the relative entropy H(x) = x 2 in (2.5), we obtain

I(t) = L -L π -π qω 2 dθ dy + π 0 q + ω 2 + dθ + 0 -π q -ω 2 -dθ and dI(t) dt ≤ 0.
Since the initial data are bounded as in (2.16), I(0) ≤ (4πL + 2π)Λ 2 < ∞. Then, since I(t) is nonnegative, it converges to a limit when t → ∞. Therefore, lim τ →∞ ∞ τ dI(t) dt dt = 0 and we find lim

τ →∞ ∞ τ D (1)
H + D

(2)

H + D (3)
H [ω, ω ± ] dt = 0.

For ω n k and ω n k ,± , this gives us lim

n k →∞ ∞ T D (1)
H + D (2)
H + D

H [ω n k , ω n k ,± ](t) dt = lim n k →∞ ∞ T +n k D (3)
H + D (1)
H + D

H [ω, ω ± ](t) dt = 0. (3)
Consequently, due to the non-negativity of each term in D

H , D

H , D

H in (2.6)-(2.8), we can extract a subsequence, still denoted by n k such that

                           ∂ω ∞ ∂θ (t, y, θ) = 0, a.e. in [T, ∞) × Ω, ∂ω ∞,± ∂θ (t, θ) = 0, a.e. in [T, ∞) × Ω ± , lim n k →∞ ω n k ,±L (t, θ) -ω n k ,± (t, θ) = 0, a.e. in [T, ∞) × Ω ± , lim n k →∞ ω n k ,±L (t, 0) -ω n k ,± (t, 0) = 0, a.e. in [T, ∞), lim n k →∞ ω n k ,±L (t, ±π) -ω n k ,± (t, ±π) = 0, a.e. in [T, ∞). (3.5a) (3.5b) (3.5c) (3.5d) (3.5e)
Second step: the limits ω ∞,± are independent of time. Since q ± are independent of time, (2.9) indicates that, for all k ∈ N, ω n k and ω n k ,+ satisfy

q + ∂ω n k ,+ ∂t = D θ (ω n k ,+ -ω n k ,+L ) ∂ 2 q + ∂θ 2 + 2 ∂q + ∂θ ∂ω n k ,+ ∂θ + q + ∂ 2 ω n k ,+ ∂θ 2 , θ ∈ (0, π). (3.6)
We choose compactly supported test function φ

+ (t, θ) ∈ L 1 ([T, ∞) × Ω + ) C ∞ ([T, ∞) × Ω + )
, multiply both sides of (3.6) by the test function φ + , integrate on [T, ∞) × (0, π), and find

∞ T π 0 q + ω n k ,+ ∂φ + ∂t dθ dt = D θ ∞ T π 0 ∂ω n k ,+ ∂θ ∂(q + φ + ) ∂θ dθ dt -2D θ ∞ T π 0 ∂q + ∂θ ∂ω n k ,+ ∂θ φ + dθ dt + D θ ∞ T π 0 (ω n k ,+L -ω n k ,+ ) ∂ 2 q + ∂θ 2 φ + dθ dt. (3.7)
Let n k → ∞ and pass to limit, from (3.5), the right hand side of (3.7) vanishes and hence, for any compactly supported

φ + (t, θ) ∈ L 1 ([T, ∞) × Ω + ) C ∞ ([T, ∞) × Ω + ), ∞ T π 0 q + ∂ω ∞,+ ∂t φ + dθ dt = - ∞ T π 0 q + ω ∞,+ ∂φ + ∂t dθ dt = 0. (3.8)
Therefore, ω ∞,+ is independent of t. Similarly, ω ∞,-is independent of t as well.

Third step: ω ∞ is independent of t and y. Since q is independent of time, from (2.13), for all k ∈ N, we have

q ∂ω n k ∂t = 2D θ ∂q ∂θ ∂ω n k ∂θ + D θ ∂ 2 ω n k ∂θ 2 q -V sin θq ∂ω n k ∂y . (3.9)
Letting n k → ∞ in (3.9) and using (3.5a), gives

q ∂ω ∞ ∂t = -V sin θq ∂ω ∞ ∂y .
Thanks to the strictly positivity of q inside Ω, we divide both sides of the above equation by q, multiply it by a test function φ(t, y) ∈ C 2 ([0, ∞) × (-L, L)) and find

∞ 0 L -L φ ∂ω ∞ ∂t dy dt = -V sin θ ∞ 0 L -L φ ∂ω ∞ ∂y dy dt. (3.10)
The right hand side of (3.10) depends on θ while the left hand side does not, which indicates that both sides are zero. Therefore, ω ∞ is independent of t and y. So far we proved ω ∞ and ω ∞,± are both constant almost everywhere in their domain of definition respectively.

Forth step: show that ω ∞ = ω ∞,± . Multiplying (3.9) by a test function φ ∈ C 2 ([0, +∞) × Ω) that is compact supported in t and θ but not in y and integrating on [T, ∞) × Ω yields

∞ T L -L π -π qω n k ∂φ ∂t dθ dy dt = ∞ T π -π V sin θ L -L qφ ∂ω n k ∂y dy dθ dt + D θ ∞ T L -L π -π ∂ω n k ∂θ ∂(qφ) ∂θ dθ dy dt -2D θ ∞ T L -L π -π ∂q ∂θ ∂ω n k ∂θ φ dθ dy dt +D θ ∞ T (ω n k ,+ -ω n k ,L ) ∂q + ∂θ φ(t, L, θ) π 0 dt + D θ ∞ T (ω n k ,--ω n k ,-L ) ∂q - ∂θ φ(t, -L, θ) 0 -π dt.
(3.11)

Here the last term on the right hand side comes from the boundary conditions for p and q in (1.5), (1.6) and (2.2). Besides, using integration by part, we find

∞ T π -π V sin θ L -L qφ ∂ω n k ∂y dy dθ dt = - ∞ T π -π V sin θ L -L ω n k ∂(qφ) ∂y dy dθ dt + ∞ T π -π V sin θ(φqω n k ) L -L dθ dt = - ∞ T π -π V sin θ L -L ω n k ∂(qφ) ∂y dy dθ dt + ∞ T π 0 V sin θφ L ω n k ,+ q +L dθ dt - ∞ T 0 -π V sin θφ -L ω n k ,-q -L dθ dt - ∞ T π 0 V sin θφ L (ω n k ,+ -ω n k ,L )q +L dθ dt + ∞ T 0 -π V sin θφ -L (ω n k ,--ω n k ,-L )q -L dθ dt.
Therefore using (3.5), when

n k → ∞, ∞ T L -L π -π ω ∞ ∂(qφ) ∂t dθ dy dt = - ∞ T L -L π -π ω ∞ V sin θ ∂(qφ) ∂y dθ dy dt + ∞ T π 0 V sin θφ L ω ∞,+ q +L dθ dt - ∞ T 0 -π V sin θφ -L ω ∞,-q -L dθ dt = -(ω ∞ -ω ∞,+ ) ∞ T L -L π 0 V sin θ ∂(qφ) ∂y dθ dy dt -(ω ∞ -ω ∞,-) ∞ T L -L 0 -π V sin θ ∂(qφ)
∂y dθ dy dt.

(3.12)

Here the last equation has used the second boundary condition in (2.2). We can choose φ(t, y, θ) that is compactly supported in (0, +∞) for t and (0, π) for θ and satisfies

ω ∞ | Ω = ω ∞,± | Ω± = C, a.e.
By the conservation of mass,

∀n k ∈ N, ω n k satisfies L -L π -π qω n k dθ dy + π 0 q + ω n k ,+ dθ + 0 -π q -ω n k ,-dθ = 0, (3.13) which indicates that L -L π -π qω ∞ dθ dy + π 0 q + ω ∞,+ dθ + 0 -π q -ω ∞,-dθ = 0. (3.14)
Therefore, the constant C has to be zero. Here C in the last inequality is some positive constant and the inequality is easy to prove by using integral by parts and Corollary 2.2. Therefore, u(t) is a Lipschitz function in time and we have shown that it converges weakly to 0, then it converges strongly. Strong convergence of ω ± in time can be derived similarly.

Numerical simulations

The steady state problem (2.1) and (2.2) is a linear eigenvalue problem and only unique under the constraint (2.3). We can simulate it by fully implicit scheme based on an entropy decreasing semi-discrete scheme for the time evolution problem (1.1). The discrete entropy estimate of the semi-discrete scheme is established in the Appendix and we observe the convergence to steady state numerically. More sophisticated methods, with exponential decay, have been proposed, not including the difficulty related to the boundary, see [START_REF] Dujardin | Coercivity, hypocoercivity, exponential time decay and simulations for discrete Fokker-Planck equations[END_REF] and the references therein.

Dimensionless parameters are used in our numerical solver and simulations. Rotational diffusion coefficient is D θ = 1, moving speed is V = 20, and distance between the two horizontal plates are 2L = 20.

Description of the finite difference scheme

We consider a uniform mesh for y and θ in the rectangular domain Ω = [-L, L] × [-π, π]. Let the mesh sizes of y and θ be respectively ∆y = L/I and ∆θ = π/J where I, J are two integers. The index sets are

V y = {1, 2, • • • , 2I}, V θ = {1, 2, • • • , 2J}, Ṽθ = {0, 1, • • • , 2J},
and

V θ+ = {J + 1, J + 2, • • • , 2J -1}, V θ-= {1, 2, • • • , J -1}, V θ0 = {0, J, 2J}.
As in Figure 1, the grids inside the computational domain are:

(y i , θ j ) = (i -I -1/2)∆y, (j -J)∆θ , i ∈ V y , j ∈ Ṽθ , (4.1)
and the nodes at the boundaries are:

(±L, θ j ), j ∈ V θ± V θ0± .
The unknowns are p ±,j ≈ p ± (t, θ j ) and

p i,j ≈ 1 ∆y yi+∆y/2
yi-∆y/2 p(t, y i , θ j ) dy.

We use upwind discretization for ∂ y p, central finite difference discretization for ∂ 2 p ∂θ 2 , and the discretizations of

y θ -L L π -π 1 2I 0 2J
Figure 1: Stencil of semi-discretized scheme. Here nodes for p(t, y i , θ j ) inside Ω are black, while nodes for p ± are blue and red respectively.

boundary conditions are chosen to preserve the total mass. The semi-discretized scheme writes

                                 dp i,j dt + V sin θ j p i,j -p i-1,j ∆y -D θ p i,j-1 + p i,j+1 -2p i,j ∆θ 2 = 0, i ∈ V y , j ∈ V θ+ , dp i,j dt + V sin θ j p i+1,j -p i,j ∆y -D θ p i,j-1 + p i,j+1 -2p i,j ∆θ 2 = 0, i ∈ V y , j ∈ V θ-, dp i,j dt -D θ p i,j-1 + p i,j+1 -2p i,j ∆θ 2 = 0, i ∈ V y \ {1, 2I}, j ∈ V θ0 , dp +,j dt -D θ p +,j-1 + p +,j+1 -2p +,j ∆θ 2 = V sin θ j p 2I,j , j ∈ V θ+ dp -,j dt -D θ p -,j-1 + p -,j+1 -2p -,j ∆θ 2 = -V sin θ j p 1,j , j ∈ V θ-. (4.2a) (4.2b) (4.2c) (4.2d) (4.2e)
From the periodic boundary conditions in (1.2) for p(t, y, θ) in θ, we use

p i,-1 = p i,2J-1 , p i,0 = p i,2J , p i,2J+1 = p i,1 , ∀i ∈ V y ;
and due to the zero boundary conditions in (1.3) and (1.4), we use

p 0,j = 0, j ∈ V θ+ , p 2I+1,j = 0, j ∈ V θ-, p +,J = p +,2J = 0, p -,0 = p -,J = 0. ( 4.3)
At the continuous level, p is integrated over the whole space, the conservation of mass can be obtained because the boundary conditions of p and p ± cancel each other. On the other hand, at the discrete level, if we replace evolution equations by boundary conditions at the four points (1, J), (1, 0), (2I, J), (2I, 2J), then the information about how density function evolves at these point are lost. Hence we prefer to discretize the equation (1.7)-(1.8) mentioned in Remark 1.1, which is equivalent to our model.

Approximate equation at the point (-L, 0) reads

dp 1,J dt -D θ p 1,J-1 + p 1,J+1 -2p 1,J ∆θ 2 = D θ p -,J -p -,J-1 ∆θ 2 ∆y . (4.4a)
Similarly, at (-L, -π), (L, 0), (-L, π), we have

dp 1,1 dt -D θ p 1,-1 + p 1,1 -2p 1,0 ∆θ 2 = -D θ p -,1 -p -,0 ∆θ 2 ∆y , (4.4b)
dp 2I,J dt -D θ p 2I,J-1 + p 2I,J+1 -2p 2I,J ∆θ 2 = -D θ p +,J+1 -p +,J ∆θ 2 ∆y , (4.4c)
dp 2I,2J dt -D θ p 2I,2J+1 + p 2I,2J-1 -2p 2I,2J ∆θ 2 = D θ p +,2J -p +,2J-1 ∆θ 2 ∆y . (4.4d)
Non-negative initial data p i,j (0) and p ±,j (0) that satisfy ∆θ∆y (i,j)∈Vy×V θ p i,j (0) + ∆θ

j∈V θ + p +,j (0) + ∆θ j∈V θ - p -,j (0) = 1. (4.5)
is used and it is easy to verify that the scheme (4.2)-(4.4) conserves the discretized total mass defined by ∆θ∆y (i,j)∈Vy×V θ p i,j (t) + ∆θ

j∈V θ + p +,j (t) + ∆θ j∈V θ - p -,j (t)
In Appendix, we proved the discrete relative entropy estimate for (4.2)-(4.5). In the subsequent part, we can see numerically that when t → ∞, the solution converges to the solution of the following discretized steady state equation:

                               V sin θ j q i,j -q i-1,j ∆y -D θ q i,j-1 + q i,j+1 -2q i,j ∆θ 2 = 0, i ∈ V y , j ∈ V θ+ , V sin θ j q i+1,j -q i,j ∆y -D θ q i,j-1 + q i,j+1 -2q i,j ∆θ 2 = 0, i ∈ V y , j ∈ V θ-, -D θ q i,j-1 + q i,j+1 -2q i,j ∆θ 2 = 0, i ∈ V y \ {1, 2I}, j ∈ V θ0+ ∪ V θ0-, -D θ q +,j-1 + q +,j+1 -2q +,j ∆θ 2 = V sin θ j q 2I,j , j ∈ V θ+ -D θ q -,j-1 + q -,j+1 -2q -,j ∆θ 2 = -V sin θ j q 1,j , j ∈ V θ-. (4.6a) (4.6b) (4.6c) (4.6d) (4.6e)
with boundary conditions q i,-1 = q i,2J-1 , q i,0 = q i,2J , q i,1 = q i,2J+1 , i ∈ V y ; q 0,j = 0, j ∈ V θ+ , q 2I+1,j = 0; j ∈ V θ-; q +,J = q +,2J = 0; q -,0 = q -,J = 0; q 1,-1 + q 1,1 -2q 1,0 + q -,1 ∆y = 0, q 1,J-1 + q 1,J+1 -2q 1,J + q -,J-1 ∆y = 0, q 2I,J-1 + q 2I,J+1 -2q 2I,J + q +,J+1 ∆y = 0, q 2I,2J-1 + q 2I,2J+1 -2q 2I,2J + q +,2J-1 ∆y = 0, (

and total mass normalization ∆θ∆y (i,j)∈Vy×V θ q i,j + ∆θ

j∈V θ + q +,j + ∆θ j∈V θ - q -,j = 1. (4.8)
We have used upwind discretization for ∂ y p, which has first order convergence w.r.t y; second order centered finite difference discretization is employed for ∂ 2 p ∂θ 2 , but due to first order discretization for derivatives in the source term (1.8) as in (4.4), the scheme has first order convergence w.r.t θ. We define the following weighted L 2 norm to measure the numerical errors:

(e, e + , e -) 2 In Table 1, numerical errors of q i,j calculated with different mesh sizes are given, where the reference solution is computed with a fine mesh 256 × 256. From Figure 2, first order convergence can be observed when I = J.

q, q + , q -) -(q ref , q ref + , q ref -) L 2 w .
The weighted L 2 norm of the numerical errors as in (4.9)

Figure 2: The convergence order of the numerical errors displayed in Table 1 when I = J.

Numerical comparisons between the SDE and Fokker-Planck models

We perform Monte-Carlo simulations for the SDE model in (0.1) with Euler-Maruyama-like scheme [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equations[END_REF]. The computational domain is Ω = {(y, θ)| -L ≤ y ≤ L, -π ≤ θ ≤ π} and we track the trajectories of K = 5000000 cells. Each cell is represented by its position y i , orientation θ i . The initial y i for all cells are uniformly distributed on (-L, L), their initial orientation θ i are uniformly distributed on [-π, π]. Let ∆t be the time step. At each step we evolve (y i , θ i ) (i = 1, • • • , K) by the following calculations:

1) Update θ i by dΘ t = √ 2D θ dW t . For each i, generate a random increment I i with normal distribution independent of previous ones. Suppose that θ i + √ 2D θ ∆tI i ∈ [2kπ -π, 2kπ + π) for some integer k, then set the new orientation θ i to be θ i + √ 2D θ ∆tI i -2kπ.

2) Update the position y i . If y i = ±L, and ± sin(θ i+1 ) ≥ 0, then y i does not change. If

|y i | < L, |y i +V sin(θ i+1 )| > L,
then the new y i is set to be sign(y i + V sin(θ i ))L. Otherwise, evolve y i by y i + V sin(θ i+1 ).

We take ∆t = 0.01 and count the number of particles inside each cell of a 200 × 200 uniform mesh and normalize them by the total number of particles. It has been found in simulations that, the probability distribution is close to steady state when T = 4. Finite difference scheme described in section 4.1 is used to compute the numerical solution of Fokker-Planck model. In practice, we take I = J = 100 such that the numerical results of two model have same resolution. Fully implicit scheme is performed on the semi-discretized equation (4.2), and we employed GMRES to solve the resulting algebraic system. The initial condition is consistent with SDE model, i.e. p 0 i,j = 1 4πL for all i ∈ V y , j ∈ V θ . The same as the SDE model, solution of the Fokker-Planck equation is close to steady state when T = 4.

The time evolutions of bulk density p(t, y, θ) by numerical simulations based on both SDE and Fokker-Planck models are shown in Figure 3. The results of the two models are close to each other qualitatively. There exist two wells when y = -L, θ ∈ (0, π) (y = L, θ ∈ (-π, 0)), since cells with positive (negative) θ tend to leave the boundary. When t is larger, the density peaks at (±L, 0), (±L, ±π), since cells switch from BCP to FSP only at these points. We can see density in most of bulk area is uniform, and decreases as t increases.

To verify the behavior of p(t, y, θ) near the boundaries, p 1,j and p 2I,j are plotted in Figure 4. As t becomes larger, the peaks at (±L, 0), (±L, ±π) become sharper. Due to the limitation of numerical resolutions, the height of peaks stop increasing when t reaches 4. It is not clear from the numerics if there exist singularities in these points. The probability density distributions of cells at BCP on both left and right boundaries are plotted in Figure 4 as well, we can see that their maximum increase with time until they reach a constant.

Conclusion

Motivated by present biophysics experiments, we have considered the stochastic process which describes self-propelled confined micro-organisms moving with constant speed, turning with rotational Brownian motion and attaching to the boundaries. Its probability distribution is singular and composed of the bulk density p(t, y, θ) and boundary densities p ± (t, θ). These functions satisfy a degenerate Fokker-Planck system for which we have established relative entropy properties, a priori estimates and long term convergence to a steady-state. We have finally compared Monte-Carlo simulations for the stochastic process with the solution of the Fokker-Planck system using a finite difference scheme. Despite the complexity and singularities of the solution, a good fit is observed.

Several questions are left opened as the regularity properties of the time dependent solution and of the steady state solution. In particular, the type of singularity at the exit points of the boundary is challenging, It would be also interesting to prove a rate of convergence to steady state and decide it is exponential or not. There are also several where equality holds if and only if ω i,j = 0 for any i ∈ V y , j ∈ V θ and ω ±,j = 0 for any j ∈ V θ,± .

Proof. As in section 2, we divide the proof in several steps. First step: estimate of relative entropy at BCP. Substitute p +,j by q +,j ω +,j +q +,j in (4.2d) and employ (4.6d), we deduce that q +,j dω +,j dt = D θ ∆θ 2 [(q +,j-1 ω +,j-1 + q +,j+1 ω +,j+1 -2q +,j ω +,j ) + (q +,j-1 + q +,j+1 -2q +,j )] + V sin θ j q 2I,j (ω 2I,j + 1)

= D θ ∆θ 2 [(q +,j-1 ω +,j-1 + q +,j+1 ω +,j+1 -2q +,j ω +,j ) -ω 2I,j (q +,j-1 + q +,j+1 -2q +,j )] = D θ ∆θ 2 [q +,j+1 (ω +,j+1 -ω +,j ) -q +,j-1 (ω +,j -ω +,j-1 ) + (ω +,j -ω 2I,j )(q +,j-1 + q +,j+1 -2q +,j )]. (A.3)
Notice that q +,j-1 H(ω +,j-1 ) + q +,j+1 H(ω +,j+1 ) -2q +,j H(ω +,j ) =q +,j+1 [H(ω +,j+1 ) -H(ω +,j )] -q +,j-1 [H(ω +,j ) -H(ω +,j-1 )] + H(ω +,j )(q +,j-1 + q +,j+1 -2q +,j ), (A.4)

q i,j dH(ω i,j ) dt - D θ ∆θ 2 [q i,j-1 H(ω i,j-1 ) + q i,j+1 H(ω i,j+1 ) -2q i,j H(ω i,j )] = - V sin θ j ∆y q i-1,j H (ω i,j )(ω i,j -ω i-1,j ) + 2D θ ∆θ 2 q i,j H(ω i,j ) + D θ ∆θ 2 q i,j+1 [H (ω i,j )(ω i,j+1 -ω i,j ) -H(ω i,j+1 )] - D θ ∆θ 2 q i,j-1 [H (ω i,j )(ω i,j -ω i,j-1 ) + H(ω i,j-1 )] = - V sin θ j ∆y [q i,j H(ω i,j ) -q i-1,j H(ω i-1,j )] - V sin θ j ∆y q i-1,j [H (ω i,j )(ω i,j -ω i-1,j ) -H(ω i,j ) + H(ω i-1,j )] - D θ ∆θ 2 q i,j+1 [H(ω i,j+1
) -H(ω i,j ) -H (ω i,j )(ω i,j+1 -ω i,j )] -D θ ∆θ 2 q i,j-1 [H (ω i,j )(ω i,j -ω i,j-1 ) -H(ω i,j ) + H(ω i,j-1 )] -V sin θ j ∆y (q i-1,j -q i,j )H(ω i,j ) -D θ ∆θ 2 (q i,j-1 + q i,j+1 -2q i,j )H(ω i,j ).

The last two terms offset according to (4.6a), so we deduce that q i,j dH(ω i,j ) dt + V sin θ j ∆y [q i,j H(ω i,j ) -q i-1,j H(ω i-1,j )] ≤ D θ ∆θ 2 [q i,j-1 H(ω i,j-1 ) + q i,j+1 H(ω i,j+1 ) -2q i,j H(ω i,j )],

(A.9) holds for i ∈ V y , j ∈ V θ+ , where equality holds if and only if ω i,j = ω i,j-1 = ω i,j+1 = ω i-1,j . Similarly, by rearranging (4.2a) and employ (4.6b), we have q i,j dH(ω i,j ) dt + V sin θ j ∆y [q i+1,j H(ω i+1,j ) -q i,j H(ω i,j )] ≤ D θ ∆θ 2 [q i,j-1 H(ω i,j-1 ) + q i,j+1 H(ω i,j+1 ) -2q i,j H(ω i,j )],

(A.10) holds for i ∈ V y , j ∈ V θ-, where equality holds if and only if ω i,j = ω i,j-1 = ω i,j+1 = ω i+1,j . For i ∈ V y \ {1, 2I}, j ∈ V θ0 , (4.2c) and (4.6c) gives q i,j dH(ω i,j ) dt ≤ D θ ∆θ 2 [q i,j-1 H(ω i,j-1 ) + q i,j+1 H(ω i,j+1 ) -2q i,j H(ω i,j )], (A.11) where equality holds if and only if ω i,j = ω i,j-1 = ω i,j+1 . Third step: special case where ∂p ∂θ jumps. On node (i, j) = (2I, J), substitute p i,j by q i,j ω i,j + q i,j in (4.4c) we obtain q 2I,J dω 2I,J dt = D θ ∆θ 2 q 2I,J-1 ω 2I,J-1 + q 2I,J+1 ω 2I,J+1 -2q 2I,J ω 2I,J + 1 ∆y q +,J+1 ω +,J+1 = D θ ∆θ 2 q 2I,J-1 (ω 2I,J-1 -ω 2I,J ) + q 2I,J+1 (ω 2I,J+1 -ω 2I,J ) + 1 ∆y q +,J+1 (ω +,J+1 -ω 2I,J ) + D θ ∆θ 2 ω 2I,J q 2I,J-1 + q 2I,J+1 -2q 2I,J + 1 ∆y q +,J+1 , (A.12) in which the last term is zero by (4.7). Multiply both sides with H (ω 2I,J ), it follows that q 2I,J dH(ω 2I,J ) dt -D θ ∆θ 2 [q 2I,J+1 H(ω 2I,J+1 ) + q 2I,J-1 H(ω 2I,J-1 ) -2q 2I,J H(ω 2I,J )] = D θ ∆θ 2 q 2I,J-1 [H (ω 2I,J )(ω 2I,J-1 -ω 2I,J ) -H(ω 2I,J-1 )] + D θ ∆θ 2 q 2I,J+1 [H (ω 2I,J )(ω 2I,J+1 -ω 2I,J ) -H(ω 2I,J+1 )] + D θ ∆θ 2 ∆y q +,J+1 H (ω 2I,J )(ω +,J+1 -ω 2I,J+1 ) + 2D θ ∆θ 2 q 2I,J H(ω 2I,J ) = D θ ∆θ 2 q 2I,J-1 [H (ω 2I,J )(ω 2I,J-1 -ω 2I,J ) -H(ω 2I,J-1 ) + H(ω 2I,J )] + D θ ∆θ 2 q 2I,J+1 [H (ω 2I,J )(ω 2I,J+1 -ω 2I,J ) -H(ω 2I,J+1 ) + H(ω 2I,J )] + D θ ∆θ 2 ∆y q +,J+1 [H (ω 2I,J )(ω +,J+1 -ω 2I,J ) -H(ω +,J+1 ) + H(ω 2I,J )] + D θ ∆θ 2 ∆y q +,J+1 H(ω +,J+1 ) -D θ ∆θ 2 q 2I,J-1 + q 2I,J+1 -2q 2I,J + 1 ∆y q +,J+1 H(ω 2I,J ).

The last term is zero according to (4.7), so we deduce that q 2I,J dH(ω 2I,J ) dt ≤ D θ ∆θ 2 [q 2I,J+1 H(ω 2I,J+1 ) + q 2I,J-1 H(ω 2I,J-1 ) -2q 2I,J H(ω 2I,J )] + D θ ∆θ 2 ∆y q +,J+1 H(ω +,J+1 ) (A. [START_REF] Figueroa-Morales | Living on the edge: transfer and traffic of E. coli in a confined flow[END_REF] where equality holds if and only if ω 2I,J-1 = ω 2I,J = ω 2I,J+1 = ω +,J+1 . Similarly, on node (i, j) = (2I, 2J), (1, 0), (1, J), we have q 2I,2J dH(ω 2I,2J ) dt ≤ D θ ∆θ 2 [q 2I,2J+1 H(ω 2I,2J+1 ) + q 2I,2J-1 H(ω 2I,2J-1 ) -2q 2I,2J H(ω 2I,2J )] + D θ ∆θ 2 ∆y q +,2J-1 H(ω +,2J-1 ) (A.14) q 1,0 dH(ω 1,0 ) dt ≤ D θ ∆θ 2 [q 1,1 H(ω 1,1 ) + q 1,-1 H(ω 1,-1 ) -2q 1,0 H(ω 1,0 )] + D θ ∆θ 2 ∆y q -,1 H(ω -,1 ) (A.15) q 1,J dH(ω 1,J ) dt ≤ D θ ∆θ 2 [q 1,J+1 H(ω 1,J+1 ) + q 1,J-1 H(ω 1,J-1 ) -2q 1,J H(ω 1,J )] + D θ ∆θ 2 ∆y q -,J-1 H(ω -,J-1 ) (A.16)

Fourth step: collection and offsets To summarize, sum up (A.9) for i ∈ V y , j ∈ V θ+ , (A.10) for i ∈ V y , j ∈ V θ-, (A.11) for i ∈ V y \ {1, 2I}, j ∈ V θ0 , together with (A.13)-(A.16), the diffusion term offsets and it follows that d dt i∈Vy j∈V θ q i,j H(ω i,j ) + j∈V θ - V sin θ j ∆y q 2I,j H(ω 2I,j ) -j∈V θ - V sin θ j ∆y q 2I,j H(ω 2I,j ) ≤ D θ ∆θ 2 ∆y [q +,J+1 H(ω +,J+1 ) + q +,2J-1 H(ω +,2J-1 ) + q -,1 H(ω -,1 ) + q -,J-1 H(ω +,J-1 )].

(A.17) Hence (A.5), (A.7) and (A.17) imply inequality (A.2) in the theorem.

Define M (t) = ∆θ∆y i∈Vy j∈V θ q i,j ω 2 i,j (t) + ∆θ j∈V θ + q +,j ω 2 +,j (t) + ∆θ j∈V θ - q -,j ω 2 -,j (t). (A.18)

V

  sin θ(qφ) L -L dθ dy dt = 0. Then (3.12) indicates that ω ∞ = ω ∞,+ = C. Similarly, we have ω ∞ = ω ∞,-. Thus there exist a constant C such that

Fifth

  step: strong convergence in time. For a smooth test function φ(y, θ) ∈ C 2 0 (Ω), we define u(t) = Ω qφω dy dθ.Multiplying both sides of (2.13) by φ(y, θ) and integrating over Ω gives Ω

Figure 3 :

 3 Figure 3: Density in the bulk p i,j at different time calculated by Fokker Plank model and SDE model.

Figure 4 :

 4 Figure 4: Density in the bulk near boundary p 1,j , p 2I,j and boundary density p -,j and p +,j (from top to bottom) at different times calculated by the Fokker Plank model (line) and SDE model (histogram). Notice that p 1,j is close to {y = -L} but not at the boundary, so the boundary condition (1.3) is not satisfied for p 1,j .

  , for |Y t | < L or |Y t | = L, and Y t Θ t ≤ 0, 0, for |Y t | = L, and Y t Θ t > 0, dΘ t = 2D θ dB t .

	(0.1)
	Here B

t is a Brownian motion, and D θ is the rotational diffusion coefficient. The process (Y t , Θ t ) stays in the set Ω ∪ Ω + ∪ Ω -.

Table 1 :

 1 (

	∆y	∆θ π/128	π/64	π/32	π/16
	L/128	5.2E-4 7.3E-4 9.4E-4 1.5E-3
	L/64		1.1E-3 1.6E-3 2.2E-3 2.5E-3
	L/32		1.9E-3 2.7E-3 3.8E-3 4.9E-3
	L/16		2.9E-3 4.1E-3 5.7E-3 7.7E-3

Choose H(x) = x 2 in (A.2) then we have dM (t)/ dt ≤ 0. Suppose the initial data of ω i,j and ω ±,j satisfies M (0) < ∞, then since M (t) decreases monotonically, and it has a lower bound 0, so M (t) converges as t goes to infinity. While the equality of (A.2) holds if and only if M (t) = 0, so we have lim t→∞ M (t) = 0, which implies lim t→∞ ω i,j (t) = lim t→∞ ω ±,j (t) = 0. We conclude that p i,j (t) converges to q i,j for any i ∈ V θ , j ∈ V y as t goes to infinity, and p ±,j (t) converges to q ±,j for any j ∈ V j,± as t goes to infinity.

B.P. has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 740623). . M. Tang is supported by NSFC 11871340 and Young Changjiang Scholar of Ministry of Education in China.

A Discrete relative entropy inequality and a priori estimate

The semi-discrete relative gap ω i,j (t) and ω ±,j (t) are defined by

Here we state the discrete relative entropy inequality which is satisfied by discrete density p i,j (t) and p ±,j (t).

Theorem A.1. For any convex function H ∈ C 2 (R), semi-discrete relative gap ω i,j (t) and ω ±,j (t) satisfy the relative entropy inequality:

where equality holds if and only if ω +,j-1 = ω +,j = ω +,j+1 = ω 2I,j .

Sum up (A.5) for j ∈ V θ+ , we have d dt

where equality holds if and only if ω +,j = ω 2I,j = 0 for any j ∈ V θ+ . Similarly, we obtain an inequality for q -,j H(ω -,j ) d dt

where equality holds if and only if ω -,j = ω 1,j = 0 for any j ∈ V θ-.

Second step: estimate of relative entropy at FSP. Substituting p i,j by q i,j ω i,j + q i,j in (4.2a) we obtain that for any i ∈ V y , j ∈ V θ+ , q i,j dω i,j dt = -V sin θ j ∆y (q i,j ω i,j -q i-1,j ω i-1,j ) + D θ ∆θ 2 (q i,j-1 ω i,j-1 + q i,j+1 ω i,j+1 -2q i,j ω i,j ) = -V sin θ j ∆y q i-1,j (ω i,j -ω i-1,j ) + D θ ∆θ 2 q i,j+1 (ω i,j+1 -ω i,j ) -D θ ∆θ 2 q i,j-1 (ω i,j -ω i,j-1 )

+ ω i,j -V sin θ j ∆y (q i,j -q i-1,j ) + D θ ∆θ 2 (q i,j-1 + q i,j+1 -2q i,j ) (A.8) By (4.6a), the last term is zero. Rearranging H (ω i,j )(A.8), we deduce