A kinetic approach to active rod dynamics in confined domains

TL;DR

This paper develops a kinetic model for active self-propelled rods in confined fluids, deriving boundary conditions and simplifying the system by neglecting inertia and diffusion, to better understand wall effects like accumulation and upstream motion.

Contribution

It provides a rigorous derivation of boundary conditions for active rods' distribution in the low inertia limit, reducing model complexity for active matter in confined domains.

Findings

Derived boundary conditions for active rods at walls.

Reduced kinetic model by neglecting inertia and diffusion.

Enabled separate analysis of wall-accumulated and bulk particles.

Abstract

The study of active matter consisting of many self-propelled (active) swimmers in an imposed flow is important for many applications. Self-propelled swimmers may represent both living and artificial ones such as bacteria and chemically driven bi-metallic nano-particles. In this work we focus on a kinetic description of active matter represented by self-propelled rods swimming in a viscous fluid confined by a wall. It is well-known that walls may significantly affect the trajectories of active rods in contrast to unbounded or periodic containers. Among such effects are accumulation at walls and upstream motion (also known as negative rheotaxis). Our first main result is the rigorous derivation of boundary conditions for the active rods' probability distribution function in the limit of vanishing inertia. Finding such a limit is important due to (i) the fact that in many examples of active matter inertia is negligible, since swimming occurs in a low Reynolds number regime, and (ii) this limit allows us to reduce the dimension - and so computational complexity - of the kinetic description. For the resulting model, we derive the system in the limit of vanishing translational diffusion which is also typically negligible for active particles. This system allows for tracking separately active particles accumulated at walls and active particles swimming in the bulk of the fluid.

Equations

$${\partial }_t{f}_{\epsilon }+\frac{1}{\epsilon }\mathbf{v}\cdot {\nabla }_{\mathbf{r}}{f}_{\epsilon }+\frac{1}{{\epsilon }^2}{\nabla }_{\mathbf{v}}\cdot \left((\epsilon \mathbf{u}-\mathbf{v})f_{\epsilon }-D_{\mathrm{tr}}{\nabla }_{\mathbf{v}}{f}_{\epsilon }\right)+$$

$$+\frac{1}{\epsilon }\omega {\partial }_{\phi }{f}_{\epsilon }+\frac{1}{{\epsilon }^2}{\partial }_{\omega }\left((\epsilon T-\omega )f_{\epsilon }-D_{\mathrm{rot}}{\partial }_{\omega }{f}_{\epsilon }\right)=0.$$

$$\mathbf{v}{f}_{\epsilon }\cdot \mathbf{n}={\mathbf{v}}^{\prime }{f}_{\epsilon }^{\prime }\cdot (-\mathbf{n}),\mathbf{r}\text{ on }\Gamma .$$

$$\left[\begin{array}[]{c}\boldsymbol{v}^{\prime}\\ \omega^{\prime}\end{array}\right]=\mathcal{C}\left[\begin{array}[]{c}\boldsymbol{v}\\ \omega\end{array}\right],\quad\mathcal{C}\in\mathbb{R}^{3\times 3},\quad|\det{\mathcal{C}}|=1.$$

$$\mathcal{D}:=\left[\begin{array}[]{ccc}D_{\mathrm{tr}}&0&0\\ 0&D_{\mathrm{tr}}&0\\ 0&0&D_{\mathrm{rot}}\end{array}\right].$$

$${\partial }_t{f}_{\epsilon }+\frac{1}{\epsilon }\mathcal{V}\cdot {\nabla }_{\mathcal{X}}{f}_{\epsilon }+\frac{1}{{\epsilon }^2}{\nabla }_{\mathcal{V}}\cdot \left((\epsilon \mathcal{U}-\mathcal{V})f_{\epsilon }-\mathcal{D}{\nabla }_{\mathcal{V}}{f}_{\epsilon }\right)=0.$$

$${\rho }_{\epsilon }(t,\mathcal{X}):={\int }_{{\mathbb{R}}^3}{f}_{\epsilon }(t,\mathcal{X},\mathcal{V})\mathrm{d}\mathcal{V},$$

$${\partial }_t\rho +{\nabla }_{\mathcal{X}}\cdot (\mathcal{U}\rho )={\nabla }_{\mathcal{X}}\cdot \mathcal{D}{\nabla }_{\mathcal{X}}\rho$$

$$D_{\mathrm{tr}}\frac{\partial \rho }{\partial \mathbf{n}}=(\mathbf{u}\cdot \mathbf{n})\rho ,\mathbf{r}\text{ on }\Gamma ,-\pi \le \phi <\pi .$$

$$\delta :=D_{\mathrm{tr}}.$$

$${\partial }_t\rho +{\nabla }_{\mathbf{r}}\cdot (\mathbf{u}\rho )+{\partial }_{\phi }(T\rho )=\delta {\Delta }_{\mathbf{r}}\rho +D_{\text{rot}}{\partial }_{\phi }^2\rho$$

$$\delta \frac{\partial \rho }{\partial \mathbf{n}}=(\mathbf{u}\cdot \mathbf{n})\rho ,\mathbf{r}\text{ on }\Gamma ,-\pi \le \phi <\pi .$$

$$\rho (t,\mathbf{r},\phi )={\psi }_{\mathrm{wall}}(t,\mathbf{r},\phi ){\delta }_{\Gamma }(\mathbf{r})+{\rho }_{\mathrm{bulk}}(t,\mathbf{r},\phi ).$$

$${\partial }_t{\rho }_{\mathrm{bulk}}+{\nabla }_{\mathbf{r}}\cdot (\mathbf{u}{\rho }_{\mathrm{bulk}})+{\partial }_{\phi }(T{\rho }_{\mathrm{bulk}})=D_{\mathrm{rot}}{\partial }_{\phi }^2{\rho }_{\mathrm{bulk}}+\sum _{i=1,2}{\chi }_i{\delta }_{\Gamma }(\mathbf{r})\delta (\phi -{\phi }_i),$$

$$\mathbf{r}\in \Omega ,-\pi \le \phi <\pi ,$$

$${\partial }_t{\psi }_{\mathrm{wall}}+{\partial }_s((\mathbf{u}\cdot \mathbf{\tau }){\psi }_{\mathrm{wall}})+{\partial }_{\phi }\left(T{\psi }_{\mathrm{wall}}\right)=D_{\mathrm{rot}}{\partial }_{\phi }^2{\psi }_{\mathrm{wall}}+(\mathbf{u}\cdot \mathbf{n}){\rho }_{\mathrm{bulk}},$$

$$\mathbf{r}\in \Gamma ,s\in [0,L],\phi \in ({\phi }_1,{\phi }_2),$$

$${\rho }_{\mathrm{bulk}}=0,\mathbf{r}\in \Gamma ,\phi \notin [{\phi }_1,{\phi }_2],{\rho }_{\mathrm{bulk}}\mathrm{is}2\pi -\mathrm{periodic}\mathrm{in}\phi \mathrm{for}\mathrm{all}\mathbf{r}\in \Omega ,$$

$${\psi }_{\mathrm{wall}}=0,s\in [0,L],\phi \in \left\{{\phi }_1,{\phi }_2\right\},$$

$${\chi }_i(s,\phi ):=(-1{)}^i\left(T{\psi }_{\mathrm{wall}}-D_{\mathrm{rot}}{\partial }_{\phi }{\psi }_{\mathrm{wall}}\right),i=1,2.$$

$$\mathbf{u}\cdot \mathbf{n}>0$$

$$J_{\epsilon }(t,\mathcal{X}):=\frac{1}{\epsilon }{\int }_{{\mathbb{R}}^3}\mathcal{V}{f}_{\epsilon }\mathrm{d}\mathcal{V}\text{and}{\mathbb{P}}_{\epsilon }(t,\mathcal{X}):={\int }_{{\mathbb{R}}^3}\mathcal{V}\otimes \mathcal{V}{f}_{\epsilon }\text{d}\mathcal{V}.$$

$${\partial }_t{\rho }_{\epsilon }+{\nabla }_{\mathcal{X}}\cdot J_{\epsilon }=0,$$

$${\epsilon }^2$$

$${\partial }_t\rho +{\nabla }_{\mathcal{X}}\cdot (\mathcal{U}\rho )={\Delta }_{\mathcal{X}}\rho .$$

$${\hat{J}}_{\epsilon }\cdot \mathbf{n}=0\text{ on }\Gamma ,$$

$$\frac{\partial \rho }{\partial \mathbf{n}}=(\mathbf{u}\cdot \mathbf{n})\rho .$$

$${\mathcal{E}}_{\epsilon }(t)$$

$$d_{\epsilon }(t,\mathcal{X},\mathcal{V})$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\mathcal{E}}_{\epsilon }(t)+\frac{1}{2{\epsilon }^2}{\int }_{{\mathbb{R}}^3}{\int }_{\Omega \times (-\pi ,\pi )}|d_{\epsilon }{|}^2\mathrm{d}\mathcal{X}\mathrm{d}\mathcal{V}<C.$$

$$\frac{{\mathbf{v}}^2}{2}+\frac{{\omega }^2}{2}+\ln f_{\epsilon }$$

$$\frac{\text{d}}{\text{d}t}\left[{\int }_{{\mathbb{R}}^3}{\int }_{\Omega \times (-\pi ,\pi )}\left\{\frac{{\mathbf{v}}^2}{2}+\frac{{\omega }^2}{2}+\ln f_{\epsilon }\right\}{f}_{\epsilon }\text{d}\mathcal{X}\text{d}\mathcal{V}\right]=$$

$$=-\frac{1}{\epsilon }{\int }_{{\mathbb{R}}^3}{\int }_{-\pi }^{\pi }{\int }_{\Gamma }\left\{\frac{{\mathbf{v}}^2}{2}+\frac{{\omega }^2}{2}+\ln f_{\epsilon }\right\}(\mathbf{v}\cdot \mathbf{n})f_{\epsilon }\text{d}{s}_{\mathbf{r}}\text{d}\phi \text{d}\mathcal{V}$$

$$-\frac{1}{{\epsilon }^2}{\int }_{{\mathbb{R}}^3}{\int }_{\Omega \times (-\pi ,\pi )}{\left|(\mathcal{V}-\epsilon \mathcal{U})+{\nabla }_{\mathcal{V}}(\ln f_{\epsilon })\right|}^2{f}_{\epsilon }\text{d}\mathcal{X}\text{d}\mathcal{V}$$

$$-\frac{1}{\epsilon }{\int }_{{\mathbb{R}}^3}{\int }_{\Omega \times (-\pi ,\pi )}\mathcal{U}\cdot \left(\mathcal{V}-\epsilon \mathcal{U}\right)f_{\epsilon }\text{d}\mathcal{X}\text{d}\mathcal{V}.$$

$$S_{\sigma >0}:=\left\{(\mathbf{v},\omega ):\sigma >0\right\}\text{ and }{S}_{\sigma <0}:=\left\{(\mathbf{v},\omega ):\sigma <0\right\},$$

Full text

1. 1 Introduction

2. 2 Main results

3. 3 Vanishing inertia limit in Fokker-Planck equation

4. 4 Boundary layer equation at wall in vanishing translational diffusion limit: derivation of (10)

5. 5 Numerical example

6. 6 Individual based model and Fokker-Planck equation for an active rod with inertia

7. 6.1 Individual based model for an active rod with inertia

8. 6.2 Fokker-Planck equation for an active rod with inertia

9. Acknowledgment

1 Introduction

Recently, active matter has attracted much attention of the scientific community (see e.g. reviews [1, 2, 3]). In general, active matter is defined as a system of many agents moving due to consumption of energy stored in the surrounding environment (e.g., chemical or food) and converting it into mechanical force which is called self-propulsion. The agents exhibiting self-propulsion are named active, as opposed to passive agents which can move only if an external field is applied. There are a vast number of examples of active matter: from suspensions of bacteria [4, 5, 6, 7, 8, 9, 10], flocks of birds [11, 12, 13], and schools of fish [14, 15] to crowds of people [16] which also meet the definition of active matter since people exhibit self-propulsion (walking). Modeling and further analysis of active matter is of great importance due to the variety of striking phenomena and promising applications (reduction of viscosity, cargo delivery for medical purposes, materials repair, etc.).

In this work, we are interested in modeling the wide class of active matter where agents are rod-shaped microswimmers, i.e., the surrounding environment is a viscous fluid and swimming occurs in the low Reynolds number regime. Examples of such microswimmers are bacteria (especially, rod-shaped B. subtilis) and active bi-metallic micro- and nanorods which swim in a viscous solution with hydrogen peroxide [17]. It was observed both theoretically and experimentally for various types of microswimmers that their trajectories are much more complex than in the case of passive swimmers which simply follow streamlines of an external field (the background flow). In particular, the following phenomena were observed in dynamics of active microswimmers: accumulation at walls (bordertaxis) and upstream motion (rheotaxis), see [18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 9, 28, 29] and references therein. Throughout the paper we use the term “active rods” for these rod-shaped active microswimmers.

There are two common mathematical approaches to describe dynamics of an active rod in a viscous fluid. The first one is based on force and torque balances for each individual active rod. This approach results in a Langevin equation (or a system of coupled Langevin equations) for unknown location, orientation and velocities, both translational and angular, of the active rod. In the second approach, which is also called a kinetic approach, the main unknown is the probability distribution function of the active rod, and the function satisfies the Fokker-Planck equation. These two approaches are directly related mathematically: roughly speaking, right hand sides of equations in the first approach are coefficients of the Fokker-Planck equation in the second approach. The second approach is more preferable if one studies statistical properties of a large number of active rods since it does not require many realizations, unlike the first direct one.

The focus of this work is on the development of a kinetic approach for active rods swimming in a container restricted by a confinement (a wall). One can formulate how an active rod behaves when it collides with the wall (a collision rule) in the first approach. On the other hand, it is not immediately clear how the collision rule translates into a boundary condition for the Fokker-Planck equation. This is because a collision rule is typically a relation between velocities before and after a collision, whereas the Fokker-Planck equation is usually written in the vanishing inertia (overdamped) limit. Thus, the active rods’ velocities are no longer variables of the unknown probability distribution function. The vanishing inertia limit of the Fokker-Planck equation is relevant for the low Reynolds number regime and important since it allows one to reduce dimension and thus drastically decrease computational complexity and even make the equation amenable for analysis.

Our first result is the rigorous derivation of this limit and, more importantly, the boundary condition for the overdamped Fokker-Planck equation in this limit. Namely, we show that the limiting probability distribution function, which depends on the active rod’s location and orientation, satisfies the no-flux condition on the wall for each given orientation. We note that similar boundary conditions were phenomenologically and independently derived in [30] to analyze the distribution of active rods inside an infinite channel.

Next, we consider the case of a small translational diffusion which is negligible in experiments for active particles and equated to zero in corresponding individual based models. By using the boundary layer multi-scale approach, we derive the kinetic system in the limit of vanishing translational diffusion. The significance of the system is that it describes explicitly the population of active rods accumulated at walls.

The structure of this paper is as follows. First, in Section 2 we formulate our two main results: on vanishing inertia and vanishing translational diffusion limits. Details of the results’ derivation are relegated to Sections 3 and 4. Next, in Section 5 we present a numerical example in which we compare the derived limiting kinetic models with Monte Carlo simulations for the corresponding individual based model. Finally, in Section 6 we provide a specific physical model of a self-propelled nano-particle swimming in a viscous flow by presenting both individual based model and the corresponding (pre-limiting) Fokker-Planck equation.

2 Main results

We start with the Fokker-Planck equation describing random dynamics of an active rod with inertia:

[TABLE]

The unknown function $f_{\varepsilon}(t,\boldsymbol{r},\boldsymbol{v},\varphi,\omega)$ is the probability distribution function of the active rod’s location ${\boldsymbol{r}}\in\Omega\subset\mathbb{R}^{2}$, translational velocity $\boldsymbol{v}\in\mathbb{R}^{2}$, orientation angle $\varphi\in[-\pi,\pi)$, and angular velocity $\omega\in\mathbb{R}$. Given functions $\boldsymbol{u}=\boldsymbol{u}(\boldsymbol{r},\varphi)$ and $T=T(\boldsymbol{r},\varphi)$ are smooth in $\Omega\times[-\pi,\pi)$ and $2\pi$-periodic in $\varphi$. Small positive parameter $\varepsilon\ll 1$ measures the effect of inertia on dynamics of an active rod. The unknown function $f_{\varepsilon}$ is $2\pi$-periodic in $\varphi$. In addition, the following boundary condition is imposed on $f_{\varepsilon}$:

[TABLE]

Here $\Gamma$ is the boundary of $\Omega$ (a wall of the container), $\boldsymbol{n}$ is the outward normal, and $f_{\varepsilon}^{\prime}=f_{\varepsilon}(t,\boldsymbol{r},\boldsymbol{v}^{\prime},\varphi,\omega^{\prime})$ where pairs $(\boldsymbol{v},\omega)$ and $(\boldsymbol{v}^{\prime},\omega^{\prime})$ represent translational and angular velocities of the active rod before and after a collision with the wall. The relation between the two pairs of velocities is given by:

[TABLE]

The specific form of matrix $\mathcal{C}$ as well as the derivation of the kinetic model (2)-(2) from the individual dynamics of an active rod are relegated to Section 6.1.

To simplify notations, denote $\mathcal{X}:=({\boldsymbol{r}},\varphi)$, ${\mathcal{V}}:=(\boldsymbol{v},\omega)$, $\mathcal{U}:=(\boldsymbol{u},T)$, and

[TABLE]

In new notations, the Fokker-Planck equation (2) is

[TABLE]

Our first main result is the reduction, for small $\varepsilon$, of the unknown function $f_{\varepsilon}(t,\mathcal{X},\mathcal{V})$ depending on a 7-dimensional variable to the unknown function $\rho(t,\mathcal{X})$ depending on a 4-dimensional variable. The result is formulated in the following theorem.

Theorem 2.1**.**

Let $\rho_{\varepsilon}$ be defined by

[TABLE]

where $f_{\varepsilon}$ solves Fokker-Planck equation (4) with boundary condition (2). Then $\rho_{\varepsilon}$ converges to $\rho$ in the distributional sense, as $\varepsilon\to 0$, where $\rho$ satisfies the following limiting equation

[TABLE]

with the boundary condition

[TABLE]

Recall that $\mathcal{U}=(\boldsymbol{u},T)$.

The proof of Theorem 2.1 is given in Section 3. This theorem means that elastic collision boundary condition (2) for $f_{\varepsilon}$ transforms into no-flux boundary condition (6) for $\rho$. The vanishing inertia limit in the Fokker-Planck equation for spherical particles (no $\varphi$) was considered in [31, 32, 33]. The main difference, besides no $\varphi$, is that in [31, 32, 33], $\boldsymbol{u}$ is not given but instead solves the Navier-Stokes equation. In principle, a system with Navier-Stokes equation is obviously more complicated; on the other hand, Navier-Stokes equation has an additional dissipation term in the energy relation. Also, the coupling term leads to a certain cancellation in the energy relation (the term $\iint\mathcal{U}\cdot(\mathcal{V}-\varepsilon\mathcal{U})f_{\varepsilon}$ in (15) from Section 3). The scaling (that is, how $\varepsilon$ is introduced in Fokker-Planck equation) in (4) is similar to [31], but the boundary conditions in [31] are simpler (periodic). In [33] the asymptotic regime of the Fokker-Planck equation is studied for the reflection boundary condition (spherical particles elastically collide with walls) as well, but due to no slip conditions for $\boldsymbol{u}$ (which imply no flow at boundary $\Gamma$), the limiting boundary conditions were not investigated.

In experimental observations, rod-like microswimmers, such as bacteria or bi-metallic particles, are more likely to spontaneously turn rather than jump to another position. In other words, random forces along the perimeter of a microswimmer caused by collisions with molecules of the fluid likely result in a significant torque whereas a net force is small. These observations imply that the translational diffusion coefficient is small, $D_{\text{tr}}\ll 1$. Hence, it is natural to study the limit $D_{\text{tr}}\to 0$ which is singular in the case of active particles, that is one cannot simply set $D_{\text{tr}}$ to zero in both Fokker-Planck equation (5) and no-flux condition (6) in order to obtain the limiting system as $D_{\text{tr}}\to 0$. This is because active particles tend to accumulate at walls and in particular they form a boundary layer of the width $\sim D_{\text{tr}}$. Moreover, from boundary conditions (6) it follows that $\nabla\rho$ may blow up near boundary $\Gamma$ in the limit $D_{\text{tr}}\to 0$; in this case for stable numerical simulation of (5)-(6) with small $D_{\text{tr}}$ one needs a very fine mesh resulting in high computational complexity of the simulations.

To investigate vanishing translational diffusion $D_{\text{tr}}$ we re-denote the translational diffusion coefficient by symbol $\delta$ which is typically used for notations of small parameters:

[TABLE]

Then problem (5)-(6) consists of the Fokker-Planck equation

[TABLE]

and the boundary condition

[TABLE]

Our second main result is the derivation of the limit in problem (7)-(8) as $\delta\to 0$. The derivation is done by formal multi-scale asymptotic expansion; we formulate the result as a conjecture since rigorous justification of the multi-scale asymptotic expansions are out of the scope of this work.

Conjecture: In the limit $\delta\to 0$ the probability distribution function $\rho(t,\boldsymbol{r},\varphi)$ has the following representation:

[TABLE]

Here $\delta_{\Gamma}(\boldsymbol{r})$ is the $\delta$-function distribution supported on $\Gamma=\left\{\boldsymbol{\gamma}(s):0\leq s\leq L\right\}$ ($s$ is arc-length parameter of curve $\Gamma$) and probability distribution functions $\rho_{\text{bulk}}$ and $\psi_{\text{wall}}$ solve the following system:

[TABLE]

Here $\boldsymbol{\tau}$ denotes the tangential vector of $\Gamma$ and angles $\varphi_{1}(s)$ and $\varphi_{2}(s)$ are introduced such that

[TABLE]

Representation (9) means that the total probability distribution function $\rho$ consists of the regular part, $\rho_{\text{bulk}}$, describing distribution of particles in the bulk, and the singular part, $\psi_{\text{wall}}\,\delta_{\Gamma}$, describing distribution of particles accumulated at wall. Derivation of system (10) is presented in Section 4. We also provide a numerical example in Section 5 in which we test the derived kinetic approaches (7)-(8) and (10) with results of Monte-Carlo simulations for the corresponding individual based model.

3 Vanishing inertia limit in Fokker-Planck equation

Proof of Theorem 2.1. In this section we take $D_{\text{tr}}=D_{\text{rot}}=1$ for the sake of simplicity. To consider the limit $\varepsilon\to 0$, introduce the mean flux (or, the mean velocity), and “the kinetic pressure”:

[TABLE]

By integration of (4) with respect to $\text{d}\mathcal{V}$ and $\varepsilon\mathcal{V}\,\text{d}\mathcal{V}$ one obtains the system for $\rho_{\varepsilon}$ and $J_{\varepsilon}$:

[TABLE]

By using arguments similar to [31], we will show that the limit of (12) is $J=\rho\,\mathcal{U}-\nabla_{\mathcal{X}}\rho$. Substitution of this formula for $J$ into the limiting version of the equation (11) (that is, equation (11) without sub-indexes $\varepsilon$), then gives:

[TABLE]

The main question is how to find the boundary condition for $\rho$. Note that from collision boundary condition (2) it follows that active rods cannot leave domain $\Omega$, that is, there is no flux through the boundary $\Gamma$:

[TABLE]

where $\hat{J}_{\varepsilon}=(J_{1},J_{2})$ (no $J_{3}$, corresponding to the flux of orientations $\varphi$). The main purpose of this section is to prove that this relation is preserved in the limit $\varepsilon\to 0$, which, taking into account the formula for the limiting flux $J=\rho\,\mathcal{U}-\nabla_{\mathcal{X}}\rho$, is equivalent to

[TABLE]

In order to prove Theorem 2.1 we will use two auxiliary propositions. In Proposition 1, the energy estimate is established. This estimate leads to a priori bounds needed to obtain that the family $\left\{f_{\varepsilon}\right\}_{\varepsilon}$ has a limit as $\varepsilon\to 0$ (Proposition 2).

First introduce the following notations:

[TABLE]

Recall that $\mathcal{X}=(\boldsymbol{r},\varphi)\in\Omega\times[-\pi,\pi)$ and $\mathcal{V}=(\boldsymbol{v},\omega)\in\mathbb{R}^{2}_{\boldsymbol{v}}\times\mathbb{R}_{\omega}=\mathbb{R}^{3}$.

Proposition 1**.**

There exists a constant $C$, independent of $\varepsilon$, such that the following estimate (the entropy inequality) holds:

[TABLE]

Proof.

Multiplication of (4) by

[TABLE]

and integration with respect to both $\mathcal{X}$ and ${\mathcal{V}}$ gives:

[TABLE]

Next we compute the boundary term in the right hand side of (15) (the one with the integral over $\Gamma$). To this end, we use the boundary condition (2). For each $\boldsymbol{r}\in\Gamma$ and $-\pi\leq\varphi<\pi$ denote

[TABLE]

where $\sigma=\sigma(\boldsymbol{r},\boldsymbol{v},\varphi,\omega)$ is given by (53) in Section 6.1. The introduced sets $S_{\sigma>0}$ and $S_{\sigma<0}$ can be understood as sets of configurations (velocities) of an active rod before and after a collision, respectively, at the given location of the boundary $\boldsymbol{r}\in\Gamma$ and the given orientation $\varphi\in[-\pi,\pi)$. Then the boundary integral can be written as follows for each $\boldsymbol{r}\in\Gamma$ and $-\pi\leq\varphi<\pi$:

[TABLE]

We claim that the two integrals in the right hand side of the equality above cancel each other. To verify this, one needs to make the substitution in the first integral $\mathcal{V}^{\prime}=\mathcal{C}\mathcal{V}$ with $\mathcal{C}$ from (55) (or equivalently (50)-(51)) and to use boundary condition (2), conservation of energy during a collision (52), and $\mathrm{d}\mathcal{V}^{\prime}=\mathrm{d}\mathcal{V}$ which follows from $|\text{det}\,\mathcal{C}|=1$ (see (3)). Hence, the boundary term in the right hand side of (15) vanishes.

Note that in the same manner one can show that the conservation of total $f_{\varepsilon}$:

[TABLE]

Indeed, by using equation (4) and integration by parts it follows that

[TABLE]

and one can show that the right hand side vanishes following the same arguments as for the boundary term in (15)

Finally, the last term in the right hand side of (15) is estimated as follows:

[TABLE]

Thus, we obtained (14) and the proposition is proved. ∎

In the standard manner (see, e.g., [31]) the entropy estimate (14) implies the following bounds:

[TABLE]

The proof of the following proposition is also standard and can be found in [31].

Proposition 2**.**

There exist such $\rho$ and $J$ that the following convergences hold as $\varepsilon\to 0$ in the distributional sense:

[TABLE]

Here $\mathbb{I}$ is identity matrix.

Now we are in position to complete the proof of Theorem 2.1.

Back to proof of Theorem 2.1.

By multiplication of (4) by a scalar test functions $\psi(t,\mathcal{X},{\mathcal{V}})$ with a finite support in $0<t<T$, $\mathcal{X}\in\Omega\times(-\pi,\pi)$ and ${\mathcal{V}}\in\mathbb{R}^{3}$, integration with respect to $t$, $\mathcal{X}$ and $\mathcal{V}$, as well as integration by parts, one obtains the following equality:

[TABLE]

where $\psi^{\prime}=\psi(t,\mathcal{X},\mathcal{V}^{\prime})$ with $\mathcal{V}^{\prime}=\mathcal{C}\mathcal{V}$, and matrix $\mathcal{C}$ and set $S_{\sigma>0}$ are defined by (55) and (16), respectively. Equality (19) can be understood as the weak formulation of the problem (4) with boundary condition (2).

Next, take a test function in (19) which is independent of $\mathcal{V}$: $\psi:=g(t,\mathcal{X})$ (for the sake of clarity, choose different symbol for the test function here: $g$ instead of $\psi$). Formally, such a test function is not admissible since it does not have a finite support in $\mathcal{V}$. On the other hand, one can use truncations $g(t,\mathcal{X})\chi_{|\mathcal{V}|<R}(\mathcal{V})$ as test functions and pass to the limit $R\to\infty$ to obtain (19) for $\psi=g(t,\mathcal{X})$. By integrating in ${\mathcal{V}}$ we obtain:

[TABLE]

Note that the boundary term in (19) vanishes for test functions independent from $\mathcal{V}$. By passing to the limit $\varepsilon\to 0$ we obtain:

[TABLE]

Take the test function of the form $\psi:=\varepsilon{v}_{i}h_{i}(t,\mathcal{X})$ for $i=1,2$:

[TABLE]

Here both the super- and sub-index $i$ stand for the coordinate number. Consider $\boldsymbol{h}=\left\{h_{i}\right\}_{i=1}^{3}$ so that $\boldsymbol{h}\cdot\boldsymbol{n}=h_{1}n_{1}+h_{2}n_{2}=0$ for all $-\pi\leq\varphi<\pi$ and $h_{3}(t,\mathcal{X})$ is arbitrary. The test function $\psi$ corresponding to $h_{3}(t,\mathcal{X})$ is $\psi=\omega h_{3}(t,\mathcal{X})$, and in this case equality (3) holds for $i=3$ with no second (boundary) term. Then taking the sum with respect to $i$ in (3) leads to that the boundary term in (3) vanishes, so that the following equality holds:

[TABLE]

Passing to the limit $\varepsilon\to 0$ and using Proposition 2 we get:

[TABLE]

This equality gives us the following relation for $J$:

[TABLE]

for all $\boldsymbol{h}$ such that $h_{1}n_{1}+h_{2}n_{2}|_{\Gamma}=0$ for all $-\pi\leq\varphi<\pi$. Finally, take any $g(t,\mathcal{X})$ in (20) such that $\left.\dfrac{\partial g}{\partial\boldsymbol{n}}\right|_{\Gamma}=0$ for all $-\pi\leq\varphi<\pi$ and $\boldsymbol{h}:=\nabla_{\mathcal{X}}g$ in (23) to express $\int_{0}^{T}\int_{\Omega\times(-\pi,\pi)}J\cdot\nabla_{\mathcal{X}}g\,\text{d}\mathcal{X}\text{d}t$ in (20). We obtain the following weak formulation for the equation for $\rho$:

[TABLE]

for all $g$ such that $\left.\dfrac{\partial g}{\partial\boldsymbol{n}}\right|_{\Gamma}=0$ for all $-\pi\leq\varphi<\pi$. By integrating by parts and using condition $\left.\dfrac{\partial g}{\partial\boldsymbol{n}}\right|_{\Gamma}=0$, (24) leads to

[TABLE]

Varying $g$ on $\Gamma$ we get our final result which is the no-flux boundary condition (13):

[TABLE]

for all $-\pi\leq\varphi<\pi$.

Thus, the proof of Theorem 2.1 is complete. ∎

4 Boundary layer equation at wall in vanishing translational diffusion limit: derivation of (10)

To describe the behavior of $\rho$ near the wall, one needs to consider the limiting behavior of $\rho$ inside the “band” $\Omega^{*}:=\left\{\boldsymbol{r}\in\Omega:\text{dist}(\boldsymbol{r},\Gamma)<c\right\}$ where $c$ is small but independent of $\delta$. Here we choose $c>0$ such that $c<\kappa_{\text{max}}^{-1}$ where $\kappa_{\text{max}}$ is the maximum curvature along $\Gamma$. If $\Gamma$ is a straight line or a segment, then curvature $\kappa$ is zero, and $c$ is an arbitrary number independent of $\delta$.

We introduce new coordinate system in band $\Omega^{*}$, related to parametrization of wall $\Gamma=\partial\Omega$. Namely, let $\Gamma=\{\boldsymbol{\gamma}(s)\colon 0\leq s\leq L\}$ be the natural parametrization of the wall $\Gamma$ (in other words, $s$ is the arc length parameter) and $\boldsymbol{n}(s)$, $\boldsymbol{\tau}(s)$ be an outward normal and tangential vectors at $\boldsymbol{\gamma}(s)\in\Gamma$, respectively. For every $\boldsymbol{r}\in\Omega^{*}$, we define

[TABLE]

where $r=\text{dist}(\boldsymbol{r},\Gamma)$ and $\boldsymbol{\gamma}(s)$ is the “projection” of $\boldsymbol{r}$ onto $\Gamma$.

In this new coordinate system we introduce the two-scale ansatz for unknown function $\rho$:

[TABLE]

Here sub-indexes ’$w$’ and ’$b$’ stand for “wall” and “bulk”, respectively. Variable $z$ denotes below the second argument of $\rho_{w}^{(k)}$, the $k$th wall (boundary layer) coefficient, i.e., $z=\delta^{-1}r$ in (27). We assume that functions $\{\rho_{w}^{(k)}(t,z,s,\varphi)\}_{k}$ vanish with all derivatives in $z$ as $z\to\infty$. In what follows, we focus on the first three terms of two-scale expansion (27) or, in other words, on terms of order $\delta^{-1}$ and $\delta^{0}$:

[TABLE]

The representations (28) and (9) are related via the following equalities:

[TABLE]

Next, we rewrite the Fokker-Planck equation (7) in the coordinate system $(r,s)$. To this end, we introduce the inverse substitution functions $R(\boldsymbol{r})$ and $S(\boldsymbol{r})$:

[TABLE]

Using chain rule and the 2D Frenet-Serret relation for the normal vector $\boldsymbol{n}^{\prime}(s)=-\kappa(s)\boldsymbol{\tau}(s)$ where $\kappa(s)$ is the curvature of $\Gamma$ at $\boldsymbol{r}=\boldsymbol{\gamma}(s)$, one obtains

[TABLE]

To compute $\nabla_{\boldsymbol{r}}\rho$, $\Delta_{\boldsymbol{r}}\rho$ and $\nabla_{\boldsymbol{r}}\cdot\boldsymbol{u}$ one can use (30) and both 2D Frenet-Serret relations, $\boldsymbol{n}^{\prime}(s)=-\kappa(s)\boldsymbol{\tau}(s)$ and $\boldsymbol{\tau}^{\prime}(s)=\kappa(s)\boldsymbol{n}(s)$:

[TABLE]

Here $u_{\text{n}}=\boldsymbol{u}\cdot\boldsymbol{n}$ and $u_{\tau}=\boldsymbol{u}\cdot\boldsymbol{\tau}$.

Then the original problem (7)-(8) converts into

[TABLE]

with boundary conditions

[TABLE]

When substituting representation (28) into equation (31), we will treat the second term in the left hand side as follows:

[TABLE]

where $u_{\text{n}}^{(k)}=1/k!\,\partial^{k}_{r}u_{\text{n}}|_{r=0}$. Note that $u_{\text{n}}^{(k)}$ is a function of $\varphi$ and $s$ for each $k=1,2,...$.

At the order $\delta^{-2}$ in (31) one has the following equality:

[TABLE]

The lowest order in the boundary conditions (8) is $\delta^{-1}$ and the corresponding equality is

[TABLE]

Combining (33) and (34) we obtain a formula for $\rho_{w}^{(-1)}$:

[TABLE]

At order $\delta^{-1}$, the Fokker-Planck equation (31) has the form

[TABLE]

Boundary conditions at the order $\delta^{0}$ looks as follows:

[TABLE]

Consider $u_{\text{n}}^{(0)}>0$. After integration (35) with respect to $z$ from [math] to $\infty$ and simplifications one obtains an equation for $\psi_{\text{wall}}=\int_{0}^{\infty}\rho_{w}^{(-1)}\,\text{d}z=B/u_{\text{n}}^{(0)}$:

[TABLE]

Equation (36) is a conservation law for the distribution of active rods $\psi_{\text{wall}}$ accumulated at wall; these active rods re-orient and move along the wall in the tangential direction. Term $u_{\text{n}}^{(0)}\rho_{b}^{(0)}$ in the left hind side of (36) accounts for particles coming from bulk. If divergence-free condition is imposed, $\nabla_{\bf r}\cdot u=0$, which is for $r=0$ has the form $u_{\text{n}}^{(1)}-\partial_{s}u_{\tau}^{(0)}-\kappa u_{\text{n}}^{(0)}=0$, then (36) can be rewritten in the form:

[TABLE]

If $u_{\text{n}}^{(0)}\leq 0$, then $\psi_{\text{wall}}=0$ and $\rho_{b}^{(0)}|_{r=0}=0$.

Next we obtain equation for $\rho_{\text{bulk}}=\rho^{(0)}_{b}$. To this end, consider equation (31) at the order $\delta^{0}$ and pass to the limit $z\to\infty$. After we rewrite the resulting equation in the original coordinate system we recover the Fokker-Planck equation for $\rho_{\text{bulk}}$

[TABLE]

Representation (9) with $\psi_{\text{wall}}$ and $\rho_{\text{bulk}}$ satisfying equations (36) and (37) respectively, is valid if $\boldsymbol{u}\cdot\boldsymbol{n}|_{\Gamma}>0$ for all $-\pi\leq\varphi<\pi$. In this case, active rods accumulate at the wall, but they cannot leave the wall. On the other hand, due to the rotational diffusion, active rods in experiments (for example those described by system (43)-(44)) may leave the wall and be re-injected into the bulk. In this case, one needs an additional boundary layer term in the representation (27) (with an additional scale different from $1$ and $\delta$) which will capture active rods at wall $\Gamma$ with $\boldsymbol{u}\cdot\boldsymbol{n}\approx 0$. It is similar to “parabolic boundary layers” in vanishing diffusion limit in elliptic equations; these boundary layers introduce two new scales $\delta^{1/3}$ and $\delta^{2/3}$, and these terms are constructed at boundary points where characteristics of the limiting hyperbolic equation are tangential to the boundary, see Sec. 2.7.5 in [34], see also [35]. In this work, the formula for $\chi$ in (10), the flux of re-injected active rods, i.e.,, with $\boldsymbol{r}\in\Gamma$, $\boldsymbol{u}\cdot\boldsymbol{n}|_{\Gamma}\approx 0$ and $\dfrac{\text{d}}{\text{d}t}\boldsymbol{u}\cdot\boldsymbol{n}|_{\Gamma}<0$, is derived from the conservation of total density:

[TABLE]

Adding terms with flux $\chi$ to the right hand side of equation (37) we derive system (10).

5 Numerical example

In this section, we provide a numerical example to illustrate the relation between a specific individual based model for an active rod and corresponding kinetic approaches discussed above.

We assume that the wall $\Gamma$ coincides with $x$-axis, i.e. $\Gamma=\left\{(x,y):y=0\right\}$, and the active rod’s probability distribution function does not depend on $x$ (see Fig. 1, left). Drag $\boldsymbol{u}$ exerted on an active rod is the sum of two components: a background shear flow $\boldsymbol{u}_{\text{BG}}=(\dot{\gamma}y,0)$ and self-propulsion $v_{\text{prop}}(\cos(\varphi),\sin\varphi)$. The vertical component $u$ of drag velocity $\boldsymbol{u}$ and torque $T$ exerted on an active rod, are defined as follows:

[TABLE]

Note that since the background shear flow $\boldsymbol{u}_{\text{BG}}$ has zero $y$-component it does not enter the formula for vertical drag $u$. Expression for $T$ follows from (45) from Section 6.1 with $T=\Phi_{\text{BG}}$. Parameters $v_{\text{prop}}=0.2$ and $\dot{\gamma}=1.0$ are self-propulsion speed and shear rate, respectively.

First, we perform Monte Carlo simulations for the individual based model for an active rod with vertical component $y(t)$ of location $\boldsymbol{r}(t)$ and orientation angle $\varphi$ swimming in $\Omega~{}=~{}\left\{y>0\right\}$:

[TABLE]

where $\zeta$ is the white noise with $\langle\zeta(t),\zeta(t^{\prime})\rangle=\delta(t-t^{\prime})$ and $D_{\text{rot}}$ is the rotational diffusion coefficient. The following “overdamped” collision-with-wall rule for $y(t)=0$ is imposed:

[TABLE]

This collision rule means that the active rod does not move if it is oriented towards the wall, i.e., downward, $\sin\varphi(t)\leq 0$. Regardless if it points downward or upward, the active rod’s orientation ${\varphi}$ is governed by the same equation as in the bulk. Note that, as it is often done in applications (see, e.g., [28]) we neglect inertia and translational diffusion in the individual based model.

Next, we compare results of Monte Carlo simulations for individual based model (38)-(39) with the initial boundary value problem derived in Theorem 2.1 from Section 3, consisting of the Fokker-Planck equation for probability distribution function $\rho(t,y,\varphi)$:

[TABLE]

and no-flux boundary condition:

[TABLE]

Translation diffusion coefficient is chosen to be small, $D_{\text{tr}}=0.05$.

Finally, we simulate (10), derived as the limit $D_{\text{tr}}\to 0$ in (40)-(41).

A sample trajectory obtained from simulating (38)-(39) for $0<t<20$ is depicted in Fig. 1, right. This trajectory is drawn in $\varphi y$-plane, and rod locations within this plane at integer moment of times are marked by red dots while the value of the corresponding moment of time is written above each dot. The trajectory demonstrates typical behavior of an active rod swimming at a wall. After collision with the wall, ($t=4$), the rod attaches to the wall (it still can swim in $x$ direction) and re-orients under background shear decreasing orientation angle $\varphi$ to $\varphi=-\pi$, and then detaches from the wall ($t\approx 5.5$). Swimming with orientation $\varphi$ close to $\pm\pi$ with background flow given by $\boldsymbol{u}_{\text{BG}}=(\dot{\gamma}y,0)$ means that the active rod swims upstream, that is, exhibits negative rheotaxis.

For Monte Carlo simulations of (38)-(39) we chose initial location $y$ and orientation angle $\varphi$ to be random with uniform distributions in intervals $0.25<y<1.25$ and $-3\pi/4<\varphi<-\pi/4$, respectively. The corresponding initial condition for both probability distribution functions $\rho^{\text{I}}$, solution of Fokker-Planck equation (40) with no flux boundary condition (41), and $\rho^{\text{II}}$, given by two-scale expansion (28), (29) with $\delta=D_{\text{tr}}=0.05$ and terms $\rho_{\text{bulk}}$ and $\psi_{\text{wall}}$ solving system (10), are

[TABLE]

The initial condition (42) is shown in Fig. 2, upper left.

Results of numerical simulations are depicted in Fig. 2 and 3. Behavior of solutions of all the three problems $-$ Monte Carlo simulations of individual based model (38)-(39), Fokker-Planck equation (40) with no flux boundary condition (41), and two-scale expansion (28), (29) with $\delta=D_{\text{tr}}=0.05$ and terms $\rho_{\text{bulk}}$ and $\psi_{\text{wall}}$ solving system (10) $-$ is qualitatively similar. Distributions of location $y$ and orientation $\varphi$ concentrate at $y\approx 0$ and $\varphi\approx-\pi$ (Fig. 2, lower row). These plots illustrate wall accumulation $y\approx 0$ and negative rheotaxis $\varphi\approx-\pi$. In this numerical example, initial condition (42) is chosen such that active rods initially point towards the wall, so the wall accumulation is somewhat enforced. Nevertheless, if domain is confined from all sides, which is a generic situation, an active rod reaches the wall with a high probability and spends a non-zero time at the wall reorienting before been detached. Therefore, the wall accumulation of active rods in confined bounded domains necessarily occurs, and the numerical example in this section investigates the situation when a population of active rods approach a wall.

We also analyzed the active rod distribution inside boundary layer $\mathcal{L}:=\{0~{}\leq~{}y~{}<~{}0.2\}$ (accumulated active rods). Comparison of angular distributions at $t=10$, obtained from the three approaches, is given in Fig. 2, upper right, and the probability of swimming inside boundary layer $\mathcal{L}$, as a function of time $0<t<20$, is depicted in Fig. 3. All the three methods show active rod accumulation increases up to a moment $3<t<5$ (depending on the value of rotational diffusion coefficient $D_{\text{rot}}$; see Fig. 3) and a peak of angular orientations forms close to $\phi=-\pi$ (see Fig. 2, upper right). However, the third method slightly underestimates probability of swimming inside the boundary layer $\mathcal{L}$ for larger values of $D_{\text{rot}}$ (specifically, for $D_{\text{rot}}=2.0$). This underlines the subtlety, described in the paragraph after (37), of the multi-scale expansion for $\rho$ with respect to $\delta$ and presence of the scales additional to $1$ and $\delta^{-1}$ in (27). Rigorous analysis of the limit $\delta\to 0$ of solution to problem (7)-(8) as well as explicit convergence estimates are left for our future work.

6 Individual based model and Fokker-Planck equation for an active rod with inertia

6.1 Individual based model for an active rod with inertia

In this section we present equations governing the motion of an active rod in the domain $\Omega\subset\mathbb{R}^{2}$ with non-empty boundary $\Gamma=\partial\Omega$. The rod is assumed to be a non-deformable one-dimensional segment of length $\ell$ swimming in a viscous fluid. At each moment of time $t$ the rod is characterized by the location of its center of mass $\boldsymbol{r}(t)$ and the orientation angle $\varphi(t)$, so that the unit vector $\boldsymbol{p}(\varphi)=(\cos\varphi,\sin\varphi)$ with $\varphi=\varphi(t)$ determines the orientation of the active rod. Equations for $\boldsymbol{r}(t)$ and $\varphi(t)$ are

[TABLE]

Here $\rho_{0}$ is the density of the rod, so $m=\rho_{0}\ell$ is its mass; $I_{\mathrm{rod}}=\rho_{0}\ell^{3}/12$ is the moment of inertia of the rod around the center of mass; $\eta_{1}$ and $\eta_{2}$ are material constants related to the background fluid viscosity. Function $\boldsymbol{u}_{\text{BG}}(\boldsymbol{r})$ is the velocity of the background flow at point $\boldsymbol{r}$. We assume that the background flow is not affected by the active rod thus $\boldsymbol{u}_{\text{BG}}$ satisfies the homogeneous Stokes equation in $\Omega$:

[TABLE]

$F_{\mathrm{thrust}}$ is the magnitude of the thrust force (the self-propulsion force); $D_{\text{tr}}$ and $D_{\text{rot}}$ are translational and rotational diffusion coefficients, respectively; $\zeta_{1}$ and $\zeta_{2}$ are uncorrelated white noises with intensities $\langle\zeta_{i}(t),\zeta_{i}(t^{\prime})\rangle=\delta(t-t^{\prime})$, $i=1,2$. Function $\Phi_{\mathrm{BG}}(\boldsymbol{r},\varphi)$ can be understood as the proportionality coefficient defining the torque due to the background flow exerted on the unit rod at location ${\bf r}$, with orientation angle $\varphi$ and zero angular velocity, and $\Phi_{\mathrm{BG}}(\boldsymbol{r},\varphi)$ is given by

[TABLE]

Equation (43) is the force balance for the active rod, and this equation reads as follows. There are three forces exerted on the active rod: (i) the viscous drag force which is proportional to the relative velocity of the active rod with respect to the background flow, (ii) the self-propulsion force directed always along the direction of the active rod, and (iii) the random (Brownian) force. Equation (44) is the torque balance for the active rod which implies that there are two torques re-orienting the active rod: viscous and random torques. Here we assume that the self-propulsion force does not affect orientation of the active rod.

The system (43)-(44) can be written in the following general form

[TABLE]

Here $u(\boldsymbol{r},\varphi,t)$ and $T(\boldsymbol{r},\varphi,t)$ are given smooth and bounded functions and $\varepsilon$ is a small parameter.

Next, we describe the rule of collision of the active rod and the wall $\Gamma$. At time of a collision, $t_{\text{coll}}$, an instantaneous force (an impulse) is exerted on the rod and this force is directed in direction normal to $\Gamma$ (see Fig. 4):

[TABLE]

where $\boldsymbol{n}$ is the outward normal vector and $\sigma$ will be determined below. One can formally add the force $F_{\text{coll}}$ to the right hand side of equation (43). This translates into the following collision rule for velocities

[TABLE]

where $\boldsymbol{v}=\dot{\boldsymbol{r}}(t_{\mathrm{coll}}-0)$ and $\boldsymbol{v}^{\prime}=\dot{\boldsymbol{r}}(t_{\mathrm{coll}}+0)$ are velocities before and after the collision, respectively. The force $F_{\mathrm{coll}}$ is exerted at the front of the active rod which touches the wall $\Gamma$. This implies that the active rod has an additional torque (in the right hand side of (44)) due to the collision, which equals to $-m\sigma\frac{\ell}{2}\,(\boldsymbol{p}\times\boldsymbol{n})\delta(t-t_{\text{coll}})$ and thus, the collision rule for angular velocities is

[TABLE]

where $\omega=\dot{\varphi}(t_{\mathrm{coll}}-0)$ and $\omega^{\prime}=\dot{\varphi}(t_{\mathrm{coll}}+0)$ are angular velocities before and after the collision, respectively, and $\text{e}_{z}=(0,0,1)$ is the unit vector orthogonal to the plane containing domain $\Omega$.

To completely describe the collision rule, the value of parameter $\sigma$ needs to be determined. To this end, we assume that the collision is perfectly elastic, that is, the kinetic energy does not change in time of the collision:

[TABLE]

Expanding the right hand side of (52) we get the formula for $\sigma$:

[TABLE]

Remark 1**.**

We can extend our consideration to imperfect elastic conditions, that is, we may assume that a part of energy is lost in time of collision:

[TABLE]

where $0\leq\epsilon\leq 1$ is coefficient of restitution which measures how elastic the collision is: if it is $1$, the collision is perfectly elastic; if $\epsilon=0$, the collision is perfectly inelastic.

Proposition 3**.**

Let $\mathcal{C}$ be the linear operator (a $3\times 3$ matrix)

[TABLE]

with $\sigma$ from (53). Then

[TABLE]

Proof.

Denote $v_{\mathrm{n}}:=\boldsymbol{v}\cdot\boldsymbol{n}$ and $v_{\tau}:=\boldsymbol{v}\cdot\boldsymbol{\tau}$ ($\boldsymbol{\tau}$ is the unit tangent vector on $\Gamma$). We also represent $\sigma$ as follows:

[TABLE]

In addition, we introduce angle $\theta$ between vectors $\boldsymbol{p}$ and $\boldsymbol{n}$. Note that $\sin\theta=(\boldsymbol{p}\times\boldsymbol{n})\cdot\mathrm{e}_{z}$.

Then the operator $\mathcal{C}$ can be represented by

[TABLE]

Finally, we compute $\mathrm{det}\,\mathcal{C}$:

[TABLE]

Thus, the proof of proposition is complete. ∎

6.2 Fokker-Planck equation for an active rod with inertia

For active rods whose motion inside domain $\Omega$ is described by (46)-(49), the Fokker-Planck equation is

[TABLE]

The unknown function $\tilde{f}_{\varepsilon}(t,\boldsymbol{r},\boldsymbol{v},\varphi,\omega)$ is the probability distribution function of the active rod’s location $\boldsymbol{r}\in\Omega$, translational velocity $\boldsymbol{v}\in\mathbb{R}^{2}$, orientation angle $\varphi\in[-\pi,\pi)$, and angular velocity $\omega\in\mathbb{R}$.

The collision rule of the active rods with the wall $\Gamma$ is given by (50)-(51). The rule translates into the following boundary conditions for $\tilde{f}_{\varepsilon}$:

[TABLE]

where $\tilde{f}_{\varepsilon}^{\prime}=\tilde{f}_{\varepsilon}(t,\boldsymbol{r},\boldsymbol{v}^{\prime},\varphi,\omega^{\prime})$ and the pair $(\boldsymbol{v}^{\prime},\omega^{\prime})$ is given by collision rule (50)-(51).

Remark 2**.**

The meaning of boundary condition (59) is as follows: the flux of incident active rods equal to the flux of reflected active rods. In other words, the flux is the same before and after collisions. Note that if one considers spherical particles with no preferred orientation (no $\varphi$ and $\omega$), then the collision rule is $\boldsymbol{v}^{\prime}=\boldsymbol{v}-2(\boldsymbol{v}\cdot\boldsymbol{n})\boldsymbol{n}$ and (59) is reduced in this case to equality of probability distribution functions $f_{\varepsilon}=f^{\prime}_{\varepsilon}$. This boundary condition is used in many works where particles are assumed to be spherical, see e.g. [31, 32, 33]. We point out here that imposing such boundary conditions (equality of probability distribution functions instead of fluxes as in (59)) for active rods leads to violation of both the mass conservation and the energy relation.

In this work we are interested in the limit $\varepsilon\to 0$. Following [31, 32, 33], to obtain a meaningful limit we first rescale $\tilde{f}_{\varepsilon}$:

[TABLE]

Then the Fokker-Planck equation for the rescaled probability distribution function $f_{\varepsilon}$ has the following form:

[TABLE]

with boundary condition (59) for $f_{\varepsilon}$.

Acknowledgment

The work of MP, LB, ER were supported by NSF DMREF grant DMS-1628411. PEJ was partially supported by NSF Grant 1614537, and NSF Grant RNMS (Ki-Net) 1107444.