Brownian Swarm Dynamics and Burgers’ Equation with Higher Order Dispersion

Abstract

The concept of ranked order probability distribution unveils natural probabilistic interpretations for the kink waves (and hence the solitons) solving higher order dispersive Burgers’ type PDEs. Thanks to this underlying structure, it is possible to propose a systematic derivation of exact solutions for PDEs with a quadratic nonlinearity of the Burgers’ type but with arbitrary dispersive orders. As illustrations, we revisit the dissipative Kotrweg de Vries, Kuramoto-Sivashinski, and Kawahara equations (involving third, fourth, and fifth order dispersion dynamics), which in this context appear to be nothing but the simplest special cases of this infinitely rich class of nonlinear evolutions.

1. Introduction

In the vast realm of nonlinear PDEs, the scalar Burgers’s (BU) Equation [1],

$${\partial }_{t}u(x,t)+u(x,t){\partial }_{x}u(x,t)+{\rho }_2{\partial }_x^{2}u(x,t)=0,$$

(1)

is doubtlessly commonly employed. Not only can BU be linearised (via the Hopf-Cole logaritmic transformation), but exact solutions in terms of kink and soliton travelling waves are also very easy to derive [1]. This makes it amazingly simple to observe the interplay between the nonlinearity and dispersion mechanisms. The higher order dispersive generalisation of BU, namely,

$${\partial }_{t}u(x,t)+u(x,t){\partial }_{x}u(x,t)+\sum _{j=2}^{n+1}{\rho }_j{\partial }_x^{j}u(x,t)=0,$$

(2)

has also sustained attention (Without lost of generality, the first oder term of the form ${\rho }_1{\partial }_{x}u(x,t)$ can always be removed from both Equations (1) and (2) via a Gallilean transformation of the form $t\mapsto t$ and $x\mapsto x-{\rho }_{1}t$.). Motivations to study Equation (2) are triggered by its numerous potential applications [2,3,4,5,6] (see review in [4]). In addition, using the so-called Tanh-method [7,8], Equation (2) also offers the possibility of displaying many exact and explicit solutions of the dynamics. In particular, a vast “zoology” of beautiful kink and soliton solutions of Equation (2) can be found in the available literature. We refrain here from proposing an extensive list of this vast corpus of contributions; rather, we focus on a small selection [9,10,11] directly relevant for our present purpose. Facing such a rich collection of solutions, one might perhaps feel a need for a unifying pattern and possibly a simple physical context helping to group many different kink solutions of Equation (2) under one footing; unveiling such a pattern is the aim of this paper. Our physical and mathematical inspiration emanates from studying Brownian swarm dynamics [12] with very large agent populations. Adopting a mean-field approach, the collective swarm dynamics can be stylised by hydrodynamic equations that possibly allow for exact solutions. Developed in Section 2, we recall (see also [13]) that, for simple mutual interactions of the type “catch the leader” or conversely “catch the laggards”, the resulting hydrodynamic evolution matches Equation (1). Such interactions can be viewed as special cases (i.e., limited to large swarms populations) of the more general class of dynamics introduced in [14,15]. In the specific context of Brownian agents where the evolution is Markovian, Equation (1) with specific boundary conditions can be alternatively interpreted as a nonlinear Fokker-Planck equation [16]. With this probabilistic interpretation, kink solutions directly describe travelling probability distributions of the form $P(\xi -vt):=P\left(x\right)$ with $\xi \in \mathbb{R}$ and v the travelling velocity. From the probability distributions $P\left(x\right)$, we can construct the class of ranked order probability (ROP) distributions $P_{k:n}\left(x\right)$ [17]. ROPs are derived by drawing independently from $P\left(x\right)$ a sample of n numbers, say $\left\{X_1,X_2,\dots ,X_n\right\}$, and then ranking them so that

$$X_{\left(1\right)}\le X_{\left(2\right)}\le \dots \le X_{\left(n\right)}.$$

The distribution $P_{k:n}\left(x\right)$ can then be expressed as [17]

$$P_{k:n}\left(x\right)=\mathrm{Prob}\left\{-\infty \le X_{\left(k\right)}\le x\right\}=\sum _{j=k}^n\left(\genfrac{}{}{0pt}{}{n}{j}\right)P^j\left(x\right){\left[1-P\left(x\right)\right]}^{n-j},k=1,2,\dots ,n.$$

(3)

Having introduced $P_{k:n}\left(x\right)$, we now can informally state the central result of the paper as follows:

Kink type probability distributions solving Equation (2) with specific coefficients ${\left\{{\rho }_j\right\}}_{j=2,\dots ,n+1}$ are nothing but ranked order distributions with a sampling size n derived from the kink type probability distribution solving Equation (1).

The paper is organised as follows: in Section 2, we introduce the Brownian swarm dynamics, its corresponding hydrodynamic picture, and the kink type probability distribution (i.e., here a logistic distribution) that emerges. Ranked order logistic distributions and some of their properties are briefly reviewed in Section 3. Section 4 contains the central result of the paper, as stated above. In Section 5, we list a collection of special cases of highly dispersive nonlinear evolutions classically encountered in mathematical physics. Namely, for $n=2$, $n=3$, and $n=4$, one recovers special cases of the dissipative Kortweg de Vries, Kuramoto-Sivashinski, and Kawahara nonlinear dynamics. Elementary but cumbersome technical details are systematically postponed to three appendices.

2. Brownian Swarms and Burgers’ Evolution

The dynamic of a swarm consists in N-interacting Brownian agents ${\mathcal{A}}_j$. $j=1,2,\dots ,N$ will be here described by a set of N stochastic differential equations (SDE) [18]:

$$\left\{\begin{array}{c}d{X}_j\left(t\right)=u\left[X_j\left(t\right),\mathbf{X}\left(t\right)\right]dt+\sigma d{W}_j\left(t\right),\hfill \\ \\ \mathbf{X}(t=0)={\mathbf{x}}_0,\hfill \end{array}\right.$$

(4)

where $\mathbf{X}\left(t\right):=\left(X_1\left(t\right),X_2\left(t\right),\dots ,X_N\left(t\right)\right)\in {\mathbb{R}}^N$, and $d{W}_j\left(t\right)$ are N-independent White Gaussian Noise (WGN) processes (i.e., formal derivatives of N-independent Brownian motions). The swarm is here homogeneous since both the drift $u\left[X_j\left(t\right),\mathbf{X}\left(t\right)\right]:\mathbb{R}\times {\mathbb{R}}^N\to \mathbb{R}$ and the noise amplitude $\sigma$ are assumed to be j-independent. In Equation (4), the WGN driving noise motivates the denomination Brownian swarms [12]. Since the drifts of the agents depend on $\mathbf{X}\left(t\right)$, they effectively stylise mutual interactions between the ${\mathcal{A}}_j$ values. To characterise the swarm’s collective evolution, we may define an empirical repartition density as

$$\mathcal{P}(x,t):=\frac{1}{N}\sum _{j=1}^N\delta (X_j\left(t\right)-x),$$

(5)

where $\delta \left(x\right)$ is the Dirac mass function. Focusing on large swarms (i.e., $N\to \infty$), it is legitimate to adopt a hydrodynamic description, write ${lim}_{N\to \infty }\mathcal{P}(x,t)\simeq P(x,t)$, and finally assume that $P(x,t)$ is ${\mathcal{C}}^2\left(\mathbb{R}\right)\times {\mathcal{C}}^1\left({\mathbb{R}}^{+}\right)$ is a probability density function, with

$$\left\{\begin{array}{c}{\int }_{\mathbb{R}}P(\xi ,t)d\xi =1,\hfill \\ \\ P(x,t=0)=p_0\left(x\right).\hfill \end{array}\right.$$

Since $P(x,t)$ is normalised density, the quantity $NP(x,t)dx$ stands therefore the instantaneous number of agents located within the infinitesimal interval $[x,x+dx]$. Under the former hypothesis, we further focus on the subclass of dynamics for which the interactions are expressible via a mean-field (MF) kernel:

$$\begin{array}{c}u\left[X_j\left(t\right),\mathbf{X}\left(t\right)\right]\equiv u\left[X_j\left(t\right),P(x,t)\right]\hfill \end{array}$$

(6)

Equation (6) describes agents mutually interacting via their own repartition density. The homogeneous character of the interactions of Equation (6) and the specific choice of WGN stochastic driving enable one to formally write the collective evolution by means of a nonlinear parabolic PDE (i.e., a nonlinear Fokker-Planck (FP) equation) [16]:

$$\left\{\begin{array}{c}{\partial }_{t}P(x,t)=-{\partial }_x\left\{u\left[x,P(x,t),t\right]p(x,t)\right\}+\frac{{\sigma }^2}{2}{\partial }_{xx}P(x,t).\hfill \\ \\ P(x,0)=p_0\left(x\right).\hfill \end{array}\right.$$

(7)

Inspired from [14,15], we now further specialise the dynamics and consider two specific types of mutual interactions [13]:

(a)

Catch the leader interactions (CLEA). In this case, at any time t, ${\mathcal{A}}_k$ at location x determines $n_k(x,t)$, which counts the number of ${\mathcal{A}}_k$ leaders. Knowing $n_k(x,t)$, ${\mathcal{A}}_k$ adjusts its drift according to the CLEA rule:

$$\left\{\begin{array}{c}{n}_k(x,t)=\sum _{j=1,j\ne k}^n\mathbb{I}\left\{x_j\left(t\right)\in \left[x,\infty \right]\right\},\hfill \\ \\ u_{\mathrm{CLEA}}\left[X_k\left(t\right),\mathbf{X}\left(t\right)\right] :=\frac{n_k(x,t)}{N}\hfill \end{array}\right.$$

(8)

where $\mathbb{I}\left\{x_j\left(t\right)\in \left[x,\infty \right]\right\}$ is the indicator function. For large swarms $N\to \infty$, the mean field (MF) interactions governing the evolution of a representative agent (we now drop the subscript k) is now written as

$$u_{\mathrm{CLEA}}\left[X_k\left(t\right),\mathbf{X}\left(t\right)\right]=\frac{n_k(x,t)}{N}=\frac{1}{N}{\sum }_{j=1,j\ne k}^n\mathbb{I}\left\{x_j\left(t\right)\in \left[x,\infty \right]\right\}\cong {\int }_x^{\infty }{P}_{\mathrm{CLEA}}(\xi ,t)d\xi .$$

(9)

Accordingly, the swarm’s evolution Equation (7) takes a special form:

$$\left\{\begin{array}{l}{\partial }_t{P}_{\mathrm{CLEA}}(x,t)=-{\partial }_x\left\{\left[{\int }_x^{+\infty }{P}_{\mathrm{CLEA}}(\xi ,t)d\xi \right]P_{\mathrm{CLEA}}(x,t)\right\}+\frac{{\sigma }^2}{2}{\partial }_x^2{P}_{\mathrm{CLEA}}(x,t),\hfill \\ P_{\mathrm{CLEA}}(x,0)=p_0\left(x\right),\\ P_{\mathrm{CLEA}}(x,t)\ge 0\mathrm{and}{\int }_{-\infty }^{+\infty }{P}_{\mathrm{CLEA}}(\xi ,t)d\xi =1.\end{array}\right.$$

(10)

(b)

Catch the laggard interactions (CLAG). Similarly, at any time t, an agent ${\mathcal{A}}_k$ at location x determines $n_k(x,t)$, which counts the number of ${\mathcal{A}}_k$ laggards. ${\mathcal{A}}_k$ then adjusts its drift according to the CLAG rule:

$$\left\{\begin{array}{c}{n}_k(x,t)=\sum _{j=1,j\ne k}^n\mathbb{I}\left\{x_j\left(t\right)\in \left[-\infty ,x\right]\right\},\hfill \\ \\ u_{\mathrm{CLAG}}\left[X_k\left(t\right),\mathbf{X}\left(t\right)\right] :=-\frac{n_k(x,t)}{N}\cong -{\int }_{-\infty }^x{P}_{\mathrm{CLAG}}(\xi ,t)d\xi ,\hfill \end{array}\right.$$

(11)

and the swarm’s evolution here is

$$\left\{\begin{array}{c}{\partial }_t{P}_{\mathrm{CLAG}}(x,t)={\partial }_x\left\{\left[{\int }_{-\infty }^x{P}_{\mathrm{CLAG}}(\xi ,t)d\xi \right]P_{\mathrm{CLAG}}(x,t)\right\}+\frac{{\sigma }^2}{2}{\partial }_x^2{P}_{\mathrm{CLAG}}(x,t),\hfill \\ \\ P_{\mathrm{CLAG}}(x,0)=p_0\left(x\right),\hfill \\ \\ P_{\mathrm{CLAG}}(x,t)\ge 0\mathrm{and}{\int }_{-\infty }^{+\infty }{P}_{\mathrm{CLAG}}(\xi ,t)d\xi =1.\hfill \end{array}\right.$$

(12)

Let us now introduce the following notations:

$$G(x,t):={\int }_x^{+\infty }{P}_{\mathrm{CLEA}}(\xi ,t)d\xi \mathrm{and}H(x,t):={\int }_{-\infty }^x{P}_{\mathrm{CLAG}}(\xi ,t)d\xi .$$

(13)

In terms of Equation (13), we now observe that both Equations (10) and (12) reduce to the Burgers’ equation, but they must fulfill different boundary conditions (BCs), namely,

$$\left\{\begin{array}{c}{P}_{\mathrm{CLEA}}(x,t)=-{\partial }_{x}G(x,t),\hfill \\ \\ {\partial }_t^{}G(x,t)=-G(x,t){\partial }_{x}G(x,t)+\frac{{\sigma }^2}{2}{\partial }_x^{2}G(x,t),\hfill \\ \\ G(-\infty ,t)=1\mathrm{and}G(+\infty ,t)=0\hfill \end{array}\right.$$

(14)

and, similarly,

$$\left\{\begin{array}{c}{P}_{\mathrm{CLAG}}(x,t)={\partial }_{x}H(x,t),\hfill \\ \\ {\partial }_t^{}H(x,t)=+H(x,t){\partial }_{x}H(x,t)+\frac{{\sigma }^2}{2}{\partial }_x^{2}H(x,t),\hfill \\ \\ H(-\infty ,t)=0\mathrm{and}H(+\infty ,t)=1.\hfill \end{array}\right.$$

(15)

Kink type travelling waves solving Equations (14) and (15) are easily derivable [1] and they respectively read

$$\left\{\begin{array}{c}G(x,t)=\frac{1}{2}-\frac{1}{2}tanh\left[\gamma \left(x-\frac{t}{2}\right)\right],\hfill \\ \\ H(x,t)=\frac{1}{2}+\frac{1}{2}tanh\left[\gamma \left(x-\frac{t}{2}\right)\right],\hfill \end{array}\right.$$

(16)

with $\gamma =\frac{1}{4{\sigma }^2}$. Writing $T:=tanh\left[\gamma \left(x-\frac{t}{2}\right)\right]$, we observe that Equations (14) and (15) enjoy the following property:

$$H\left(T\right)=1-G\left(T\right)=G(-T).$$

(17)

Note finally that, in this specific propagating mode, we have $P_{\mathrm{CLEA}}\left(T\right)=P_{\mathrm{CLAG}}\left(T\right)$, which is compatible with Equation (13).

3. Ranked Order Logistic Distribution

The specific BCs in Equations (14) and (15), implying that both $G\left(T\right)$ and $H\left(T\right)$ are probability distributions, allow us to now introduce sets of ROPs. Focusing on $G\left(T\right)$, we obtain ${\mathcal{G}}_n\left(T\right):={\left\{G_{k:n}\left(T\right)\right\}}_{k=1,2,\dots ,n}$:

$$\left\{\begin{array}{c}G\left(T\right)=\frac{1}{2}(1-T),\hfill \\ G_{k:n}\left(T\right):=\sum _{j=k}^n\left(\genfrac{}{}{0pt}{}{n}{j}\right){\left[\frac{1}{2}(1-T)\right]}^j{\left[\frac{1}{2}(1+T)\right]}^{n-j},k=1,2,\dots ,n,\hfill \\ T:=T(x,t)=\tanh \left[\frac{1}{4{\sigma }^2}(x-\frac{\mathrm{t}}{2})\right]\hfill \\ G_{k:n}(-1)=1\hspace{1em}\hspace{1em}\mathrm{and}\hspace{1em}\hspace{1em}{G}_{k:n}(+1)=0\hfill \end{array}\right.$$

(18)

and similarly for the set ${\mathcal{H}}_n\left(T\right):={\left\{H_{k:n}\left(T\right)\right\}}_{k=1,2,\dots ,n}$. The elements of the sets ${\mathcal{G}}_n\left(T\right)$ and ${\mathcal{H}}_n\left(T\right)$ are T-polynomials, all with degree n. From the definition, we have the following property:

Lemma 1.

$$G_{k:n}(-T)=1-G_{n-k+1:n}\left(T\right)=H_{k:n}\left(T\right).$$

(19)

Proof. 

$$G_{k:n}\left(T\right)=\sum _{j=k}^n\left(\genfrac{}{}{0pt}{}{n}{j}\right){\left[\frac{1}{2}(1-T)\right]}^j{\left[\frac{1}{2}(1+T)\right]}^{n-j}.$$

The identity ${(a+b)}^n={\sum }_{j=0}^n\left(\genfrac{}{}{0pt}{}{n}{j}\right)a^j{b}^{n-j}$ enables one to write

$$\begin{array}{cc}{G}_{k:n}(-T)\hfill & =\sum _{j=k}^n\left(\genfrac{}{}{0pt}{}{n}{j}\right){\left[\frac{1+T}{2}\right]}^j{\left[\frac{1-T}{2}\right]}^{n-j}=H_{k:n}\left(T\right)=\hfill \\ & {\sum }_{m=0}^{m=n-k}\left(\genfrac{}{}{0pt}{}{n}{n-m}\right){\left[\frac{1-T}{2}\right]}^m{\left[\frac{1+T}{2}\right]}^{n-m}=\hfill \\ & \frac{{\left[1+T+1-T\right]}^n}{2}-{\sum }_{m=n-k+1}^n\left(\genfrac{}{}{0pt}{}{n}{m}\right){\left[\frac{1-T}{2}\right]}^m{\left[\frac{1+T}{2}\right]}^{n-m}=1-G_{n-k+1:n}\left(T\right).\hfill \end{array}$$

□

Logistic Distribution

The logistic distribution $\mathcal{L}\left(z\right)$, discussed for example in [17] (see Chapter 5), reads as

$$\mathcal{L}\left(2z\right):=\frac{1}{1+e^{2z}}=\frac{e^{-z}}{\left[2cosh\left(z\right)\right]}=\frac{1}{2}[1-tanh\left(z\right)].$$

(20)

Using Equation (20), we can directly express the solutions of the Burgers’ equation, Equation (14), as $G\left(T\right)=\mathcal{L}\left[2\gamma (x-t/2)\right]$ (and similarly for Equation (15)). The associated ranked order logistic distribution is denoted as ${\mathcal{L}}_{k:n}\left[2\gamma (x-t/2)\right]$, and its moment generating function $M_{k:n}\left(u\right)$ is calculated in [17] (see Chapter 5) and reads

$$\begin{array}{c}{M}_{k:n}\left(u\right):=\mathbb{E}\left\{e^{u\left[x-\frac{t}{2}\right]}\right\}={\int }_{\mathbb{R}}{e}^{u\left[x-\frac{t}{2}\right]}d{G}_{k:n}\left(T\right)={\int }_{\mathbb{R}}{e}^{u\left[x-\frac{t}{2}\right]}d{\mathcal{L}}_{k:n}\left[2\gamma \left(x-\frac{t}{2}\right)\right]=\hfill \\ \\ {\int }_{\mathbb{R}}{e}^{2\gamma u\zeta }d{\mathcal{L}}_{k:n}\left(\zeta \right)=\frac{\Gamma \left(k+\frac{u}{2\gamma }\right)\Gamma \left(n-k+1-\frac{u}{2\gamma }\right)}{\Gamma \left(k\right)\Gamma (n-k+1)},\left(\gamma =\frac{1}{4{\sigma }^2}\right).\hfill \end{array}$$

(21)

In particular, for the mean and the variance, we have [17] (see Chapter 5)

$$\begin{array}{c}{\mu }_{k:n}\left(t\right)=\frac{t}{2}+\frac{1}{2\gamma }\left[\Psi \left(k\right)-\Psi (n-k+1)\right],\hfill \\ \\ {\sigma }_{k:n}^2=\frac{1}{4{\gamma }^2}\left[{\Psi }^{\left(1\right)}\left(k\right)+{\Psi }^{\left(1\right)})(n-k+1)\right],\hfill \end{array}$$

(22)

where $\Psi \left(z\right):=\frac{d}{dz}log[\Gamma \left(z\right)]$ is the digamma function, and ${\Psi }^{\left(1\right)}\left(z\right)=\frac{d^2}{d{z}^2}\Psi \left(z\right)$.

From Equation (22), we have

$$\begin{array}{c}{\mu }_{k:n}\left(t\right)=t-{\mu }_{n-k+1:n}\left(t\right),\hfill \\ {\sigma }_{k:n}^2={\sigma }_{n-k+1:n}^2.\hfill \end{array}$$

(23)

Remark 1.

Consistent with basic intuition, we observe that the higher the sampling size n is, the smaller the resulting variance is. From Equation (23), we also observe that, for odd sampling n, we have

$$\begin{array}{c}{\mu }_{\frac{1}{2}(n+1):n}\left(t\right)=\frac{1}{2}t\Rightarrow d{G}_{\frac{1}{2}(n+1):n}\left[\gamma (x-t/2)\right]=:\hfill \\ \\ g_{\frac{1}{2}(n+1):n}\left[\gamma (x-t/2)\right]dx=g_{\frac{1}{2}(n+1):n}[-\gamma (x-t/2)]dx,\forall n\mathrm{odd},\hfill \end{array}$$

(24)

with $g_{\frac{1}{2}(n+1):n}(x,t)$ being the corresponding ranked order probability densities. Hence, for odd sampling size n, the middle position is characterised as $g_{\frac{1}{2}(n+1):n}(x,t)$, which is described by a symmetric probability density (here, a symmetric soliton). Conversely, for arbitrary n and $k\le n$, the corresponding probability densities $g_{k:n}(x,t)$ deriving from Equations (14) and (15) propagate as skew solitons.

Remark 2.

The average distance $\Delta \left(n\right)$ separating the laggard and the leader positions reads as

$$\Delta \left(n\right):={\mu }_{n:n}\left(t\right)-{\mu }_{1:n}\left(t\right)=\frac{1}{\gamma }\left[\Psi \left(n\right)-\Psi \left(1\right)\right],n=1,2,\dots .$$

(25)

Note that, since we have a kink type evolution (i.e., a stationary regime), $\Delta \left(n\right)$ is necessarily time-independent. It is monotonously increasing with the sampling size n, and, for large samplings n, we have $\Psi \left(n\right)\simeq ln\left(n\right)$, implying that $\Delta \left(n\right)\simeq \frac{1}{\gamma }ln\left(n\right)$.

4. Nonlinear Evolution Equations Solved by Ranked Order Distributions

We now raise questions regarding the dynamics of the ranked order distributions $G_{k:n}\left(T\right)$ derived from the Burgers’ kink type solutions, Equations (14) and (15). Specifically, we are looking for the sets of PDEs describing the evolution of $G_{k:n}\left(T\right)$ (respectively $H_{k:n}\left(T\right)$) for $k=1,2,\dots ,n$. Since we already know that both $G_{k:n}\left(T\right)$ and $H_{k:n}\left(T\right)$ are T-polynomials of degree n, it is natural to invoke the well known Tanh-method [7,8], commonly used to derive solutions of nonlinear PDEs. This enables us to assert the following:

Proposition 1.

The set ${\mathcal{G}}_n\left(T\right)$ consisting of the n ranked order distributions ${\left\{G_{k:n}\left(T\right)\right\}}_{k=1,\dots ,n}$ with $T:=tanh\left[\frac{1}{4{\sigma }^2}\left(x-\frac{t}{2}\right)\right]$ solves a set of n distinct nth-order dispersion PDEs:

$$\begin{array}{c}{\partial }_t{G}_{k:n}\left(T\right)+G_{k:n}\left(T\right)\left[{\partial }_x{G}_{k:n}\left(T\right)\right]+\sum _{j=2}^{n+1}\left[{\rho }_{j,\left\{k:n\right\}}\left(\sigma \right)\right]{\partial }_x^j{G}_{k:n}\left(T\right)=0,n=1,2,3,\dots .\hfill \end{array}$$

(26)

The sets of coefficients ${\left\{{\rho }_{j,\left\{k:n\right\}}\left(\sigma \right)\right\}}_{j=2,\dots ,n+1}$ are solutions of n distinct sets of $(2n+2)$ nonlinear algebraic relations.

Proof. 

Consider a high order dispersive PDE of the type Equation (26):

$$\begin{array}{c}\left\{\begin{array}{c}{\partial }_t{\mathcal{P}}_n\left(T\right)+{\mathcal{P}}_n\left(T\right)[{\partial }_x{\mathcal{P}}_n\left(T\right)+\sum _{j=2}^{n+1}\left[{\rho }_j\right]{\partial }_x^j{\mathcal{P}}_n\left(T\right)=0,n=1,2,3,\dots ,\hfill \\ \\ T:=T(x,t):=tanh\left[\gamma \left(x-\frac{t}{2}\right)\right]:=tanh\left(\zeta \right),\hfill \\ \\ {\mathcal{P}}_n(-\infty )=1\mathrm{and}{\mathcal{P}}_n(+\infty )=0.\hfill \end{array}\right.\hfill \end{array}$$

(27)

Assume that ${\mathcal{P}}_n\left(T\right)$ is an nth degree T-polynomial solution of Equation (27). As a function of $\zeta =\gamma (x-\frac{t}{2})$ itself, Equation (27) is the ODE:

$$-\frac{\gamma }{2}{\partial }_{\zeta }{\mathcal{P}}_n\left(T\right)+\gamma {\mathcal{P}}_n\left(T\right)\left[{\partial }_{\zeta }{\mathcal{P}}_n\left(T\right)\right]+\sum _{j=2}^{n+1}\left[{\rho }_j{\gamma }^j\right]{\partial }_{\zeta }^j{\mathcal{P}}_n\left(T\right)=0,n=1,2,3,\dots ,$$

(28)

After integrating once this ODE with respect to $\zeta$ (and with zero integration constant) and imposing the boundary condition ${\mathcal{P}}_n(-\infty )=1$ and ${\mathcal{P}}_n(\infty )=0$, we immediately verify that the kink’s travelling velocity is indeed $\frac{1}{2}$.

The polynomial ${\mathcal{P}}_n\left(T\right)$ is defined via a set of $n+1$ coefficients $a_j$ for $j=0,1,\dots ,n$ and the single parameter $\gamma$. Therefore, the set $\Omega :=\left\{\gamma ,{\left\{a_j\right\}}_{j=0,1,\dots ,n},{\left\{{\rho }_j\right\}}_{j=2,\dots ,n+1}\right\}$ contains $(2n+1)$ parameters. The Tanh-method [7,8] consists in introducing an $n^{\mathrm{th}}$ degree T-polynomial ${\mathcal{P}}_n\left(T\right)$ into the evolution Equation (27) and in successively balancing all identical $T^j$ contributions. This leads to a system of $2n+2$ nonlinear algebraic relations that connect the $2n+1$ parameters of the set $\Omega$. Since we have more relations than parameters, there is a priori no guarantee that, in general, the Tanh-method [7] actually provides a solution.

However, the ranked order set ${\mathcal{G}}_n\left(T\right):={\left\{G_{k:n}\left(T\right)\right\}}_{k=1,2\dots ,n}$ is a very specific subclass of dynamics enjoying an extra symmetry structure. As seen before, ${\mathcal{G}}_n\left(T\right)$ consists of n distinct $n^{\mathrm{th}}$ degree T-polynomials. Therefore, for ${\mathcal{G}}_n\left(T\right)$, we actually have n times sets of the $\Omega$-type, thus leading to $n\times (2n+1)=2n^2+n$ coefficients. Applying the Tanh-method to the set ${\mathcal{G}}_n\left(T\right)$, we shall therefore obtain $n\times (2n+2)=2n^2+2n$ algebraic relations. However, for the specific set ${\mathcal{G}}_n\left(T\right)$, Equation (19) imposes n extra symmetry relations. Taking into account this extra symmetry, we conclude that, for ${\mathcal{G}}_n\left(T\right)$, we effectively end with $(2n^2+2n)-n=2n^2+n$ nonlinear relations to determine $2n^2+n$ parameters. Hence, for ${\mathcal{G}}_n\left(T\right)$, the number of relations now matches the number of parameters. Due to the nonlinearity, this system of relations still does not necessarily possess a solution in general. However, focusing on ${\mathcal{G}}_n\left(T\right)$, we consider the set of coefficients ${\left\{a_{j,k:n}\right\}}_{j=0,1,\dots ,n}$, implicitly defined as

$${\mathcal{P}}_n\left(T\right)=\sum _{j=0}^n{a}_{j,k:n}{T}^j:=G_{k:n}\left(T\right)=\frac{1}{2^n}\sum _{j=k}^n\left(\genfrac{}{}{0pt}{}{n}{j}\right){(1-T)}^j{(1+T)}^{n-j}.$$

(29)

Going back to Equation (28), let us now adopt a similar notation and write ${\rho }_j\mapsto {\rho }_{j;k:n}$ to notify the specific evolution of interest. Apply the Tanh-method, by introducing Equation (29) into Equation (28) and use the specific coefficients ${\left\{a_{j,k:n}\right\}}_{j=0,1,\dots ,n}$. The successive balancing of the $T^j$ contributions produces only linear relations between the rescaled coefficients ${\widehat{\rho }}_{j,k:n}:={\rho }_{j,k:n}{\gamma }^{j-1}$. Linearity is due to the fact that, in Equation (27), only a linear superposition of dispersive terms occurs. Summarising, for ${\mathcal{G}}_n\left(T\right)$, the number of algebraic relations matches the number of parameters, and the choice of coefficients specified by Equation (29) produces linear algebraic sub-systems connecting the remaining unknowns ${\widehat{\rho }}_{j,k:n}$. Hence, one can conclude that the Tanh-method always yields a kink type solution for the specific subclass of ${\mathcal{G}}_n\left(T\right)$ dynamics. □

Remark 3.

Using Proposition 1, we see that the travelling waves ${\partial }_x{G}_{k:n}\left(T\right)$ are a superposition of a bell and a kink profile waves and hence are skew solitons. The skewness which is due to the kink profile component results exclusively from the underlying ranked order mechanism. Several explicit illustrations based on the method exposed in the proof of Proposition 1 will be found in Section 5 and the corresponding appendices.

Corollary 1.

(Antisymmetry between the dispersion coefficients)

The ranked order travelling kink $G_{k:n}\left(T\right)$ obeys the higher order dispersive Burgers’ equation of Equation (26), in which the dispersive coefficients ${\rho }_{j,k:n}\left(\sigma \right)$ satisfy an antisymmetry relation:

$${\rho }_{j,k:n}={(-1)}^j{\rho }_{j,n-k+1:n}$$

(30)

Proof. 

(Corollary 1)

Define $G_{k:n}\left(T\right)=1-\phi (-T)$. Therefore,

$$\left\{\begin{array}{c}{\partial }_t{G}_{k:n}\left(T\right)+G_{k:n}\left(T\right)\left[{\partial }_x{G}_{k:n}\left(T\right)\right]+\sum _{j=2}^{n+1}\left[{\rho }_{j,k:n}\right]{\partial }_x^j{G}_{k:n}\left(T\right)=0,\hfill \\ \\ {\partial }_t\phi (-T)+{\partial }_x\phi (-T)-\phi \left(T\right)\left[{\partial }_x\phi (-T)\right]+\sum _{j=2}^{n+1}\left[{\rho }_{j,k:n}\right]{\partial }_x^j\phi (-T)=0,\hfill \end{array}\right.$$

(31)

Introduce the change of variables:

$$\left\{\begin{array}{c}t\mapsto \tau =t,{\partial }_t\mapsto {\partial }_{\tau }+{\partial }_z\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\Rightarrow \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\\ x\mapsto z=-x+t,{\partial }_x\mapsto -{\partial }_z\hfill \end{array}\right.$$

(32)

Equation (32) implies $T\mapsto -T$. With Equation (32), the second line of Equation (31) becomes

$${\partial }_{\tau }\phi \left(T\right)+\phi \left(T\right)\left[{\partial }_z\phi \left(T\right)\right]+\sum _{j=2}^{n+1}\left[{(-1)}^j{\rho }_{j,k:n}\right]{\partial }_z^j\phi \left(T\right)=0.$$

(33)

In particular, for the choice $\phi \left(T\right)=G_{n-k+1:n}\left(T\right)$, Equation (33) implies

$${\partial }_{\tau }{G}_{n-k+1:n}\left(T\right)+G_{n-k+1:n}\left(T\right)\left[{\partial }_z{G}_{n-k+1:n}\left(T\right)\right]+\sum _{j=2}^{n+1}\left[{(-1)}^j{\rho }_{j,k:n}\right]{\partial }_z^j{G}_{n-k+1:n}\left(T\right)=0.$$

From the definition of $G_{n-k+1:n}\left(T\right)$, we can alternatively write

$${\partial }_{\tau }{G}_{n-k+1:n}\left(T\right)+G_{n-k+1:n}\left(T\right)\left[{\partial }_z{G}_{n-k+1:n}\left(T\right)\right]+\sum _{j=2}^{n+1}\left[{\rho }_{j,n-k+1:n}\right]{\partial }_z^j{G}_{n-k+1:n}\left(T\right)=0.$$

Therefore, we conclude that ${(-1)}^j{\rho }_{j,k:n}={\rho }_{j,n-k+1:n}$.

□

Corollary 2. 

(Parity relation)

For n odd, we have ${\rho }_{j,\frac{n+1}{2}:n}\left(\sigma \right)=0$ for all j odd.

Proof. 

(Corollary 2).

Taking $k=\frac{n+1}{2}$ in Equation (30) of Corollary 1, we have ${\rho }_{j,\frac{n+1}{2}:n}={(-1)}^j{\rho }_{j,\frac{n+1}{2}:n}$. This implies that ${\rho }_{j,\frac{n+1}{2}:n}=0$ for any odd j. □

While the explicit dispersion amplitudes ${\left\{{\rho }_{j,k:n}\right\}}_{j=2,\dots ,n+1}$ are generally cumbersome to write down, the Tanh-method yields a simple form for the highest dispersion amplitude ${\rho }_{n+1,k:n}$. Specifically, we obtain the following:

Corollary 3.

(Recurrence relations for the highest dispersion coefficients)

$${\rho }_{n+1,k:n}\left(\sigma \right)={(-1)}^{n+1}\left[\sum _{j=k}^n{(-1)}^j\left(\genfrac{}{}{0pt}{}{n}{j}\right)\right]\frac{2^n\Gamma \left(n\right)}{\Gamma \left(2n\right)}{\sigma }^{2n}.$$

(34)

Proof. 

(Corollary 3)

Denote ${\varphi }_n^{\left(l\right)}\left(T\right)$ for $l=1,2,3$ three unspecified T-polynomials of degree n. From the definition of the ranked order distributions, we can write

$$G_{k:n}\left(T\right)=\left[\sum _{j=k}^n{(-1)}^j\left(\genfrac{}{}{0pt}{}{n}{j}\right)\right]\frac{T^n}{2^n}+{\varphi }_{n-1}^{\left(1\right)}\left(T\right),$$

which implies

$$\left\{\begin{array}{l}{G}_{k:n}\left(T\right)\left[{\partial }_x{G}_{k:n}\left(T\right)\right]=-n{\left[{\sum }_{j=k}^n{(-1)}^j\left(\genfrac{}{}{0pt}{}{n}{j}\right)\right]}^2\gamma \frac{T^{2n+1}}{2^{2n}}+{\varphi }_{2n}^{\left(2\right)}\left(T\right),\left(\gamma =\frac{1}{4{\sigma }^2}\right),\\ {\partial }_x^{n+1}{G}_{k:n}\left(T\right)={(-1)}^{n+1}\frac{\Gamma \left(2n\right)}{\Gamma (n-1)}\left[{\sum }_{j=k}^n{(-1)}^j\left(\genfrac{}{}{0pt}{}{n}{j}\right)\right]{\gamma }^{n+1}\frac{T^{2n+1}}{2^n}+{\varphi }_{2n}^{\left(3\right)}\left(T\right).\end{array}\right.$$

According to the Tanh-method, we introduce $G_{k:n}\left(T\right)$ into Equation (26) and balance the contributions of degree $(2n+1)$. This leads to

$$-n{\left[\sum _{j=k}^n{(-1)}^j\left(\genfrac{}{}{0pt}{}{n}{j}\right)\right]}^2\gamma \frac{T^{2n+1}}{2^{2n}}+{\rho }_{n+1,k:n}\left(\sigma \right){(-1)}^{n+1}\frac{\Gamma \left(2n\right)}{\Gamma (n-1)}\left[\sum _{j=k}^n{(-1)}^j\left(\genfrac{}{}{0pt}{}{n}{j}\right)\right]{\gamma }^{n+1}\frac{T^{2n+1}}{2^n}=0,$$

from which Equation (34) follows directly. □

5. Illustrations

The lower order sampling cases $n=2$, $n=3$, and $n=4$ coincide with well studied nonlinear high order dispersive evolutions of mathematical physics. Below we list simple illustrations and postpone all details of calculations and technicalities to Appendix A, Appendix B and Appendix C.

5.1. Dissipative Kortweg de Vries Dynamics (Case n = 2 ⇒ Third Order Dispersion)

Here the equation of interest reads as

$$\left\{\begin{array}{l}{\partial }_t{G}_{1:2}\left(T\right)+G_{1:2}\left(T\right)\left[{\partial }_x{G}_{1:2}\left(T\right)\right]+{\rho }_{2,\left\{1:2\right\}}\left(\sigma \right){\partial }_x^2\left[G_{1:2}\left(T\right)\right]+{\rho }_{3,\left\{1:2\right\}}\left(\sigma \right){\partial }_x^3\left[G_{1:2}\left(T\right)\right]=0,\hfill \\ G_{1:2}\left(T\right)=\frac{1}{4}\left[3-2T-T^2\right],T=tanh\left[\frac{1}{4{\sigma }^2}(x-\frac{t}{2})\right],\hfill \\ {\rho }_{2,\left\{1:2\right\}}\left(\sigma \right)=-\frac{5}{6}{\sigma }^2,\\ {\rho }_{3,\left\{1:2\right\}}\left(\sigma \right)=(2-1)\frac{2^2\Gamma \left(2\right)}{\Gamma \left(4\right)}{\sigma }^4=\frac{{\sigma }^4}{3}.\end{array}\right.$$

(35)

Using Equation (30) of Corollary 1, we also immediately have

$$\left\{\begin{array}{l}{\partial }_t{G}_{2:2}\left(T\right)+G_{2:2}\left(T\right)\left[{\partial }_x{G}_{2:2}\left(T\right)\right]{\rho }_{2,\left\{2:2\right\}}\left(\sigma \right){\partial }_x^2\left[G_{2:2}\left(T\right)\right]+{\rho }_{3,\left\{2:2\right\}}\left(\sigma \right){\partial }_x^3\left[G_{2:2}\left(T\right)\right]=0,\hfill \\ G_{2:2}\left(T\right)=\frac{1}{4}\left[1-2T+T^2\right]T=tanh\left[\frac{1}{4{\sigma }^2}(x-\frac{t}{2})\right],\\ {\rho }_{2,\left\{2:2\right\}}\left(\sigma \right)=\frac{5}{6}{\sigma }^2,\\ {\rho }_{3,\left\{2:2\right\}}\left(\sigma \right)=-\frac{2^2\Gamma \left(2\right)}{\Gamma \left(4\right)}{\sigma }^4=-\frac{{\sigma }^4}{3}.\end{array}\right.$$

(36)

The solution $G_{1:2}\left(T\right)$ has been previously derived in [10]. The dissipative Kortweg de Vries pde also referred as the Burgers-Kortweg de Vries equation is among the simplest evolution where nonlinearity, dispersion and dissipation coexist and an extensive list of physically relevant contexts is given in [10].

5.2. Kuramoto-Sivashinsky (KS) Dynamics (Case n = 3 ⇒ Fourth Order Dispersion)

$$\left\{\begin{array}{c}{\partial }_t{G}_{2:3}\left(T\right)+G_{2:3}\left(T\right)\left[{\partial }_x{G}_{2:3}(x,t)\right]+\sum _{j=2}^4{\rho }_{j,\left\{2:3\right\}}\left(\sigma \right){\partial }_x^j\left[G_{2:3}\left(T\right)\right]=0,\hfill \\ G_{2:3}(x,t)=\frac{1}{4}\left[T^3-3T+2\right],T=tanh\left[\frac{1}{4{\sigma }^2}(x-\frac{t}{2})\right],\hfill \\ {\rho }_{2,\left\{2:3\right\}}\left(\sigma \right)=-\frac{19}{3}{\sigma }^2,\hfill \\ {\rho }_{3,\left\{2:3\right\}}\left(\sigma \right)=0,\hfill \\ {\rho }_{4,\left\{2:3\right\}}\left(\sigma \right)=2\left[\frac{2^3\Gamma \left(3\right)}{\Gamma \left(6\right)}\right]{\sigma }^6=\frac{2}{15}{\sigma }^6.\hfill \end{array}\right.$$

(37)

and

$$\left\{\begin{array}{c}{\partial }_t{G}_{3:3}\left(T\right)+G_{3:3}\left(T\right)\left[{\partial }_x{G}_{3:3}\left(T\right)\right]+\sum _{j=2}^4{\rho }_{j,\left\{3:3\right\}}\left(\sigma \right){\partial }_x^j\left[G_{3:3}\left(T\right)\right]=0,\hfill \\ \\ G_{3:3}\left(T\right)={\left[\frac{1-T}{2}\right]}^3,T=tanh\left[\frac{1}{4{\sigma }^2}(x-\frac{t}{2})\right],\hfill \\ \\ {\rho }_{2,\left\{3:3\right\}}\left(\sigma \right)=-\frac{47}{60}{\sigma }^2,\hfill \\ \\ {\rho }_{3,\left\{3:3\right\}}\left(\sigma \right)=-\frac{2}{5}{\sigma }^4,\hfill \\ \\ {\rho }_{4,\left\{3:3\right\}}\left(\sigma \right)=-\frac{2^3\Gamma \left(3\right)}{\Gamma \left(6\right)}{\sigma }^6=-\frac{1}{15}{\sigma }^6.\hfill \end{array}\right.$$

(38)

Observe that ${\rho }_{3,\left\{2:3\right\}}\left(\sigma \right)=0$ in accordance with the parity constraint expressed in Corollary 2. Finally, the use of Equation (30) in Corollary 1 yields immediately ${\left\{{\rho }_{j,\left\{1:3\right\}}\right\}}_{j=1,2,3}$ for the dispersion amplitudes of $G_{1:3}\left(T\right)$.

Strictly speaking the KS system does not includes a third order derivatives. Whenever one is interested in evolutions which simultaneously include nonlinearity, dispersion, dissipation and instability, the simplest pde possible includes all derivatives up to fourth order as in Equations (37) and (38). This generalised version is known as the Kuramto-Sivanshinsky-Benney (KSB) equation [19]. It was initially introduced to describe the nonlinear evolution of a fluid flowing on an inclined plane [20]. The instability present in KSB is due to the fact that in this case after integrating once Equation (28), one ends with third order nonlinear ode’s for which chaotic evolutions exist [19].

5.3. The Kawahara Fifth Order Dispersive Dynamics (Case n = 4 ⇒ Fifth Order Dispersion)

Since the full exact expressions become rapidly cumbersome, we limit here to the $G_{4:4}\left(T\right)$ kink type solution. The solution $G_{4:4}\left(T\right)$ has been derived in [9].

$$\left\{\begin{array}{l}{\partial }_t{G}_{4:4}\left(T\right)+G_{4:4}\left(T\right)\left[{\partial }_x{G}_{4:4}\left(T\right)\right]+\sum _{j=2}^5{\rho }_{j,\left\{4:4\right\}}\left(\sigma \right){\partial }_x^j\left[G_{4:4}\left(T\right)\right]=0,\hfill \\ G_{4:4}\left(T\right)={\left[\frac{1-T}{2}\right]}^4,T:=tanh\left[\frac{1}{4{\sigma }^2}\left(x-\frac{t}{2}\right)\right],\hfill \\ {\rho }_{2,\left\{4:4\right\}}\left(\sigma \right)=-\frac{29}{40}{\sigma }^2,\\ {\rho }_{3,\left\{4:4\right\}}\left(\sigma \right)=-\frac{179}{480}{\sigma }^4,\\ {\rho }_{4,\left\{4:4\right\}}\left(\sigma \right)=-\frac{11}{120}{\sigma }^6,\\ {\rho }_{5,\left\{4:4\right\}}\left(\sigma \right)=-\frac{2^4\Gamma \left(4\right)}{\Gamma \left(8\right)}{\sigma }^8=-\frac{1}{120}{\sigma }^8,\end{array}\right.$$

(39)

Similar results follow for

$$\left\{\begin{array}{c}{G}_{3:4}\left(T\right)=5-12T+6T^2+4T^3-3T^4,\hfill \\ \\ G_{2:4}\left(T\right)=11-12T-6T^2+4T^3+3T^4=1-G_{3:4}(-T),\hfill \\ \\ G_{1:4}\left(T\right)=15-4T-6T^2-4T^3-T^4=1-G_{4:4}(-T).\hfill \end{array}\right.$$

In particular, by using Equation (30) of Corollary 1, we can read the coefficients defining the evolution of $G_{1:4}\left(T\right)$. In addition, from Corollary 3, we directly obtain the higher order dispersive coefficients:

$$\left\{\begin{array}{l}{\rho }_{5,\left\{1:4\right\}}\left(\sigma \right)=\frac{2^4\Gamma \left(4\right)}{\Gamma \left(8\right)}{\sigma }^8=\frac{1}{120}{\sigma }^8=-{\rho }_{5,\left\{4:4\right\}}\left(\sigma \right),\hfill \\ {\rho }_{5,\left\{2:4\right\}}\left(\sigma \right)=3\frac{2^4\Gamma \left(4\right)}{\Gamma \left(8\right)}{\sigma }^8=\frac{1}{30}{\sigma }^8=-{\rho }_{5,\left\{3:4\right\}}\left(\sigma \right),\hfill \\ {\rho }_{5,\left\{3:4\right\}}\left(\sigma \right)=-3\frac{2^4\Gamma \left(4\right)}{\Gamma \left(8\right)}{\sigma }^8=-\frac{1}{30}{\sigma }^8,\\ {\rho }_{5,\left\{4:4\right\}}\left(\sigma \right)=-\frac{2^4\Gamma \left(4\right)}{\Gamma \left(8\right)}{\sigma }^8=-\frac{1}{120}{\sigma }^8.\end{array}\right.$$

The Kawahara fifth order dispersive dynamics is relevant to model surface and internal waves [6] and for nonlinear waves in a viscoelastic tube [9]. Here Equation (28) reduces to a nonlinear fourth order ode for which a foison of different behaviours exist and in particular including chaotic ones.

6. Conclusions and Perspectives

The rich modelling platform offered by the Burgers’ pde is here used in the context of swarms of Brownian agents. Adopting an hydrodynamic picture the agents’ spatial probability distributions obey to a Burger’s equation. The ranked order interactions of the type “catch the leader” or “catch the laggard” generate the Burgers’ quadratic nonlinearity and impose specific boundary conditions. It results that the swarm evolution solving the Burgers’ eq. exhibits a $tanh$-kink type traveling wave. The joint presence of ranked order interactions and probability distributions trigger natural questions regarding the evolution of ranked order distributions (ROD) based on the nominal kink solution of the Burgers’ eq. In a traveling reference frame, this Burger’s kink wave is a stationary and so will necessarily be all the derived ROD’s. Moreover as shown here, the ROD are also traveling kinks but solving Burger’s dynamics with higher orders dispersion terms (HODBU). This last amazing property emanates from the rather exceptional properties of the nominal $tanh$-kink (logistic type) distribution. Indeed the $tanh$ functions are solutions of Riccati equations, a property at the cornerstone of the $tanh$-method developed to solved nonlinear pde’s and in particular the HODBU. In the light of these observations, one may rise questions about potential feasibility of this approach either to other stochastic Markovian agents, (Cauchy agents for examples) or to vectorial Brownian agents. In this last case, higher orders Burgers’ eqs. will be involved and obviously the concept of spatial order for these higher dimension cases has to be suitably redefined. Such investigations remain so far fully open. In addition, let us point out some byproducts and other perspectives:

(a)

The HODBU evolutions possess a direct probabilistic interpretation of and hence their associated PDEs enjoy the property of positivity conservation.

(b)

A physically intuitive and particularly simple interpretation is immediately available for the kink type solutions.

(c)

The HODBU kinks are generally skewed and the origin of the skewness is clearly understood from the underlying construction of the ranked order distributions.

(d)

The unveiled ranked order structure opens imagination to write down further nonlinear evolution. For example the corresponding PDEs for joint ranked order distributions as defined in [17].

(e)

The intimate relation with swarm dynamics opens possibilities for applications. In the domain of mean-field games for example, the HODBU kink type solutions can be interpreted as the quasi-ergodic states of games jointly solving a Fokker-Planck and a Hamilton-Jacobi-Belmann system of PDEs [21].