# Strongly Nonlinear Stochastic Processes in Physics and the Life Sciences

Authors: T. D. Frank
Journal: ISRN Mathematical Physics
 (2013)

---

## Abstract

Strongly nonlinear stochastic processes can be found in many applications in physics and the life sciences. In particular, in physics, strongly nonlinear stochastic processes play an important role in understanding nonlinear Markov diffusion processes and have frequently been used to describe order-disorder phase transitions of equilibrium and nonequilibrium systems. However, diffusion processes represent only one class of strongly nonlinear stochastic processes out of four fundamental classes of time-discrete and time-continuous processes evolving on discrete and continuous state spaces. Moreover, strongly nonlinear stochastic processes appear both as Markov and non-Markovian processes. In this paper the full spectrum of strongly nonlinear stochastic processes is presented. Not only are processes presented that are defined by nonlinear diffusion and nonlinear Fokker-Planck equations but also processes are discussed that are defined by nonlinear Markov chains, nonlinear master equations, and strongly nonlinear stochastic iterative maps. Markovian as well as non-Markovian processes are considered. Applications range from classical fields of physics such as astrophysics, accelerator physics, order-disorder phase transitions of liquids, material physics of porous media, quantum mechanical descriptions, and synchronization phenomena in equilibrium and nonequilibrium systems to problems in mathematics, engineering sciences, biology, psychology, social sciences, finance, and economics.


---

## Body

## 1. Introduction

### 1.1. McKean and a New Class of Stochastic Processes

In two seminal papers, McKean introduced a new type of stochastic processes [1, 2]. The stochastic processes studied by McKean satisfy drift-diffusion equations for probability densities that are nonlinear with respect to their probability densities. While nonlinear drift-diffusion equations have frequently been used to describe how concentration fields evolve in time [3], the models considered by McKean describe stochastic processes.  A stochastic process not only describes how expectation values (such as the mean value) evolve in time, but also, provides a complete  description of all possible correlations between two or more than two time points [4–6]. Consequently, the study of stochastic processes goes beyond the study of concentration fields or the evolution of single-time point expectation values such as the mean. The observation by McKean that nonlinear drift-diffusion equations can describe stochastic processes opened a new avenue of research.


### 1.2. From McKean’s Stochastic Processes to Strongly Nonlinear Stochastic Processes

Due to the key role that the nonlinearity has for the processes considered by McKean and due to the fact that the processes satisfy the Markov property, processes of this type have been called “nonlinear Markov processes” (see, e.g., [7]). Moreover, in this context the term “nonlinear Fokker-Planck equations” has also been introduced because the drift-diffusion equations of nonlinear Markov processes exhibit formally the mathematical structure of Fokker-Planck equations. However, the phrases “nonlinear Markov processes” and “nonlinear Fokker-Planck equations” are misleading because they refer both to processes that are subjected to nonlinear forces and to processes that are described by drift-diffusion equations that are nonlinear with respect to their probability densities. In search for a more specific label, two alternative terms have been proposed in the literature. Wehner and Wolfer [8] suggested the term “truly nonlinear Fokker-Planck equations,” whereas Frank [9] coined the term “strongly nonlinear Fokker-Planck equations.” Only recently, it has become clear that the class of stochastic processes addressed by McKean is not limited to be described by drift-diffusion equations of the Fokker-Planck type. Therefore, the research field inspired by the work of McKean is about nonlinear stochastic processes in general. To give processes of this kind a more specific name, we may adopt the suggestions by Wehner and Wolfer or Frank and say that there is a novel class of processes called “truly nonlinear stochastic processes” or “strongly nonlinear stochastic processes.” In what follows, we will use the term “strongly” that should be considered as a synonym to “truly” just to avoid repeating the phrase “truly or strongly” over and over again.For the purpose of this paper, the following unifying definition of a strongly nonlinear process will be used:A strongly nonlinear stochastic process is a stochastic process for which the evolution of realizations depends on at least one expectation value of the process or on the appropriately defined probability density of that process.Note that for strongly nonlinear processes evolving on discrete state spaces the realizations are required to depend on at least one expectation value or on the occupation probability rather than the probability density. However, as we will demonstrate in Section2, for our purposes the occupation probability of a discrete-state process can be derived from an appropriately defined probability density. Therefore, the definition as given above holds both for processes evolving on continuous and discrete state spaces.In this review of strongly nonlinear stochastic processes a distinction will be made between four classes of processes. These classes arise naturally if it is taken into consideration that in many disciplines time is not necessarily considered to be a continuous variable. Rather, time might be considered as a discrete variable such that processes evolve in discrete time steps. For example, a process describing the evolution of a stock price on a day-to-day basis is typically modelled by means of a time-discrete process [10, 11]. In short, time can be a discrete or continuous variable. Likewise, the state space of a process might be discrete or continuous. The four possible combinations of time and space characteristics give rise to four classes of stochastic processes. We will address these classes in separate sections, see Table 1.Table 1
Four classes of strongly nonlinear stochastic processes and sections in which they will be addressed.
Time discrete
Time continuous
Space discrete
Section3
Section5
Space continuous
Section4
Section6Before we proceed in Sections3, 4, 5, and 6 with the systematic presentation outlined in Table 1, we will present in Section 2 benchmark examples of processes for each of the four classes. For the examples discussed in Section 2 a unifying mathematical definition irrespective of class membership will be presented in the final subsection of Section 2. Mean field theory and the mathematical literature will be addressed in Sections 7 and 8. As it turns out, all examples reviewed in Sections 2 to 6 exhibit the Markov property. Therefore, Section 9 is devoted to the review of some more recent studies on non-Markovian strongly nonlinear stochastic processes. In Section 10 some key issues of strongly nonlinear stochastic processes will be highlighted.


## 1.1. McKean and a New Class of Stochastic Processes

In two seminal papers, McKean introduced a new type of stochastic processes [1, 2]. The stochastic processes studied by McKean satisfy drift-diffusion equations for probability densities that are nonlinear with respect to their probability densities. While nonlinear drift-diffusion equations have frequently been used to describe how concentration fields evolve in time [3], the models considered by McKean describe stochastic processes.  A stochastic process not only describes how expectation values (such as the mean value) evolve in time, but also, provides a complete  description of all possible correlations between two or more than two time points [4–6]. Consequently, the study of stochastic processes goes beyond the study of concentration fields or the evolution of single-time point expectation values such as the mean. The observation by McKean that nonlinear drift-diffusion equations can describe stochastic processes opened a new avenue of research.


## 1.2. From McKean’s Stochastic Processes to Strongly Nonlinear Stochastic Processes

Due to the key role that the nonlinearity has for the processes considered by McKean and due to the fact that the processes satisfy the Markov property, processes of this type have been called “nonlinear Markov processes” (see, e.g., [7]). Moreover, in this context the term “nonlinear Fokker-Planck equations” has also been introduced because the drift-diffusion equations of nonlinear Markov processes exhibit formally the mathematical structure of Fokker-Planck equations. However, the phrases “nonlinear Markov processes” and “nonlinear Fokker-Planck equations” are misleading because they refer both to processes that are subjected to nonlinear forces and to processes that are described by drift-diffusion equations that are nonlinear with respect to their probability densities. In search for a more specific label, two alternative terms have been proposed in the literature. Wehner and Wolfer [8] suggested the term “truly nonlinear Fokker-Planck equations,” whereas Frank [9] coined the term “strongly nonlinear Fokker-Planck equations.” Only recently, it has become clear that the class of stochastic processes addressed by McKean is not limited to be described by drift-diffusion equations of the Fokker-Planck type. Therefore, the research field inspired by the work of McKean is about nonlinear stochastic processes in general. To give processes of this kind a more specific name, we may adopt the suggestions by Wehner and Wolfer or Frank and say that there is a novel class of processes called “truly nonlinear stochastic processes” or “strongly nonlinear stochastic processes.” In what follows, we will use the term “strongly” that should be considered as a synonym to “truly” just to avoid repeating the phrase “truly or strongly” over and over again.For the purpose of this paper, the following unifying definition of a strongly nonlinear process will be used:A strongly nonlinear stochastic process is a stochastic process for which the evolution of realizations depends on at least one expectation value of the process or on the appropriately defined probability density of that process.Note that for strongly nonlinear processes evolving on discrete state spaces the realizations are required to depend on at least one expectation value or on the occupation probability rather than the probability density. However, as we will demonstrate in Section2, for our purposes the occupation probability of a discrete-state process can be derived from an appropriately defined probability density. Therefore, the definition as given above holds both for processes evolving on continuous and discrete state spaces.In this review of strongly nonlinear stochastic processes a distinction will be made between four classes of processes. These classes arise naturally if it is taken into consideration that in many disciplines time is not necessarily considered to be a continuous variable. Rather, time might be considered as a discrete variable such that processes evolve in discrete time steps. For example, a process describing the evolution of a stock price on a day-to-day basis is typically modelled by means of a time-discrete process [10, 11]. In short, time can be a discrete or continuous variable. Likewise, the state space of a process might be discrete or continuous. The four possible combinations of time and space characteristics give rise to four classes of stochastic processes. We will address these classes in separate sections, see Table 1.Table 1
Four classes of strongly nonlinear stochastic processes and sections in which they will be addressed.
Time discrete
Time continuous
Space discrete
Section3
Section5
Space continuous
Section4
Section6Before we proceed in Sections3, 4, 5, and 6 with the systematic presentation outlined in Table 1, we will present in Section 2 benchmark examples of processes for each of the four classes. For the examples discussed in Section 2 a unifying mathematical definition irrespective of class membership will be presented in the final subsection of Section 2. Mean field theory and the mathematical literature will be addressed in Sections 7 and 8. As it turns out, all examples reviewed in Sections 2 to 6 exhibit the Markov property. Therefore, Section 9 is devoted to the review of some more recent studies on non-Markovian strongly nonlinear stochastic processes. In Section 10 some key issues of strongly nonlinear stochastic processes will be highlighted.


## 2. Benchmark Examples of Strongly Nonlinear Processes and on a General Stochastic Evolution Equation for a Broad Class of Strongly Nonlinear Markov Processes

### 2.1. Some Notations

At the present stage, it is useful to introduce some notations in order to give some benchmark examples for each of the four classes indicated in Table1.Time-discrete processes evolve on discrete time stepst=1,2,3,…. In contrast, time-continuous processes are described by a continuous time variable t defined on some interval. For our purposes this interval will go from zero to infinity, meaning that we will consider initial value problems starting at an initial time t=0 and extending in time infinitely. Note that the symbol t will be used to denote both discrete-valued and continuous time variables.Processes evolving in discrete, finite state spaces exhibit a number ofN states labeled by k=1,2,…,N. Let X denote the state variable of a state discrete process, which implies that X assumes an integer between 1 and N. Moreover, let Xr denote the rth realization of the process. A fundamental (although not complete) stochastic description of processes evolving on finite, discrete state spaces is given by the occupation probabilities p(k). The occupation probability p(k) defines the probability that the state k is occupied. It can be computed from a number of R realizations in the limiting case of R→∞ like
(1)p(k)=limR→∞1Rδ(Xr,k).
The function δ(·) is the Kronecker function, which equals 1 for Xr=k and equals 0 otherwise. The probabilities p(k) defined by (1) satisfy the normalization condition
(2)∑k=1Np(k)=1.Processes evolving on finite dimensional, continuous state spaces can be described by a state vectorX with coefficients X1,…,XN. The elements or coordinates Xk are defined on certain domains of definition Ωk. For example, a domain of definition may correspond to the real line. An important measure characterizing stochastic processes evolving on continuous state spaces is the probability density P defined as the ensemble average over the Dirac delta function δ(·). Note that the same symbol δ(·) will be used for the Kronecker function and the Dirac delta function in order to stress the analogies between the classes listed in Table 1. Mathematically speaking, the probability density P can be defined by
(3)P(x)=limR→∞1Rδ(X1r-x1)⋯δ(XNr-xN).
In (3) the variable x is the state vector with coordinates x1,…,xN, whereas the elements  X1r,…,XNr are the elements of the rth vector-valued realization Xr of the process. The probability density satisfies the normalization condition
(4)∫ΩtotP(x)∏k=1Ndxk=1,Ωtot=∏k=1NΩk,
where Ωk are the domains of definitions mentioned previously.Let us return to processes defined on a discrete state space. In addition, we putN=1 and assume that Ωtot=Ω1 denotes the real half-line from 0 to infinity. Formally, the probability density P(x) assumes the form
(5)P(x)=∑k=1Np(k)δ(x-k),
where p(k) are the occupation probabilities of the states k defined in (1). In turn, the occupation probabilities can be calculated from P(x) like
(6)p(k)=∫k-εk+εP(x)dx,
where ε is a positive number smaller than 1 (e.g., ε=0.1). In view of (6) the occupation probabilities p(k) may be considered as functions of an appropriately defined probability density P, as stated in Section 1. With these notations at hand, some benchmark examples of strongly nonlinear stochastic processes can be given.


### 2.2. Strongly Nonlinear Markov Processes Described by Nonlinear Markov Chains

In general, the occupation probabilitiesp(k) of time-discrete stochastic processes evolving on discrete state spaces constitute a time-dependent vector p(t) with coefficients p(1,t),…,p(N,t). For the important special case of nonlinear Markov processes described by nonlinear Markov chains, the occupation probabilities satisfy [12]
(7)p(k,t+1)=∑i=1NT(i⟶k,t,p(t))p(i,t).
Here, T describes the transition probability matrix of the process. The matrix is an N times N matrix with coefficients T(i→k). The coefficients T(i→k) describe the conditional probabilities that transitions from state i to state k occur given that the process is in the state i. The conditional probabilities T(i→k) may depend explicitly on time t such that T(i→k,t). A nonlinear Markov process is characterized by the property that at least one conditional probability depends at least on one expectation value computed from the process at time step t or on one occupation probability p(k,t). In this case we have T(i→k,t,p(t)) as indicated in (7). A complete description of the stochastic process can be derived from (7) as shown in Frank [12]. In view of the definition of a strongly nonlinear stochastic process given in Section 1.2 that is based on realizations of stochastic processes, it is worth while to consider an alternative definition of the process. To this end, the stochastic map f(·,ε) can be defined that maps a state X(t) to a state X(t+1) like
(8)X(t+1)=f(X(t),t,p(t),ε(t)).
The mapping f(·,ε) depends on a random variable ε(t) that is uniformly distributed on the interval [0,1] at any time step t. Across time, ε(t) is statistically independent. The stochastic map f(·,ε) describes the jumps from X(t) to X(t+1) satisfying the transition probabilities T(i→k). In order to derive an explicit definition of f(·,ε), we first note that the conditional probabilities T(i→k) are by definition real numbers between 0 and 1 that add up to unity. Next, it us useful to consider the cumulative conditional probabilities S(i,k) defined for i=1,2,…N and k=0,1,…,N by
(9)S(i,k)=∑j=1kT(i⟶j,t,p(t)),
for k>0. For k=0 we put S(i,0)=0. For k>0 the cumulative measures S(i,k) describe the probabilities that transitions from state i to one of the states 1,2,…,k occur at time t given that the process has an occupation probability vector p at time t. The cumulative measures represent a set of monotonically increasing points on the interval [0,1] such that the interval between two adjacent points S(i,k-1) and S(i,k) equals T(i→k). Consequently, the point spectrum S(i,k) in combination with the random variable ε(t) can be used to construct realizations of the process under consideration. If X(t)=i holds and if ε(t) falls into the interval [S(i,k-1),S(i,k)), then we put X(t+1)=k. If we do so, then transitions from i to k happen with the probability T(i→k) if the process is in state i at t and distributed like p(t) at t—as required. A concise mathematical formulation of this idea can be obtained with the help of the indicator function I(i→k,ε) defined by
(10)I(i⟶k,t,p,ε)={1forε∈[S(i,k-1),S(i,k))0otherwise.
Then, the map f(·,ε) can be defined as follows:
(11)f(X(t)=i,t,p(t),ε(t))=∑k=1Nk·I(i⟶k,t,p(t),ε(t)).
As indicated in (11), the mapping depends on the probability vector p(t) because the indicator function I(·) depends on p(t) via the conditional probabilities S and T. Therefore, each realization Xr(t) of the process defined by (8)–(11) depends on the occupation probabilities p(k,t) of the process. The occupation probabilities in turn can be computed from the realizations; see (1). Therefore, the mathematical model given by (1), and (8)–(11) provides a closed (and complete) description for the strongly nonlinear stochastic process described by the nonlinear Markov chain (7).


### 2.3. Strongly Nonlinear Markov Processes Described by Stochastic Iterative Maps

The realization of a time-discrete stochastic process on anN-dimensional, continuous state space is described by a time-dependent state vector X(t) with time-dependent components X1(t),…,XN(t). From (3) it follows that at any time step t, the vector X(t) is distributed according to the time-dependent probability density:
(12)P(x,t)=limR→∞1Rδ(X1r(t)-x1)⋯δ(XNr(t)-xN).
Strongly nonlinear Markov processes defined by stochastic iterative maps with nonparametric white noise satisfy equations of the general form
(13)X(t+1)=h(X(t),t,θX(t)[P(x,t)])+g(X(t),t,θX(t)[P(x,t)])·ε(t).
A simplified version of (13) was introduced earlier in Frank [13]. In (13) the vector-valued variable h(·) is called the drift function and describes the deterministic evolution of the process. Second, the term g(·)·ε(t) in (13) describes the vector-valued fluctuating force of the stochastic process. The random vector ε(t) with coefficients ε1(t),…,εN(t) is an N-dimensional normally distributed white noise process. That is, each coefficient εk(t) is normally distributed at any time step t and statistically independent in time. The random coefficients εk(t) are statistically independent among each other. The variable g is an N times N matrix with coefficients gjk. The coefficients gjk with k=1,…,N for a fixed index j weight the contributions that the random components ε1(t),…,εN(t) have on the evolution of the process in the direction (or along the coordinate) xj. In particular, for gjk(·,t)=0 the random coefficient εk(t) does not occur in the evolution equation of Xj at time t. Therefore, g may be considered as weight matrix. Importantly, the coefficients gjk measure the strength of the fluctuating force g(·)·ε(t) of the stochastic process described by (13). The reason for this is that the random vector ε(t) is normalized at any time step t. Most importantly, in the case of a strongly nonlinear Markov process, the drift function h(·) or the fluctuating force g(·)·ε(t) or both depend on at least one expectation value or on the probability density P as indicated in (13). The operator θ maps the probability density P to a real number or a set of real numbers. This mapping may depend on the value of X. For example, the probability density may be evaluated at the value of the process like
(14)θX(t)[P(x,t)]=P(x,t)|x=X(t).
In this context, it should be noted that the formal dependencies h(θX[P]) and g(θX[P]) include the dependencies on expectation values (such as the mean) as special cases. For example,
(15)θX(t)[P(x,t)]=∫ΩtotxP(x,t)∏k=1Ndxk=〈X(t)〉,
where 〈·〉 is the operator symbol for an ensemble average such that 〈X〉 denotes the mean of the random vector X. Note that (12) and (13) provide a closed description of a stochastic process. Moreover, from the structure of (13) it follows immediately that processes defined by (12) and (13) satisfy the definition of strongly nonlinear stochastic processes presented in Section 1.2. An illustrative example of a stochastic iterative map (13) is the autoregressive model of order 1 subjected to a so-called mean field force hMF proportional to the mean 〈X(t)〉 of the process. For N=1, the model reads [13]
(16)X(t+1)=AX(t)+B〈X(t)〉+gε(t),〈X(t)〉=limR→∞1R∑r=1RXr(t),
with g≥0, where the parameters A and B are scalars. We will return to model (16) in Section 4.


### 2.4. Strongly Nonlinear Markov Processes Described by Nonlinear Master Equations

For time-continuous stochastic processes evolving on discrete state spaces the occupation probabilitiesp(k) are functions of time t, where t is a continuous variable. Likewise, the probability vector p(t) becomes a time-dependent variable with coefficients p(1,t),…,p(N,t) just as in Section 2.2 but with t being a continuous rather than a discrete variable. Time-continuous, strongly nonlinear Markov processes on discrete state spaces can be described by means of nonlinear master equations of the form (see, e.g., [14])
(17)p(k,t)=p(k,t′)+∫t′t[∑i=1Nw(i⟶k,s,p(s))·p(i,s)-p(k,s)∑i=1Nw(k⟶i,s,p(s))]ds,
for t>t′ and k=1,…,N with w(k→k)=0 for all states k. Here, the variables w(i→k)≥0 are transition rates that describe the rate of transitions (number of transitions per unit time) from states i to states k. We will use the convention w(k→k)=0 for all states k. Note that as opposed to the conditional probabilities T(·) occurring in a Markov chain, the rates w(·) of the master equation are not normalized. Transition rates may depend explicitly on time such that w(i→k,t). For strongly nonlinear stochastic processes transition rates do not only depend on the states i and k and on time but also depend on an expectation value or the occupation probabilities p(k,t). Note that (17) is the integral form of the master equation. If the integral is evaluated for a small time interval Δt=t-t′→0, then the integral version exhibits the structure of a Markov chain (7) with
(18)T(i⟶k,t,Δt,p(t))=Δt·w(i⟶k,t,p(t)),i≠k,T(k⟶k,t,Δt,p(t))=1-Δt∑j=1Nw(k⟶j,t,p(t)),
which can be written in a more concise form like
(19)T(i⟶k,t,Δt,p(t))=Δtw(i⟶k,t,p(t))+δ(i,k)·[1-Δt∑j=1Nw(i⟶j,t,p(t))].
Note again that we will use the convention that the diagonal elements w(k→k) of the transition rates equal zero, which implies that for the diagonal elements T(k→k) only the second term in (19) matters as stated in (18). For the sake of completeness, note that the differential version of the master equation (17) reads
(20)ddtp(k,t)=∑i=1Nw(i⟶k,t,p(t))p(i,t)-p(k,t)∑i=1Nw(k⟶i,t,p(t)).
Just as in Section 2.2, in order to describe the four classes of stochastic processes listed in Table 1 on an equal footing, we consider an alternative description of the process defined by (20). Accordingly, we introduce the flow function Φ(t,ε) that maps initial states X(0) defined on the integer set 1,…,N to states X(t):
(21)X(t)=Φ(X(0),ε).
The variable ε is a random variable that is uniformly distributed on the interval [0,1] at every time point t and is statistically independently distributed across time. In view of the fact that the time-discrete version of the master equation exhibits the structure of a Markov chain, trajectories defined by the flow function can be regarded as time-continuous counterparts to the trajectories of Markov chains as defined by (1), and (8)–(11). This analogy can be exploited to construct iteratively the flow function Φ in the limiting case Δt→0 like
(22)Φ:X(t+Δt)=f(X(t),t,Δt,p(t),ε(t)),
with
(23)f(X(t)=i,t,Δt,p(t),ε(t))=∑k=1Nk·I(i⟶k,t,Δt,p(t),ε(t)).
The iterative map f(·) involves the indicator function
(24)I(i⟶k,ε,t,Δt,p)={1forε∈[S(i,k-1),S(i,k))0otherwise
and the cumulative transition probabilities
(25)S(i,k)=∑j=1kT(i⟶j,t,Δt,p(t)).
Equations (1), (18), and (22)–(25) provide a closed approximate description of the trajectories of interest. By definition, the description becomes exact in the limiting case of infinitely small time intervals Δt. Note that the conditional (transition) probabilities T(·) defined by (18) are in the interval [0,1] provided that Δt is chosen sufficiently small. In the limiting case Δt→0 this is always the case assuming that the rates w(·) are bounded from above. Note also that the probabilities T(·) defined by (18) are normalized to 1. Furthermore, from (18), and (22)–(25) it follows that f(·)=i holds for Δt=0. That is, the iterative map (22) describes transitions from states X(t)=i to states X(t+Δt) in the small time interval Δt. Most importantly, according to (22), the evolution of realizations depends on the occupation probability vector p(t) at every time point t. Therefore, processes defined by (1), (18), and (22)–(25) belong to the class of strongly nonlinear stochastic processes (see the definition in Section 1.2).


### 2.5. Strongly Nonlinear Markov Diffusion Processes Defined by Strongly Nonlinear Langevin Equations and Strongly Nonlinear Fokker-Planck Equations

Realizations of time-continuous stochastic processes evolving onN-dimensional, continuous state spaces are described by time-dependent state vectors X(t) with time-dependent components X1(t),…,XN(t). At every time point t, the vector X(t) is distributed according to the probability density P(x,t) defined by (12), where the variable x is the state vector with coordinates x1,…,xN. Strongly nonlinear Markov diffusion processes have realizations (or trajectories) defined in the limiting case of infinitely small time intervals Δt→0 by strongly nonlinear stochastic difference equations [9] of the form
(26)X(t+Δt)=X(t)+Δth(X(t),t,θX(t)[P(x,t)])+2Δtg(X(t),t,θX(t)[P(x,t)])·ε(t).
Carrying out the limiting procedure, (26) can be written in a more concise way in terms of a so-called Ito-Langevin equation:
(27)ddtX(t)=h(X(t),t,θX(t)[P(x,t)])+g(X(t),t,θX(t)[P(x,t)])·Γ(t).
Just as in the case of the stochastic iterative maps discussed in Section 2.3, (26) and (27) exhibit two different components. First, the vector-valued function h(·) occurring in (26) and (27) is the drift function introduced in Section 2.3 that accounts for the deterministic evolution of a process. Second, the expressions g·ε and g·Γ, respectively, describe a vector-valued fluctuating force. The variable ε(t) is the N-dimensional random variable introduced in Section 2.3. The function Γ(t) in the Ito-Langevin equation (27) is the counterpart to the random vector ε(t) in (26) and corresponds to a so-called vector-valued Langevin force [4, 6]. More precisely, each component Γk(t) is a Langevin force. We will consider Langevin forces normalized with respect to a magnitude of 2, like
(28)limR→∞1RΓjr(t)Γkr(s)=2δ(j,k)δ(t-s),
where Γkr(t) are realizations of the kth component. The variable g is the N times N weight matrix introduced in Section 2.3.The stochastic process defined by the Ito-Langevin equation (27) can equivalently be expressed by means of a strongly nonlinear Fokker-Planck equation:
(29)∂∂tP(x,t)=-∑k=1N∂∂xk[hk(x,t,θx[P(·,t)])P(x,t)]+∑i=1,j=1,k=1N∂2∂xi∂xj[gik(·)gjk(·)P(x,t)].
The evolution equation for P(x,t) is nonlinear for P(x,t). In contrast, the corresponding evolution equation for the conditional probability density T(x,t∣y,s) with t>s is linear with respect to T(x,t∣y,s) and reads
(30)∂∂tT(x,t∣y,s)=-∑k=1N∂∂xk[hk(x,t,θx[P(·,t)])T(x,t∣y,s)]+∑i=1,j=1,k=1N∂2∂xi∂xj[gik(·)gjk(·)T(x,t∣y,s)].
For details, see Frank [9]. A complete description of the stochastic process can be obtained either from the realizations defined by (26) and (27) or from the evolution equations (29) and (30) for P(x,t) and T(x,t∣y,s). The expression
(31)Dij(x,t,θx[P(·,t)])=∑k=1Ngik(·)gjk(·)
occurring in (29) and (30) is the diffusion coefficient of the process.A fundamental example of a strongly nonlinear Markov diffusion process is a generalized Ornstein-Uhlenbeck process involving a mean field forcehMF proportional to the mean 〈X(t)〉 of the process [9, 15, 16]. For N=1, the stochastic difference equation reads
(32)X(t+Δt)=X(t)-Δt[AX(t)+B〈X(t)〉]+2Δtgε(t),〈X(t)〉=limR→∞1R∑r=1RXr(t),
with g≥0. Typically, the case A, B>0 is considered. For (32) the strongly nonlinear Langevin equation reads
(33)ddtX(t)=-[AX(t)+B〈X(t)〉]+gΓ(t)〈X(t)〉=limR→∞1R∑r=1RXr(t).
Finally, the corresponding strongly nonlinear Fokker-Planck equation reads
(34)∂∂tP(x,t)=∂∂x[(Ax+BM(t))P(x,t)]+D∂2∂x2P(x,t),∂∂tT(x,t∣y,s)=∂∂x[(Ax+BM(t))T(x,t∣y,s)]+D∂2∂x2T(x,t∣y,s),M(t)=∫-∞∞xP(x,t)dx,D=g2.
The generalized Ornstein-Uhlenbeck process in its form as stochastic difference equation (32) corresponds to a special case of the generalized autoregressive model of order 1 defined by (16). Model (34) has been studied in the context of the Shimizu-Yamada model. We will return to this issue later.


### 2.6. On a General Stochastic Evolution Equation for Strongly Nonlinear Markov Processes

In the previous sections we considered some benchmark examples of strongly nonlinear stochastic processes belonging to the four classes of processes defined in Table1. All processes considered in Sections 2.2 to 2.5 satisfy the Markov property. A more general, formal description of strongly nonlinear Markov processes that holds for all four classes reads
(35)X(t+ξ)=f(X(t),t,ξ,θX(t)[P(x,t)],ε(t)).
The additional variable ξ defines the characteristic jumping interval of the process under consideration and can be used to distinguish between time-discrete and time-continuous processes. For ξ=1 we are dealing with time-discrete stochastic processes. In the case of processes defined on discrete state spaces P(x,t) is related to the occupation probability p(k,t) by means of (5) and (6) and (35) reduces to the special case of (8) of the Markov chain model. For processes defined on continuous state spaces and ξ=1 the general form (35) becomes (13) in the case of stochastic processes driven by nonparametric white noise. In the limiting case ξ→0 we are dealing with time-continuous stochastic processes. In this case, the function f(·) possesses the property f(·)=X(t) for ξ=0. For processes that evolve on discrete state spaces and satisfy a nonlinear master equation it follows that (35) becomes the iterative map (22) provided that the argument P(x,t) in f(·) is expressed in terms of the occupation probability p(k,t); see (6). In contrast, for processes evolving on continuous state spaces and satisfying strongly nonlinear Ito-Langevin equations, the general form (35) reduces to the stochastic difference equation (26) when we consider the limiting case ξ→0.


## 2.1. Some Notations

At the present stage, it is useful to introduce some notations in order to give some benchmark examples for each of the four classes indicated in Table1.Time-discrete processes evolve on discrete time stepst=1,2,3,…. In contrast, time-continuous processes are described by a continuous time variable t defined on some interval. For our purposes this interval will go from zero to infinity, meaning that we will consider initial value problems starting at an initial time t=0 and extending in time infinitely. Note that the symbol t will be used to denote both discrete-valued and continuous time variables.Processes evolving in discrete, finite state spaces exhibit a number ofN states labeled by k=1,2,…,N. Let X denote the state variable of a state discrete process, which implies that X assumes an integer between 1 and N. Moreover, let Xr denote the rth realization of the process. A fundamental (although not complete) stochastic description of processes evolving on finite, discrete state spaces is given by the occupation probabilities p(k). The occupation probability p(k) defines the probability that the state k is occupied. It can be computed from a number of R realizations in the limiting case of R→∞ like
(1)p(k)=limR→∞1Rδ(Xr,k).
The function δ(·) is the Kronecker function, which equals 1 for Xr=k and equals 0 otherwise. The probabilities p(k) defined by (1) satisfy the normalization condition
(2)∑k=1Np(k)=1.Processes evolving on finite dimensional, continuous state spaces can be described by a state vectorX with coefficients X1,…,XN. The elements or coordinates Xk are defined on certain domains of definition Ωk. For example, a domain of definition may correspond to the real line. An important measure characterizing stochastic processes evolving on continuous state spaces is the probability density P defined as the ensemble average over the Dirac delta function δ(·). Note that the same symbol δ(·) will be used for the Kronecker function and the Dirac delta function in order to stress the analogies between the classes listed in Table 1. Mathematically speaking, the probability density P can be defined by
(3)P(x)=limR→∞1Rδ(X1r-x1)⋯δ(XNr-xN).
In (3) the variable x is the state vector with coordinates x1,…,xN, whereas the elements  X1r,…,XNr are the elements of the rth vector-valued realization Xr of the process. The probability density satisfies the normalization condition
(4)∫ΩtotP(x)∏k=1Ndxk=1,Ωtot=∏k=1NΩk,
where Ωk are the domains of definitions mentioned previously.Let us return to processes defined on a discrete state space. In addition, we putN=1 and assume that Ωtot=Ω1 denotes the real half-line from 0 to infinity. Formally, the probability density P(x) assumes the form
(5)P(x)=∑k=1Np(k)δ(x-k),
where p(k) are the occupation probabilities of the states k defined in (1). In turn, the occupation probabilities can be calculated from P(x) like
(6)p(k)=∫k-εk+εP(x)dx,
where ε is a positive number smaller than 1 (e.g., ε=0.1). In view of (6) the occupation probabilities p(k) may be considered as functions of an appropriately defined probability density P, as stated in Section 1. With these notations at hand, some benchmark examples of strongly nonlinear stochastic processes can be given.


## 2.2. Strongly Nonlinear Markov Processes Described by Nonlinear Markov Chains

In general, the occupation probabilitiesp(k) of time-discrete stochastic processes evolving on discrete state spaces constitute a time-dependent vector p(t) with coefficients p(1,t),…,p(N,t). For the important special case of nonlinear Markov processes described by nonlinear Markov chains, the occupation probabilities satisfy [12]
(7)p(k,t+1)=∑i=1NT(i⟶k,t,p(t))p(i,t).
Here, T describes the transition probability matrix of the process. The matrix is an N times N matrix with coefficients T(i→k). The coefficients T(i→k) describe the conditional probabilities that transitions from state i to state k occur given that the process is in the state i. The conditional probabilities T(i→k) may depend explicitly on time t such that T(i→k,t). A nonlinear Markov process is characterized by the property that at least one conditional probability depends at least on one expectation value computed from the process at time step t or on one occupation probability p(k,t). In this case we have T(i→k,t,p(t)) as indicated in (7). A complete description of the stochastic process can be derived from (7) as shown in Frank [12]. In view of the definition of a strongly nonlinear stochastic process given in Section 1.2 that is based on realizations of stochastic processes, it is worth while to consider an alternative definition of the process. To this end, the stochastic map f(·,ε) can be defined that maps a state X(t) to a state X(t+1) like
(8)X(t+1)=f(X(t),t,p(t),ε(t)).
The mapping f(·,ε) depends on a random variable ε(t) that is uniformly distributed on the interval [0,1] at any time step t. Across time, ε(t) is statistically independent. The stochastic map f(·,ε) describes the jumps from X(t) to X(t+1) satisfying the transition probabilities T(i→k). In order to derive an explicit definition of f(·,ε), we first note that the conditional probabilities T(i→k) are by definition real numbers between 0 and 1 that add up to unity. Next, it us useful to consider the cumulative conditional probabilities S(i,k) defined for i=1,2,…N and k=0,1,…,N by
(9)S(i,k)=∑j=1kT(i⟶j,t,p(t)),
for k>0. For k=0 we put S(i,0)=0. For k>0 the cumulative measures S(i,k) describe the probabilities that transitions from state i to one of the states 1,2,…,k occur at time t given that the process has an occupation probability vector p at time t. The cumulative measures represent a set of monotonically increasing points on the interval [0,1] such that the interval between two adjacent points S(i,k-1) and S(i,k) equals T(i→k). Consequently, the point spectrum S(i,k) in combination with the random variable ε(t) can be used to construct realizations of the process under consideration. If X(t)=i holds and if ε(t) falls into the interval [S(i,k-1),S(i,k)), then we put X(t+1)=k. If we do so, then transitions from i to k happen with the probability T(i→k) if the process is in state i at t and distributed like p(t) at t—as required. A concise mathematical formulation of this idea can be obtained with the help of the indicator function I(i→k,ε) defined by
(10)I(i⟶k,t,p,ε)={1forε∈[S(i,k-1),S(i,k))0otherwise.
Then, the map f(·,ε) can be defined as follows:
(11)f(X(t)=i,t,p(t),ε(t))=∑k=1Nk·I(i⟶k,t,p(t),ε(t)).
As indicated in (11), the mapping depends on the probability vector p(t) because the indicator function I(·) depends on p(t) via the conditional probabilities S and T. Therefore, each realization Xr(t) of the process defined by (8)–(11) depends on the occupation probabilities p(k,t) of the process. The occupation probabilities in turn can be computed from the realizations; see (1). Therefore, the mathematical model given by (1), and (8)–(11) provides a closed (and complete) description for the strongly nonlinear stochastic process described by the nonlinear Markov chain (7).


## 2.3. Strongly Nonlinear Markov Processes Described by Stochastic Iterative Maps

The realization of a time-discrete stochastic process on anN-dimensional, continuous state space is described by a time-dependent state vector X(t) with time-dependent components X1(t),…,XN(t). From (3) it follows that at any time step t, the vector X(t) is distributed according to the time-dependent probability density:
(12)P(x,t)=limR→∞1Rδ(X1r(t)-x1)⋯δ(XNr(t)-xN).
Strongly nonlinear Markov processes defined by stochastic iterative maps with nonparametric white noise satisfy equations of the general form
(13)X(t+1)=h(X(t),t,θX(t)[P(x,t)])+g(X(t),t,θX(t)[P(x,t)])·ε(t).
A simplified version of (13) was introduced earlier in Frank [13]. In (13) the vector-valued variable h(·) is called the drift function and describes the deterministic evolution of the process. Second, the term g(·)·ε(t) in (13) describes the vector-valued fluctuating force of the stochastic process. The random vector ε(t) with coefficients ε1(t),…,εN(t) is an N-dimensional normally distributed white noise process. That is, each coefficient εk(t) is normally distributed at any time step t and statistically independent in time. The random coefficients εk(t) are statistically independent among each other. The variable g is an N times N matrix with coefficients gjk. The coefficients gjk with k=1,…,N for a fixed index j weight the contributions that the random components ε1(t),…,εN(t) have on the evolution of the process in the direction (or along the coordinate) xj. In particular, for gjk(·,t)=0 the random coefficient εk(t) does not occur in the evolution equation of Xj at time t. Therefore, g may be considered as weight matrix. Importantly, the coefficients gjk measure the strength of the fluctuating force g(·)·ε(t) of the stochastic process described by (13). The reason for this is that the random vector ε(t) is normalized at any time step t. Most importantly, in the case of a strongly nonlinear Markov process, the drift function h(·) or the fluctuating force g(·)·ε(t) or both depend on at least one expectation value or on the probability density P as indicated in (13). The operator θ maps the probability density P to a real number or a set of real numbers. This mapping may depend on the value of X. For example, the probability density may be evaluated at the value of the process like
(14)θX(t)[P(x,t)]=P(x,t)|x=X(t).
In this context, it should be noted that the formal dependencies h(θX[P]) and g(θX[P]) include the dependencies on expectation values (such as the mean) as special cases. For example,
(15)θX(t)[P(x,t)]=∫ΩtotxP(x,t)∏k=1Ndxk=〈X(t)〉,
where 〈·〉 is the operator symbol for an ensemble average such that 〈X〉 denotes the mean of the random vector X. Note that (12) and (13) provide a closed description of a stochastic process. Moreover, from the structure of (13) it follows immediately that processes defined by (12) and (13) satisfy the definition of strongly nonlinear stochastic processes presented in Section 1.2. An illustrative example of a stochastic iterative map (13) is the autoregressive model of order 1 subjected to a so-called mean field force hMF proportional to the mean 〈X(t)〉 of the process. For N=1, the model reads [13]
(16)X(t+1)=AX(t)+B〈X(t)〉+gε(t),〈X(t)〉=limR→∞1R∑r=1RXr(t),
with g≥0, where the parameters A and B are scalars. We will return to model (16) in Section 4.


## 2.4. Strongly Nonlinear Markov Processes Described by Nonlinear Master Equations

For time-continuous stochastic processes evolving on discrete state spaces the occupation probabilitiesp(k) are functions of time t, where t is a continuous variable. Likewise, the probability vector p(t) becomes a time-dependent variable with coefficients p(1,t),…,p(N,t) just as in Section 2.2 but with t being a continuous rather than a discrete variable. Time-continuous, strongly nonlinear Markov processes on discrete state spaces can be described by means of nonlinear master equations of the form (see, e.g., [14])
(17)p(k,t)=p(k,t′)+∫t′t[∑i=1Nw(i⟶k,s,p(s))·p(i,s)-p(k,s)∑i=1Nw(k⟶i,s,p(s))]ds,
for t>t′ and k=1,…,N with w(k→k)=0 for all states k. Here, the variables w(i→k)≥0 are transition rates that describe the rate of transitions (number of transitions per unit time) from states i to states k. We will use the convention w(k→k)=0 for all states k. Note that as opposed to the conditional probabilities T(·) occurring in a Markov chain, the rates w(·) of the master equation are not normalized. Transition rates may depend explicitly on time such that w(i→k,t). For strongly nonlinear stochastic processes transition rates do not only depend on the states i and k and on time but also depend on an expectation value or the occupation probabilities p(k,t). Note that (17) is the integral form of the master equation. If the integral is evaluated for a small time interval Δt=t-t′→0, then the integral version exhibits the structure of a Markov chain (7) with
(18)T(i⟶k,t,Δt,p(t))=Δt·w(i⟶k,t,p(t)),i≠k,T(k⟶k,t,Δt,p(t))=1-Δt∑j=1Nw(k⟶j,t,p(t)),
which can be written in a more concise form like
(19)T(i⟶k,t,Δt,p(t))=Δtw(i⟶k,t,p(t))+δ(i,k)·[1-Δt∑j=1Nw(i⟶j,t,p(t))].
Note again that we will use the convention that the diagonal elements w(k→k) of the transition rates equal zero, which implies that for the diagonal elements T(k→k) only the second term in (19) matters as stated in (18). For the sake of completeness, note that the differential version of the master equation (17) reads
(20)ddtp(k,t)=∑i=1Nw(i⟶k,t,p(t))p(i,t)-p(k,t)∑i=1Nw(k⟶i,t,p(t)).
Just as in Section 2.2, in order to describe the four classes of stochastic processes listed in Table 1 on an equal footing, we consider an alternative description of the process defined by (20). Accordingly, we introduce the flow function Φ(t,ε) that maps initial states X(0) defined on the integer set 1,…,N to states X(t):
(21)X(t)=Φ(X(0),ε).
The variable ε is a random variable that is uniformly distributed on the interval [0,1] at every time point t and is statistically independently distributed across time. In view of the fact that the time-discrete version of the master equation exhibits the structure of a Markov chain, trajectories defined by the flow function can be regarded as time-continuous counterparts to the trajectories of Markov chains as defined by (1), and (8)–(11). This analogy can be exploited to construct iteratively the flow function Φ in the limiting case Δt→0 like
(22)Φ:X(t+Δt)=f(X(t),t,Δt,p(t),ε(t)),
with
(23)f(X(t)=i,t,Δt,p(t),ε(t))=∑k=1Nk·I(i⟶k,t,Δt,p(t),ε(t)).
The iterative map f(·) involves the indicator function
(24)I(i⟶k,ε,t,Δt,p)={1forε∈[S(i,k-1),S(i,k))0otherwise
and the cumulative transition probabilities
(25)S(i,k)=∑j=1kT(i⟶j,t,Δt,p(t)).
Equations (1), (18), and (22)–(25) provide a closed approximate description of the trajectories of interest. By definition, the description becomes exact in the limiting case of infinitely small time intervals Δt. Note that the conditional (transition) probabilities T(·) defined by (18) are in the interval [0,1] provided that Δt is chosen sufficiently small. In the limiting case Δt→0 this is always the case assuming that the rates w(·) are bounded from above. Note also that the probabilities T(·) defined by (18) are normalized to 1. Furthermore, from (18), and (22)–(25) it follows that f(·)=i holds for Δt=0. That is, the iterative map (22) describes transitions from states X(t)=i to states X(t+Δt) in the small time interval Δt. Most importantly, according to (22), the evolution of realizations depends on the occupation probability vector p(t) at every time point t. Therefore, processes defined by (1), (18), and (22)–(25) belong to the class of strongly nonlinear stochastic processes (see the definition in Section 1.2).


## 2.5. Strongly Nonlinear Markov Diffusion Processes Defined by Strongly Nonlinear Langevin Equations and Strongly Nonlinear Fokker-Planck Equations

Realizations of time-continuous stochastic processes evolving onN-dimensional, continuous state spaces are described by time-dependent state vectors X(t) with time-dependent components X1(t),…,XN(t). At every time point t, the vector X(t) is distributed according to the probability density P(x,t) defined by (12), where the variable x is the state vector with coordinates x1,…,xN. Strongly nonlinear Markov diffusion processes have realizations (or trajectories) defined in the limiting case of infinitely small time intervals Δt→0 by strongly nonlinear stochastic difference equations [9] of the form
(26)X(t+Δt)=X(t)+Δth(X(t),t,θX(t)[P(x,t)])+2Δtg(X(t),t,θX(t)[P(x,t)])·ε(t).
Carrying out the limiting procedure, (26) can be written in a more concise way in terms of a so-called Ito-Langevin equation:
(27)ddtX(t)=h(X(t),t,θX(t)[P(x,t)])+g(X(t),t,θX(t)[P(x,t)])·Γ(t).
Just as in the case of the stochastic iterative maps discussed in Section 2.3, (26) and (27) exhibit two different components. First, the vector-valued function h(·) occurring in (26) and (27) is the drift function introduced in Section 2.3 that accounts for the deterministic evolution of a process. Second, the expressions g·ε and g·Γ, respectively, describe a vector-valued fluctuating force. The variable ε(t) is the N-dimensional random variable introduced in Section 2.3. The function Γ(t) in the Ito-Langevin equation (27) is the counterpart to the random vector ε(t) in (26) and corresponds to a so-called vector-valued Langevin force [4, 6]. More precisely, each component Γk(t) is a Langevin force. We will consider Langevin forces normalized with respect to a magnitude of 2, like
(28)limR→∞1RΓjr(t)Γkr(s)=2δ(j,k)δ(t-s),
where Γkr(t) are realizations of the kth component. The variable g is the N times N weight matrix introduced in Section 2.3.The stochastic process defined by the Ito-Langevin equation (27) can equivalently be expressed by means of a strongly nonlinear Fokker-Planck equation:
(29)∂∂tP(x,t)=-∑k=1N∂∂xk[hk(x,t,θx[P(·,t)])P(x,t)]+∑i=1,j=1,k=1N∂2∂xi∂xj[gik(·)gjk(·)P(x,t)].
The evolution equation for P(x,t) is nonlinear for P(x,t). In contrast, the corresponding evolution equation for the conditional probability density T(x,t∣y,s) with t>s is linear with respect to T(x,t∣y,s) and reads
(30)∂∂tT(x,t∣y,s)=-∑k=1N∂∂xk[hk(x,t,θx[P(·,t)])T(x,t∣y,s)]+∑i=1,j=1,k=1N∂2∂xi∂xj[gik(·)gjk(·)T(x,t∣y,s)].
For details, see Frank [9]. A complete description of the stochastic process can be obtained either from the realizations defined by (26) and (27) or from the evolution equations (29) and (30) for P(x,t) and T(x,t∣y,s). The expression
(31)Dij(x,t,θx[P(·,t)])=∑k=1Ngik(·)gjk(·)
occurring in (29) and (30) is the diffusion coefficient of the process.A fundamental example of a strongly nonlinear Markov diffusion process is a generalized Ornstein-Uhlenbeck process involving a mean field forcehMF proportional to the mean 〈X(t)〉 of the process [9, 15, 16]. For N=1, the stochastic difference equation reads
(32)X(t+Δt)=X(t)-Δt[AX(t)+B〈X(t)〉]+2Δtgε(t),〈X(t)〉=limR→∞1R∑r=1RXr(t),
with g≥0. Typically, the case A, B>0 is considered. For (32) the strongly nonlinear Langevin equation reads
(33)ddtX(t)=-[AX(t)+B〈X(t)〉]+gΓ(t)〈X(t)〉=limR→∞1R∑r=1RXr(t).
Finally, the corresponding strongly nonlinear Fokker-Planck equation reads
(34)∂∂tP(x,t)=∂∂x[(Ax+BM(t))P(x,t)]+D∂2∂x2P(x,t),∂∂tT(x,t∣y,s)=∂∂x[(Ax+BM(t))T(x,t∣y,s)]+D∂2∂x2T(x,t∣y,s),M(t)=∫-∞∞xP(x,t)dx,D=g2.
The generalized Ornstein-Uhlenbeck process in its form as stochastic difference equation (32) corresponds to a special case of the generalized autoregressive model of order 1 defined by (16). Model (34) has been studied in the context of the Shimizu-Yamada model. We will return to this issue later.


## 2.6. On a General Stochastic Evolution Equation for Strongly Nonlinear Markov Processes

In the previous sections we considered some benchmark examples of strongly nonlinear stochastic processes belonging to the four classes of processes defined in Table1. All processes considered in Sections 2.2 to 2.5 satisfy the Markov property. A more general, formal description of strongly nonlinear Markov processes that holds for all four classes reads
(35)X(t+ξ)=f(X(t),t,ξ,θX(t)[P(x,t)],ε(t)).
The additional variable ξ defines the characteristic jumping interval of the process under consideration and can be used to distinguish between time-discrete and time-continuous processes. For ξ=1 we are dealing with time-discrete stochastic processes. In the case of processes defined on discrete state spaces P(x,t) is related to the occupation probability p(k,t) by means of (5) and (6) and (35) reduces to the special case of (8) of the Markov chain model. For processes defined on continuous state spaces and ξ=1 the general form (35) becomes (13) in the case of stochastic processes driven by nonparametric white noise. In the limiting case ξ→0 we are dealing with time-continuous stochastic processes. In this case, the function f(·) possesses the property f(·)=X(t) for ξ=0. For processes that evolve on discrete state spaces and satisfy a nonlinear master equation it follows that (35) becomes the iterative map (22) provided that the argument P(x,t) in f(·) is expressed in terms of the occupation probability p(k,t); see (6). In contrast, for processes evolving on continuous state spaces and satisfying strongly nonlinear Ito-Langevin equations, the general form (35) reduces to the stochastic difference equation (26) when we consider the limiting case ξ→0.


## 3. Time-Discrete Strongly Nonlinear Stochastic Processes on Discrete State Spaces

### 3.1. The Strongly Nonlinear Two-State Markov Chain Model

A fundamental strongly nonlinear stochastic process is the two-state nonlinear Markov chain process [17]. For the sake of convenience, the states k=1 and k=2 will be labeled with “A” and “B.” Transitions from A to A, A to B, B to A, and B to B occur with conditional probabilities T(A→A), T(A→B), T(B→A), and T(B→B), respectively. Accordingly, the nonlinear Markov chain equation (7) reads
(36)p(A,t+1)=T(A⟶A,t,p(t))p(A,t)+T(B⟶A,t,p(t))p(B,t)p(B,t+1)=T(A⟶B,t,p(t))p(A,t)+T(B⟶B,t,p(t))p(B,t).
Taking the normalization condition p(A)+p(B)=1 into consideration, (36) can equivalently be expressed by
(37)p(A,t+1)=[T(A⟶A,t,p(A,t))-T(B⟶A,t,p(A,t))]p(A,t)+T(B⟶A,t,p(A,t)),
which is a closed nonlinear equation in p(A,t). Model (37) involves two functions only: the conditional probabilities T(A→A) and T(B→A) which can be chosen independently from each other and in the case of a nonlinear Markov process depend on p(A,t) as indicated in (37). The remaining conditional probabilities can be calculated from T(A→B)=1-T(A→A) and T(B→B)=1-T(B→A). Model (37) has been examined in physics, social sciences, and psychology as will be discussed next.


### 3.2. Chaos in Probability: Time-Discrete Case

Following Frank [18], let us put T(A→A)=T(B→A)=T0. Then, (37) becomes
(38)p(A,t+1)=T0(t,p(A,t)).
Although (38) describes a stochastic process, formally (38) exhibits the structure of a nonlinear deterministic map subjected to the constraint that the function T0 must be in the interval [0,1] at any time step t and for any occupation probability p(A). Iterative maps in general are known to exhibit chaos [19]. Therefore, (38) suggests that time-discrete strongly nonlinear Markov processes exhibit chaos as well. Such strongly nonlinear chaotic processes exhibit chaos in the evolution of the occupation probabilities and the conditional probabilities. For example, in Frank [18] the logistic function T0(z)=α·z(1-z) [19, 20] was considered for which (38) reads
(39)p(A,t+1)=αp(A,t)[1-p(A,t)].
For α in [0,4] the function T0(·) is in the interval [0,1] such that (39) maps the occupation probability p(A) from the interval [0,1] to the interval [0,1] as required. For parameters α≥αcrit≈3.570 [20] the map exhibits chaotic solutions. That is, p(A,t) evolves in a chaotic fashion. However, as stated previously, not only the occupation probability exhibits chaos, but the conditional probabilities T(A→A)=T(B→A)=T0 exhibit a chaotic behavior as well which was demonstrated explicitly by numerical simulations [18]. Chaos in strongly nonlinear time-discrete Markov processes is not limited to the two-state models. In fact, the two-state model can be generalized to an N-state model. In this case, (38) becomes [18]
(40)p(k,t+1)=T0(t,p(t)),
for k=1,…,N-1. The occupation probability p(N,t) is computed from the normalization condition. Consequently, T0 depends only on the occupation probabilities p(1,t),…,p(N-1,t). For example, the two-dimensional chaotic Baker map [20] was described by means of (40) using N=3 [18].


### 3.3. Social Sciences and Group Decision Making: Time-Discrete Case

Group decision making has frequently been described by nonlinear stochastic processes  [21–24]. The reason for this is that under certain circumstances individuals are likely to be affected by the opinion that is presented by a group as a whole rather than by the opinion of individuals. This so-called group opinion in turn is produced by the opinions and decisions of the individual group members. Consequently, the decision making process of an individual is affected by the sample average obtained from the opinions of the others. If the groups can be considered as being relatively large, the sample average may be replaced by the ensemble average. The decision making behavior of an individual is then considered as a realization of a stochastic process that is affected by the ensemble average of the process. In the context of voting behavior, the two-state model introduced in Section 3.1 was applied [17]. The members of the US House of Representatives had the choice to vote for or against passing the so-called bailout bill (US House of Representatives, USA, 2008). The bill failed to pass in the September 29 vote. In a second vote, on October 3, the bill finally passed. That is, a critical number of “No” voters changed their opinions between September 29 and October 3. The states k=1 and k=2 were labeled by “Y” and “N” representing “Yes” and “No” voters. According to model (37), the conditional probabilities T(N→Y) and T(Y→Y) describing a change in the opinion from “No” to “Yes” and a consistent “Yes” opinion were considered. Using the notation of (37), the voting model proposed in Frank [18] reads
(41)p(Y,t+1)=[T(Y⟶Y,t,p(Y,t))-T(N⟶Y,t,p(Y,t))]p(Y,t)+T(N⟶Y,t,p(Y,t)).
The data available in the literature suggested that T(Y→Y)=1. For the T(N→Y) the function T(N→Y)=A-B·p(Y,t) was used to describe the impact of the group opinion on the voting behavior. This relationship assumes that the pressure to change from a “No” voter to a “Yes” voter was high immediately subsequent to the September 29 votes, when the probability p(Y,t) at t = Sept 29 was not high enough to make the bill pass. In this context, it should be mentioned that the failure of the bill to pass on September 29 triggered a dramatic breakdown of the US stock market. That is, the members of the House of Representatives who had voted with “No” were not only subjected to political pressure but also experienced essential pressure from the economy. For T(Y→Y)=1 and T(N→Y)=A-Bp(Y,t) the stationary occupation probability p(Y,st) is given by p(Y,st)=A/B such that (41) eventually reads
(42)p(Y,t+1)=p(Y,t)+B[p(Y,st)-p(Y,t)][1-p(Y,t)].
The stationary value p(Y,st) was estimated from the data. A detailed model-based analysis of the data showed that the members of the two main parties, the Democratic and Republican members, were differently affected by the group pressure.


### 3.4. Psychology and Human Memory

A pioneering experiment in human memory research is the free recall experiment in which participants are asked to recall items of a particular category [25, 26]. For example, they are asked to recall animal names. Typically, participants are instructed to recall as many items as possible and to do so as fast as possible and without repetition. The total number of recalled items is by definition a monotonically increasing function of time. Strikingly, free recall involves clusters in which item subcategories are recalled more or less as a whole. For example, when recalling animal names a participant may recall a number of insect names in succession. From a theoretical point of view, at any given time point during the free recall process the ability to recall items in clusters depends on the size of the pool of items available for recalling. In other words, if there is a large pool of items that have not yet been recalled, then there is a high chance that a whole item subcategory can be recalled. In line with this argument, a stochastic model for free recall dynamics was developed in which the recall process depends on the size of the pool of items available for recall [27]. More precisely, time was parceled in discrete time intervals in which recall attempts are executed. The aim was to develop a model for the probability that a single item has been or has not yet been recalled at a given time step t. To this end, the two-state model of Section 3.1 was applied. Accordingly, the states k=1 and k=2 correspond to being recalled or being not yet been recalled and are labeled by “R” and “N,” respectively. The two relevant conditional probabilities are T(R→R) and T(N→R). By the design of the experiment, we have T(R→R)=1. Moreover, as argued previously, if the pool of items that have not yet recalled is large, that is, if the probability p(N) that an item has not yet been recalled is high, then the conditional probability T(N→R) should be high as well. Likewise, if the probability p(R) that an item has been recalled is high, then T(N→R) should be relatively low. Therefore, it was suggested to put T(N→R)=B·(1-p(R,t)) with B>0. This function implies that the occupation probability p(R,t) converges to the stationary fixed p(R,st)=1 at which T(N→R)=0. The full model can be obtained by substituting T(R→R) and T(N→R) into (37) and reads
(43)p(R,t+1)=p(R,t)+B[1-p(R,t)]2.
Note that from a mathematical point of view, this model is a special case of model (42) for voting behavior (hint: put p(Y,st)=1 in (42)). Since the model describes the single item recall probability, the model must be scaled to the number of items recalled by any given participant. To this end, a parameter C describing the maximal number of items for recall was introduced. The unknown parameters B and C were estimated for free recall data collected in a study by Rhodes and Turvey [28], and the sample mean values B=0.009 (SD = 0.005) and C=307 (SD = 138) were found for a sample of eight participants [27]. Note that (43) predicts that the fixed point p(R,st)=1 will never be reached in a finite number of time steps. Consequently, the model predicts that the total number of recalled items M at the end of the experiment is smaller than the number of available items for recall. This prediction was confirmed by the data. Moreover, the model-based analysis revealed a statistically significant positive correlation between M and C for the participant sample studied in Frank and Rhodes [27].


### 3.5. Biology and Evolutionary Population Dynamics

A primary concern of evolutionary population dynamics is the understanding of the origin of life. Among the many models that have been proposed in this context, the studies by Crutchfield and Görnerup [29] and Görnerup and Crutchfield [30] are of interest for our review. In these studies a number of N macromolecules with different functional properties were considered. These macromolecules compete with each other for survival during evolutionary time scales. However, not only competition but also cooperative behavior is taken into account. Roughly speaking, it is assumed that in the early stages of the development of life, macromolecules of different types existed and could change their types due to interactions with other macromolecules of the same or different types. Eventually, only a subset of functionally different types survived the selection process. Let k=1,…,N denote the functionally different types of macromolecules and p(k,t) the occupation probability of the kth type. Then, at every discrete reaction time step t the macromolecules can change their types such that transitions from type i to type k occur. This so-called metamodel can be formulated in terms of a nonlinear Markov chain model [31]. Accordingly, the conditional (transition) probability T(i→k) describes the probability of a transition from type i to type k given that a macromolecule is of type i at time step t. As suggested by the studies of Crutchfield and Görnerup [29] and Görnerup and Crutchfield [30] interactions between two different types of macromolecules can be modeled by assuming that the evolution of occupation probabilities is affected by terms of the form p(j,t)p(i,t). More precisely, in general, a macromolecule of type i may interact with a macromolecule of type j and turn into a macromolecule of type k. It can be shown that from the metamodel it follows that T(i→k) assumes the form [31]
(44)T(i⟶k)=∑j=1Ng(ijk)p(j,t).
In (44) the parameters g(ijk)≥0 are weight coefficients that describe the relevance of interactions between macromolecules of different types. Note that the conditional probabilities T(i→k) satisfy the normalization condition that the sum over all probabilities with index k equals 1. This normalization can be guaranteed by requiring that the sum over all weight coefficients g with index k equals 1. Substituting (44) into (7), the metamodel for evolutionary population dynamics can be cast into a nonlinear Markov chain model and reads
(45)p(k,t+1)=∑i,j=1Ng(ijk)p(j,t)p(i,t).
It has been shown that the nonlinear Markov chain model (45) is a useful tool for investigating evolutionary population dynamics. First of all, trajectories that describe the evolution of individual macromolecules were computed numerically using (1), and (8)–(11); see Section 2.2. That is, the model allows for generating temporal sequences of type changes for individual macromolecules. Second, it was shown that the degree to which a type k is attractive can be quantified. In particular, for evolutionary selection process that become stationary, the attractors of the surviving macromolecule types can be analyzed and distinctions between weakly and strongly attractive types can be made (e.g., see Figure 3 in [31]).


## 3.1. The Strongly Nonlinear Two-State Markov Chain Model

A fundamental strongly nonlinear stochastic process is the two-state nonlinear Markov chain process [17]. For the sake of convenience, the states k=1 and k=2 will be labeled with “A” and “B.” Transitions from A to A, A to B, B to A, and B to B occur with conditional probabilities T(A→A), T(A→B), T(B→A), and T(B→B), respectively. Accordingly, the nonlinear Markov chain equation (7) reads
(36)p(A,t+1)=T(A⟶A,t,p(t))p(A,t)+T(B⟶A,t,p(t))p(B,t)p(B,t+1)=T(A⟶B,t,p(t))p(A,t)+T(B⟶B,t,p(t))p(B,t).
Taking the normalization condition p(A)+p(B)=1 into consideration, (36) can equivalently be expressed by
(37)p(A,t+1)=[T(A⟶A,t,p(A,t))-T(B⟶A,t,p(A,t))]p(A,t)+T(B⟶A,t,p(A,t)),
which is a closed nonlinear equation in p(A,t). Model (37) involves two functions only: the conditional probabilities T(A→A) and T(B→A) which can be chosen independently from each other and in the case of a nonlinear Markov process depend on p(A,t) as indicated in (37). The remaining conditional probabilities can be calculated from T(A→B)=1-T(A→A) and T(B→B)=1-T(B→A). Model (37) has been examined in physics, social sciences, and psychology as will be discussed next.


## 3.2. Chaos in Probability: Time-Discrete Case

Following Frank [18], let us put T(A→A)=T(B→A)=T0. Then, (37) becomes
(38)p(A,t+1)=T0(t,p(A,t)).
Although (38) describes a stochastic process, formally (38) exhibits the structure of a nonlinear deterministic map subjected to the constraint that the function T0 must be in the interval [0,1] at any time step t and for any occupation probability p(A). Iterative maps in general are known to exhibit chaos [19]. Therefore, (38) suggests that time-discrete strongly nonlinear Markov processes exhibit chaos as well. Such strongly nonlinear chaotic processes exhibit chaos in the evolution of the occupation probabilities and the conditional probabilities. For example, in Frank [18] the logistic function T0(z)=α·z(1-z) [19, 20] was considered for which (38) reads
(39)p(A,t+1)=αp(A,t)[1-p(A,t)].
For α in [0,4] the function T0(·) is in the interval [0,1] such that (39) maps the occupation probability p(A) from the interval [0,1] to the interval [0,1] as required. For parameters α≥αcrit≈3.570 [20] the map exhibits chaotic solutions. That is, p(A,t) evolves in a chaotic fashion. However, as stated previously, not only the occupation probability exhibits chaos, but the conditional probabilities T(A→A)=T(B→A)=T0 exhibit a chaotic behavior as well which was demonstrated explicitly by numerical simulations [18]. Chaos in strongly nonlinear time-discrete Markov processes is not limited to the two-state models. In fact, the two-state model can be generalized to an N-state model. In this case, (38) becomes [18]
(40)p(k,t+1)=T0(t,p(t)),
for k=1,…,N-1. The occupation probability p(N,t) is computed from the normalization condition. Consequently, T0 depends only on the occupation probabilities p(1,t),…,p(N-1,t). For example, the two-dimensional chaotic Baker map [20] was described by means of (40) using N=3 [18].


## 3.3. Social Sciences and Group Decision Making: Time-Discrete Case

Group decision making has frequently been described by nonlinear stochastic processes  [21–24]. The reason for this is that under certain circumstances individuals are likely to be affected by the opinion that is presented by a group as a whole rather than by the opinion of individuals. This so-called group opinion in turn is produced by the opinions and decisions of the individual group members. Consequently, the decision making process of an individual is affected by the sample average obtained from the opinions of the others. If the groups can be considered as being relatively large, the sample average may be replaced by the ensemble average. The decision making behavior of an individual is then considered as a realization of a stochastic process that is affected by the ensemble average of the process. In the context of voting behavior, the two-state model introduced in Section 3.1 was applied [17]. The members of the US House of Representatives had the choice to vote for or against passing the so-called bailout bill (US House of Representatives, USA, 2008). The bill failed to pass in the September 29 vote. In a second vote, on October 3, the bill finally passed. That is, a critical number of “No” voters changed their opinions between September 29 and October 3. The states k=1 and k=2 were labeled by “Y” and “N” representing “Yes” and “No” voters. According to model (37), the conditional probabilities T(N→Y) and T(Y→Y) describing a change in the opinion from “No” to “Yes” and a consistent “Yes” opinion were considered. Using the notation of (37), the voting model proposed in Frank [18] reads
(41)p(Y,t+1)=[T(Y⟶Y,t,p(Y,t))-T(N⟶Y,t,p(Y,t))]p(Y,t)+T(N⟶Y,t,p(Y,t)).
The data available in the literature suggested that T(Y→Y)=1. For the T(N→Y) the function T(N→Y)=A-B·p(Y,t) was used to describe the impact of the group opinion on the voting behavior. This relationship assumes that the pressure to change from a “No” voter to a “Yes” voter was high immediately subsequent to the September 29 votes, when the probability p(Y,t) at t = Sept 29 was not high enough to make the bill pass. In this context, it should be mentioned that the failure of the bill to pass on September 29 triggered a dramatic breakdown of the US stock market. That is, the members of the House of Representatives who had voted with “No” were not only subjected to political pressure but also experienced essential pressure from the economy. For T(Y→Y)=1 and T(N→Y)=A-Bp(Y,t) the stationary occupation probability p(Y,st) is given by p(Y,st)=A/B such that (41) eventually reads
(42)p(Y,t+1)=p(Y,t)+B[p(Y,st)-p(Y,t)][1-p(Y,t)].
The stationary value p(Y,st) was estimated from the data. A detailed model-based analysis of the data showed that the members of the two main parties, the Democratic and Republican members, were differently affected by the group pressure.


## 3.4. Psychology and Human Memory

A pioneering experiment in human memory research is the free recall experiment in which participants are asked to recall items of a particular category [25, 26]. For example, they are asked to recall animal names. Typically, participants are instructed to recall as many items as possible and to do so as fast as possible and without repetition. The total number of recalled items is by definition a monotonically increasing function of time. Strikingly, free recall involves clusters in which item subcategories are recalled more or less as a whole. For example, when recalling animal names a participant may recall a number of insect names in succession. From a theoretical point of view, at any given time point during the free recall process the ability to recall items in clusters depends on the size of the pool of items available for recalling. In other words, if there is a large pool of items that have not yet been recalled, then there is a high chance that a whole item subcategory can be recalled. In line with this argument, a stochastic model for free recall dynamics was developed in which the recall process depends on the size of the pool of items available for recall [27]. More precisely, time was parceled in discrete time intervals in which recall attempts are executed. The aim was to develop a model for the probability that a single item has been or has not yet been recalled at a given time step t. To this end, the two-state model of Section 3.1 was applied. Accordingly, the states k=1 and k=2 correspond to being recalled or being not yet been recalled and are labeled by “R” and “N,” respectively. The two relevant conditional probabilities are T(R→R) and T(N→R). By the design of the experiment, we have T(R→R)=1. Moreover, as argued previously, if the pool of items that have not yet recalled is large, that is, if the probability p(N) that an item has not yet been recalled is high, then the conditional probability T(N→R) should be high as well. Likewise, if the probability p(R) that an item has been recalled is high, then T(N→R) should be relatively low. Therefore, it was suggested to put T(N→R)=B·(1-p(R,t)) with B>0. This function implies that the occupation probability p(R,t) converges to the stationary fixed p(R,st)=1 at which T(N→R)=0. The full model can be obtained by substituting T(R→R) and T(N→R) into (37) and reads
(43)p(R,t+1)=p(R,t)+B[1-p(R,t)]2.
Note that from a mathematical point of view, this model is a special case of model (42) for voting behavior (hint: put p(Y,st)=1 in (42)). Since the model describes the single item recall probability, the model must be scaled to the number of items recalled by any given participant. To this end, a parameter C describing the maximal number of items for recall was introduced. The unknown parameters B and C were estimated for free recall data collected in a study by Rhodes and Turvey [28], and the sample mean values B=0.009 (SD = 0.005) and C=307 (SD = 138) were found for a sample of eight participants [27]. Note that (43) predicts that the fixed point p(R,st)=1 will never be reached in a finite number of time steps. Consequently, the model predicts that the total number of recalled items M at the end of the experiment is smaller than the number of available items for recall. This prediction was confirmed by the data. Moreover, the model-based analysis revealed a statistically significant positive correlation between M and C for the participant sample studied in Frank and Rhodes [27].


## 3.5. Biology and Evolutionary Population Dynamics

A primary concern of evolutionary population dynamics is the understanding of the origin of life. Among the many models that have been proposed in this context, the studies by Crutchfield and Görnerup [29] and Görnerup and Crutchfield [30] are of interest for our review. In these studies a number of N macromolecules with different functional properties were considered. These macromolecules compete with each other for survival during evolutionary time scales. However, not only competition but also cooperative behavior is taken into account. Roughly speaking, it is assumed that in the early stages of the development of life, macromolecules of different types existed and could change their types due to interactions with other macromolecules of the same or different types. Eventually, only a subset of functionally different types survived the selection process. Let k=1,…,N denote the functionally different types of macromolecules and p(k,t) the occupation probability of the kth type. Then, at every discrete reaction time step t the macromolecules can change their types such that transitions from type i to type k occur. This so-called metamodel can be formulated in terms of a nonlinear Markov chain model [31]. Accordingly, the conditional (transition) probability T(i→k) describes the probability of a transition from type i to type k given that a macromolecule is of type i at time step t. As suggested by the studies of Crutchfield and Görnerup [29] and Görnerup and Crutchfield [30] interactions between two different types of macromolecules can be modeled by assuming that the evolution of occupation probabilities is affected by terms of the form p(j,t)p(i,t). More precisely, in general, a macromolecule of type i may interact with a macromolecule of type j and turn into a macromolecule of type k. It can be shown that from the metamodel it follows that T(i→k) assumes the form [31]
(44)T(i⟶k)=∑j=1Ng(ijk)p(j,t).
In (44) the parameters g(ijk)≥0 are weight coefficients that describe the relevance of interactions between macromolecules of different types. Note that the conditional probabilities T(i→k) satisfy the normalization condition that the sum over all probabilities with index k equals 1. This normalization can be guaranteed by requiring that the sum over all weight coefficients g with index k equals 1. Substituting (44) into (7), the metamodel for evolutionary population dynamics can be cast into a nonlinear Markov chain model and reads
(45)p(k,t+1)=∑i,j=1Ng(ijk)p(j,t)p(i,t).
It has been shown that the nonlinear Markov chain model (45) is a useful tool for investigating evolutionary population dynamics. First of all, trajectories that describe the evolution of individual macromolecules were computed numerically using (1), and (8)–(11); see Section 2.2. That is, the model allows for generating temporal sequences of type changes for individual macromolecules. Second, it was shown that the degree to which a type k is attractive can be quantified. In particular, for evolutionary selection process that become stationary, the attractors of the surviving macromolecule types can be analyzed and distinctions between weakly and strongly attractive types can be made (e.g., see Figure 3 in [31]).


## 4. Time-Discrete Strongly Nonlinear Stochastic Processes on Continuous State Spaces

### 4.1. The Generalized Autoregressive Model of Order 1 as a Strongly Nonlinear Markov Process

The generalized autoregressive model of order 1 defined by (16) was studied in Frank [13] as time-discrete version of the time-continuous strongly nonlinear Shimizu-Yamada model; see Section 6. Let us discuss this generalized autoregressive model in more detail. To make this discussion more self-consistent, (16) is presented here once again
(46)X(t+1)=AX(t)+B〈X(t)〉+gε(t)
as (46). Since the random variable ε exhibits a normal distribution at any time step t, the conditional probability density can be calculated that describes the distribution of X at time t+1 given the value of X=y at the previous time step t and the probability density P(x,t) of X at time step t. Assuming that the variance of ε equals 2, the solution reads
(47)T(x,t+1∣y,t;θx[P(·,t)])=T(x,t+1∣y,t;〈X(t)〉)=1g4πexp(-[x+Ay+B〈X(t)〉]24g2).
As indicated in (47) the conditional probability depends only on the first moments (i.e., the mean value) of the probability density P(x,t). By means of the conditional probability density P(x,t+1) can be computed from P(x,t) like
(48)P(x,t+1)=∫-∞∞T(x,t+1∣y,t;〈X(t)〉)P(y,t)dy.
From (46) it follows that the mean value M1(t)=〈X(t)〉 evolves like
(49)M1(t+1)=(A+B)M1(t)⟹M1(t)=M1(1)[A+B](t-1).
For 0<A+B<1 the mean value is an exponentially decaying function. For -1<A+B<0 the mean value alternate between negative and positive values but decays in amplitude monotonically. In summary, for -1<A+B<1 the mean value converges to zero in the limiting case t→∞. A detailed calculation shows that the second moment M2=〈X2(t)〉 evolves like
(50)M2(t)=AM2(t-1)+ABM1(t-1)+BM1(t)M1(t-1)+2g2.
For -1<A+B<1 and -1<A<1 the first moment converges to zero which implies that the second moment converges to a stationary value given by M2(st)=2g2/(1-A2) such that the variance of the stationary processes is determined by Var(X)=M2(st). Let us substitute the variable u=g·ε (which is a normally distributed random variable with variance 2g2) into (40). Then, for -1<A+B<1 and -1<A<1 the stationary variance of X is given by
(51)Var(X)=Var(u)1-A2.
Not surprisingly, (51) is just the expression for the variance of the ordinary stationary autoregressive processes of order 1. Moreover, for normally distributed initial values X(1), the probability density P(x,t) satisfies a normal distribution for all times t. This follows from (47) and (48) or alternatively from the linearity of (46). Consequently, for -1<A+B<1 and -1<A<1 the probability density P(x,t) converges in this case to a normal distribution with zero mean value and variance given by (51). In the stationary case, the correlation function Corr(s)=〈X(t)X(t-s)〉/Var(X) for s≥0 has exactly the same form as the correlation function of the ordinary autoregressive process of order 1 because the mean value M1 equals zero and drops out of (46):
(52)Corr(s+1)=A·Corr(s).
In conclusion, for -1<A+B<1 and -1<A<1 the generalized autoregressive process (46) behaves in the stationary case exactly as the ordinary autoregressive process (i.e., the process with B=0). The two processes however exhibit different transient characteristic properties.


### 4.2. The Generalized Deterministic Autoregressive Process of Order 1 and the Sensitivity of Strongly Nonlinear Processes to Initial Conditions

In the deterministic case, we putg=0 in (13) such that
(53)X(t+1)=h(X(t),t,θX(t)[P(x,t)]).
Let us consider the case N=1 and assume that the drift function h(·) depends only linearly on the state X(t). If, in addition, h(·) does not explicitly depend on time and θ[·] is a map (functional or linear operator) that maps P(x,t) to a real number, we arrive at the model
(54)X(t+1)=h(θ[P(x,t)])X(t)=L[P(x,t)]X(t).
In (54) we simplified the notation by introducing the operator L(·) that maps P(x,t) to a real number (e.g., L calculates expectation values of X and then uses them as arguments in a function). Equation (54) may be considered as a generalized deterministic autoregressive model of order 1. Model (54) was studied in detail in Frank [32]. A striking feature of the model is that the stationary probability density P(x,st) if it exists depends crucially on the initial probability density P(x,t=1). For sake of readability, let us put u(x)=P(x,t=1). Then, it can be shown that P(x,st) if it exists is a rescaled function of u(·) like [32]
(55)P(x,st)=1zu(xz),
with the scaling parameter
(56)z=lims→∞∏t=1sL[P(x,t)].
In other words, the strongly nonlinear process (54) might exhibit infinitely many stationary probability densities. This was demonstrated more explicitly for a specific form of L(·) which will be discussed next.


### 4.3. Economics and the Emergence of Stationary Volatilities as a Group Dynamics Process

In economics the standard deviation or the variance of a price of an asset is a measure for how risky the assist is. That is, the variability measures reflect the beliefs of traders how risky it is to hold the asset. Zero variability reflects a risk-free asset. In this context, the variability of a price is also called the price volatility. Volatility is observed by traders and informs individual traders how risky the market as a whole considers a particular asset. On the other hand, the market is composed of individual traders. An autoregressive model of the form (54) can be used to model such a circular relationship. To this end, we assume that the operator L(·) calculates the variance of X. In this case, the evolution of X is determined by the variance of X, on the one hand. On the other hand, the variance by definition is calculated from the realization of X. Let us consider the following special case of (54) (for details see [32]):
(57)X(t+1)=[1-γ(Var(X)-β)]X(t),
with γ,β>0. It can be shown that for γ<2/β the process exhibits stable but not asymptotically stationary probability densities. The shape of these probability densities depends on the initial distributions as discussed in Section 4.2. However, the variance of all stable stationary distributions equals β. This is intuitively clear, when we realize that for Var(X)=β, (57) becomes the identity X(t+1)=X(t). Model (57) also exhibits unstable stationary probability densities. One of these unstable solutions is the Dirac delta function P(x)=δ(x). It can be shown that any process converges to a stationary process with variance β provided that the initial probability density of the process has a variance Var(X) smaller than a certain critical value. For example, if the Dirac delta function is perturbed slightly such that the variance becomes finite but is still close to zero, then the process will not converge back to the Dirac delta function. Rather, the process will converge to a stationary process with variance equal to β.We distinguished previously between stable and asymptotically stable probability densities. Stable but not asymptotically stable stationary solutions have also been called marginally stable stationary solutions. Recently, research has been focused on the importance of marginally stabile stationary solutions for biochemical and biological systems [33–35]. In the context of strongly nonlinear processes, stable but not asymptotically stable probability densities or alternatively marginally stable probability densities refer to processes that exhibit a family of stationary probability densities with two properties [9]. First, by an infinitely small perturbation a particular stationary probability density can be transformed into another stationary probability density of the family of probability densities. Second, the family of probability densities as a whole acts as an attractor. That is, for any perturbation that does not transform a member of the family into another member of the family, the perturbed stationary probability density eventually converges to one of the members of the family of marginally stable probability densities. A simple example can be given for model (57). As mentioned previously, for the marginally stable probability densities with Var(X)=β, (57) reads X(t+1)=X(t). A shift of the coordinate system is consistent with the identity X(t+1)=X(t) and does not affect the variance. Therefore, any stationary probability density with Var(X)=β can be shifted along the x-axis by an infinitely small amount to give rise to another stationary probability density. Clearly, this kind of perturbation will not decay.In closing the discussion on model (54), let us return to the interpretation of the model as a model for volatility dynamics. The model predicts that the emergence of a stationary volatility in the market is not guaranteed but depends on market parameters (here the inequality γ<2/β must be satisfied) and the initial variance. Furthermore, the model predicts that the price itself is not affected whether or not a stable, stationary volatility measure emerges.


### 4.4. Phase Synchronization of Inhomogeneous Ensemble of Globally Coupled Oscillatory Units

Phase synchronization of coupled oscillatory subunits has frequently been studied by means of time-continuous strongly nonlinear stochastic model, see Section6. However, also time-discrete models have been considered. Phase synchronization is a special case of synchronization when the amplitude dynamics is neglected. Accordingly, each oscillatory unit is described by a single real number representing a phase X(t). Let us envision each oscillatory unit as an oscillator evolving along a closed orbit in an appropriately defined state space. Then the phase X(t) is the coordinate along that closed orbit. The phase X(t) is a periodic variable. The oscillators are de-synchronized if X has a uniform distribution on the domain of definition. If X exhibits a nonuniform distribution (in particular, a unimodal distribution), then the oscillators are partially synchronized. If X=A for all oscillators, then there is complete synchronization.Let us consider an all-to-all coupling between the oscillators. In this case, the evolution of the phaseX(t) of an individual oscillator depends on the evolution of all other oscillator phases. If the set of oscillators is relatively large, the limit of infinitely many oscillators may be considered as a good approximation and the sample of oscillators may be replaced by a statistical ensemble. In this case, the phase dynamics of an individual oscillator depends on the expectation value of the statistical ensemble of oscillators. Moreover, any individual oscillator represents a realization of the ensemble. In summary, a closed description of the phase dynamics in terms of a strongly nonlinear stochastic process can be obtained.In the deterministic case, we putg=0 in (13) such that (13) becomes (53) mentioned in Section 4.2. Without loss of generality, the period can be put equal to unity such that X is defined in the interval [0,1], and X=0 and X=1 indicate the same location on the closed orbit. In line with seminal works by Kuramoto [36], Daido [37] proposed an explicit form of (53) that can be written for individual realizations like
(58)Xr(t+1)=Xr(t)+Ur-κq2πsin(2πXr(t)),q=limR→∞1R∑r=1Rcos(2πXr(t)).
The variable q is the order parameter of the oscillator system. For desynchronized oscillators (i.e., for a uniformly distributed X) we have q=0. If q differs from zero, the oscillators are partially synchronized. In (58) the parameter κ≥0 describes the coupling strength between the oscillators. More precisely, it is the coupling between individual oscillators and the mean field force defined by hMF(X)=q·sin(2πX)/(2π) that acts on an individual oscillator and is generated by all oscillators. In (58) the parameter U denotes the phase velocity for the uncoupled oscillator. As indicated, U is taken from a distribution and differs among the oscillators. Therefore, model (58) describes an inhomogeneous ensemble of globally coupled oscillatory units.Model (58) is a useful model for studying the emergence of phase synchronization from the perspective of nonequilibrium phase transitions. At a critical value of the coupling parameter κ, a nonequilibrium phase transition occurs and the probability density P(x) bifurcates from a uniform distribution to a nonuniform distribution indicating that the ensemble makes a transition from a desynchronized state to a partially synchronized state [37]. For Lorentzian distributed phase velocities U, an analytical expression of the critical coupling parameter can be found [37]. Similar aspects of time-discrete models for synchronizations have been studied in Daido [38, 39] and by Teramae and Kuramoto [40].Following more recent work [32], the conditional probability density of the phase synchronization model (58) in terms of the so-called Frobenius-Perron operator involving the Dirac delta function δ(·) can be determined and reads
(59)T(x,t+1∣X(t)=y,θ[P(·,t)],U)=δ(mod[x-y-U+κq2πsin(2πy),1])q=∫01cos(2πx)P(x,t)dx.
The modulus-1 operator may be eliminated by casting (59) into the form
(60)T(x,t+1∣X(t)=y,θ[P(·,t)],U)=∑j=-∞∞δ(x-y-U+κq2πsin(2πy)+j).
The conditional probability density holds for a particular unit with phase velocity U. Let ρ(U) describe the probability density of phase velocities. Then, P(x,t+1) can at least formally be computed from P(x,t) and ρ(U) via the integral equation
(61)P(x,t+1)=∫y=1,U=-∞0,∞T(x,t+1∣X(t)=y,θ[P(·,t)],U)×P(y,t)ρ(U)dydU.
Accordingly, stationary solutions P(x,st) satisfy
(62)P(x,st)=∫y=1,U=-∞0,∞T(x,t+1∣X(t)=y,θ[P(·,st)],U)×P(y,st)ρ(U)dydU.
The uniform distribution P(x)=1 satisfies (62) for any coupling parameter κ≥0 and consequently is a stationary probability density. However, for oscillator systems with Lorentzian velocity distributions ρ(U), the uniform distribution is an unstable stationary probability density if the coupling parameter exceeds the critical value which can be shown by numerical simulations [37].


## 4.1. The Generalized Autoregressive Model of Order 1 as a Strongly Nonlinear Markov Process

The generalized autoregressive model of order 1 defined by (16) was studied in Frank [13] as time-discrete version of the time-continuous strongly nonlinear Shimizu-Yamada model; see Section 6. Let us discuss this generalized autoregressive model in more detail. To make this discussion more self-consistent, (16) is presented here once again
(46)X(t+1)=AX(t)+B〈X(t)〉+gε(t)
as (46). Since the random variable ε exhibits a normal distribution at any time step t, the conditional probability density can be calculated that describes the distribution of X at time t+1 given the value of X=y at the previous time step t and the probability density P(x,t) of X at time step t. Assuming that the variance of ε equals 2, the solution reads
(47)T(x,t+1∣y,t;θx[P(·,t)])=T(x,t+1∣y,t;〈X(t)〉)=1g4πexp(-[x+Ay+B〈X(t)〉]24g2).
As indicated in (47) the conditional probability depends only on the first moments (i.e., the mean value) of the probability density P(x,t). By means of the conditional probability density P(x,t+1) can be computed from P(x,t) like
(48)P(x,t+1)=∫-∞∞T(x,t+1∣y,t;〈X(t)〉)P(y,t)dy.
From (46) it follows that the mean value M1(t)=〈X(t)〉 evolves like
(49)M1(t+1)=(A+B)M1(t)⟹M1(t)=M1(1)[A+B](t-1).
For 0<A+B<1 the mean value is an exponentially decaying function. For -1<A+B<0 the mean value alternate between negative and positive values but decays in amplitude monotonically. In summary, for -1<A+B<1 the mean value converges to zero in the limiting case t→∞. A detailed calculation shows that the second moment M2=〈X2(t)〉 evolves like
(50)M2(t)=AM2(t-1)+ABM1(t-1)+BM1(t)M1(t-1)+2g2.
For -1<A+B<1 and -1<A<1 the first moment converges to zero which implies that the second moment converges to a stationary value given by M2(st)=2g2/(1-A2) such that the variance of the stationary processes is determined by Var(X)=M2(st). Let us substitute the variable u=g·ε (which is a normally distributed random variable with variance 2g2) into (40). Then, for -1<A+B<1 and -1<A<1 the stationary variance of X is given by
(51)Var(X)=Var(u)1-A2.
Not surprisingly, (51) is just the expression for the variance of the ordinary stationary autoregressive processes of order 1. Moreover, for normally distributed initial values X(1), the probability density P(x,t) satisfies a normal distribution for all times t. This follows from (47) and (48) or alternatively from the linearity of (46). Consequently, for -1<A+B<1 and -1<A<1 the probability density P(x,t) converges in this case to a normal distribution with zero mean value and variance given by (51). In the stationary case, the correlation function Corr(s)=〈X(t)X(t-s)〉/Var(X) for s≥0 has exactly the same form as the correlation function of the ordinary autoregressive process of order 1 because the mean value M1 equals zero and drops out of (46):
(52)Corr(s+1)=A·Corr(s).
In conclusion, for -1<A+B<1 and -1<A<1 the generalized autoregressive process (46) behaves in the stationary case exactly as the ordinary autoregressive process (i.e., the process with B=0). The two processes however exhibit different transient characteristic properties.


## 4.2. The Generalized Deterministic Autoregressive Process of Order 1 and the Sensitivity of Strongly Nonlinear Processes to Initial Conditions

In the deterministic case, we putg=0 in (13) such that
(53)X(t+1)=h(X(t),t,θX(t)[P(x,t)]).
Let us consider the case N=1 and assume that the drift function h(·) depends only linearly on the state X(t). If, in addition, h(·) does not explicitly depend on time and θ[·] is a map (functional or linear operator) that maps P(x,t) to a real number, we arrive at the model
(54)X(t+1)=h(θ[P(x,t)])X(t)=L[P(x,t)]X(t).
In (54) we simplified the notation by introducing the operator L(·) that maps P(x,t) to a real number (e.g., L calculates expectation values of X and then uses them as arguments in a function). Equation (54) may be considered as a generalized deterministic autoregressive model of order 1. Model (54) was studied in detail in Frank [32]. A striking feature of the model is that the stationary probability density P(x,st) if it exists depends crucially on the initial probability density P(x,t=1). For sake of readability, let us put u(x)=P(x,t=1). Then, it can be shown that P(x,st) if it exists is a rescaled function of u(·) like [32]
(55)P(x,st)=1zu(xz),
with the scaling parameter
(56)z=lims→∞∏t=1sL[P(x,t)].
In other words, the strongly nonlinear process (54) might exhibit infinitely many stationary probability densities. This was demonstrated more explicitly for a specific form of L(·) which will be discussed next.


## 4.3. Economics and the Emergence of Stationary Volatilities as a Group Dynamics Process

In economics the standard deviation or the variance of a price of an asset is a measure for how risky the assist is. That is, the variability measures reflect the beliefs of traders how risky it is to hold the asset. Zero variability reflects a risk-free asset. In this context, the variability of a price is also called the price volatility. Volatility is observed by traders and informs individual traders how risky the market as a whole considers a particular asset. On the other hand, the market is composed of individual traders. An autoregressive model of the form (54) can be used to model such a circular relationship. To this end, we assume that the operator L(·) calculates the variance of X. In this case, the evolution of X is determined by the variance of X, on the one hand. On the other hand, the variance by definition is calculated from the realization of X. Let us consider the following special case of (54) (for details see [32]):
(57)X(t+1)=[1-γ(Var(X)-β)]X(t),
with γ,β>0. It can be shown that for γ<2/β the process exhibits stable but not asymptotically stationary probability densities. The shape of these probability densities depends on the initial distributions as discussed in Section 4.2. However, the variance of all stable stationary distributions equals β. This is intuitively clear, when we realize that for Var(X)=β, (57) becomes the identity X(t+1)=X(t). Model (57) also exhibits unstable stationary probability densities. One of these unstable solutions is the Dirac delta function P(x)=δ(x). It can be shown that any process converges to a stationary process with variance β provided that the initial probability density of the process has a variance Var(X) smaller than a certain critical value. For example, if the Dirac delta function is perturbed slightly such that the variance becomes finite but is still close to zero, then the process will not converge back to the Dirac delta function. Rather, the process will converge to a stationary process with variance equal to β.We distinguished previously between stable and asymptotically stable probability densities. Stable but not asymptotically stable stationary solutions have also been called marginally stable stationary solutions. Recently, research has been focused on the importance of marginally stabile stationary solutions for biochemical and biological systems [33–35]. In the context of strongly nonlinear processes, stable but not asymptotically stable probability densities or alternatively marginally stable probability densities refer to processes that exhibit a family of stationary probability densities with two properties [9]. First, by an infinitely small perturbation a particular stationary probability density can be transformed into another stationary probability density of the family of probability densities. Second, the family of probability densities as a whole acts as an attractor. That is, for any perturbation that does not transform a member of the family into another member of the family, the perturbed stationary probability density eventually converges to one of the members of the family of marginally stable probability densities. A simple example can be given for model (57). As mentioned previously, for the marginally stable probability densities with Var(X)=β, (57) reads X(t+1)=X(t). A shift of the coordinate system is consistent with the identity X(t+1)=X(t) and does not affect the variance. Therefore, any stationary probability density with Var(X)=β can be shifted along the x-axis by an infinitely small amount to give rise to another stationary probability density. Clearly, this kind of perturbation will not decay.In closing the discussion on model (54), let us return to the interpretation of the model as a model for volatility dynamics. The model predicts that the emergence of a stationary volatility in the market is not guaranteed but depends on market parameters (here the inequality γ<2/β must be satisfied) and the initial variance. Furthermore, the model predicts that the price itself is not affected whether or not a stable, stationary volatility measure emerges.


## 4.4. Phase Synchronization of Inhomogeneous Ensemble of Globally Coupled Oscillatory Units

Phase synchronization of coupled oscillatory subunits has frequently been studied by means of time-continuous strongly nonlinear stochastic model, see Section6. However, also time-discrete models have been considered. Phase synchronization is a special case of synchronization when the amplitude dynamics is neglected. Accordingly, each oscillatory unit is described by a single real number representing a phase X(t). Let us envision each oscillatory unit as an oscillator evolving along a closed orbit in an appropriately defined state space. Then the phase X(t) is the coordinate along that closed orbit. The phase X(t) is a periodic variable. The oscillators are de-synchronized if X has a uniform distribution on the domain of definition. If X exhibits a nonuniform distribution (in particular, a unimodal distribution), then the oscillators are partially synchronized. If X=A for all oscillators, then there is complete synchronization.Let us consider an all-to-all coupling between the oscillators. In this case, the evolution of the phaseX(t) of an individual oscillator depends on the evolution of all other oscillator phases. If the set of oscillators is relatively large, the limit of infinitely many oscillators may be considered as a good approximation and the sample of oscillators may be replaced by a statistical ensemble. In this case, the phase dynamics of an individual oscillator depends on the expectation value of the statistical ensemble of oscillators. Moreover, any individual oscillator represents a realization of the ensemble. In summary, a closed description of the phase dynamics in terms of a strongly nonlinear stochastic process can be obtained.In the deterministic case, we putg=0 in (13) such that (13) becomes (53) mentioned in Section 4.2. Without loss of generality, the period can be put equal to unity such that X is defined in the interval [0,1], and X=0 and X=1 indicate the same location on the closed orbit. In line with seminal works by Kuramoto [36], Daido [37] proposed an explicit form of (53) that can be written for individual realizations like
(58)Xr(t+1)=Xr(t)+Ur-κq2πsin(2πXr(t)),q=limR→∞1R∑r=1Rcos(2πXr(t)).
The variable q is the order parameter of the oscillator system. For desynchronized oscillators (i.e., for a uniformly distributed X) we have q=0. If q differs from zero, the oscillators are partially synchronized. In (58) the parameter κ≥0 describes the coupling strength between the oscillators. More precisely, it is the coupling between individual oscillators and the mean field force defined by hMF(X)=q·sin(2πX)/(2π) that acts on an individual oscillator and is generated by all oscillators. In (58) the parameter U denotes the phase velocity for the uncoupled oscillator. As indicated, U is taken from a distribution and differs among the oscillators. Therefore, model (58) describes an inhomogeneous ensemble of globally coupled oscillatory units.Model (58) is a useful model for studying the emergence of phase synchronization from the perspective of nonequilibrium phase transitions. At a critical value of the coupling parameter κ, a nonequilibrium phase transition occurs and the probability density P(x) bifurcates from a uniform distribution to a nonuniform distribution indicating that the ensemble makes a transition from a desynchronized state to a partially synchronized state [37]. For Lorentzian distributed phase velocities U, an analytical expression of the critical coupling parameter can be found [37]. Similar aspects of time-discrete models for synchronizations have been studied in Daido [38, 39] and by Teramae and Kuramoto [40].Following more recent work [32], the conditional probability density of the phase synchronization model (58) in terms of the so-called Frobenius-Perron operator involving the Dirac delta function δ(·) can be determined and reads
(59)T(x,t+1∣X(t)=y,θ[P(·,t)],U)=δ(mod[x-y-U+κq2πsin(2πy),1])q=∫01cos(2πx)P(x,t)dx.
The modulus-1 operator may be eliminated by casting (59) into the form
(60)T(x,t+1∣X(t)=y,θ[P(·,t)],U)=∑j=-∞∞δ(x-y-U+κq2πsin(2πy)+j).
The conditional probability density holds for a particular unit with phase velocity U. Let ρ(U) describe the probability density of phase velocities. Then, P(x,t+1) can at least formally be computed from P(x,t) and ρ(U) via the integral equation
(61)P(x,t+1)=∫y=1,U=-∞0,∞T(x,t+1∣X(t)=y,θ[P(·,t)],U)×P(y,t)ρ(U)dydU.
Accordingly, stationary solutions P(x,st) satisfy
(62)P(x,st)=∫y=1,U=-∞0,∞T(x,t+1∣X(t)=y,θ[P(·,st)],U)×P(y,st)ρ(U)dydU.
The uniform distribution P(x)=1 satisfies (62) for any coupling parameter κ≥0 and consequently is a stationary probability density. However, for oscillator systems with Lorentzian velocity distributions ρ(U), the uniform distribution is an unstable stationary probability density if the coupling parameter exceeds the critical value which can be shown by numerical simulations [37].


## 5. Time-Continuous Strongly Nonlinear Stochastic Processes on Discrete State Spaces

Strongly nonlinear stochastic processes evolving on discrete state spaces but exhibiting a continuous time variable have been studied in the context of nonlinear master equations as introduced in Section2.4; that is, in the context of Markov processes. Frequently, nonlinear master equations have been studied as a tool to derive nonlinear Fokker-Planck equations (see, e.g., [14, 41–44]). However, nonlinear master equations have also been studied in their own merit (see, e.g., [45]).


### 5.1. Nonlinear Master Equations as a Tool for Identifying the Generic Forms of Nonlinear Fokker-Planck Equations

As mentioned in Section2.4, transition rates of master equations may depend on occupation probabilities. Curado and Nobre [14] considered only transition rates between neighboring states k. In addition, the explicit time dependency of rates was neglected. In this case, (20) reads
(63)ddtp(k,t)=w(k-1⟶k,p(t))p(k-1,t)+w(k+1⟶k,p(t))p(k+1,t)-p(k,t)[w(k⟶k+1,p(t))+w(k⟶k-1,p(t))].
Note that there are only two types of transition rates. First, there are transition rates that describe the rate of transitions from a level to the next higher level (“upward rates”). Second, there are the transition rates of transitions from a level to the next lower level (“down rates”). Using w(k,up,p)=w(k→k+1,p) and w(k,down,p)=w(k→k-1,p), we can re-write (63) as
(64)ddtp(k,t)=w(k-1,up,p(t))p(k-1,t)+w(k+1,down,p(t))p(k+1,t)-p(k,t)[w(k,up,p(t))+w(k,down,p(t))].
Furthermore, Curado and Nobre [14] assumed that the states k represent positions on a one-dimensional grid with spacing Δ. They assumed that the transition rates depend on the spacing Δ but also on the occupation probabilities p(k,t). More precisely, the model involves the following up- and downrates:
(65)w(k,up,p(t))=AΔ2[p(k,t)]q-1,w(k,down,p(t))=-1Δh(kΔ)+AΔ2[p(k,t)]q-1.
In the limiting case of an infinitely small grid (i.e., for Δ→0), the master equation (64) with (65) behaves like the nonlinear Fokker-Planck equation of the form (29). The diffusion term of that Fokker-Planck equation is proportional to the probability density Pq. This model, in turn, is a standard model for diffusion in porous media [46, 47]; see Section 6.


### 5.2. Strongly Discrete-Valued Nonlinear Stochastic Processes Maximizing Generalized Entropy Measures

A key property of stochastic processes described by ordinary master equations such as (20) with rates w(i→k) independent of the occupation probabilities p(k) is that the processes approach stationarity in the long time limit [48]. More precisely, processes evolve such that the Kullback distance measure [49] is a monotonically decreasing function of time and approaches a stationary value. The stationarity of the Kullback distance measure indicates that the process under consideration has become stationary. The Kullback distance measure is defined on the basis of the Boltzmann-Gibbs-Shannon entropy, which is a fundamental entropy measure not only for equilibrium systems [50] but also for nonequilibrium systems [51]. It was suggested to exploit this link between the Boltzmann-Gibbs-Shannon entropy, the Kullback distance measure, and the ordinary master equation to define nonlinear master equations based on generalized entropy and generalized Kullback distance measures [9, 45]. Accordingly, entropy measures S of the form
(66)S(p)=∑k=1Ns(p(k))
are considered that involve a kernel function s(·). Furthermore, it is assumed that the kernel s(·) satisfies certain properties. The second derivative of the kernel s(z) with respect to z is negative for any argument z in the interval [0,1] which implies that the entropy S(p) is concave [52]. In addition, the first derivative of the kernel s(z) with respect to z in the limiting case z→0 goes to plus infinity (i.e., exhibits a singularity). This property is important for the master equation that will be defined in the following and guarantees that occupation probabilities p(k) solving the master equation do not become negative. The following master equation has been proposed [45]:
(67)ddtp(k,t)=∑i=1NA(i⟶k,t)F[p(i,t)]-F[p(k,t)]∑i=1NA(k⟶i,t)F(z)=exp{-1-ds(z)dz}.
The nonlinearity in the master equation (67) is conceptually different from the nonlinearity in the master equation (20). While in (20) it is assumed that the transition rates depend on occupation probabilities, in (67) the transition rates A(i→k) are independent of the occupation probabilities. However, transitions are not proportional to the occupation probabilities p(k) as in (20). Rather, they are proportional to the terms F(p(k)), whose explicit form depends on the entropy kernel s(·). In view of this property, transitions are entropy driven. It can be shown that stochastic processes satisfying (67) converge to occupation probabilities that maximize the entropy measures S. Therefore, we may say that the transitions of processes defined by (67) are proportional to an entropy drive F(·) that maximizes the entropy S under consideration. Note that for the special case of the Boltzmann-Gibbs-Shannon entropy the function F reduces to the identity F(z)=z, which implies that (67) includes the ordinary, linear master equation as special case.In closing our considerations on the nonlinear master equation (67), it is useful to note that formally (67) can be cast into the form of the nonlinear master equation (20) irrespective of any conceptual differences. If we put
(68)w(i⟶k,t,p(i,t))p(i,t)=A(i⟶k,t)F[p(i,t)]⟹w(i⟶k,t,p(i,t))=A(i⟶k,t)F[p(i,t)]p(i,t)
for p(i,t)>0, then model (67) assumes the form of (20). Note that in the limiting case p→0, we have F→0 because of the aforementioned assumption that ds(z)/dz→∞ holds for z→0. This, in turn, implies that w(i→k,t,p(i,t)) for p(i,t)=0 can be defined as the product of A(i→k,t) and the limiting case “0/0,” where the ratio “0/0” has to evaluated for any given entropy kernel s(z). The relation (68) is useful because it can be used to compute realizations of the nonlinear master equation (67) from (1), (18), and (22)–(25).


### 5.3. Physiological, Partitioned (Discrete-Valued) Signals and Strongly Nonlinear Stochastic Processes Exhibiting Stationary Cut-Off Distributions

It has been argued that any probability density with tails that extend to infinity fails to capture the boundedness of physiological signals [53]. In other words, while probability densities with tails that extend to infinity such as the normal distribution are crucially important for developing theory and for modeling and are frequently good approximations, they imply that signals can assume arbitrarily large values in magnitude. This is not plausible. Physiological signals such as muscle force produced by human and animals, limb velocities, limb accelerations, blood flow velocities, heart rates, electroencephalographic recorded brain signals, pulse rates of neurons, and dendritic currents are bounded and cannot exceed certain maximum values. Taking an information theoretical point of view [51], stochastic processes describing physiological signals may maximize information or entropy measures. However, they cannot maximize the Boltzmann-Gibbs-Shannon measure under ordinary constraints because this measure yields normal distributions or other distributions with tails that extend to infinity. In contrast, physiological signals may maximize generalized information and entropy measures that exhibit the so-called cut-off distributions as maximum entropy distributions. This argument was developed for stochastic processes defined on continuous state spaces [53] but holds for processes evolving on discrete-state spaces as well. In particular, signals from sensory systems are known to be partitioned in a certain way. Differences between two magnitudes can only be detected if they exceed certain thresholds. For this reason, modeling of sensory physiological signals and physiological signals in general may assume that the state space of consideration is of discrete nature.In Frank [45] the coordination of rhythmic movements of two limbs has been addressed exploiting the master equation model (67). It has been assumed that the physiological relevant variable, namely, the relative phase between the limbs, is appropriately described on a partitioned (i.e., discrete-value) state space. Furthermore, it has been assumed that the coordination dynamics corresponds to a stochastic process maximizing an entropy measure that exhibits a maximum entropy cut-off distribution. Out of all possible entropy measures, the entropy measure, nowadays called Tsallis entropy, has been used for that purpose [54–56]. The model can be considered as a modified version of the seminal coordination model by Haken et al. [57]. The benefits of the nonlinear master equation model (67) relative to the original Haken-Kelso-Bunz model are twofold. First, model (67) offers a consistent approach that addresses physiological data within the framework of cut-off distributions. Second, the model predicts that coordination is bistable in a distributional sense. That is, under certain conditions the model exhibits two stationary distributions p(k,st). This feature is in line with experimental observations. For human coordination models with continuous state spaces similar attempts have been made to deal with bistability in a distributional sense [58, 59].


## 5.1. Nonlinear Master Equations as a Tool for Identifying the Generic Forms of Nonlinear Fokker-Planck Equations

As mentioned in Section2.4, transition rates of master equations may depend on occupation probabilities. Curado and Nobre [14] considered only transition rates between neighboring states k. In addition, the explicit time dependency of rates was neglected. In this case, (20) reads
(63)ddtp(k,t)=w(k-1⟶k,p(t))p(k-1,t)+w(k+1⟶k,p(t))p(k+1,t)-p(k,t)[w(k⟶k+1,p(t))+w(k⟶k-1,p(t))].
Note that there are only two types of transition rates. First, there are transition rates that describe the rate of transitions from a level to the next higher level (“upward rates”). Second, there are the transition rates of transitions from a level to the next lower level (“down rates”). Using w(k,up,p)=w(k→k+1,p) and w(k,down,p)=w(k→k-1,p), we can re-write (63) as
(64)ddtp(k,t)=w(k-1,up,p(t))p(k-1,t)+w(k+1,down,p(t))p(k+1,t)-p(k,t)[w(k,up,p(t))+w(k,down,p(t))].
Furthermore, Curado and Nobre [14] assumed that the states k represent positions on a one-dimensional grid with spacing Δ. They assumed that the transition rates depend on the spacing Δ but also on the occupation probabilities p(k,t). More precisely, the model involves the following up- and downrates:
(65)w(k,up,p(t))=AΔ2[p(k,t)]q-1,w(k,down,p(t))=-1Δh(kΔ)+AΔ2[p(k,t)]q-1.
In the limiting case of an infinitely small grid (i.e., for Δ→0), the master equation (64) with (65) behaves like the nonlinear Fokker-Planck equation of the form (29). The diffusion term of that Fokker-Planck equation is proportional to the probability density Pq. This model, in turn, is a standard model for diffusion in porous media [46, 47]; see Section 6.


## 5.2. Strongly Discrete-Valued Nonlinear Stochastic Processes Maximizing Generalized Entropy Measures

A key property of stochastic processes described by ordinary master equations such as (20) with rates w(i→k) independent of the occupation probabilities p(k) is that the processes approach stationarity in the long time limit [48]. More precisely, processes evolve such that the Kullback distance measure [49] is a monotonically decreasing function of time and approaches a stationary value. The stationarity of the Kullback distance measure indicates that the process under consideration has become stationary. The Kullback distance measure is defined on the basis of the Boltzmann-Gibbs-Shannon entropy, which is a fundamental entropy measure not only for equilibrium systems [50] but also for nonequilibrium systems [51]. It was suggested to exploit this link between the Boltzmann-Gibbs-Shannon entropy, the Kullback distance measure, and the ordinary master equation to define nonlinear master equations based on generalized entropy and generalized Kullback distance measures [9, 45]. Accordingly, entropy measures S of the form
(66)S(p)=∑k=1Ns(p(k))
are considered that involve a kernel function s(·). Furthermore, it is assumed that the kernel s(·) satisfies certain properties. The second derivative of the kernel s(z) with respect to z is negative for any argument z in the interval [0,1] which implies that the entropy S(p) is concave [52]. In addition, the first derivative of the kernel s(z) with respect to z in the limiting case z→0 goes to plus infinity (i.e., exhibits a singularity). This property is important for the master equation that will be defined in the following and guarantees that occupation probabilities p(k) solving the master equation do not become negative. The following master equation has been proposed [45]:
(67)ddtp(k,t)=∑i=1NA(i⟶k,t)F[p(i,t)]-F[p(k,t)]∑i=1NA(k⟶i,t)F(z)=exp{-1-ds(z)dz}.
The nonlinearity in the master equation (67) is conceptually different from the nonlinearity in the master equation (20). While in (20) it is assumed that the transition rates depend on occupation probabilities, in (67) the transition rates A(i→k) are independent of the occupation probabilities. However, transitions are not proportional to the occupation probabilities p(k) as in (20). Rather, they are proportional to the terms F(p(k)), whose explicit form depends on the entropy kernel s(·). In view of this property, transitions are entropy driven. It can be shown that stochastic processes satisfying (67) converge to occupation probabilities that maximize the entropy measures S. Therefore, we may say that the transitions of processes defined by (67) are proportional to an entropy drive F(·) that maximizes the entropy S under consideration. Note that for the special case of the Boltzmann-Gibbs-Shannon entropy the function F reduces to the identity F(z)=z, which implies that (67) includes the ordinary, linear master equation as special case.In closing our considerations on the nonlinear master equation (67), it is useful to note that formally (67) can be cast into the form of the nonlinear master equation (20) irrespective of any conceptual differences. If we put
(68)w(i⟶k,t,p(i,t))p(i,t)=A(i⟶k,t)F[p(i,t)]⟹w(i⟶k,t,p(i,t))=A(i⟶k,t)F[p(i,t)]p(i,t)
for p(i,t)>0, then model (67) assumes the form of (20). Note that in the limiting case p→0, we have F→0 because of the aforementioned assumption that ds(z)/dz→∞ holds for z→0. This, in turn, implies that w(i→k,t,p(i,t)) for p(i,t)=0 can be defined as the product of A(i→k,t) and the limiting case “0/0,” where the ratio “0/0” has to evaluated for any given entropy kernel s(z). The relation (68) is useful because it can be used to compute realizations of the nonlinear master equation (67) from (1), (18), and (22)–(25).


## 5.3. Physiological, Partitioned (Discrete-Valued) Signals and Strongly Nonlinear Stochastic Processes Exhibiting Stationary Cut-Off Distributions

It has been argued that any probability density with tails that extend to infinity fails to capture the boundedness of physiological signals [53]. In other words, while probability densities with tails that extend to infinity such as the normal distribution are crucially important for developing theory and for modeling and are frequently good approximations, they imply that signals can assume arbitrarily large values in magnitude. This is not plausible. Physiological signals such as muscle force produced by human and animals, limb velocities, limb accelerations, blood flow velocities, heart rates, electroencephalographic recorded brain signals, pulse rates of neurons, and dendritic currents are bounded and cannot exceed certain maximum values. Taking an information theoretical point of view [51], stochastic processes describing physiological signals may maximize information or entropy measures. However, they cannot maximize the Boltzmann-Gibbs-Shannon measure under ordinary constraints because this measure yields normal distributions or other distributions with tails that extend to infinity. In contrast, physiological signals may maximize generalized information and entropy measures that exhibit the so-called cut-off distributions as maximum entropy distributions. This argument was developed for stochastic processes defined on continuous state spaces [53] but holds for processes evolving on discrete-state spaces as well. In particular, signals from sensory systems are known to be partitioned in a certain way. Differences between two magnitudes can only be detected if they exceed certain thresholds. For this reason, modeling of sensory physiological signals and physiological signals in general may assume that the state space of consideration is of discrete nature.In Frank [45] the coordination of rhythmic movements of two limbs has been addressed exploiting the master equation model (67). It has been assumed that the physiological relevant variable, namely, the relative phase between the limbs, is appropriately described on a partitioned (i.e., discrete-value) state space. Furthermore, it has been assumed that the coordination dynamics corresponds to a stochastic process maximizing an entropy measure that exhibits a maximum entropy cut-off distribution. Out of all possible entropy measures, the entropy measure, nowadays called Tsallis entropy, has been used for that purpose [54–56]. The model can be considered as a modified version of the seminal coordination model by Haken et al. [57]. The benefits of the nonlinear master equation model (67) relative to the original Haken-Kelso-Bunz model are twofold. First, model (67) offers a consistent approach that addresses physiological data within the framework of cut-off distributions. Second, the model predicts that coordination is bistable in a distributional sense. That is, under certain conditions the model exhibits two stationary distributions p(k,st). This feature is in line with experimental observations. For human coordination models with continuous state spaces similar attempts have been made to deal with bistability in a distributional sense [58, 59].


## 6. Time-Continuous Strongly Nonlinear Stochastic Processes on Continuous State Spaces

There is a relatively large literature on time-continuous strongly nonlinear stochastic processes defined on continuous state spaces. Most of these studies are concerned with the Markov diffusion processes. Reviews can be found in Frank [9, 60]. Strongly nonlinear phase synchronization models based on the Kuramoto model [36] are reviewed in Acebrón et al. [61]. In this section, some of applications in physics and the life sciences will be highlighted without going into much detail.


### 6.1. Chaos in Probability: Time-Continuous Case

In Section3.2 strongly nonlinear time-discrete stochastic processes were considered that exhibit occupation probabilities that evolve in a chaotic fashion. As an example, the logistic map model defined by (39) was discussed. A key feature of this model and of similar models is that the chaotic behavior arises from the coupling of the dynamics of the realizations to the occupation probabilities. In particular, without this coupling equation (39) reads p(A,t+1)=T0=α/4 with fixed parameters α taken from the interval [0,4] and clearly does not exhibit a chaotic solution. Chaos is due to the fact that transition probability T0=α/4 is replaced by a transition probability that depends on the process occupation probabilities like T0=α·p(A,t)(1-p(A,t)). That is, probability measures of strongly nonlinear processes can exhibit chaotic behavior due to the very nature of strongly nonlinear processes, which is the influence of a probability measure on realizations from which the probability measure is calculated. Such a circular causality structure may be realized in many-body systems composed of interacting subsystems. In such systems, chaos may emerge due to the coupling of the single subsystem dynamics to the relative frequency with which states are occupied by subsystems. Chaos emerges even if the isolated subsystem dynamics is not chaotic [18]. Consequently, chaos may be considered as a self-organization phenomenon. Recall that self-organizing phenomena in general are emerging properties of many-body systems that do not exist on the level of the subsystems [62]. In our context, the many-body system as a whole behaves qualitatively differently from the individual, isolated parts. The catchphrases “the whole is different from the sum of its parts” or “more is different” apply [63].This feature of strongly nonlinear stochastic processes has also been demonstrated for time-continuous models involving continuous state spaces [64–67]. Subsystems described by ordinary stochastic differential equations that do not exhibit chaotic solutions have been coupled together and the limiting case of a many-body system composed of infinitely many subsystems has been considered. The limiting case models were described by strongly nonlinear stochastic processes in the sense defined in Section 1.2. In particular, in these studies the dynamics of realization was coupled to certain expectation values of the respective processes. It was found that the expectation values under appropriate conditions exhibit a chaotic dynamics.


### 6.2. Plasma Physics and the Landau Form of Strongly Nonlinear Fokker-Planck Equations

Plasma physics is concerned with the evolution of many-particle systems composed of charged particles. The particles can interact with each other in various ways. Two interaction forms are particle-particle collisions and Coulomb interaction. The Landau form of the Fokker-Planck equation accounts for these two interactions and describes the evolution of the particle density in the three-dimensional velocity space. That is, letv=(v1,v2,v3) denote the velocity vector of a particle. Let ρ(v,t) denote the particle mass density at time t. Assuming that particles are neither created nor destroyed, the total mass M is invariant over time. The particle probability density P(v,t)=ρ(v,t)/M can be introduced. According to the Landau form, the probability density P(v,t) satisfies a strongly nonlinear Fokker-Planck equation of the form (29). More precisely, the Landau form reads [68–74]
(69)∂∂tP(v,t)=-∑k=13∂∂vk[hk(v,θv[P(·,t)])P(v,t)]+∑j=1,k=13∂2∂vj∂vk[Djk(·)P(v,t)],hk(v,θv[P(·,t)])=A∂∂vk∫ΩP(v′,t)|v-v′|d3v′,Djk(·)=B∂2∂vj∂vk∫Ω|v-v′|P(v′,t)d3v′.
Stationary solutions P(v,st) of (69) have been studied, for example, by Soler et al. [75]. In addition, several studies have been concerned with the derivation of transient solutions P(v,t) using different (semi)numerical approaches [74, 76–79].


### 6.3. Astrophysics: Stellar Cut-Off Mass Distributions Defined by Strongly Nonlinear Stochastic Models

Astrophysical mass distributions may exhibit relative sharp boundaries. More precisely, particle distributions in the stellar space may be spatially bounded due to gravity. The gravitational forces are produced by the mass distributions themselves. That is, the dynamics of an individual particle of a mass distribution is affected by the particle distribution as a whole. In turn, the distribution is composed of its individual particles. In view of this circular causality, it does not come as a surprise that strongly nonlinear stochastic models have been investigated that describe cut-off distributions of stellar particles [80–84]. The models typically assume the form of the strongly nonlinear Fokker-Planck equation (29).Two more comments may be appropriate. First, in Section5.3 we addressed cut-off distributions in the context of physiological signals. However, the motivation in Section 5.3 for considering cut-off distributions was quite different from the argument used in astrophysical applications. Second, in addition to the aforementioned studies on strongly nonlinear processes describing astrophysical cut-off mass distribution, strongly nonlinear stochastic processes satisfying the Landau form (69) have also been used for studying astrophysical problems [85, 86].


### 6.4. Accelerator Physics: Shortening of Bunch-Particle Distribution of Particle Beams

Particles in accelerator beams can travel in localized cluster of particles. In a longitudinal beam, the so-called bunch-particle distribution is a mass density that describes the distribution of particles with respect to the center of the cluster. Letx1 denote relative position of a particle with respect to the bunch center. Typically, the particle dynamics also depends on the particle energy. That is, beam particle dynamics involves a second coordinate x2 which is a rescaled energy measure (rather than a velocity or momentum variable). Dividing the mass density with the total mass, the bunch-particle distribution becomes a probability density P(x,t) of the two-dimensional vector x=(x1,x2). Particles in accelerator beams are charged particles that are accelerated by means of electromagnetic forces. Although particles may interact with each other by means of Coulomb forces, the crucial force that determines the shape and length of a cluster is the so-called wake-field [87, 88]. The wake-field is an electromagnetic field generated by the whole bunch of particles that travels ahead of the bunch but is reflected at the accelerator walls and then acts back on the particles of the bunch. The wake-field can cause a shortening of the particle bunch. That is, under the impact of the wake-field the cluster is smaller in size than it would be theoretically in the absence of the wake-field. The understanding of bunch shortening is an important step towards an understanding of beam instabilities that should be avoided.The dynamics of bunch particles is typically described by strongly nonlinear Markov diffusion processes. The reason for this is that the dynamics of individual particles is affected by the wake-field which can be calculated from an appropriate expectation value of the bunch particle probability densityP(x,t). As a result, in various studies it has been proposed that P(x,t) satisfies a strongly nonlinear Fokker-Planck equation of the form (29) (see, e.g., [89–98]). A useful model, in particular for the purpose of theory development, is the Haissinski model [95, 96, 99–104] for which analytical solutions of the reduced density P(x1) can be derived. In particular, according to the Haissinski model, the dynamics of single particles satisfies a strongly nonlinear Langevin equation (27), which reads [100]
(70)ddtX1(t)=X2(t)ddtX2(t)=h(X(t),θX1(t)[P(x,t)])+gΓ2(t)h=-X1(t)-γX2(t)-∂∂X1∫ΩUW(X1-x1′)P(x1′,x2,t)dx1′dx2,
where γ>0 describes a phenomenological damping constant. For the Haissinski model the wake-field function UW assumes the form of a Dirac delta function: UW(z)=κ·δ(z). The parameter κ measures the strength of the impact of the wake-field. For κ<0 the width of the reduced probability density P(x1) becomes smaller when increasing the magnitude |κ| of the impact of the wake-field. In doing so, model (70) provides a dynamic account for the squeezing effect of wake-fields on particle clusters traveling in accelerator beams.


### 6.5. Physics of Liquid Crystal Phase Transitions from the Perspective of Strongly Nonlinear Stochastic Processes

Liquid crystals typically consist of rod-shape macromolecules that interact with each other by means of weak Van-der-Waals forces. Due to this interaction, molecules can align themselves such that the molecules point with their primary axes more or less in the same direction. The physical properties of a crystal with aligned molecules are qualitatively different from the properties that the crystal exhibits when the macromolecules are isotropically (randomly) distributed. Therefore, liquid crystals exhibit two phases: an ordered phase with molecules that are aligned at least to some degree, which is called the nematic phase, and a disorder phase with molecules pointing in random directions, which is called the isotrope phase. A phase transition from the isotropic to the nematic phase occurs when the temperature is decreased below a critical temperature [105–108]. The degree of orientational order can be captured by the Maier-Saupe order parameter [102, 105–111]. For the isotropic phase the order parameter equals zero. In the nematic phase, the Maier-Saupe order parameter is larger than zero.Stochastic models describing the emergence of orientational order (i.e., isotropic-nematic phase transitions of liquid crystals) exploit the notion that the order parameter itself expresses the magnitude of an overall force that acts on individual molecules and tends to align them with the mean orientation of the liquid crystal macromolecules. The model is based conceptually on the notions of circular causality and self-organization: individual molecules are affected by the order parameter and at the same time constitute the order parameter. Hess [112] as well as Doi and Edwards [113] have proposed a strongly nonlinear stochastic process describing the rotational Brownian motion of the molecules under the impact of the (self-generated) order parameter. The Hess-Doi-Edwards model exhibits the structure of a strongly nonlinear Fokker-Planck equation (29) and can be cast into the form of a strongly nonlinear Langevin equation (27) (see, e.g., [114]). Exploiting the concept of strongly nonlinear stochastic processes, the martingale for the Hess-Doi-Edwards model can be defined [60].Transient solutions of the Hess-Doi-Edwards model may be computed from approximate coupled moment equations [115]. Using Lyapunov’s direct method [62, 116], it can be shown that the ordered state exhibits an asymptotically stable probability density [114]. Importantly, since stochastic models provide a dynamic account, it is possible to study response functions of liquid crystals. Transient responses describing the relaxation to equilibrium [117, 118] as well as correlation functions and mean square displacement functions in the stationary case [114] can be determined at least semi-numerically.Beyond the application of the concept of strongly nonlinear stochastic processes to the Hess-Doi-Edwards Fokker-Planck equation, in general, strongly nonlinear stochastic models have turned out to be useful models in liquid polymer physics (see, e.g., [119–121]).


### 6.6. Classical Descriptions of Quantum Systems by means of Strongly Nonlinear Stochastic Processes

The relaxation towards equilibrium of quantum mechanical Fermi and Bose systems can be modeled using a classical approach by means generalized Boltzmann equations that exhibit appropriately modified collision integrals [122–125]. Applying a diffusion approximation to these collision integrals, such quantum mechanical Boltzmann equations reduce to Fokker-Planck equations that are nonlinear with respect to their probability densities [122–124]. In addition to this approach via the Boltzmann equation, alternative derivations of classical quantum mechanical Fokker-Planck equations that are nonlinear with respect to probability densities have been proposed [126–133]. A key property of these models is that they exhibit stationary probability densities that correspond to Fermi or Bose distributions. Interpreting these models in terms of strongly nonlinear Fokker-Planck equations of form (29) allows us not only to consider the evolution of probability densities but also to examine predictions of semiclassical stochastic processes for quantum systems. For example, solutions of trajectories computed from the corresponding strongly nonlinear Langevin equations of form (27) have been studied [126, 127]. In particular, the relaxation of the free electron gas towards its stationary Fermi energy distribution can be determined by computing numerically individual electron trajectories in the relevant energy space [9, 126]. Moreover, martingales for quantum processes within this classical Langevin approach can be defined [60].Finally, note that classical stochastic descriptions of Fermi and Bose systems do not necessarily involve strongly nonlinear stochastic processes. It has been shown that ordinary Fokker-Planck equations with appropriately defined drift functionsh(·) can also exhibit Fermi and Bose distributions as stationary probability densities (see, e.g., [134]).


### 6.7. Material Physics of Porous Media

Gas particles and fluids diffusing through porous media satisfy a nonlinear diffusion equation, frequently called the porous media equation. In one spatial dimension, the porous medium equation for the particle mass densityρ(x,t) reads [3, 9, 46, 47, 135]
(71)∂∂tρ(x,t)=A∂2∂x2[ρ(x,t)]q,
with A>0. Dividing the mass density ρ(x,t) by the total mass M, (71) becomes an evolution equation for a probability density P(x,t) normalized to unity
(72)∂∂tP(x,t)=D∂2∂x2[P(x,t)]q,
with D=A·M(q-1), and assumes the form of the nonlinear Fokker-Planck equation (29). The connection between the porous medium equation and nonlinear master equations has been addressed in Section 5.1 (see the aforementioned). Equations (71) and (72) exhibit exact time-dependent solutions, frequently called Barenblatt-Pattle solutions [136, 137].The parameterq captures the nonlinearity of the equation. For q=1 (linear case) the porous medium equation reduces to the ordinary diffusion equations and the variance of the diffusion process (or the second moment of P(x,t) assuming a zero first moment) increases linearly with time t. In contrast, for q>1/3 and q≠1 (nonlinear case) the variance increases as a nonlinear function of t like tr with r=2/(1+q) as can be shown by various methods [9, 138–140].Interpreting (72) in terms of a strongly nonlinear Fokker-Planck equation (29), trajectories of realizations can be computed from the corresponding strongly nonlinear Langevin equation [138]. For a generalized version of (72) involving a drift force such a Langevin equation approach has been exploited to derive analytical expressions of certain autocorrelation functions [141].The porous medium equation in its form as a nonlinear Fokker-Planck equation (72) plays an important role in the theory development of nonextensive thermostatistics. As it turns out, a particular generalization of (72), the so-called Plastino-Plastino model, exhibits stationary solutions that maximize the nonextensive entropy proposed by Tsallis. We will return to this issue later.


### 6.8. Material Physics and the Liouville Dynamics of Many-Body Systems with Interacting Subsystems

The particle distribution of a many-body system in the overdamped, deterministic case satisfies a Liouville equation. For many-body systems with interacting subsystems the Liouville equation involves an integral term describing the interaction force on a single subsystem. This interaction integral may be evaluated using a diffusion approximation. In doing so, the Liouville equation assumes the form of the nonlinear Fokker-Planck equation (29) when considering the probability density rather than the subsystem density. The nonlinearity occurs in the diffusion term. This approach has been developed in a series of studies by Andrade et al. [142] and Ribeiro et al. [143, 144]. For example, the distribution functions of interacting particles in dusty plasmas, interacting vortices, and interacting polymer chains in complex fluids have been studied within this framework.Note that as mentioned previously the dynamics of the subsystems is deterministic. Nevertheless, the effective model for the subsystem probability density involves a diffusion term. That is, according to the effective, approximative model in terms of a nonlinear Fokker-Planck equation, due to the interaction between the subsystems the many-body system as a whole generates a fluctuating force that acts on the individual subsystems of the many-body system.


### 6.9. Nonextensive Physical Systems and the Plastino-Plastino Model

Tsallis (Tsallis, [54, 55]; see also Abe and Okamoto [56]) introduced in physics a generalized form of the Boltzmann-Gibbs-Shannon entropy, nowadays frequently called the Tsallis entropy. The Tsallis entropy exhibits a parameter q. For q=1 the entropy reduces to the Boltzmann-Gibbs-Shannon entropy capturing properties of extensive systems. For q≠1 the Tsallis entropy describes nonextensive systems.  A. R Plastino and A. Plastino [145] proposed a nonlinear Fokker-Planck equation of the form (29) that may be written like
(73)∂∂tP(x,t)=-∂∂x[h(x)P(x,t)]+D∂2∂x2[P(x,t)]q.
The model describes the probability density P(x,t) of a particle living in a one-dimensional continuous state space given by the whole real line. The term involving a force function h(x) in (73) is linear with respect to the probability density P(x,t). The diffusion term of the model involves the parameter q, which using the notation suggested by (73) corresponds to the parameter q of the Tsallis entropy (see Frank [9]). Note that in the original model a slightly different notation was proposed [145]. A key feature of the Plastino-Plastino model (73) is that stationary probability densities P(x,st) of (73) maximize the Tsallis entropy. For q=1, that is, in the special case in which the Tsallis entropy reduces to the Boltzmann-Gibbs-Shannon entropy, the model is linear with respect to the probability density P(x,t). For q≠1, that is, for the general case that the Tsallis entropy describes a nonextensive system, the diffusion term is nonlinear with respect to P(x,t). In view of these observations, the Plastino-Plastino model establishes a connection between nonextensivity and nonlinearity.The Plastino-Plastino model (73) has been examined in various studies. In particular, it has frequently been pointed out that (73) may be considered as a generalized version of the porous medium model (72) for gas and fluid flows that are subjected to an additional driving force captured by the force function h(x).Exact analytical solutions for the time-dependent probability densityP(x,t) are available in the case of a linear drift force [140, 145, 146]. Trajectories describing the stochastic dynamics of individual particles can be computed when interpreting (73) as a strongly nonlinear Fokker-Planck equation by means of the corresponding strongly nonlinear Langevin equation [138, 141, 147]. In this case, second order statistics measures such as autocorrelation functions can be determined [141] and martingales can be found [60]. Exact time-dependent solutions have been derived for generalized versions of model (73) that account for source terms [148–150]. The Kramers escape rate has been determined for processes described by the Plastino-Plastino model involving an appropriately defined drift force [151]. The multivariate generalization of (73) with [139, 152, 153] and without [129, 154] time-dependent coefficients has been considered. The Plastino-Plastino model (73) with explicit state-dependent diffusion coefficients D(x) has been studied [153, 155]. The model has been generalized for the Sharma-Mittal entropy that involves two parameters and includes the Tsallis entropy as a special case [156]. Interestingly, H-theorems have been found that show that transient solutions P(x,t) of (73) converge to stationary probabilities densities P(x,st) provided that h(x) is chosen appropriately [122, 157–164].


### 6.10. Oscillatory Coupled Nonequilibrium Systems: Canonical-Dissipative Approach

Coupled oscillators with short-range interactions can be studied within the framework of strongly nonlinear stochastic processes [165]. In this study isolated oscillators were modeled using the canonical-dissipative approach [166–168] that allows to exploit the concepts of Hamiltonian and Newtonian dynamics, on the one hand, and to address nonequilibrium systems such as self-excited limit cycle oscillators, on the other hand.


### 6.11. Phase Synchronization and Kuramoto Model: Time-Continuous Case

Synchronization is a phenomenon studied in various disciplines ranging from physics to the social sciences [169]. In the context of strongly nonlinear stochastic processes, we are primarily concerned with the synchronization of a large ensemble of oscillatory subunits. Synchronization can be studied both in the phase space of the oscillators (e.g., the space spanned by position and momentum variables) and in reduced spaces. An example for the latter case was discussed already in Section 4.4: phase synchronization. In fact, we will focus in this section on the phenomenon of phase synchronization.A benchmark model that describes the emergence of phase synchronization in a population of phase oscillators is the Kuramoto model [36]. Accordingly, each oscillator is described by a phase X that evolves on a continuous time scale t and represents a periodic variable. Each oscillator is coupled to all remaining oscillators and the limiting case of a group of infinitely many oscillators is considered. The oscillators as a whole generate an order parameter, which is a complex-valued variable. The complex-valued order parameter has an amplitude, the so-called cluster amplitude A, and an angle, the so-called cluster phase α. Both cluster phase α and cluster amplitude A are measures defined in circular statistics (see, e.g., [36, 170, 171]). Mathematically speaking, for an oscillator phase X defined on the interval [0,2π] the complex order parameter q is defined by the expectation value
(74)q(θ[P(x,t)])=limR→∞1R∑r=1Rexp(iXr(t))=∫x=02πexp(ix)P(x,t)dx=A(t)exp(iα(t)),
where i is the imaginary unit (i.e., the square root of −1) and the cluster phase α is chosen such that A≥0 in any case. Note that both A and α are time-dependent variables.The order parameterq induces a force that acts on the individual phase oscillators. That is, the oscillators interact with each other indirectly via the self-generated order parameter. The population of phase oscillators can be regarded as a statistical ensemble. Accordingly, the order parameter can be computed as an expectation value from the probability density P(x,t), that is, the probability density of the phase of an arbitrarily chosen phase oscillator. Since the order parameter acts back on the realizations (individual phase oscillators), the model satisfies the definition of a strongly nonlinear stochastic process provided in Section 1.2. The strongly nonlinear Langevin equation of the Kuramoto model reads
(75)ddtX(t)=-κA(θ[P(x,t)])sin(X(t)-α(θ[P(x,t)]))+gΓ(t)
and assumes the form (27). In (75), the parameter κ≥0 measures the strength of the force induced by the order parameter q.The uniform probability densityP=1/(2π) describing a desynchronized oscillator population is a stationary solution for any parameter κ because for P=1/(2π) the cluster amplitude A equals zero. However, only for κ<2g2 the uniform distribution is asymptotically stable. For κ>2g2 the uniform distribution is unstable meaning that for any initial condition that slightly deviates from the uniform probability density P=1/(2π) the strongly nonlinear stochastic process will not converge to a stationary process with uniform probability density. Rather, for κ>2g2 the Kuramoto model exhibits stable, stationary unimodal distributions [9, 36]. These unimodal probability densities indicate that the oscillator population is partially synchronized. Moreover, the order parameter amplitude A is larger than zero for these unimodal stationary probability densities. The unimodal probability densities are stable but not asymptotically stable (see, e.g., [9]). A shift of such a stationary probability density along the x-axis transforms a member of the family of stable stationary probability densities into another member of that family. This feature becomes intuitively clear when we realize that the form of (75) is invariant against transformation X→X+s and α→α+s, where s is an arbitrarily small real number. We may use the term “marginally” stable to characterize the stability of these stationary probability densities (see also Section 4.3).The Kuramoto model has been reviewed in Strogatz [172] and Acebrón et al. [61]. In particular, there exists an H-theorem of the Kuramoto model showing that transient probability densities converge to stationary ones [173]. This H-theorem is related to a free energy approach to the Kuramoto model (Frank [174]; see also Frank [9]).A key question in the context of the Kuramoto model concerns the bifurcation from the desynchronized state to the synchronized state for oscillator populations with inhomogeneous eigenfrequencies. In this case, (75) may be written for realizations Xr of the process like
(76)ddtXr(t)=Ur-κA(θ[P(x,t)])sin(Xr(t)-α(θ[P(x,t)]))+gΓr(t),
where Ur is the phase velocity of the rth realization of the process. If we restrict our considerations to symmetric distributions of velocities and to symmetric initial probability densities, then the symmetry property remains conserved during the evolution of the process and the cluster phase can assume only one of the two values α=0 or α=π. Moreover, the order parameter q becomes a real-valued variable. It can be shown that in this case (76) reduces to
(77)ddtXr(t)=Ur-κqsin(Xr(t))+gΓr(t)q=limR→∞1R∑r=1Rcos(Xr(t)).
Comparing (77) with the Daido model (58) discussed in Section 4.4, it becomes at this stage clear that the Daido model (58) can be considered as a time-discrete counterpart of the time-continuous, fully symmetric Kuramoto model (77) with distributed eigenfrequencies. Note that the Daido model exhibits phases defined on the interval [0,1], whereas in the context of the Kuramoto model typically phases defined on the interval [0,2π] are considered.As mentioned previously, reviews for the Kuramoto model are available in the literature and the reader may consult these works for more details ([61, 172]; see also Tass [175] and Section 7.5 in [9]).


### 6.12. Biology: Population Dispersion and Population Dynamics

The spatial distribution of insects forming swarms may be described in terms of cut-off distributions. For example, the insect densityρ(x)=A·[{1-B·|x|}+]2 has been proposed [176], where the coordinate x measures the location of an insect with respect to the center of the swarm and A,B>0. The operator {·}+ corresponds to the identify for positive arguments and equals zero for negative arguments (i.e., {z}+=z for z>0 and {z}+=0 otherwise). Normalizing ρ by the total mass M, a stochastic model for the probability density P(x)=ρ(x)/M needs to explain the emergence of a stationary insect cut-off probability density. In fact, to this end, a model similar to the Plastino-Plastino model (73) with an appropriate drift term was proposed [176, 177]:
(78)∂∂tP(x,t)=-∂∂x[γsgn(x)P(x,t)]+D∂∂x{[P(x,t)]μ∂∂xP(x,t)},
with γ>0. The stationary probability P(x,st)=C·[{1-B·|x|}+]2 is obtained for μ=0.5. However, exponents different from μ=0.5 have also been proposed [178, 179]. Moreover, descriptions of population dynamics in terms of nonlinear Fokker-Planck equations alternative to (78) have been considered in the literature as well ([81, 168]; see also the Shigesada-Teramoto model in Okubo and Levin, [176], Section 5.6).


### 6.13. Social Sciences and Group Decision Making: Time-Continuous Case

Group decision making has already been addressed in Sections3.3 and 4.3 in the context of time-discrete models. In a series of studies, Boster and coworkers developed the so-called linear discrepancy model of group decision making and compared predictions with experimental data [180–182]. In particular, the model accounts for a key phenomenon in the social sciences: the risky shift. A risky shift occurs when people arrive at riskier decisions when discussing a given topic in a group than when making their own decisions individually. According to the linear discrepancy model due to discussions the opinion in favor or against a particular issue shifts by an amount that is proportional to the opinion that an individual holds and the opinion that is presented in the group by another individual. The linear discrepancy model predicts that a risky shift can be observed when the initial distribution (i.e., the prediscussion distribution) of opinions is skewed. Accordingly, during the discussion session, the opinion distribution tends to become more symmetric, which is accompanied with a shift of the mean value. Some of the mathematical properties of the linear discrepancy model have been studied in more detail within the framework of a strongly nonlinear stochastic process [183]. Interestingly, the stationary, symmetric probability densities describing the distribution of opinions among the group members are only stable and not asymptotically stable. In particular, a shift along the opinion axis transforms one member of the family of stable, stationary probability density into another member. In this regard, the group decision making model exhibits the same feature as the Kuramoto model.


### 6.14. Finance and Option Pricing

A stochastic option pricing model has been developed on the basis of the Plastino-Plastino model (73). In particular, the option pricing model exhibits a diffusion term similar to the one of the Plastino-Plastino model [184–186]. A particular benefit of the model is that it features non-Gaussian distributions. This is due to the nonlinearity of the diffusion term (or the nonextensivity of the Tsallis entropy measure involved in the definition of the process). In fact, non-Gaussian distributions frequently fit financial data better than Gaussian distributions.


### 6.15. Economics and Modeling the Dynamics of Target Leverage Ratios

The total capital (firm value) of a company is typically composed of borrowed money and equity. The leverage of a company or the leverage ratio is the ratio of borrowed money relative to equity. A company needs to decide what leverage ratio the company should have. The decision is not necessarily made once and for all. Rather, the desired leverage ratio (i.e., the target leverage ratio) may be updated dynamically depending on the internal and external circumstances. For this reason, the leverage ratios of companies frequently correspond to dynamic variables subjected to fluctuations.Lo and Hui [187] proposed a strongly nonlinear stochastic model for the dynamics of the leverage ratio. The model is based on a model developed earlier by Hui et al. [188] that involved an explicitly time-dependent function. In the strongly nonlinear stochastic model, this function is eliminated, and in this sense the strongly nonlinear stochastic model can be regarded as being closed. In order to achieve this closure, Lo and Hui [187] assumed that the leverage dynamics of a company is affected by the leverage ratios of similar companies. Within the framework of strongly nonlinear stochastic processes this implies that the leverage ratio dynamics of companies depends on an expectation value of the leverage ratio process. Formally, the model by Lo and Hui is closely related to the Shimizu-Yamada model that will be discussed later.


### 6.16. Engineering Sciences: Generators for Noise Sources

In engineering sciences and related disciplines a common problem is to simulate numerically signals of noise sources. These signals can then be used to test the response of engineered systems to random signals emerging from noise sources. A characteristic feature of noise sources is that they are stationary. In view of this property, strongly nonlinear stochastic processes can be defined which for any initial probability density are stationary [147, 189, 190]. For example, the strongly nonlinear Langevin equation
(79)ddtX(t)=-AX(t)+g(X(t),θ[P(x,t)])Γ(t)g=BP(y,t)[C-∫-∞yP(z,t)dz]|y=X(t),
with A,B,C>0, corresponds to a nonlinear Fokker-Planck equation of the form (29) whose right-hand side is by definition equal to zero [9, 147]. Consequently, the Langevin equation defines for any initial probability density P(x,t=0) a stochastic process with a stationary probability density P(x). The parameters A, B, C can be chosen in order to select a desired correlation time of the process.


### 6.17. Further Benchmark Models

#### 6.17.1. The Shimizu-Yamada Model

The Shimizu-Yamada model is motivated by research on muscle contraction (Shimizu [191], and Shimizu and Yamada [192]; see also Section 5.5.4 in [9]). The Shimizu-Yamada model corresponds to the generalized Ornstein-Uhlenbeck process that is defined by (33) and involves a mean-field force proportional to the mean value of the process. The Shimizu-Yamada model is highly relevant for the theory development of strongly nonlinear Markov diffusion processes because it can be solved exactly. That is, exact transient solutions can be found and exact solutions for the conditional probability density T(x,t∣y,s,〈X(s)〉) are available. The conditional probability density, in turn, can be used to calculate all correlation functions of interest both in the transient and in the stationary cases. For details, see Frank [9, 193]. Since the model exhibits exact solutions, it has also been used to test the accuracy of numerical algorithms for solving nonlinear Markov diffusion processes  [16].


#### 6.17.2. The Desai-Zwanzig Model

The Desai-Zwanzig model involves a mean-field force proportional to the mean value of the process just as the Shimizu-Yamada model. However, the individual, isolated subsystems are assumed to perform an overdamped motion in a double-well potential [194, 195]. The strongly nonlinear Fokker-Planck equation of the Desai-Zwanzig model reads
(80)∂∂tP(x,t)=-∂∂x[h(x,θ[P(x,t)])P(x,t)]+D∂2∂x2P(x,t)h(x,θ[P(x,t)])=Ax-Bx3-κ(x-M(t))M(t)=∫-∞∞xP(x,t)dx
for a stochastic process defined on the whole real line and A,B,κ>0. The term Ax-Bx3in (80) describes the potential force derived from the aforementioned double-well potential. The term -κ(x-M) describes the mean-field force proportional to the mean value of the process. The force attracts the realizations of the process towards the mean value of the process. The parameter κ measures the strength of the mean-field force. The model is symmetric with respect to x=0. In view of this observation it is intuitively clear that the model exhibits a symmetric stationary probability density with M=0 for all parameters κ. It can be shown that for parameters κ smaller than a certain critical value this symmetric probability density is asymptotically stable and is the only stationary probability density. However, when κ exceeds the critical value, the symmetric stationary probability density becomes unstable. Two asymptotically stable, stationary probability densities emerge with mean values M=z and M=-z, where z>0 [9, 194, 196]. The mean value M has typically been considered as order parameter of a many-body system described by the strongly nonlinear stochastic process (80). For M=0 the many-body system exhibits a disordered state. For M≠0 the system exhibits some order. Therefore, the Desai-Zwanzig model describes an order-disorder phase transition.The special importance of the Desai-Zwanzig model is its feature that provides a fundamental stochastic model for an order-disorder phase transition. The Desai-Zwanzig model and the Kuramoto model may be viewed as the two key models in this regard. Both models are able to capture order-disorder phase transitions. While the Kuramoto model is a prototype system for many-body systems with state variables subjected to periodic boundary conditions, the Desai-Zwanzig model is a prototype system for many-body systems featuring state variables satisfying natural boundary conditions. Moreover, the Desai-Zwanzig model has played a crucial role in the development of the H-theorem for strongly nonlinear Fokker-Planck equations (we will return to this issue later).Various studies have addressed at least to some extent the Desai-Zwanzig model. This can be shown by interpreting the Desai-Zwanzig model in the context of the free energy approach to nonlinear Fokker-Planck equations [9]. Accordingly, the Desai-Zwanzig model does not only involve the mean value M of the process but also the process variance K=〈X2〉-M2. In order to see this, we need to take advantage of variational calculus. Let δK/δP denote the variational derivative of K with respect to the probability density P(x). A detailed calculation shows that
(81)M=∫-∞∞xP(x)dx,K=∫-∞∞x2P(x)dx-M2⟹δKδP=x2-2xM⟹12∂∂xδKδP=x-M.
Exploiting this result, the free energy F can be defined like
(82)F[P]=∫-∞∞V(x)P(x)dx+κ2K[P]-DS[P],S[P]=-∫-∞∞P(x)ln(P(x))dx,
where S denotes the Boltzmann-Gibbs-Shannon entropy of the probability density P(x). The potential V(x) denotes the double-well potential V(x)=-Ax2/2+Bx4/4. The strongly nonlinear Fokker-Planck equation (80) can then equivalently be expressed by [9, 194, 196]
(83)∂∂tP(x,t)=∂∂x[P(x,t)∂∂xδFδP],
where δF/δP is the variational derivative of the free energy F with respect to the probability density P(x,t). According to (82) and (83) the Desai-Zwanzig model involves the variance K in the free energy. It has been suggested to call strongly nonlinear stochastic processes “variance models” or “K models” provided that they involve the variational derivatives of the variance (i.e., they involve a coupling to the process mean value). A minireview of the Desai-Zwanzig models and related “variance models” or “K models” can be found in Frank [9]. Finally, note that the Shimizu-Yamada model discussed in Section 6.17.1, that is, the generalized Ornstein-Uhlenbeck model (33), belongs to the class of “variance models” as well.


### 6.18. Shiino’s H-Theorem for Nonlinear Fokker-Planck Equations

As mentioned in Section6.17.2, the Desai-Zwanzig model played a crucial role in the theory development of the H-theorem for strongly nonlinear Fokker-Planck equations. While it was well known that stochastic processes described by ordinary Fokker-Planck equations under appropriate circumstances converge to stationary processes in the long time limit, the situation in this regard was not clear for strongly nonlinear Markov diffusion processes described by strongly nonlinear Fokker-Planck equations. In physics, the proof of convergence to the stationary case is typically given in terms of a so-called H-theorem. In a seminal work, Shiino [196] provided an H-theorem for the Desai-Zwanzig model. This H-theorem was generalized and discussed in various studies [122, 157–164, 173, 197–201]; see also our comments in Section 6.8 and 6.10 for H-theorems related to the Plastino-Plastino model and the Kuramoto model, respectively. A selection of H-theorems can be found in Frank [9]. Finally, note that one can show that not only the single time point probability density P(x,t) of a strongly nonlinear stochastic process converges to a stationary probability density. But, the whole process becomes stationary as well [202].In the context of Shiino’s work on the H-theorem, we would like to point out that the study by Shiino [196] also established a novel approach to conduct a stability analysis for nonlinear Fokker-Planck equation. A minireview of the stability analysis by Shiino can be found in Frank [9].


## 6.1. Chaos in Probability: Time-Continuous Case

In Section3.2 strongly nonlinear time-discrete stochastic processes were considered that exhibit occupation probabilities that evolve in a chaotic fashion. As an example, the logistic map model defined by (39) was discussed. A key feature of this model and of similar models is that the chaotic behavior arises from the coupling of the dynamics of the realizations to the occupation probabilities. In particular, without this coupling equation (39) reads p(A,t+1)=T0=α/4 with fixed parameters α taken from the interval [0,4] and clearly does not exhibit a chaotic solution. Chaos is due to the fact that transition probability T0=α/4 is replaced by a transition probability that depends on the process occupation probabilities like T0=α·p(A,t)(1-p(A,t)). That is, probability measures of strongly nonlinear processes can exhibit chaotic behavior due to the very nature of strongly nonlinear processes, which is the influence of a probability measure on realizations from which the probability measure is calculated. Such a circular causality structure may be realized in many-body systems composed of interacting subsystems. In such systems, chaos may emerge due to the coupling of the single subsystem dynamics to the relative frequency with which states are occupied by subsystems. Chaos emerges even if the isolated subsystem dynamics is not chaotic [18]. Consequently, chaos may be considered as a self-organization phenomenon. Recall that self-organizing phenomena in general are emerging properties of many-body systems that do not exist on the level of the subsystems [62]. In our context, the many-body system as a whole behaves qualitatively differently from the individual, isolated parts. The catchphrases “the whole is different from the sum of its parts” or “more is different” apply [63].This feature of strongly nonlinear stochastic processes has also been demonstrated for time-continuous models involving continuous state spaces [64–67]. Subsystems described by ordinary stochastic differential equations that do not exhibit chaotic solutions have been coupled together and the limiting case of a many-body system composed of infinitely many subsystems has been considered. The limiting case models were described by strongly nonlinear stochastic processes in the sense defined in Section 1.2. In particular, in these studies the dynamics of realization was coupled to certain expectation values of the respective processes. It was found that the expectation values under appropriate conditions exhibit a chaotic dynamics.


## 6.2. Plasma Physics and the Landau Form of Strongly Nonlinear Fokker-Planck Equations

Plasma physics is concerned with the evolution of many-particle systems composed of charged particles. The particles can interact with each other in various ways. Two interaction forms are particle-particle collisions and Coulomb interaction. The Landau form of the Fokker-Planck equation accounts for these two interactions and describes the evolution of the particle density in the three-dimensional velocity space. That is, letv=(v1,v2,v3) denote the velocity vector of a particle. Let ρ(v,t) denote the particle mass density at time t. Assuming that particles are neither created nor destroyed, the total mass M is invariant over time. The particle probability density P(v,t)=ρ(v,t)/M can be introduced. According to the Landau form, the probability density P(v,t) satisfies a strongly nonlinear Fokker-Planck equation of the form (29). More precisely, the Landau form reads [68–74]
(69)∂∂tP(v,t)=-∑k=13∂∂vk[hk(v,θv[P(·,t)])P(v,t)]+∑j=1,k=13∂2∂vj∂vk[Djk(·)P(v,t)],hk(v,θv[P(·,t)])=A∂∂vk∫ΩP(v′,t)|v-v′|d3v′,Djk(·)=B∂2∂vj∂vk∫Ω|v-v′|P(v′,t)d3v′.
Stationary solutions P(v,st) of (69) have been studied, for example, by Soler et al. [75]. In addition, several studies have been concerned with the derivation of transient solutions P(v,t) using different (semi)numerical approaches [74, 76–79].


## 6.3. Astrophysics: Stellar Cut-Off Mass Distributions Defined by Strongly Nonlinear Stochastic Models

Astrophysical mass distributions may exhibit relative sharp boundaries. More precisely, particle distributions in the stellar space may be spatially bounded due to gravity. The gravitational forces are produced by the mass distributions themselves. That is, the dynamics of an individual particle of a mass distribution is affected by the particle distribution as a whole. In turn, the distribution is composed of its individual particles. In view of this circular causality, it does not come as a surprise that strongly nonlinear stochastic models have been investigated that describe cut-off distributions of stellar particles [80–84]. The models typically assume the form of the strongly nonlinear Fokker-Planck equation (29).Two more comments may be appropriate. First, in Section5.3 we addressed cut-off distributions in the context of physiological signals. However, the motivation in Section 5.3 for considering cut-off distributions was quite different from the argument used in astrophysical applications. Second, in addition to the aforementioned studies on strongly nonlinear processes describing astrophysical cut-off mass distribution, strongly nonlinear stochastic processes satisfying the Landau form (69) have also been used for studying astrophysical problems [85, 86].


## 6.4. Accelerator Physics: Shortening of Bunch-Particle Distribution of Particle Beams

Particles in accelerator beams can travel in localized cluster of particles. In a longitudinal beam, the so-called bunch-particle distribution is a mass density that describes the distribution of particles with respect to the center of the cluster. Letx1 denote relative position of a particle with respect to the bunch center. Typically, the particle dynamics also depends on the particle energy. That is, beam particle dynamics involves a second coordinate x2 which is a rescaled energy measure (rather than a velocity or momentum variable). Dividing the mass density with the total mass, the bunch-particle distribution becomes a probability density P(x,t) of the two-dimensional vector x=(x1,x2). Particles in accelerator beams are charged particles that are accelerated by means of electromagnetic forces. Although particles may interact with each other by means of Coulomb forces, the crucial force that determines the shape and length of a cluster is the so-called wake-field [87, 88]. The wake-field is an electromagnetic field generated by the whole bunch of particles that travels ahead of the bunch but is reflected at the accelerator walls and then acts back on the particles of the bunch. The wake-field can cause a shortening of the particle bunch. That is, under the impact of the wake-field the cluster is smaller in size than it would be theoretically in the absence of the wake-field. The understanding of bunch shortening is an important step towards an understanding of beam instabilities that should be avoided.The dynamics of bunch particles is typically described by strongly nonlinear Markov diffusion processes. The reason for this is that the dynamics of individual particles is affected by the wake-field which can be calculated from an appropriate expectation value of the bunch particle probability densityP(x,t). As a result, in various studies it has been proposed that P(x,t) satisfies a strongly nonlinear Fokker-Planck equation of the form (29) (see, e.g., [89–98]). A useful model, in particular for the purpose of theory development, is the Haissinski model [95, 96, 99–104] for which analytical solutions of the reduced density P(x1) can be derived. In particular, according to the Haissinski model, the dynamics of single particles satisfies a strongly nonlinear Langevin equation (27), which reads [100]
(70)ddtX1(t)=X2(t)ddtX2(t)=h(X(t),θX1(t)[P(x,t)])+gΓ2(t)h=-X1(t)-γX2(t)-∂∂X1∫ΩUW(X1-x1′)P(x1′,x2,t)dx1′dx2,
where γ>0 describes a phenomenological damping constant. For the Haissinski model the wake-field function UW assumes the form of a Dirac delta function: UW(z)=κ·δ(z). The parameter κ measures the strength of the impact of the wake-field. For κ<0 the width of the reduced probability density P(x1) becomes smaller when increasing the magnitude |κ| of the impact of the wake-field. In doing so, model (70) provides a dynamic account for the squeezing effect of wake-fields on particle clusters traveling in accelerator beams.


## 6.5. Physics of Liquid Crystal Phase Transitions from the Perspective of Strongly Nonlinear Stochastic Processes

Liquid crystals typically consist of rod-shape macromolecules that interact with each other by means of weak Van-der-Waals forces. Due to this interaction, molecules can align themselves such that the molecules point with their primary axes more or less in the same direction. The physical properties of a crystal with aligned molecules are qualitatively different from the properties that the crystal exhibits when the macromolecules are isotropically (randomly) distributed. Therefore, liquid crystals exhibit two phases: an ordered phase with molecules that are aligned at least to some degree, which is called the nematic phase, and a disorder phase with molecules pointing in random directions, which is called the isotrope phase. A phase transition from the isotropic to the nematic phase occurs when the temperature is decreased below a critical temperature [105–108]. The degree of orientational order can be captured by the Maier-Saupe order parameter [102, 105–111]. For the isotropic phase the order parameter equals zero. In the nematic phase, the Maier-Saupe order parameter is larger than zero.Stochastic models describing the emergence of orientational order (i.e., isotropic-nematic phase transitions of liquid crystals) exploit the notion that the order parameter itself expresses the magnitude of an overall force that acts on individual molecules and tends to align them with the mean orientation of the liquid crystal macromolecules. The model is based conceptually on the notions of circular causality and self-organization: individual molecules are affected by the order parameter and at the same time constitute the order parameter. Hess [112] as well as Doi and Edwards [113] have proposed a strongly nonlinear stochastic process describing the rotational Brownian motion of the molecules under the impact of the (self-generated) order parameter. The Hess-Doi-Edwards model exhibits the structure of a strongly nonlinear Fokker-Planck equation (29) and can be cast into the form of a strongly nonlinear Langevin equation (27) (see, e.g., [114]). Exploiting the concept of strongly nonlinear stochastic processes, the martingale for the Hess-Doi-Edwards model can be defined [60].Transient solutions of the Hess-Doi-Edwards model may be computed from approximate coupled moment equations [115]. Using Lyapunov’s direct method [62, 116], it can be shown that the ordered state exhibits an asymptotically stable probability density [114]. Importantly, since stochastic models provide a dynamic account, it is possible to study response functions of liquid crystals. Transient responses describing the relaxation to equilibrium [117, 118] as well as correlation functions and mean square displacement functions in the stationary case [114] can be determined at least semi-numerically.Beyond the application of the concept of strongly nonlinear stochastic processes to the Hess-Doi-Edwards Fokker-Planck equation, in general, strongly nonlinear stochastic models have turned out to be useful models in liquid polymer physics (see, e.g., [119–121]).


## 6.6. Classical Descriptions of Quantum Systems by means of Strongly Nonlinear Stochastic Processes

The relaxation towards equilibrium of quantum mechanical Fermi and Bose systems can be modeled using a classical approach by means generalized Boltzmann equations that exhibit appropriately modified collision integrals [122–125]. Applying a diffusion approximation to these collision integrals, such quantum mechanical Boltzmann equations reduce to Fokker-Planck equations that are nonlinear with respect to their probability densities [122–124]. In addition to this approach via the Boltzmann equation, alternative derivations of classical quantum mechanical Fokker-Planck equations that are nonlinear with respect to probability densities have been proposed [126–133]. A key property of these models is that they exhibit stationary probability densities that correspond to Fermi or Bose distributions. Interpreting these models in terms of strongly nonlinear Fokker-Planck equations of form (29) allows us not only to consider the evolution of probability densities but also to examine predictions of semiclassical stochastic processes for quantum systems. For example, solutions of trajectories computed from the corresponding strongly nonlinear Langevin equations of form (27) have been studied [126, 127]. In particular, the relaxation of the free electron gas towards its stationary Fermi energy distribution can be determined by computing numerically individual electron trajectories in the relevant energy space [9, 126]. Moreover, martingales for quantum processes within this classical Langevin approach can be defined [60].Finally, note that classical stochastic descriptions of Fermi and Bose systems do not necessarily involve strongly nonlinear stochastic processes. It has been shown that ordinary Fokker-Planck equations with appropriately defined drift functionsh(·) can also exhibit Fermi and Bose distributions as stationary probability densities (see, e.g., [134]).


## 6.7. Material Physics of Porous Media

Gas particles and fluids diffusing through porous media satisfy a nonlinear diffusion equation, frequently called the porous media equation. In one spatial dimension, the porous medium equation for the particle mass densityρ(x,t) reads [3, 9, 46, 47, 135]
(71)∂∂tρ(x,t)=A∂2∂x2[ρ(x,t)]q,
with A>0. Dividing the mass density ρ(x,t) by the total mass M, (71) becomes an evolution equation for a probability density P(x,t) normalized to unity
(72)∂∂tP(x,t)=D∂2∂x2[P(x,t)]q,
with D=A·M(q-1), and assumes the form of the nonlinear Fokker-Planck equation (29). The connection between the porous medium equation and nonlinear master equations has been addressed in Section 5.1 (see the aforementioned). Equations (71) and (72) exhibit exact time-dependent solutions, frequently called Barenblatt-Pattle solutions [136, 137].The parameterq captures the nonlinearity of the equation. For q=1 (linear case) the porous medium equation reduces to the ordinary diffusion equations and the variance of the diffusion process (or the second moment of P(x,t) assuming a zero first moment) increases linearly with time t. In contrast, for q>1/3 and q≠1 (nonlinear case) the variance increases as a nonlinear function of t like tr with r=2/(1+q) as can be shown by various methods [9, 138–140].Interpreting (72) in terms of a strongly nonlinear Fokker-Planck equation (29), trajectories of realizations can be computed from the corresponding strongly nonlinear Langevin equation [138]. For a generalized version of (72) involving a drift force such a Langevin equation approach has been exploited to derive analytical expressions of certain autocorrelation functions [141].The porous medium equation in its form as a nonlinear Fokker-Planck equation (72) plays an important role in the theory development of nonextensive thermostatistics. As it turns out, a particular generalization of (72), the so-called Plastino-Plastino model, exhibits stationary solutions that maximize the nonextensive entropy proposed by Tsallis. We will return to this issue later.


## 6.8. Material Physics and the Liouville Dynamics of Many-Body Systems with Interacting Subsystems

The particle distribution of a many-body system in the overdamped, deterministic case satisfies a Liouville equation. For many-body systems with interacting subsystems the Liouville equation involves an integral term describing the interaction force on a single subsystem. This interaction integral may be evaluated using a diffusion approximation. In doing so, the Liouville equation assumes the form of the nonlinear Fokker-Planck equation (29) when considering the probability density rather than the subsystem density. The nonlinearity occurs in the diffusion term. This approach has been developed in a series of studies by Andrade et al. [142] and Ribeiro et al. [143, 144]. For example, the distribution functions of interacting particles in dusty plasmas, interacting vortices, and interacting polymer chains in complex fluids have been studied within this framework.Note that as mentioned previously the dynamics of the subsystems is deterministic. Nevertheless, the effective model for the subsystem probability density involves a diffusion term. That is, according to the effective, approximative model in terms of a nonlinear Fokker-Planck equation, due to the interaction between the subsystems the many-body system as a whole generates a fluctuating force that acts on the individual subsystems of the many-body system.


## 6.9. Nonextensive Physical Systems and the Plastino-Plastino Model

Tsallis (Tsallis, [54, 55]; see also Abe and Okamoto [56]) introduced in physics a generalized form of the Boltzmann-Gibbs-Shannon entropy, nowadays frequently called the Tsallis entropy. The Tsallis entropy exhibits a parameter q. For q=1 the entropy reduces to the Boltzmann-Gibbs-Shannon entropy capturing properties of extensive systems. For q≠1 the Tsallis entropy describes nonextensive systems.  A. R Plastino and A. Plastino [145] proposed a nonlinear Fokker-Planck equation of the form (29) that may be written like
(73)∂∂tP(x,t)=-∂∂x[h(x)P(x,t)]+D∂2∂x2[P(x,t)]q.
The model describes the probability density P(x,t) of a particle living in a one-dimensional continuous state space given by the whole real line. The term involving a force function h(x) in (73) is linear with respect to the probability density P(x,t). The diffusion term of the model involves the parameter q, which using the notation suggested by (73) corresponds to the parameter q of the Tsallis entropy (see Frank [9]). Note that in the original model a slightly different notation was proposed [145]. A key feature of the Plastino-Plastino model (73) is that stationary probability densities P(x,st) of (73) maximize the Tsallis entropy. For q=1, that is, in the special case in which the Tsallis entropy reduces to the Boltzmann-Gibbs-Shannon entropy, the model is linear with respect to the probability density P(x,t). For q≠1, that is, for the general case that the Tsallis entropy describes a nonextensive system, the diffusion term is nonlinear with respect to P(x,t). In view of these observations, the Plastino-Plastino model establishes a connection between nonextensivity and nonlinearity.The Plastino-Plastino model (73) has been examined in various studies. In particular, it has frequently been pointed out that (73) may be considered as a generalized version of the porous medium model (72) for gas and fluid flows that are subjected to an additional driving force captured by the force function h(x).Exact analytical solutions for the time-dependent probability densityP(x,t) are available in the case of a linear drift force [140, 145, 146]. Trajectories describing the stochastic dynamics of individual particles can be computed when interpreting (73) as a strongly nonlinear Fokker-Planck equation by means of the corresponding strongly nonlinear Langevin equation [138, 141, 147]. In this case, second order statistics measures such as autocorrelation functions can be determined [141] and martingales can be found [60]. Exact time-dependent solutions have been derived for generalized versions of model (73) that account for source terms [148–150]. The Kramers escape rate has been determined for processes described by the Plastino-Plastino model involving an appropriately defined drift force [151]. The multivariate generalization of (73) with [139, 152, 153] and without [129, 154] time-dependent coefficients has been considered. The Plastino-Plastino model (73) with explicit state-dependent diffusion coefficients D(x) has been studied [153, 155]. The model has been generalized for the Sharma-Mittal entropy that involves two parameters and includes the Tsallis entropy as a special case [156]. Interestingly, H-theorems have been found that show that transient solutions P(x,t) of (73) converge to stationary probabilities densities P(x,st) provided that h(x) is chosen appropriately [122, 157–164].


## 6.10. Oscillatory Coupled Nonequilibrium Systems: Canonical-Dissipative Approach

Coupled oscillators with short-range interactions can be studied within the framework of strongly nonlinear stochastic processes [165]. In this study isolated oscillators were modeled using the canonical-dissipative approach [166–168] that allows to exploit the concepts of Hamiltonian and Newtonian dynamics, on the one hand, and to address nonequilibrium systems such as self-excited limit cycle oscillators, on the other hand.


## 6.11. Phase Synchronization and Kuramoto Model: Time-Continuous Case

Synchronization is a phenomenon studied in various disciplines ranging from physics to the social sciences [169]. In the context of strongly nonlinear stochastic processes, we are primarily concerned with the synchronization of a large ensemble of oscillatory subunits. Synchronization can be studied both in the phase space of the oscillators (e.g., the space spanned by position and momentum variables) and in reduced spaces. An example for the latter case was discussed already in Section 4.4: phase synchronization. In fact, we will focus in this section on the phenomenon of phase synchronization.A benchmark model that describes the emergence of phase synchronization in a population of phase oscillators is the Kuramoto model [36]. Accordingly, each oscillator is described by a phase X that evolves on a continuous time scale t and represents a periodic variable. Each oscillator is coupled to all remaining oscillators and the limiting case of a group of infinitely many oscillators is considered. The oscillators as a whole generate an order parameter, which is a complex-valued variable. The complex-valued order parameter has an amplitude, the so-called cluster amplitude A, and an angle, the so-called cluster phase α. Both cluster phase α and cluster amplitude A are measures defined in circular statistics (see, e.g., [36, 170, 171]). Mathematically speaking, for an oscillator phase X defined on the interval [0,2π] the complex order parameter q is defined by the expectation value
(74)q(θ[P(x,t)])=limR→∞1R∑r=1Rexp(iXr(t))=∫x=02πexp(ix)P(x,t)dx=A(t)exp(iα(t)),
where i is the imaginary unit (i.e., the square root of −1) and the cluster phase α is chosen such that A≥0 in any case. Note that both A and α are time-dependent variables.The order parameterq induces a force that acts on the individual phase oscillators. That is, the oscillators interact with each other indirectly via the self-generated order parameter. The population of phase oscillators can be regarded as a statistical ensemble. Accordingly, the order parameter can be computed as an expectation value from the probability density P(x,t), that is, the probability density of the phase of an arbitrarily chosen phase oscillator. Since the order parameter acts back on the realizations (individual phase oscillators), the model satisfies the definition of a strongly nonlinear stochastic process provided in Section 1.2. The strongly nonlinear Langevin equation of the Kuramoto model reads
(75)ddtX(t)=-κA(θ[P(x,t)])sin(X(t)-α(θ[P(x,t)]))+gΓ(t)
and assumes the form (27). In (75), the parameter κ≥0 measures the strength of the force induced by the order parameter q.The uniform probability densityP=1/(2π) describing a desynchronized oscillator population is a stationary solution for any parameter κ because for P=1/(2π) the cluster amplitude A equals zero. However, only for κ<2g2 the uniform distribution is asymptotically stable. For κ>2g2 the uniform distribution is unstable meaning that for any initial condition that slightly deviates from the uniform probability density P=1/(2π) the strongly nonlinear stochastic process will not converge to a stationary process with uniform probability density. Rather, for κ>2g2 the Kuramoto model exhibits stable, stationary unimodal distributions [9, 36]. These unimodal probability densities indicate that the oscillator population is partially synchronized. Moreover, the order parameter amplitude A is larger than zero for these unimodal stationary probability densities. The unimodal probability densities are stable but not asymptotically stable (see, e.g., [9]). A shift of such a stationary probability density along the x-axis transforms a member of the family of stable stationary probability densities into another member of that family. This feature becomes intuitively clear when we realize that the form of (75) is invariant against transformation X→X+s and α→α+s, where s is an arbitrarily small real number. We may use the term “marginally” stable to characterize the stability of these stationary probability densities (see also Section 4.3).The Kuramoto model has been reviewed in Strogatz [172] and Acebrón et al. [61]. In particular, there exists an H-theorem of the Kuramoto model showing that transient probability densities converge to stationary ones [173]. This H-theorem is related to a free energy approach to the Kuramoto model (Frank [174]; see also Frank [9]).A key question in the context of the Kuramoto model concerns the bifurcation from the desynchronized state to the synchronized state for oscillator populations with inhomogeneous eigenfrequencies. In this case, (75) may be written for realizations Xr of the process like
(76)ddtXr(t)=Ur-κA(θ[P(x,t)])sin(Xr(t)-α(θ[P(x,t)]))+gΓr(t),
where Ur is the phase velocity of the rth realization of the process. If we restrict our considerations to symmetric distributions of velocities and to symmetric initial probability densities, then the symmetry property remains conserved during the evolution of the process and the cluster phase can assume only one of the two values α=0 or α=π. Moreover, the order parameter q becomes a real-valued variable. It can be shown that in this case (76) reduces to
(77)ddtXr(t)=Ur-κqsin(Xr(t))+gΓr(t)q=limR→∞1R∑r=1Rcos(Xr(t)).
Comparing (77) with the Daido model (58) discussed in Section 4.4, it becomes at this stage clear that the Daido model (58) can be considered as a time-discrete counterpart of the time-continuous, fully symmetric Kuramoto model (77) with distributed eigenfrequencies. Note that the Daido model exhibits phases defined on the interval [0,1], whereas in the context of the Kuramoto model typically phases defined on the interval [0,2π] are considered.As mentioned previously, reviews for the Kuramoto model are available in the literature and the reader may consult these works for more details ([61, 172]; see also Tass [175] and Section 7.5 in [9]).


## 6.12. Biology: Population Dispersion and Population Dynamics

The spatial distribution of insects forming swarms may be described in terms of cut-off distributions. For example, the insect densityρ(x)=A·[{1-B·|x|}+]2 has been proposed [176], where the coordinate x measures the location of an insect with respect to the center of the swarm and A,B>0. The operator {·}+ corresponds to the identify for positive arguments and equals zero for negative arguments (i.e., {z}+=z for z>0 and {z}+=0 otherwise). Normalizing ρ by the total mass M, a stochastic model for the probability density P(x)=ρ(x)/M needs to explain the emergence of a stationary insect cut-off probability density. In fact, to this end, a model similar to the Plastino-Plastino model (73) with an appropriate drift term was proposed [176, 177]:
(78)∂∂tP(x,t)=-∂∂x[γsgn(x)P(x,t)]+D∂∂x{[P(x,t)]μ∂∂xP(x,t)},
with γ>0. The stationary probability P(x,st)=C·[{1-B·|x|}+]2 is obtained for μ=0.5. However, exponents different from μ=0.5 have also been proposed [178, 179]. Moreover, descriptions of population dynamics in terms of nonlinear Fokker-Planck equations alternative to (78) have been considered in the literature as well ([81, 168]; see also the Shigesada-Teramoto model in Okubo and Levin, [176], Section 5.6).


## 6.13. Social Sciences and Group Decision Making: Time-Continuous Case

Group decision making has already been addressed in Sections3.3 and 4.3 in the context of time-discrete models. In a series of studies, Boster and coworkers developed the so-called linear discrepancy model of group decision making and compared predictions with experimental data [180–182]. In particular, the model accounts for a key phenomenon in the social sciences: the risky shift. A risky shift occurs when people arrive at riskier decisions when discussing a given topic in a group than when making their own decisions individually. According to the linear discrepancy model due to discussions the opinion in favor or against a particular issue shifts by an amount that is proportional to the opinion that an individual holds and the opinion that is presented in the group by another individual. The linear discrepancy model predicts that a risky shift can be observed when the initial distribution (i.e., the prediscussion distribution) of opinions is skewed. Accordingly, during the discussion session, the opinion distribution tends to become more symmetric, which is accompanied with a shift of the mean value. Some of the mathematical properties of the linear discrepancy model have been studied in more detail within the framework of a strongly nonlinear stochastic process [183]. Interestingly, the stationary, symmetric probability densities describing the distribution of opinions among the group members are only stable and not asymptotically stable. In particular, a shift along the opinion axis transforms one member of the family of stable, stationary probability density into another member. In this regard, the group decision making model exhibits the same feature as the Kuramoto model.


## 6.14. Finance and Option Pricing

A stochastic option pricing model has been developed on the basis of the Plastino-Plastino model (73). In particular, the option pricing model exhibits a diffusion term similar to the one of the Plastino-Plastino model [184–186]. A particular benefit of the model is that it features non-Gaussian distributions. This is due to the nonlinearity of the diffusion term (or the nonextensivity of the Tsallis entropy measure involved in the definition of the process). In fact, non-Gaussian distributions frequently fit financial data better than Gaussian distributions.


## 6.15. Economics and Modeling the Dynamics of Target Leverage Ratios

The total capital (firm value) of a company is typically composed of borrowed money and equity. The leverage of a company or the leverage ratio is the ratio of borrowed money relative to equity. A company needs to decide what leverage ratio the company should have. The decision is not necessarily made once and for all. Rather, the desired leverage ratio (i.e., the target leverage ratio) may be updated dynamically depending on the internal and external circumstances. For this reason, the leverage ratios of companies frequently correspond to dynamic variables subjected to fluctuations.Lo and Hui [187] proposed a strongly nonlinear stochastic model for the dynamics of the leverage ratio. The model is based on a model developed earlier by Hui et al. [188] that involved an explicitly time-dependent function. In the strongly nonlinear stochastic model, this function is eliminated, and in this sense the strongly nonlinear stochastic model can be regarded as being closed. In order to achieve this closure, Lo and Hui [187] assumed that the leverage dynamics of a company is affected by the leverage ratios of similar companies. Within the framework of strongly nonlinear stochastic processes this implies that the leverage ratio dynamics of companies depends on an expectation value of the leverage ratio process. Formally, the model by Lo and Hui is closely related to the Shimizu-Yamada model that will be discussed later.


## 6.16. Engineering Sciences: Generators for Noise Sources

In engineering sciences and related disciplines a common problem is to simulate numerically signals of noise sources. These signals can then be used to test the response of engineered systems to random signals emerging from noise sources. A characteristic feature of noise sources is that they are stationary. In view of this property, strongly nonlinear stochastic processes can be defined which for any initial probability density are stationary [147, 189, 190]. For example, the strongly nonlinear Langevin equation
(79)ddtX(t)=-AX(t)+g(X(t),θ[P(x,t)])Γ(t)g=BP(y,t)[C-∫-∞yP(z,t)dz]|y=X(t),
with A,B,C>0, corresponds to a nonlinear Fokker-Planck equation of the form (29) whose right-hand side is by definition equal to zero [9, 147]. Consequently, the Langevin equation defines for any initial probability density P(x,t=0) a stochastic process with a stationary probability density P(x). The parameters A, B, C can be chosen in order to select a desired correlation time of the process.


## 6.17. Further Benchmark Models

### 6.17.1. The Shimizu-Yamada Model

The Shimizu-Yamada model is motivated by research on muscle contraction (Shimizu [191], and Shimizu and Yamada [192]; see also Section 5.5.4 in [9]). The Shimizu-Yamada model corresponds to the generalized Ornstein-Uhlenbeck process that is defined by (33) and involves a mean-field force proportional to the mean value of the process. The Shimizu-Yamada model is highly relevant for the theory development of strongly nonlinear Markov diffusion processes because it can be solved exactly. That is, exact transient solutions can be found and exact solutions for the conditional probability density T(x,t∣y,s,〈X(s)〉) are available. The conditional probability density, in turn, can be used to calculate all correlation functions of interest both in the transient and in the stationary cases. For details, see Frank [9, 193]. Since the model exhibits exact solutions, it has also been used to test the accuracy of numerical algorithms for solving nonlinear Markov diffusion processes  [16].


### 6.17.2. The Desai-Zwanzig Model

The Desai-Zwanzig model involves a mean-field force proportional to the mean value of the process just as the Shimizu-Yamada model. However, the individual, isolated subsystems are assumed to perform an overdamped motion in a double-well potential [194, 195]. The strongly nonlinear Fokker-Planck equation of the Desai-Zwanzig model reads
(80)∂∂tP(x,t)=-∂∂x[h(x,θ[P(x,t)])P(x,t)]+D∂2∂x2P(x,t)h(x,θ[P(x,t)])=Ax-Bx3-κ(x-M(t))M(t)=∫-∞∞xP(x,t)dx
for a stochastic process defined on the whole real line and A,B,κ>0. The term Ax-Bx3in (80) describes the potential force derived from the aforementioned double-well potential. The term -κ(x-M) describes the mean-field force proportional to the mean value of the process. The force attracts the realizations of the process towards the mean value of the process. The parameter κ measures the strength of the mean-field force. The model is symmetric with respect to x=0. In view of this observation it is intuitively clear that the model exhibits a symmetric stationary probability density with M=0 for all parameters κ. It can be shown that for parameters κ smaller than a certain critical value this symmetric probability density is asymptotically stable and is the only stationary probability density. However, when κ exceeds the critical value, the symmetric stationary probability density becomes unstable. Two asymptotically stable, stationary probability densities emerge with mean values M=z and M=-z, where z>0 [9, 194, 196]. The mean value M has typically been considered as order parameter of a many-body system described by the strongly nonlinear stochastic process (80). For M=0 the many-body system exhibits a disordered state. For M≠0 the system exhibits some order. Therefore, the Desai-Zwanzig model describes an order-disorder phase transition.The special importance of the Desai-Zwanzig model is its feature that provides a fundamental stochastic model for an order-disorder phase transition. The Desai-Zwanzig model and the Kuramoto model may be viewed as the two key models in this regard. Both models are able to capture order-disorder phase transitions. While the Kuramoto model is a prototype system for many-body systems with state variables subjected to periodic boundary conditions, the Desai-Zwanzig model is a prototype system for many-body systems featuring state variables satisfying natural boundary conditions. Moreover, the Desai-Zwanzig model has played a crucial role in the development of the H-theorem for strongly nonlinear Fokker-Planck equations (we will return to this issue later).Various studies have addressed at least to some extent the Desai-Zwanzig model. This can be shown by interpreting the Desai-Zwanzig model in the context of the free energy approach to nonlinear Fokker-Planck equations [9]. Accordingly, the Desai-Zwanzig model does not only involve the mean value M of the process but also the process variance K=〈X2〉-M2. In order to see this, we need to take advantage of variational calculus. Let δK/δP denote the variational derivative of K with respect to the probability density P(x). A detailed calculation shows that
(81)M=∫-∞∞xP(x)dx,K=∫-∞∞x2P(x)dx-M2⟹δKδP=x2-2xM⟹12∂∂xδKδP=x-M.
Exploiting this result, the free energy F can be defined like
(82)F[P]=∫-∞∞V(x)P(x)dx+κ2K[P]-DS[P],S[P]=-∫-∞∞P(x)ln(P(x))dx,
where S denotes the Boltzmann-Gibbs-Shannon entropy of the probability density P(x). The potential V(x) denotes the double-well potential V(x)=-Ax2/2+Bx4/4. The strongly nonlinear Fokker-Planck equation (80) can then equivalently be expressed by [9, 194, 196]
(83)∂∂tP(x,t)=∂∂x[P(x,t)∂∂xδFδP],
where δF/δP is the variational derivative of the free energy F with respect to the probability density P(x,t). According to (82) and (83) the Desai-Zwanzig model involves the variance K in the free energy. It has been suggested to call strongly nonlinear stochastic processes “variance models” or “K models” provided that they involve the variational derivatives of the variance (i.e., they involve a coupling to the process mean value). A minireview of the Desai-Zwanzig models and related “variance models” or “K models” can be found in Frank [9]. Finally, note that the Shimizu-Yamada model discussed in Section 6.17.1, that is, the generalized Ornstein-Uhlenbeck model (33), belongs to the class of “variance models” as well.


## 6.17.1. The Shimizu-Yamada Model

The Shimizu-Yamada model is motivated by research on muscle contraction (Shimizu [191], and Shimizu and Yamada [192]; see also Section 5.5.4 in [9]). The Shimizu-Yamada model corresponds to the generalized Ornstein-Uhlenbeck process that is defined by (33) and involves a mean-field force proportional to the mean value of the process. The Shimizu-Yamada model is highly relevant for the theory development of strongly nonlinear Markov diffusion processes because it can be solved exactly. That is, exact transient solutions can be found and exact solutions for the conditional probability density T(x,t∣y,s,〈X(s)〉) are available. The conditional probability density, in turn, can be used to calculate all correlation functions of interest both in the transient and in the stationary cases. For details, see Frank [9, 193]. Since the model exhibits exact solutions, it has also been used to test the accuracy of numerical algorithms for solving nonlinear Markov diffusion processes  [16].


## 6.17.2. The Desai-Zwanzig Model

The Desai-Zwanzig model involves a mean-field force proportional to the mean value of the process just as the Shimizu-Yamada model. However, the individual, isolated subsystems are assumed to perform an overdamped motion in a double-well potential [194, 195]. The strongly nonlinear Fokker-Planck equation of the Desai-Zwanzig model reads
(80)∂∂tP(x,t)=-∂∂x[h(x,θ[P(x,t)])P(x,t)]+D∂2∂x2P(x,t)h(x,θ[P(x,t)])=Ax-Bx3-κ(x-M(t))M(t)=∫-∞∞xP(x,t)dx
for a stochastic process defined on the whole real line and A,B,κ>0. The term Ax-Bx3in (80) describes the potential force derived from the aforementioned double-well potential. The term -κ(x-M) describes the mean-field force proportional to the mean value of the process. The force attracts the realizations of the process towards the mean value of the process. The parameter κ measures the strength of the mean-field force. The model is symmetric with respect to x=0. In view of this observation it is intuitively clear that the model exhibits a symmetric stationary probability density with M=0 for all parameters κ. It can be shown that for parameters κ smaller than a certain critical value this symmetric probability density is asymptotically stable and is the only stationary probability density. However, when κ exceeds the critical value, the symmetric stationary probability density becomes unstable. Two asymptotically stable, stationary probability densities emerge with mean values M=z and M=-z, where z>0 [9, 194, 196]. The mean value M has typically been considered as order parameter of a many-body system described by the strongly nonlinear stochastic process (80). For M=0 the many-body system exhibits a disordered state. For M≠0 the system exhibits some order. Therefore, the Desai-Zwanzig model describes an order-disorder phase transition.The special importance of the Desai-Zwanzig model is its feature that provides a fundamental stochastic model for an order-disorder phase transition. The Desai-Zwanzig model and the Kuramoto model may be viewed as the two key models in this regard. Both models are able to capture order-disorder phase transitions. While the Kuramoto model is a prototype system for many-body systems with state variables subjected to periodic boundary conditions, the Desai-Zwanzig model is a prototype system for many-body systems featuring state variables satisfying natural boundary conditions. Moreover, the Desai-Zwanzig model has played a crucial role in the development of the H-theorem for strongly nonlinear Fokker-Planck equations (we will return to this issue later).Various studies have addressed at least to some extent the Desai-Zwanzig model. This can be shown by interpreting the Desai-Zwanzig model in the context of the free energy approach to nonlinear Fokker-Planck equations [9]. Accordingly, the Desai-Zwanzig model does not only involve the mean value M of the process but also the process variance K=〈X2〉-M2. In order to see this, we need to take advantage of variational calculus. Let δK/δP denote the variational derivative of K with respect to the probability density P(x). A detailed calculation shows that
(81)M=∫-∞∞xP(x)dx,K=∫-∞∞x2P(x)dx-M2⟹δKδP=x2-2xM⟹12∂∂xδKδP=x-M.
Exploiting this result, the free energy F can be defined like
(82)F[P]=∫-∞∞V(x)P(x)dx+κ2K[P]-DS[P],S[P]=-∫-∞∞P(x)ln(P(x))dx,
where S denotes the Boltzmann-Gibbs-Shannon entropy of the probability density P(x). The potential V(x) denotes the double-well potential V(x)=-Ax2/2+Bx4/4. The strongly nonlinear Fokker-Planck equation (80) can then equivalently be expressed by [9, 194, 196]
(83)∂∂tP(x,t)=∂∂x[P(x,t)∂∂xδFδP],
where δF/δP is the variational derivative of the free energy F with respect to the probability density P(x,t). According to (82) and (83) the Desai-Zwanzig model involves the variance K in the free energy. It has been suggested to call strongly nonlinear stochastic processes “variance models” or “K models” provided that they involve the variational derivatives of the variance (i.e., they involve a coupling to the process mean value). A minireview of the Desai-Zwanzig models and related “variance models” or “K models” can be found in Frank [9]. Finally, note that the Shimizu-Yamada model discussed in Section 6.17.1, that is, the generalized Ornstein-Uhlenbeck model (33), belongs to the class of “variance models” as well.


## 6.18. Shiino’s H-Theorem for Nonlinear Fokker-Planck Equations

As mentioned in Section6.17.2, the Desai-Zwanzig model played a crucial role in the theory development of the H-theorem for strongly nonlinear Fokker-Planck equations. While it was well known that stochastic processes described by ordinary Fokker-Planck equations under appropriate circumstances converge to stationary processes in the long time limit, the situation in this regard was not clear for strongly nonlinear Markov diffusion processes described by strongly nonlinear Fokker-Planck equations. In physics, the proof of convergence to the stationary case is typically given in terms of a so-called H-theorem. In a seminal work, Shiino [196] provided an H-theorem for the Desai-Zwanzig model. This H-theorem was generalized and discussed in various studies [122, 157–164, 173, 197–201]; see also our comments in Section 6.8 and 6.10 for H-theorems related to the Plastino-Plastino model and the Kuramoto model, respectively. A selection of H-theorems can be found in Frank [9]. Finally, note that one can show that not only the single time point probability density P(x,t) of a strongly nonlinear stochastic process converges to a stationary probability density. But, the whole process becomes stationary as well [202].In the context of Shiino’s work on the H-theorem, we would like to point out that the study by Shiino [196] also established a novel approach to conduct a stability analysis for nonlinear Fokker-Planck equation. A minireview of the stability analysis by Shiino can be found in Frank [9].


## 7. Mean Field Theory and Strongly Nonlinear Stochastic Processes

While from a mathematical point of view strongly nonlinear stochastic processes are well defined, the question arises to what extent strongly nonlinear stochastic processes describe physical systems. In other words, we may ask what kind of physical systems give rise to strongly nonlinear stochastic processes. An important class of physical systems that has been discussed in this context is the class of many-body systems that can be treated at least approximately with the help of mean-field theory (see, e.g., [36, 175, 195, 203, 204]). The departure point is a many-body system composed of R subsystems. The subsystems are coupled with each other by means of an all-to-all coupling which involves an average over a function W of the state variables of the subsystems. The limiting case R→∞ (the so-called thermodynamic limit) is considered next. In this limiting case, it is assumed that (i) the subsystems become statistically independent of each other and (ii) the average over W becomes an expectation value of W of an arbitrarily chosen representative subsystem. The function W frequently represents a component of a force or corresponds to a force itself. This force is called a mean-field force. A key problem in this regard is to provide a mathematical rigorous proof that in the limiting case R→∞ the many-body system composed of R subsystems can be regarded as a statistical ensemble of realizations. That is, it is often a challenge to prove the convergence of an ordinary multivariate stochastic process to a strongly nonlinear stochastic process. In particular, the mathematical literature is concerned with this problem (see the following).An alternative bottom-up approach to strongly nonlinear stochastic processes uses as departure point a many-body system in the limiting caseR→∞. Again, the subsystems are assumed to be coupled by means of an average over a function W of the subsystem state variables. Furthermore, it is assumed that the subsystems of the many-body system are initially statistically independent and identically distributed. It can then be shown with relatively little effort that if an arbitrary finite subset of size L is considered, then the subsystems of this subset remain statistically independent at all times. The implication is that in the limiting case L→∞ the subset converges to a statistical ensemble and the many-body system can be described by means of a strongly nonlinear stochastic process [9, 202].


## 8. Strongly Nonlinear Stochastic Processes in the Mathematical Literature: Mean-Field Approach, H-Theorems, and Martingales

In the mathematical literature, nonlinear Markov processes have been discussed in the monograph by Kolokoltsov [205]. Besides the seminal work by McKean that opened the avenue to the novel research field on strongly nonlinear stochastic processes (see Sections 1.1 and 1.2), the mathematical literature on strongly nonlinear stochastic processes has traditionally been concerned with the convergence of a multivariate stochastic process to a strongly nonlinear stochastic process in the limiting case in which the system size goes to infinity [7, 195, 206–210]. That is, the argument advocated by mean-field theory presented in Section 7 has been examined in those studies from a mathematical perspective. A second line of inquiry concerns the convergence of processes towards stationarity. That is, the issue of the existence of H-theorems discussed in Section 6.17 has also attracted the attention of researchers in mathematics [197, 211–214]. Furthermore, the existence of martingales for strongly nonlinear Markov processes has been pointed out in the mathematical literature [7, 208, 209, 215–220] and has been taken from there to physics and the life sciences [60].Strongly nonlinear stochastic processes have also been approached from two perspectives slightly different from the mean-field theory approach. Dawson et al. [221] approached strongly nonlinear stochastic processes from the perspective of measure-valued processes, whereas Stroock [222] provided a geometrical approach to the notion of strongly nonlinear Markov processes.Finally, nonlinear Fokker-Planck equations have frequently been addressed in the mathematical literature in the context of semilinear and nonlinear parabolic partial differential equations. In this context, the reader may consult the books by Friedman [223, 224] and Logan [225] and the works by Milstein and colleagues [226–228] on numerical solution methods for nonlinear parabolic partial differential equations.


## 9. Strongly Nonlinear Non-Markovian Processes

The examples reviewed in the previous sections belong to the class of Markov processes. The definition of strongly nonlinear stochastic processes given in Section1.2 however does not imply that strongly nonlinear stochastic processes satisfy the Markov property. In fact, non-Markovian strongly nonlinear stochastic processes can be defined with little effort. For example, while the generalized autoregressive model of order 1 defined by (16) satisfies the Markov property, the generalized autoregressive model of order p defined by
(84)X(t)=∑k=1p(AkX(t-k)+BkM(t-k))+gε(t),M(t)=limR→∞1R∑r=1RXr(t),
with Ap≠0 or Bp≠0, describes a strongly nonlinear stochastic non-Markovian process for p>1. Likewise, strongly nonlinear stochastic non-Markovian processes have been studied evolving on continuous time scales and continuous state spaces [229–231], in particular, in the context of phase synchronization and the Kuramoto model [232] and in the context of muscle contraction and the Shimizu-Yamada model [233].


## 10. Coda

Let us highlight some aspects of strongly nonlinear stochastic processes.Any closed description of an evolution equation for the realization of a process provides a complete and well-defined description of the process because all quantities of interest can at least numerically be computed from (an infinitely larger) ensemble of realizations. Therefore, in Section1.2 strongly nonlinear stochastic processes have been defined on the basis of evolution equations for realizations of stochastic processes. In Section 2 closed descriptions of such evolution equations have been provided. In view of these considerations, we conclude that strongly nonlinear stochastic processes are well defined. That is, it is advocated here that strongly nonlinear stochastic processes are processes worth being studied in their own merit irrespective of underlying physical models or derivations based, for example, on generalized thermostatistics (Section 6.8) or mean-field theory (Sections 7 and 8).The term “strongly nonlinear” should be considered as a placeholder for a variety of other phrases that have been used in the literature to label processes of the kind discussed in this paper. In particular, it should be considered as a synonym to the phrase “truly nonlinear.” Using the phrase “strongly nonlinear,” we may also regard “weakly nonlinear” processes. Such processes exhibit stochastic evolution equations for realizations that do not depend on expectation values or probability densities but feature other types of nonlinearities, for example, nonlinear functions of the relevant state variables. The stochastic Gompertz model and the stochastic logistic model [234] of population growth are two examples of processes that may be called “weakly nonlinear.”From our considerations in Sections2 to 6 it becomes clear that the concept of strongly nonlinear stochastic processes is a fundamental one and does not depend on the explicit representation of time and space. In particular, strongly nonlinear stochastic processes are not restricted to the class of processes originally considered by McKean. That is, strongly nonlinear stochastic processes are not limited to the class of diffusion processes evolving on continuous state spaces and on a continuous time scale.Likewise, as it obvious from our discussion in Section9, the concept of strongly nonlinear stochastic processes and the concept of Markov processes are two independent concepts. A strongly nonlinear stochastic process may or may not exhibit the Markov property. In this context, it should be mentioned that textbook definitions of the Markov property frequently give rise to some confusion because such textbook definitions typically make use of the notion of a “state variable.” However, the key issue concerning the Markov property is the clear separation between future, presence, and past. Accordingly, a process is a Markov process if the future depends only on the present characteristic properties of the process and is independent of the characteristic properties that the process exhibited in the past. Subsequently, we may ask the question about what are the characteristic properties. These properties include the state variables of the process but may also include the occupation probabilities for a process defined on discrete states or the probability density for a process defined on a continuous state space. That is, in the case of a Markov process the present state as well as the present probability density determines entirely the future of a realization of a process. The knowledge about states prior to the present state as well as the knowledge about the probability density at times prior to the present time is irrelevant. This notion is at the heart of the definition of a Markov process [13, 60] and applies to the processes reviewed in Sections 2 to 6. In this context, the fact that in the mathematical literature Markov diffusion processes have been frequently related to martingales (see Section 8) should be highlighted because martingales nicely illustrate the dependency of the future of a process on the present properties of the process and the independency of the future of the process on its past properties.



---

*Source: 149169-2013-03-28.xml*