Information Geometry of Non-Equilibrium Processes in a Bistable System with a Cubic Damping

Abstract

A probabilistic description is essential for understanding the dynamics of stochastic systems far from equilibrium, given uncertainty inherent in the systems. To compare different Probability Density Functions (PDFs), it is extremely useful to quantify the difference among different PDFs by assigning an appropriate metric to probability such that the distance increases with the difference between the two PDFs. This metric structure then provides a key link between stochastic systems and information geometry. For a non-equilibrium process, we define an infinitesimal distance at any time by comparing two PDFs at times infinitesimally apart and sum these distances in time. The total distance along the trajectory of the system quantifies the total number of different states that the system undergoes in time and is called the information length. By using this concept, we investigate the information geometry of non-equilibrium processes involved in disorder-order transitions between the critical and subcritical states in a bistable system. Specifically, we compute time-dependent PDFs, information length, the rate of change in information length, entropy change and Fisher information in disorder-to-order and order-to-disorder transitions and discuss similarities and disparities between the two transitions. In particular, we show that the total information length in order-to-disorder transition is much larger than that in disorder-to-order transition and elucidate the link to the drastically different evolution of entropy in both transitions. We also provide the comparison of the results with those in the case of the transition between the subcritical and supercritical states and discuss implications for fitness.

1. Introduction

The spontaneous emergence of order out of disorder is one of the most fascinating phenomena in nature and laboratory experiments, attracting ever-increasing interest. Important examples include phase transition/critical phenomena in cosmology, elementary particle theory, condensed matter, chemistry, biology and social-economic movement [1,2,3,4,5]. Order is usually quantified by a non-zero value of a macroscopic observable (global mode). While triggered by an external parameter (such as temperature) or spontaneously, a macroscopic observable often does not simply evolve passively, but undergoes an indispensable interaction with fluctuations (microscopic variables). The self-regulation between macroscopic and microscopic variables leads to a dynamical equilibrium (self-organisation), which involves fluctuations as an essential part [6]. There has been accumulating evidence for relevance and important role of self-organisation in different systems such as shear flows or vortices in fluids or plasmas, pattern formation in chemical oscillators, homoeostasis in biosystems and even traffic flows [4,7,8,9,10,11,12,13,14,15,16]. In particular, self-organised shear (zonal) flows are now believed to play a crucial role in stabilising laboratory plasmas, beneficial for extracting fusion energy [7]. Due to large fluctuations involved in order-disorder transition or self-organising systems, it is essential to use statistical tools to describe these systems.

The aim of this paper is to understand order-disorder transition from the perspective of information change associated with transition and uncover a geometrical structure in a statistical space, which can be utilised to understand ever-increasing experimental/observational data. To this end, we investigate a bistable stochastic system that is often invoked as a canonical model of self-regulating systems, e.g., in electric circuits [17], in various cellular processes such as cycles, differentiation and apoptosis, regulation of heart, brain, etc. [18,19,20,21,22,23]. In this model, we calculate time-dependent Probability Density Function (PDF) and the total number of statistically different states that the system undergoes in time. The latter is defined by the dimensionless information length [24,25,26,27,28] (see Appendix A):

$$\mathcal{L}\left(t\right)={\int }_0^t\frac{d{t}_1}{\tau \left(t_1\right)}={\int }_0^{t}d{t}_1\sqrt{\int dx\frac{1}{p(x,t_1)}{\left[\frac{\partial p(x,t_1)}{\partial t_1}\right]}^2},$$

(1)

where $p(x,t)$ is a time-dependent PDF for a stochastic variable x. In Equation (1), $\tau \left(t\right)$ is the time-varying “time-unit”:

$$\begin{array}{ccc}\hfill \mathcal{E}\left(t\right)\equiv \frac{1}{{\left[\tau \left(t\right)\right]}^2}& =& \int dx\frac{1}{p(x,t)}{\left[\frac{\partial p(x,t)}{\partial t}\right]}^2.\hfill \end{array}$$

(2)

$\tau \left(t\right)$ in Equation (2) has dimensions of time and quantifies the correlation time over which the (dimensionless) information changes, thereby serving as the time unit for information change. Note that in equilibrium where $\frac{\partial p}{\partial t}=0$, $\tau \to \infty$. Measuring the total elapsed time in units of $\tau$ between time $t=0$ and t gives the information length in Equation (1). The latter thus establishes a distance between the initial and final PDFs in the statistical space.

We note that our information length is based on Fisher information (cf. [29]) and is a generalisation of statistical distance mainly used in equilibrium or near-equilibrium systems [30,31,32,33,34,35,36,37,38,39] to non-equilibrium systems [24,25,26,27,28]. In particular, the linear geometry of a linear Ornstein-Uhlenbeck (O-U) process was captured by the linear relation ${\mathcal{L}}_{\infty }\propto x_0$, while the power-law geometry of a nonlinear (cubic) stochastic process was revealed by power-law scalings ${\mathcal{L}}_{\infty }\propto x_0^n$ ($n\sim$ 1.5–1.9) [28], where $x_0$ is the peak position of an initial narrow PDF. Furthermore, interesting geodesic solutions were found [27] by time-periodic modulation of the model parameters in an O-U process, which by itself does not support a geodesic solution without the modulation of parameters. As a geodesic is a unique path along which a system undergoes the minimum number of changes in the statistical states given the initial and final conditions, it is beneficial to a system where the change is costly. It is thus important to elucidate the key characteristics of a stochastic system that permits or facilitates the existence of a geodesic in general (without the modulation of model parameters). Finally, we emphasise that Equation (1) can be applied to any data; [25] constructed time-dependent PDFs from MIDI files of music and elucidated the information change in music by $\mathcal{L}$ and $\tau$.

In this paper, in order to gain a key insight, we focus on a zero-dimensional (0D) model, which has one control parameter, and propose an on-quenching experiment by a sudden change of a control parameter from the critical to subcritical and from the subcritical to critical values to trigger disorder-to-order and order-to-disorder transitions, respectively. A pair of disorder-to-order and order-to-disorder transitions with the suitable choice of initial conditions then provides a simple model in which a continuous switch between ordered and disordered states can be studied in great detail. Each pair of disorder-to-order and order-to-disorder transitions models a burst (e.g., a burst in the gene expression consisting of a pair of induction and repression). Since an initial condition represents the “resting” state of a stochastic system in between the two bursts, it is important to understand the effect of different initial conditions on information change. In particular, we aim to elucidate what might be an optimal initial “resting” state that minimizes the information change, sustaining a robust geodesic solution. To this end, we provide detailed comparison of on-quenching processes in this paper with off-quenching processes, where the control parameter changes between subcritical and supercritical, reported in [40]. While some aspects are similar between the two processes, there are also important differences that will be presented and discussed here. The remainder of this paper is organized as follows. Section 2 presents our model, and Section 3 provides details on the time evolution of PDFs. We discuss information length in Section 4 and differential entropy and Fisher information in Section 5. We conclude in Section 6. Appendices contain the derivation of equations used in the main text.

2. Models

We consider the following 0D Ginzburg-Landau model [41] for a stochastic variable x:

$$\frac{dx}{dt}=F\left(x\right)+\xi =-\lambda x-\mu x^3+\xi .$$

(3)

Here, $F\left(x\right)=-\lambda x-\mu x^3$ is a deterministic force; $\xi$ is a white noise with the following statistical property:

$$\langle \xi \left(t\right)\xi \left(t^{\prime }\right)\rangle =2D\delta (t-t^{\prime }),$$

(4)

where D is the strength of the forcing, and the angular brackets denote the average over $\xi$. In our study, $\lambda$ is a control parameter. In the numerical computations, we will fix the value of $\mu$ ($=1$) with no loss of generality, while keeping track of $\mu$ in analytical calculations for clarity. $\lambda$ represents the deviation from the critical value (e.g., $\lambda \propto T-T_c$ for the temperature T where $T_c$ is the critical temperature). That is, our system is subcritical for $\lambda <0$, supercritical for $\lambda >0$ and critical at $\lambda =0$. Equation (3) is the extension of our recent work [28,42] to a cyclic transition. The Fokker–Planck equation [43,44] corresponding to Equations (3) and (4) is as follows:

$$\frac{\partial }{\partial t}p(x,t)=\frac{\partial }{\partial x}\left[-F\left(x\right)+D\frac{\partial }{\partial x}\right]p(x,t).$$

(5)

In this paper, we consider the transition between the critical state and subcritical state by changing $\lambda$ between zero and $-\gamma$ ($\gamma >0$). Here, we note that we are using $\gamma >0$ to explicitly represent a positive growth rate in the subcritical state $\lambda <0$. To this end, we induce the phase transition from the initial state $\lambda =0$ by a sudden change in $\lambda$. This mimics the “on-quenching” experiment where the quenching occurs at the critical state; disorder-to-order transition (order-to-disorder transition) represents the transition from the critical to subcritical (the subcritical to critical) states. This thus contrasts to the case of the transition from the supercritical to subcritical state studied in [40] where the quenching occurs off the critical state.

Specifically, our model is described as follows:

• Forward Process (FP): $\lambda =-\gamma <0$: at $t=0$, a unimodal PDF with a peak at $x=0$, which evolves into a bimodal PDF with peaks at $x=\pm \sqrt{\gamma /\mu }\ne 0$ as $t\to \infty$;
• Backward Process (BP): $\lambda =0$: at $t=0$, a bimodal PDF with peaks at $x=\pm \sqrt{\gamma /\mu }\ne 0$, which evolves into a unimodal PDF with a peak at $x=0$ as $t\to \infty$.

FP and BP have the following stationary distributions $p_F\left(x\right)$ and $p_B\left(x\right)$, respectively:

$$\begin{array}{ccc}\hfill p_F\left(x\right)& \propto & e^{-\frac{\mu }{4D}{(x^2-\frac{\gamma }{\mu })}^2},\hfill \end{array}$$

(6)

$$\begin{array}{ccc}\hfill p_B\left(x\right)& \propto & e^{-\frac{\mu }{4D}{x}^4}.\hfill \end{array}$$

(7)

That is, FP has a stationary bimodal distribution peaked at $x=\pm \sqrt{\frac{\gamma }{\mu }}$, which can be approximated as a double Gaussian for small D as follows:

$$\begin{array}{ccc}\hfill p_F\left(x\right)& \sim & \frac{\sqrt{{\beta }_F}}{2\sqrt{\pi }}\left[e^{-{\beta }_F{\left(x+\sqrt{\frac{\gamma }{\mu }}\right)}^2}+e^{-{\beta }_F{\left(x-\sqrt{\frac{\gamma }{\mu }}\right)}^2}\right],\hfill \end{array}$$

(8)

where ${\beta }_F=\frac{\gamma }{D}$. Equation (8) represents the sum of the two Gaussians (double Gaussian) with the peak at $\pm \sqrt{\frac{\gamma }{\mu }}$ and variance ${\sigma }_F=\langle x^2\rangle -{\langle x\rangle }^2=\frac{D}{2\gamma }=\frac{1}{2{\beta }_F}$. In comparison, BP has a unimodal quartic exponential PDF centred around $x=0$ in equilibrium. To model a cyclic transition between ordered and disordered states, we use $p_B\left(x\right)$ as an initial distribution for FP and $p_F\left(x\right)$ for BP, respectively. Consequently, the initial PDFs in both FP and BP are strongly out of equilibrium. We investigate time-dependent PDFs and information length during this transient relaxation. In particular, we are interested in how $\mathcal{L}$ depends on the deviation from the critical value ($\gamma$) and the strength of the stochastic noise (D).

Table 1. Summary of Forward Process (FP), Backward Process (BP): FP and BP have equilibrium PDFs (in the limit of $t\to \infty$), $p_F$ and $p_B$, respectively. FP and BP start with the initial PDFs, $p(x,0)=p_B$ and $p(x,0)=p_F$, respectively and reach the equilibrium PDFs, $p(x,t\to \infty )=p_F$ and $p_B$, respectively. $p_F$ and $p_B$ have equilibrium variances, $D/2\gamma$ and $2\sqrt{D/\mu }\Gamma \left(\frac{3}{4}\right)/\Gamma \left(\frac{1}{4}\right)$, respectively. ${\sigma }_0$ and ${\sigma }_F$ are the initial and final variances at $t=0$ and $t\to \infty$, respectively, for each process.

Case	$\mathit{\lambda }$	$\mathit{p}(\mathit{x},0)$	$\mathit{p}(\mathit{x},\mathit{t}\to \mathit{\infty })$	${\mathit{\sigma }}_0$	${\mathit{\sigma }}_{\mathit{F}}$
FP	$-\gamma$	$p_B\left(x\right)$	$p_F\left(x\right)$	$2\sqrt{D/\mu }\Gamma \left(\frac{3}{4}\right)/\Gamma \left(\frac{1}{4}\right)$	$D/2\gamma$
BP	0	$p_F\left(x\right)$	$p_B\left(x\right)$	$D/2\gamma$	$2\sqrt{D/\mu }\Gamma \left(\frac{3}{4}\right)/\Gamma \left(\frac{1}{4}\right)$

summarizes the value of $\lambda$ in Equation (3) and initial conditions for FP and BP, together with the variance $\sigma =\langle x^2\rangle -{\langle x\rangle }^2$ of the initial and final PDFs, where the angular brackets denote the average over the stochastic noise $\xi$. We note that for small D, the equilibrium variance for FP, $D/2\gamma$, is much smaller than that for BP, $2\sqrt{D/\mu }\Gamma \left(\frac{3}{4}\right)/\Gamma \left(\frac{1}{4}\right)$.

3. Time-Evolution of PDFs

To solve Equation (5) numerically, we note first that $\gamma$, D and t can always be rescaled such that $\mu =1$, thereby reducing the number of parameters that need to be explored to only $\gamma$ and D. Any numerical solution also requires x to be restricted to a finite interval, which can always be rescaled to $x\in [-1,1]$ without any loss of generality. If $\gamma$ and D are chosen such that $p_F\left(x\right)$ is restricted well away from the boundaries $\left|x\right|=1$, then this finite interval in x is an excellent match to the analytically more convenient infinite extent. Taking $\gamma \le 0.7$ and $D\le {10}^{-3}$ ensures that the bimodal peaks at $x=\pm \sqrt{\gamma }$ are still sufficiently far from the boundaries, and sufficiently narrow, that $p=0$ can simply be imposed as the boundary condition at $x=\pm 1$. The details of the numerical implementation then involve second-order finite-differencing in x and t, using up to $O\left({10}^6\right)$ grid-points in space, and time-steps as small as $O\left({10}^{-7}\right)$. This spatial resolution allows D to be reduced down to ${10}^{-7}$ while still fully resolving the bimodal peaks of width $O\left(D^{1/2}\right)$.

3.1. Overall Comparison of FP and BP

One of the most significant differences between FP and BP is the time scale on which the process evolves and settles in. There are many diagnostic quantities that could be used to quantify this, but a useful one is simply the ratio $\sqrt{\langle x^4\rangle }/\langle x^2\rangle$. This can be evaluated analytically for the two end states, yielding 1.48 for $p_B\left(x\right)$ and 1 for $p_F\left(x\right)$ (taking $D\ll 1$). The evolution in time must therefore be that FP is $1.48\to 1$ and BP is $1\to 1.48$, and the question is on what time scales this happens.

As shown in Figure 1a, for FP the dependence on D is such that every reduction in D by a factor of 100 shifts the curves by a constant amount in t. That is, the ratio does not deviate significantly from 1.48 until a time $c\ln D^{-1}$ has elapsed, but except for this shift the three curves are essentially identical. The numerically determined value of c is 0.355 and is very close to the factor of $\frac{1}{4\gamma }=0.357$ in $t_2$ in Equation (23), discussed later. It is interesting that this value of c is exactly half of that in [40], but otherwise, the scaling with D is the same. The reason for the faster adjustment in this case is because the initial condition already starts out much broader, with a width $O\left(D^{1/4}\right)$ here, as opposed to $O\left(D^{1/2}\right)$ in [40]. We recall that $x=0$ is an unstable equilibrium point when $\xi =0$, and the instability slowly builds up due to $\xi$ (e.g., see [45,46]) and a finite width of the initial PDF until $t\propto O\left(\right|\ln D\left|\right)$. If the initial condition is already broader, then it is not surprising that the instability can develop sooner. However, this factor of half in the settling time does not mean that the evolution of PDFs in on- and off-quenching processes has any similarity. In fact, we will show that they are quite different and that the off-quenching cannot simply be made up of the two phases where effectively $\lambda =\gamma \to 0$ and $\lambda =0\to -\gamma$ (as for the on-quenching case here).

Figure 1b shows that for BP the dependence on D is very different, with every reduction in D by a factor of 100 shifting the curves by a factor of 10; that is, time scales as $D^{-1/2}$. The backward process is initially driven by the movement of the two peaks towards $x=0$ before diffusion becomes crucial in forming a single peak at $x=0$; it is this final diffusive adjustment process that requires an $O\left(D^{-1/2}\right)$ time to happen. This is in sharp contrast with the results in [40], where FP and BP both exhibited the same $c\ln D^{-1}$ scaling (and even with the same value of c). In comparison with the off-quenching in [40], the on-quenching considered here thus has a forward process that is faster by a factor of two, but a backward process that is much slower, with a completely different $D^{-1/2}$ as opposed to $\ln D^{-1}$ scaling.

3.2. PDF of Forward Process

Figure 2 shows the structure of $p(x,t)$ for FP. The particular times are chosen to take the shift $c\ln D^{-1}$ into account; that is, different values of D are shown at the times where they have the same ratios $\sqrt{\langle x^4\rangle }/\langle x^2\rangle$. The results are seen to be identical for the three different values of D. The initial condition obviously does depend on D (as indicated also by the red curves), but once a certain broadening has occurred, in this time frame $c\ln D^{-1}$, the subsequent evolution is independent of D, relying only on the instability process (as measured by $\gamma$). It is only in the last stages of the evolution (not shown in Figure 2), when the solution settles in to the final bimodal structure, that diffusion plays a role again and determines the $O\left(D^{1/2}\right)$ width of the peaks. Comparing with results in [40], it is interesting to note also that here a double-peak emerges essentially immediately, whereas in [40] a finite time had to elapse before the central peak split into two separate peaks. The reason for this is that here the initial condition is at a critical state with a much broader profile than the Gaussian profile at the supercritical state considered in [40].

To understand these results better, it is of value to perform analytical analysis in the limiting cases. To this end, we transform the nonlinear term in Equation (3) into a linear (anti)damping term [47] by seeking a variable y such that Equation (3) becomes $dy/dt=\gamma y+\xi F\left(y\right)$ where $F\left(y\right)$ is a function of y. We then easily show that $dy/dx=\gamma y/(\gamma x-\mu x^3)$ which has the solution $x=y/\sqrt{1+\alpha y^2}$ ($\alpha =\frac{\mu }{\gamma }$). Specifically, y satisfies:

$$\frac{dy}{dt}=\gamma y+\xi {(1+\alpha y^2)}^{\frac{3}{2}}.$$

(9)

Equation (9) provides a convenient way of computing the PDF of x through y by approximating $\xi {(1+\alpha y^2)}^{\frac{3}{2}}\sim \xi$ for small y [45]. Thus, to leading order y is a Gaussian process, simply given by the Ornstein-Uhlenbeck process [44] with a negative damping. The transition probability is thus Gaussian:

$$p(y,t;y_0,0)=\sqrt{\frac{{\beta }_1}{\pi }}{e}^{-{\beta }_1{\left(y-y_0\right)}^2},$$

(10)

where $y_0=\frac{x_0}{\sqrt{1-\alpha x_0^2}}$ and $\frac{1}{{\beta }_1\left(t\right)}=\frac{D}{\gamma }(e^{2\gamma t}-1)$.

By using the conservation of the probability $p(x,t)dx=p(y,t)dy$ and Equation (10), we obtain the transition probability of x as follows:

$$p(x,t;x_0,0)=\frac{1}{{(1-\alpha x^2)}^{\frac{3}{2}}}\sqrt{\frac{{\beta }_1}{\pi }}\exp \left[-{\beta }_1{\left(\frac{x}{\sqrt{1-\alpha x^2}}-\frac{x_0}{\sqrt{1-\alpha x_0^2}}{e}^{\gamma t}\right)}^2\right],$$

(11)

which recovers the previous results [45,46] when $\alpha =1$ and $x_0=0$.

To obtain a (marginal) PDF of x, we recall that FP has an initial PDF given by a quartic exponential:

$$p(x_0,0)=M{\beta }_0^{\frac{1}{4}}{e}^{-{\beta }_0{x}_0^4}.$$

(12)

Here, ${\beta }_0=\frac{\mu }{4D}$; $M=2\Gamma {\left(\frac{1}{4}\right)}^{-1}$ is a normalisation constant where $\Gamma \left(x\right)$ is the Gamma function. Equations (10) and (12) then give:

$$p(x,t)=M{\beta }_0^{\frac{1}{4}}\sqrt{\frac{{\beta }_1}{\pi }}{\int }_{-\infty }^{\infty }d{x}_0\frac{1}{{(1-\alpha x^2)}^{\frac{3}{2}}}\exp \left[-{\beta }_1{\left(\frac{x}{\sqrt{1-\alpha x^2}}-\frac{x_0}{\sqrt{1-\alpha x_0^2}}{e}^{\gamma t}\right)}^2\right]e^{-{\beta }_0{x}_0^4}.$$

(13)

Here, $\frac{1}{{\beta }_1\left(t\right)}=\frac{D}{\gamma }(e^{2\gamma t}-1)$.

We now show that the initial quartic exponential PDF (12) undergoes roughly two stages of the time evolution: stage (i) driven by diffusion/advection with the continuous movement of the PDF peak from $x=0$ towards $\sqrt{\frac{\gamma }{\mu }}$, and then stage (ii) of settling into an equilibrium PDF with the adjustment of the PDF shape. To this end, we examine the behaviour of $p(x,t)$ in Equation (13) for a sufficiently large ${\beta }_1$ such that:

$${\beta }_1=\frac{\gamma }{D(e^{2\gamma t}-1)}\gg 1.$$

(14)

Equation (14) will later be shown to be valid in stage (i) (e.g., for $t<t_2$ where $t_2$ is defined in Equation (23)). By using Equation (14), we can approximate the exponential function in Equation (13):

$$\begin{array}{c}\sqrt{\frac{{\beta }_1}{\pi }}\exp \left[-{\beta }_1{\left(\frac{x}{\sqrt{1-\alpha x^2}}-\frac{x_0}{\sqrt{1-\alpha x_0^2}}{e}^{\gamma t}\right)}^2\right]\sim \delta \left(\frac{x}{\sqrt{1-\alpha x^2}}-\frac{x_0}{\sqrt{1-\alpha x_0^2}}{e}^{\gamma t}\right)\\ =e^{-\gamma t}{(1-\alpha x^2)}^{\frac{3}{2}}{(1-\overline{\alpha }{x}^2)}^{-\frac{3}{2}}\delta \left(x_0-\frac{x{e}^{-\gamma t}}{\sqrt{1-\overline{\alpha }{x}^2}}\right),\end{array}$$

(15)

where:

$$\overline{\alpha }\equiv \alpha (1-e^{-2\gamma t}).$$

(16)

Then, using Equation (15) in Equation (13) gives us:

$$\begin{array}{ccc}\hfill p(x,t)& \sim & M{\beta }_0^{\frac{1}{4}}\frac{e^{-\gamma t}}{{(1-\overline{\alpha }{x}^2)}^{\frac{3}{2}}}{e}^{-{\beta }_0{e}^{-4\gamma t}{\left(\frac{x^2}{1-\overline{\alpha }{x}^2}\right)}^2}\equiv M{\beta }_0^{\frac{1}{4}}{e}^{-\gamma t}{e}^{-\psi },\hfill \end{array}$$

(17)

where $\psi$ is defined by

$$\psi ={\beta }_0{e}^{-4\gamma t}{\left(\frac{x^2}{1-\overline{\alpha }{x}^2}\right)}^2+\frac{3}{2}\ln (1-\overline{\alpha }{x}^2).$$

(18)

By using $\overline{\alpha }=0$ at $t=0$, we can easily show that Equation (17) matches the initial PDF $p(x,0)$ at $t=0$ in Equation (12). From $\frac{\partial \psi }{\partial x}{|}_{x=x_{1,2}}=0$, we find the value of $x=x_{1,2}$ where the PDF takes its local maximum or minimum:

$$x_1=0,4{\beta }_0{e}^{-4\gamma t}{x}_2^2=3\overline{\alpha }{(1-\overline{\alpha }{x}_2^2)}^2.$$

(19)

Since ${\partial }_{xx}\ln p(x,t)=-{\partial }_{xx}\psi >0$ at $x=0$, $x=x_1=0$ is a local minimum of $p(x,t)$ for all $t>0$. On the other hand, two values of $x_2$ (note $\overline{\alpha }{x}_2^2<1$) represent the location of the local maximum in $x>0$ and $x<0$, respectively. Thus, $x=x_1=0$, the (global) maximum at $t=0$, becomes a local minimum for any infinitesimal time $t>0$, two peaks forming at $x_2$. For instance, for $\overline{\alpha }{x}_2^2\ll 1$, we can show that:

$$x_2^2\sim \frac{3D}{\gamma }\left(e^{4\gamma t}-e^{2\gamma t}\right).$$

(20)

This reveals the diffusive nature of the peak movement from $x=0$ due to instability $\gamma$ towards the equilibrium value $\pm \sqrt{\frac{\gamma }{\mu }}$. In the limit of a very small time $t\ll \frac{1}{4\gamma }$, Equation (20) gives $x_2\sim \pm \sqrt{6Dt}$ by using $e^{2\gamma t}=1+2\gamma t+\dots$ and ${\beta }_0=\frac{\mu }{4D}$, showing that the initial movement of the two peaks is via random walk. Our numerical solutions confirmed the predicted scaling of $x_2\propto D^{1/2}$ in Equation (20), as well as $x_2\propto \sqrt{Dt}$ for small time, followed by almost exponential increase (no figure is shown here).

To examine the evolution in more detail, we consider the characteristic time $t_2$ where the width of the PDF in Equation (17) becomes comparable to the peak position $x_2$ in Equation (20). To estimate the PDF width R in Equation (17), we use the variance of the quartic exponential function (e.g., [42]):

$$R^2=\frac{\Gamma \left(\frac{3}{4}\right)}{\Gamma \left(\frac{1}{4}\right)}\sqrt{\frac{1}{{\beta }_0{e}^{-4\gamma t}}},$$

(21)

for $\overline{\alpha }{x}^2\ll 1$. By using ${\beta }_0=\mu /4D$, we simplify Equation (21) as:

$$R^2\sim 3e^{2\gamma t}{\left(\frac{4D}{\mu }\right)}^{1/2}.$$

(22)

By equating R in Equation (22) and $x_2$ in Equation (20) and using $e^{4\gamma t}-e^{2\gamma t}\sim e^{4\gamma t}$ for $\gamma t\gg 1$, we find the characteristic time $t_2$ as follows:

$$t_2\sim \frac{1}{4\gamma }\ln \left(\frac{4{\gamma }^2}{\mu D}\right),$$

(23)

where ${\beta }_0=\frac{\mu }{4D}$ was used. For $\gamma =0.7$ and $D={10}^{-3},{10}^{-5},{10}^{-7}$, $t_2=2.7,4.35,6.0$. Notably, the value of $t_2$ will later be shown to be very close to the other time scales $t_m$ signifying order formation. $t_2$ in Equation (23) marks the time when the peak position becomes comparable to the rms (Gaussian) fluctuation. For $t\gg t_2$, PDF in settling into the final equilibrium PDF is approximated by the Gaussian [45]. We have confirmed this prediction from our numerical solutions (as discussed in more detail later). Finally, we have checked that Equation (14) is valid for $t_2$ given in Equation (23).

3.3. PDF of Backward Process

We recall that BP starts with a bimodal PDF $p_F$, which has two peaks at $\pm \sqrt{\gamma /\mu }$, which is the final equilibrium PDF of FP. For sufficiently small D, the distance between these two peaks is much larger than the width of the PDFs and are thus well separated so that PDF is approximated as the sum of the two Gaussian (double Gaussian) PDFs given by Equation (8). The latter evolve almost independently in $x>0$ and $x<0$, respectively, until $t\sim O\left(1\right)$ when PDFs undergo significant change in the shape with large fluctuation. Since the Gaussian evolution is completely determined by mean value and variance, we now compute mean value and variance in $x>0$ or $x<0$ separately by taking advantage of small fluctuations compared to mean value. Specifically, we let $x=z+\delta x$ where $z=\langle x\rangle$ is the mean component averaged over $\xi$ and the initial PDF in $x>0$ (or $x<0$) while $\delta x$ is the fluctuation $\langle \delta x\rangle =0$. This gives us:

$$\begin{array}{ccc}\hfill \frac{dz}{dt}& =& -\mu z^3-3\mu \langle {\left(\delta x\right)}^2\rangle \sim -\mu z^3,\hfill \end{array}$$

(24)

$$\begin{array}{ccc}\hfill \frac{d\sigma }{dt}& =& -6\mu z^2\sigma +2D.\hfill \end{array}$$

(25)

The solutions to Equations (24) and (25) are as follows (c.f. [48]):

$$\begin{array}{ccc}z\left(t\right)& =& \frac{z_0}{{(1+2\mu z_0^{2}t)}^{1/2}},\\ \sigma \left(t\right)& =& \frac{{\sigma }_0}{G^{\prime }\left(t\right)}+2D\frac{G\left(t\right)}{G^{\prime }\left(t\right)},\end{array}$$

(26)

where $G^{\prime }\left(t\right)={(1+2\mu z_0^{2}t)}^3$ and $G={\int }_0^t{G}^{\prime }\left(t\right)dt$, $z_0=\sqrt{\frac{\gamma }{\mu }}$ and ${\sigma }_0=\frac{D}{2\gamma }$ (see Equation (8) and Table 1). Equations (25) and (26) will be used for computing $\mathcal{L}$ in the next section.

We show $p(x,t)$ for BP in Figure 3. We choose the particular times again to take the $D^{-1/2}$ scaling into account, and show results at the times where they have the same ratios $\sqrt{\langle x^4\rangle }/\langle x^2\rangle$. The initial evolution (not shown in Figure 3) consists simply of a motion toward the origin, with the width of the peaks remaining $O\left(D^{1/2}\right)$. The scaling of this movement is found to be consistent with the prediction in Equation (26).

Once the peaks get within a distance $D^{1/4}$ of the origin they start to sense the presence of the potential well, and diffusion starts to collapse them to a single peak. As seen in Figure 3, if x is rescaled as $D^{1/4}$, and p correspondingly rescaled as $D^{-1/4}$, then the results again look the same for all three values of D. This final diffusive adjustment to the single central peak is very slow though, resulting in the $D^{-1/2}$ scaling in time.

3.4. Energy Diagnostics

We now elucidate the role of the linear growth term (positive feedback) and cubic damping (negative feedback) in FP in energy balance and geodesic. To this end, we multiply Equation (3) by x and take the average over $\xi$ and initial condition to obtain the following equation:

$$\frac{1}{2}\frac{d\langle x^2\rangle }{dt}=\gamma \langle x^2\rangle -\mu \langle x^4\rangle +D.$$

(27)

Here, the last term D, representing the rate of energy injection by $\xi$, was calculated as $\langle \xi \left(t\right)x\left(t\right)\rangle =\langle \xi \left(t\right){\int }_0^{t}d{t}_1[\gamma x\left(t_1\right)-\mu x{\left(t_1\right)}^3+\xi \left(t_1\right)]\rangle =D$ (also confirmed by the numerical calculations). The middle term $\gamma \langle x^2\rangle -\mu \langle x^4\rangle \equiv H$ represents the energy into the system or environment, depending on the sign. When $H>0$, the energy goes into the system, contributing to the increase in $\langle x^2\rangle$; when $H<0$, the energy is dissipated in the system, increasing heat in the environment.

Figure 4a shows $H=\gamma \langle x^2\rangle -\mu \langle x^4\rangle$ for $D={10}^{-3},{10}^{-5},{10}^{-7}$. Unlike $\langle x^2\rangle$ and $\langle x^4\rangle$, which each increase monotonically in time, H reaches a peak at some time $t=t_m$, and then decreases to the negative value $-D$ in settling in to the equilibrium PDF. The maximum H signifies when the positive feedback by the linear growth rate most dominates over the negative feedback by the nonlinear damping. It is notable that the times $t_m=2.6,4.25,5.9$ in Figure 4a for the maximum H are similar to the times $t_2=2.7,4.35,6.0$ given in Equation (23), with both exhibiting the same $c\ln D^{-1}$ ($c=1/4\gamma$) scaling. $t_m$ will also be shown to be very close to the time for the maximum entropy in Section 5.

Physically, $t_m\sim t_2$ signifies the start of order formation. Another diagnostic for the latter is $\Delta =\sqrt{\langle x^4\rangle }-\langle x^2\rangle$, also shown in Figure 4b, where similar non-monotonic behaviour is prominent, with $\Delta$ peaking at the same times as H. This large fluctuation $\Delta$ signifies the phase transition from disordered to ordered states due to the development of the two peaks, which occurs on time scales increasing with $c\ln D^{-1}$ as discussed above.

For BP, $H=-\langle x^4\rangle$ and $\Delta$ in Figure 4c,d are monotonic during the return to the disordered state. The monotonic evolution of H and $\Delta$ for BP is also reflected in the evolution of the differential entropy in Section 5.

4. Information Length

We calculate information length in Equation (1) and explore geometric structure during phase transition. Figure 5 and Figure 6, for FP and BP respectively, show how $\mathcal{E}$ and $\mathcal{L}$ evolve in time, as well as how the total $\mathcal{L}(t\to \infty )={\mathcal{L}}_{\infty }$ depends on $\gamma$. Since FP and BP switch between $\lambda =0$ and $-\gamma \ne 0$, $\gamma$ in Figure 5 and Figure 6 always refer to the non-zero value. We are especially interested also in comparing the on-quenching results computed here with the previous off-quenching results from [40], shown as the dashed lines in Figure 5 and Figure 6.

4.1. Forward Process

For FP, it is useful to consider times less than or greater than $t_2$ in Equation (23) separately, by approximating the time-dependent PDF as a quartic exponential and Gaussian in $t<t_2$ and $t>t_2$, respectively. First, for $t<t_2$, by ignoring the contribution from the mean value $\langle y\rangle =z$ compared with that from the variance, we obtain $\tau \left(t\right)$ in Equation (2) (see Appendix B):

$$\begin{array}{ccc}\hfill \frac{1}{{\left[\tau \left(t\right)\right]}^2}& \sim & \frac{1}{4\beta {\left(t\right)}^2}{\left(\frac{d\beta }{dt}\right)}^2.\hfill \end{array}$$

(28)

Equations (2), (23) and (28) then give $\mathcal{L}$ between the time $t=0$ and $t\sim t_2$:

$$\mathcal{L}\left(t_2\right)\sim \frac{1}{2}\left|\ln \left(\frac{{\beta }_0{e}^{-4\gamma t_2}}{{\beta }_0}\right)\right|\sim \frac{1}{2}\left|\ln \left(\frac{3\mu D}{{\gamma }^2}\right)\right|=\frac{1}{2}\ln \left(\frac{{\gamma }^2}{3\mu D}\right).$$

(29)

On the other hand, during the time between $t\sim t_2$ and $t\to \infty$, the PDF settles into the double Gaussians so that we can estimate the total $\mathcal{L}$ between $t\sim t_2$ and $t\to \infty$ by using Equation (A13) (see Appendix C) as:

$$\begin{array}{ccc}\hfill {\mathcal{L}}_{\infty }-\mathcal{L}(t\sim t_2)& \sim & \frac{1}{\sqrt{2}}\left|\ln \left(\frac{\sigma \left(t_2\right)}{\sigma (t\to \infty )}\right)\right|\sim \frac{1}{\sqrt{2}}\ln \left[\frac{2{\gamma }^2\Gamma \left(\frac{3}{4}\right)}{\sqrt{3}\mu D\Gamma \left(\frac{1}{4}\right)}\right],\hfill \end{array}$$

(30)

where $\sigma \left(t_2\right)=\frac{\Gamma \left(\frac{3}{4}\right)}{\Gamma \left(\frac{1}{4}\right)}\frac{1}{\sqrt{2{\beta }_0{e}^{-4\gamma t_2}}}=\frac{\Gamma \left(\frac{3}{4}\right)}{\Gamma \left(\frac{1}{4}\right)}\sqrt{\frac{4{\gamma }^2}{3{\mu }^2}}$ for FP (see above), $\alpha =\frac{\mu }{\gamma }$ and $\sigma (t\to \infty )=\frac{D}{2\gamma }$ were used. Equations (29) and (30) have the same dependence on D, $\mu$ and $\gamma$. The sum of Equations (29) and (30) gives the total:

$${\mathcal{L}}_{\infty }\sim -1.2+\frac{\sqrt{2}+1}{2}\ln \left(\frac{{\gamma }^2}{D}\right),$$

(31)

when $\mu =1$ and numerical values for the $\Gamma$ functions are inserted.

Figure 5 shows the numerically computed $\mathcal{E}$ and $\mathcal{L}$ for FP. We see how $\mathcal{E}$ starts out essentially constant, corresponding to a geodesic solution [27]. This constant plateau continues until the $O\left(\right|\ln D\left|\right)$ equilibration time scale previously also seen in Figure 1. After this time $\mathcal{E}$ decreases exponentially. Comparing $\mathcal{E}$ here with the previous off-quenching results, we notice three differences: (a) the previous initial adjustment before the plateau regime is absent here, and the curves are essentially flat from the initial condition onward; (b) the plateau here is higher than before; (c) the equilibration and, hence, the exponential decrease in $\mathcal{E}$, happen sooner.

Turning to $\mathcal{L}$ next, the combination that the plateau is higher, but ends sooner, has the interesting consequence that initially $\mathcal{L}$ is greater than in the off-quenching case, but the final values ${\mathcal{L}}_{\infty }$ are always lower. Figure 5c shows the variation of ${\mathcal{L}}_{\infty }$ with $\gamma$, and the same pattern persists throughout; ${\mathcal{L}}_{\infty }$ is consistently $\sim 1$ less than before, with the resulting best-fit formula:

$${\mathcal{L}}_{\infty }\approx 2.1\ln \gamma -1.05\ln D.$$

(32)

The coefficients of $\ln \gamma$ and $\ln D$ are both in generally good agreement with the analytic predictions from Equation (31), which has $2.4$ and $-1.2$. The constant terms, zero versus $-1.2$, match less well, but this term is also strongly affected by the best fit to the $\ln D$ term, since, e.g., $|\ln {10}^{-7}|=16$ is already as large as the largest ${\mathcal{L}}_{\infty }$ values. (Note finally that the deviation from straight lines for large D and small $\gamma$ has the same origin as before in [40]: the “initial” and “final” states are then so broad (large D) and so close to each other (small $\gamma$) that they overlap, causing the dynamics to be different, but also not very interesting in this regime.)

4.2. Backward Process

Figure 6 shows corresponding results for BP. $\mathcal{E}$ now starts off lower than in the off-quenching case, but the final equilibration is much slower, again as seen previously in Figure 1. The result of the initially smaller $\mathcal{E}$ is that for small times $\mathcal{L}$ is a factor of two less than in the off-quenching case. See also Equation (33) below, which confirms this analytically. Because the equilibration is so slow though, there is an additional contribution to ${\mathcal{L}}_{\infty }$ that is not present before. Curiously, this seems to result in the final ${\mathcal{L}}_{\infty }$ values always being a factor of 1.5 less than in the off-quenching case. The precise origin of this particular factor, or indeed why it is always the same, independent of D, is not fully understood. As seen in Figure 6c, the results are summarized by the formula ${\mathcal{L}}_{\infty }\approx 0.9\gamma D^{-1/2}$.

To quantify this scaling, we use Equation (A13) with Equation (A24) and $z\sim 0$:

$$\begin{array}{c}{\mathcal{L}}_{\infty }\sim {\int }_0^{\infty }dt\frac{1}{\sigma }\frac{dz}{dt}\hfill \\ \sim c_1\left(z_0\right){\int }_0^{\infty }\frac{dt}{\sqrt{{\sigma }_0+2DG}}\hfill \\ \sim \mu z_0^3{\int }_0^{\frac{1}{2\gamma }}\frac{dt}{\sqrt{{\sigma }_0+2Dt}}\sim \frac{\sqrt{3}-1}{\sqrt{2}}\frac{\gamma }{\sqrt{\mu D}},\hfill \end{array}$$

(33)

where ${\sigma }_0=\frac{D}{2\gamma }$ for BP (see Table 1) was used. The variation with $\gamma$ and D is exactly as in the numerical results, whereas the constant factor is an under-estimate, 0.5 versus 0.9. Given that Equation (33) only represents the early-time contribution to $\mathcal{L}$ though, we would expect the true ${\mathcal{L}}_{\infty }$ to be larger.

5. Differential Entropy and Fisher Information

Entropy is most commonly used to describe complexity. In a continuous system, it is given by the (Gibbs) differential entropy (e.g., see [49]) defined by:

$$S(t)=-\int dxp(x,t)\ln p(x,t).$$

(34)

Here, we use units in which the Boltzmann constant $K_B=1$. Unlike the usual entropy, the absolute value of the differential entropy does not have a physical meaning, only the difference between two values of the differential entropy being meaningful.

To elucidate the difference in S between the critical and subcritical states, we use equilibrium PDFs of FP and BP ($p_F$ and $p_B$ in Equations (7) and (8), respectively) and quantify the difference between $S(t=0)$ and $S(t\to \infty )$ in FP and BP. For the equilibrium of FP $p_F$ in Equation (8), we can show that the entropy Equation (34) takes the following form [49]:

$$\begin{array}{ccc}\hfill S_F& =& \frac{1}{2}\left[1+\ln \frac{\pi }{{\beta }_F}\right]+2{\beta }_F{x}_0^2\left[1-erf\left(\sqrt{{\beta }_F}{x}_0\right)\right]-\sqrt{\frac{{\beta }_F}{\pi }}2x_0{e}^{-{\beta }_F{x}_0^2}+\Delta .\hfill \end{array}$$

(35)

Here, $erf\left(x\right)=\frac{2}{\sqrt{\pi }}{\int }_0^{x}du\exp (-u^2)$ is the error function; ${\beta }_F=\gamma /D$; $\Delta$ is a function of ${\beta }_F$ and $x_0$, taking the value $0\le \Delta \le \ln 2$. For a sufficiently narrow PDF with ${\beta }_F{x}_0^2\gg 1$, $\Delta$ takes the maximum value $\ln 2$ (see [49]). Since in this limit ${\beta }_F{x}_0^2\gg 1$, $erf\left(\sqrt{{\beta }_F}{x}_0\right)\to 1$, Equation (35) is simplified as:

$$S_F\sim \frac{1}{2}\left[1+\ln \frac{\pi }{{\beta }_F}\right]+\ln 2=\frac{1}{2}\left[1+\ln \frac{\pi D}{\gamma }\right]+\ln 2.$$

(36)

For small values of D as used in our numerical computations, $S_F$ is negative, signifying a strongly localised PDF.

For BF, for simplicity we use the final equilibrium $p_B$ in Equation (7) or (12) and ${\beta }_B=\frac{\gamma }{2D}$ to obtain the differential entropy $S_B$:

$$S_B=-{\int }_{-\infty }^{\infty }dx{p}_B\ln p_B=\frac{1}{4}\left[\ln \frac{\Gamma \left(\frac{1}{4}\right)}{2}+1+\ln \frac{4D}{\mu }\right].$$

(37)

For small D, $S_F\sim \frac{1}{2}\ln \frac{D\pi }{\gamma }$ while $S_B\sim \ln \frac{4D}{\mu }$. Thus, the difference in differential entropy between $p_F$ and $p_B$ is:

$$\Delta S=S_F-S_B=\frac{1}{4}\ln \frac{D\mu \pi }{{\gamma }^2},$$

(38)

which is negative for small D. That is, the quartic exponential PDF at the critical state has much larger entropy than the bimodal PDF at the subcritical state.

Figure 7 shows the time evolution of (Gibbs) differential entropy defined by Equation (34). As theoretically predicted above, we see much larger difference between the initial and final states compared with the off-quenching case [40]. It is interesting to observe that for the forward process, S takes its maximum values at times 2.4, 4.0, and 5.5, very close to where H took its maximum values, and both broadly following the $t_2$ scaling from Equation (23).

To complement $S\left(t\right)$, we also show in Figure 7 the Fisher information defined by:

$$I\left(t\right)=\int \frac{1}{p}{\left[\frac{\partial p(x,t)}{\partial x}\right]}^{2}dx.$$

(39)

As the Fisher information measures the degree of “order”, increasing as the PDF develops large gradients, it shows the opposite tendency to $S\left(t\right)$, which increases with the degree of “disorder”. In particular, the Fisher information $I\left(t\right)$ in Figure 7b for FP takes the minimum value around $t_m$ where the entropy $S\left(t\right)$ is maximum and starts increasing beyond $t>t_m$. These results thus confirm that $t_m$ marks the start of the formation of order, as noted previously.

6. Conclusions

We investigated information geometry associated with order-to-disorder and disorder-to-order transitions in a 0D Ginzburg–Landau model where the formation (disappearance) of an ordered state is modelled by the transition from a unimodal (bimodal) to bimodal (unimodal) PDF of a stochastic variable x. Our 0D model permitted us to perform a detailed statistical analysis. We considered on-critical quenching with a pair of forward and backward processes FP and BP for disorder-to-order (critical to subcritical) and order-to-disorder (subcritical to critical) transitions, respectively by selecting the initial PDF of FP/BP the same as the final equilibrium PDF of BP/FP. A pair of disorder-to-order and order-to-disorder transitions models a burst, for example, in the gene expression consisting of a pair of induction and repression (e.g., see [50]). In such bistable systems, a continuous switching between ordered and disordered states is often observed, the transition occurring in bursts interspersed by a quiescent period (e.g., see [50]). For our cyclic order-disorder transition, an initial condition represents the “resting” state between the two bursts. We thus paid particular attention to the effect of initial conditions on information change by comparing on-quenching and off-quenching cases.

We showed that FP and BP exhibit strikingly different evolution of time-dependent PDFs during transient relaxation due to non-equilibrium initial PDFs. In particular, FP driven by instability undergoes the broadening of the PDF with a large increase in (anomalous) fluctuations before the transition to the ordered state accompanied by narrowing the PDF width/decreased fluctuation. This large fluctuation essentially facilitates the existence of a geodesic solution in FP. This geodesic solution is a result of the self-regulation between the positive feedback ($\gamma x$) and the negative feedback ($-\mu x^3$), which regulate each other, minimising the information change. In a biological context, this minimal geodesic path could be understood in terms of “fitness” in the growth phase (e.g., gene expression). This suggests that the predator-prey type self-regulation with a nonlinear interaction facilitates a geodesic. In comparison, BP is mainly driven by the macroscopic motion due to the movement of the PDF peak, with much less prominent appearance of a geodesic solution. Specifically, the information length $\mathcal{L}$ was found to be much larger in BP than in FP, scaling as $0.9\gamma D^{-1/2}$ for BP, but only $2.1\ln \left(\gamma D^{-1/2}\right)$ for FP, where D is the strength of an additive stochastic noise with a short correlation time. These results demonstrate a great advantage of $\mathcal{L}$ in revealing different physical processes (diffusion/advection) and the different role of diffusion D in transition.

To elucidate the importance of the initial condition between two bursts in cyclic transition, we summarise the striking differences between on-quenching and off-quenching as follows: (i) for FP, double-peaks emerge essentially immediately in on-quenching compared to their appearance only after a finite time in off-quenching; (ii) for FP, the on-quenching has a equilibration time shorter by a factor of two and information length ${\mathcal{L}}_{\infty }$ slightly less than in off-quenching; (iii) for BP, the equilibration time is much longer in on-quenching than in off-quenching, because the final state is at critical; (iv) for BP, the information length ${\mathcal{L}}_{\infty }$ is nevertheless reduced by a factor of 1.5 than in off-quenching. It is worth noting that from the perspective of a system’s “fitness”, the result (ii) could be advantageous when adjusting to a changing environment is costly, and thus, the minimum total change (measured by ${\mathcal{L}}_{\infty }$) and the minimum equilibration time are beneficial (see below). We highlight that ${\mathcal{L}}_{\infty }$ is a “Lagrangian” measure that quantifies the total change in information content in the system over time. We discuss this further in the following.

We note that our control parameter models the effect of environment (e.g., the temperature of the heat bath, etc.), and thus, a sudden change in the control parameter represents a sudden change in the environment. The time-evolution of PDFs occurs in order for the system to reach a new equilibrium state as the equilibrium state is optimal for the given new parameter (for the new environment). On the other hand, the smaller information length represents the smaller number of different states that a system undergoes to reach this new equilibrium state. Intriguingly, these seem to be closely related to the novel concept in microbial metabolism that states evolve under the trade-off between two principles: optimality under one given condition and minimal adjustment between conditions [51]. That is, when an environment changes, the initial state (optimal for the old environment) should change to the new optimal state (the final equilibrium) by undergoing time-evolution. Additionally, the smaller the information length, the less change in the system in adjustment. Thus, our results suggest that the initial “critical” state would be more advantageous for the system in changing environment.

Finally, in future work, it is planned to extend this model to more realistic cases (e.g., 1D or 2D models, a system of coupled equations, etc.).