Transition to Turbulence in Driven Active Matter

TL;DR

This paper investigates a Lorenz-like model for driven active matter, revealing a transition to chaos via period-doubling bifurcations and identifying a unique route to turbulence in such systems.

Contribution

It introduces the first Lorenz-like model demonstrating a sequence of period-doubling bifurcations leading to turbulence in active matter.

Findings

Passage to chaos through period-doubling bifurcations.

Coexistence of strange attractor and fixed points beyond homoclinic point.

Disappearance of strange attractor at Hopf point, leading to limit cycle.

Abstract

A Lorenz-like model was set up recently, to study the hydrodynamic instabilities in a driven active matter system. This Lorenz model differs from the standard one in that all three equations contain non-linear terms. The additional non-linear term comes from the active matter contribution to the stress tensor. In this work, we investigate the non-linear properties of this Lorenz model both analytically and numerically. The significant feature of the model is the passage to chaos through a complete set of period-doubling bifurcations above the Hopf point for inverse Schmidt numbers above a critical value. Interestingly enough, at these Schmidt numbers a strange attractor and stable fixed points coexist beyond the homoclinic point. At the Hopf point, the strange attractor disappears leaving a high-period periodic orbit. This periodic state becomes the expected limit cycle through a set of bifurcations and then undergoes a sequence of period-doubling bifurcations leading to the formation of a strange attractor. This is the first situation where a Lorenz-like model has shown a set of consecutive period-doubling bifurcations in a physically relevant transition to turbulence.

Equations

$$\dot{X}=\sigma (-X+Y)$$

$$\dot{Y}=-Y+rX-XZ$$

$$\dot{Z}=XY-bZ$$

$${\partial }_t{u}_i+u_j{\partial }_j{u}_i$$

$${\partial }_i{u}_i$$

$${\Sigma }_{ij}=-\zeta \left({\partial }_i\varphi {\partial }_j\varphi -\frac{{\delta }_{ij}}{3}{\left(\nabla \varphi \right)}^2\right)$$

$${\partial }_t\varphi (r,t)+u_j{\partial }_j\varphi (r,t)=D{\nabla }^2\varphi (r,t)$$

$${\varphi }_0(r)={\varphi }_{00}+z\frac{\Delta {\varphi }_0}{L}$$

$$u_0(r)=0$$

$${\partial }_t\delta u_i-\nu {\nabla }^2\delta u_i=-{\partial }_i\delta p-\zeta \frac{\Delta {\varphi }_0}{L}\left(\frac{2}{3}{\partial }_z{\partial }_i\delta \varphi +{\nabla }^2\delta \varphi \right)$$

$${\nabla }^2\delta p=-\frac{5}{3}\zeta \frac{\Delta {\varphi }_0}{L}{\partial }_z{\nabla }^2\delta \varphi$$

$${\nabla }^2\left({\partial }_t\delta u_i-\nu {\nabla }^2\delta u_i\right)=\zeta \frac{\Delta {\varphi }_0}{L}\left({\partial }_z{\partial }_i{\nabla }^2\delta \varphi -{\nabla }^4\delta \varphi {\delta }_{i^3}\right)$$

$$({\partial }_t-D{\nabla }^2)\delta \varphi =-\frac{\Delta {\varphi }_0}{L}\delta w$$

$$({\partial }_t\delta w-\nu {\nabla }^2\delta w)=-\zeta \frac{\Delta {\varphi }_0}{L}{\nabla }_{\perp }^2\delta \varphi$$

$$({\partial }_t-D{\nabla }^2)({\partial }_t-\nu {\nabla }^2)\delta w=\zeta {\left(\frac{\Delta {\varphi }_0}{L}\right)}^2{\nabla }_{\perp }^2\delta w$$

$$\delta w(r,t)=\sum _n{a}_n{e}^{i{k}_{\perp }\cdot r_{\perp }}\sin \left(\frac{n\pi z}{L}\right)e^{{\lambda }_n(k_{\perp })t}$$

$$\left[{\lambda }_n+D\left(\frac{n^2{\pi }^2}{L^2}+k_{\perp }^2\right)\right]\left[{\lambda }_n+\nu \left(\frac{n^2{\pi }^2}{L^2}+k_{\perp }^2\right)\right]$$

$$=-\zeta {\left(\frac{\Delta {\varphi }_0}{L}\right)}^2{k}_{\perp }^2$$

$$|\zeta |{\left(\frac{\Delta {\varphi }_0}{L}\right)}^2=D\nu {\left(\frac{\frac{n^2{\pi }^2}{L^2}+k_{\perp }^2}{k_{\perp }^2}\right)}^2$$

$$|\zeta |{\left(\frac{\Delta {\varphi }_0}{L}\right)}^2\frac{L^2}{D\nu }=\frac{({\pi }^2+k_{\perp }^2{L}^2{)}^2}{L^2{k}_{\perp }^2}$$

$$N_c=4{\pi }^2$$

$$u_z=A(t)\cos \pi x\sin \pi z$$

$$u_x=-A(t)\sin \pi x\cos \pi z$$

$$\delta \varphi =\varphi -{\varphi }_0(r)=B(t)\cos \pi x\sin \pi z+C(t)\sin 2\pi z$$

$$\dot{X}=\sigma (-X+rY+rYZ)$$

$$\dot{Y}=-XZ+X-Y$$

$$\dot{Z}=-2Z+XY$$

$$V(X,Y,Z)=X^2+2\sigma r{Y}^2+\sigma r(Z-3{)}^2$$

$$\dot{V}=2\left[-\sigma X^2-2\sigma r{Y}^2-2\sigma r(Z-\frac{3}{2}{)}^2+\frac{9}{2}\sigma r\right]$$

$$X=Y=Z=0$$

$$X_0=\pm \sqrt{2}{\left[(r-1)\pm \sqrt{(r-1{)}^2+(r-1)}\right]}^{\frac{1}{2}}$$

$$Y_0=\frac{X_0}{1+\frac{X_0^2}{2}}$$

$$Z_0=\frac{X_0{Y}_0}{2}$$

$$\dot{\delta X}=\sigma (-\delta X+r\delta Y)$$

$$\dot{\delta Y}=\delta X-\delta Y$$

$$\dot{\delta Z}=-2\delta Z$$

$$\left|\begin{array}{cc}\lambda +\sigma & -\sigma r\\ -1& \lambda +1\end{array}\right|=0$$

Full text

I Introduction

The route to turbulence in a fluid as some control parameter is varied (e.g. increasing the mean flow velocity, or increasing the density gradient if the instabilities are in a stratified fluid) was first conjectured by Landau [1]. In Landau’s scenario, a steady motion becomes a periodic motion and then the periodic motion is destabilized to yield a quasi-periodic motion with two incommensurable frequencies and subsequently, the number of incommensurate frequencies increases with increasing control parameter and eventually the motion becomes a-periodic. This sequence of events was shown to be improbable by Ruelle, Takens, and Newhouse [2, 3] who argued that any more than three quasi-periodic frequencies would lead to an a-periodic state. In a parallel development, to study the a-periodic flows in a fluid stratified by heating from below, Lorenz [4] decided to use a Galerkin truncation procedure to reduce the coupled (velocity and local temperature) partial differential equations to a set of ordinary differential equations. Because of the nature of the fluid (having viscosity and thermal diffusivity), the system had to be dissipative. The minimum number of modes in such a truncation would have to be three to bypass the Poincare-Bendixon theorem [5]. Choosing the three most relevant modes Lorenz obtained a set of three coupled ordinary differential equations now known as the Lorenz model. The control parameter ($r$) of the model was the temperature gradient while there were two other system parameters. One of them ($b$) describes the geometry of the setup ($b=\frac{8}{3}$ for Lorenz) and the other is the Prandtl number of the fluid. The Lorenz model took the form

[TABLE]

The modes $X$ and $Y$ stand for the convective roll patterns set up by the velocity and temperature fields and $Z$ is the amount of convective heat transfer across the fluid. For very small gradients ($r<1$) the conduction state ($X=Y=Z=0$) is stable. It becomes unstable and steady convection sets in for $r>1$. The steady-state is destabilised by a Hopf bifurcation at $r=r_{c}=\sigma\frac{\sigma+b+3}{\sigma-b-1}$ but a periodic state is not observed for $r>r_{c}$ since the bifurcation is backward. Lorenz found that the a-periodic behavior sets in almost immediately above $r_{c}$ and further, the trajectories although confined in a region, are very sensitive to the initial conditions. The set of points to which all trajectories are attracted is called a strange attractor. The signature of fluid turbulence in the truncated model of Eqs. (1.1a)—(1.1c) was the extreme sensitivity of trajectories to initial conditions—a state of affairs described as chaos. Thus the Lorenz model described a new scenario in which the onset of turbulence followed the sequence: trivial steady-state $\to$ convective steady-state $\to$ unstable periodic state via a Hopf bifurcation $\to$ turbulence.

More than a decade later a new paradigm was discovered in the study of low dimensional return maps [6]—[8]. This was an infinite sequence of bifurcations leading to periodic orbits with higher and higher time periods. If the first Hopf bifurcation produces a periodic state of time period $T$, then the successive states are characterized by $2T,2^{2}T,\dots,2^{n}T,\dots$ with new periods appearing with only a slight change in the control parameter. Soon, the system becomes a-periodic at a critical value of the control parameter and also shows sensitive dependence on initial conditions. An expected pathway to turbulence could then be steady-state $\to$ periodic state of period $T$ $\to$ $2T$ $\to$ $2^{2}T$ $\to$ $\dots$ $\to$ a-periodic state. This sequence, or more correctly, some part of it has been seen in some damped driven Duffing oscillators [9], but never in any physically motivated Lorenz-like model of a system governed by non-linear partial differential equations. In real experiments, one of the earliest observations of a few period doublings is that of Libchaber et al [10].

A slightly different kind of hydrodynamical problem is posed by “active matter” (bacteria swimming in a fluid) [11]—[13]. Here, because of its own energy source, active matter can exert an additional stress in the Navier-Stokes equation. The usual stress tensor for the velocity dynamics of an incompressible fluid is $T_{ij}=-p\delta_{ij}+\eta(v_{i,j}+v_{j,i})$ where $v_{i,j}=\partial_{j}v_{i}$, $p$ is the pressure and $\eta$ is the shear viscosity. The additional contribution to $T_{ij}$ can take different forms [14]—[17], depending on what is being studied and a particular form [18] of this extra term is reminiscent of model H in the different universality classes of dynamic critical phenomena [19]. In Ref. [20] the driven form of active matter was studied where a gradient in the density of active particles was considered in one particular direction. To study the effect of non-linearity in this problem a Lorenz-like model was set up. The crucial difference with the usual Lorenz model of Eqs. (1.1a)—(1.1c) is the existence of an extra non-linear term in Eqn. (1.1a). In Ref. [20] we studied the properties of the steady-state and showed how they were significantly different from that of the standard Lorenz model.

In this work, we concentrate on the nature of the Hopf bifurcation. We find that this version of the Lorenz model has the following series of bifurcations for small Schmidt numbers:

Trivial steady-state $\to$ convective steady-state with two stable fixed points $\to$ a Hopf bifurcation to two symmetrically situated periodic orbits of period $T$ about the previously stable fixed points $\to$ a sequence of period-halving bifurcations leading to one periodic orbit encircling both the unstable fixed points $\to$ a sequence of period-doubling bifurcations leading to chaos.

As far as we know, this is the first situation where a Lorenz-like model has shown a set of consecutive period-doubling bifurcations in a physically relevant scenario. A complete summary of our principal results is shown in Fig. 1.

The layout of the paper is as follows: In Sec. II we study the convective instability of the trivial state in the governing partial differential equations and then set up the relevant Lorenz model and study the non-trivial steady state. In Sec. III, the Hopf bifurcation is examined and the system reduces to the canonical form near the bifurcation point. The linear part is shown to lead to a growing mode as the control parameter is increased. The non-linear part shows a forward bifurcation for $\sigma\gg 1$ but a backward bifurcation for small $\sigma$. In the process, we extend the Krylov-Bogoliubov technique [21, 22] beyond its usual habitat of two-dimensional dynamical systems, the details of which are given in Appendix A. In Sec. IV we present detailed numerical results for the model which show the existence of a complete sequence of period-doubling bifurcations leading to chaos as shown in Fig. 1. We examine closely the region $r_{0}<r<r_{c}$, $r_{0}$ being the point where the homoclinic bifurcation occurs and $r_{c}$ a point slightly above the Hopf bifurcation point $r_{H}$, and find some rather striking features. We conclude with a summary in Sec. V.

II The Lorenz Model for Active Matter

We consider wet active matter in it’s simplest form, by ignoring the orientational degrees of freedom. This model was proposed by Tiribocchi et al. [18] and is the active fluid version of Model H in the Halperin-Hohenberg classification scheme.

The variables of interest are:

The scalar concentration field $\phi(\mathbb{r},t)$ of the active particles (also called swimmers) which is proportional to their density, and 2. 2.

The vector velocity field $\mathbb{u}(\mathbb{r},t)$ which represents the velocity of the fluid medium in which the active matter moves about.

The fluid velocity follows an incompressible Navier-Stokes equation which takes the form

[TABLE]

where, $p(\mathbb{r},t)$ is the pressure field, $\nu$ is the kinematic viscosity and $\Sigma_{ij}$ is the stress tensor generated by the active particles. To the lowest non-trivial order allowed by symmetry considerations, it has the form

[TABLE]

where $\zeta$ is a constant which can be termed the activity coefficient. This additional term in the Navier-Stokes equation is the non-linear Burnett term. The concentration field $\phi(\mathbb{r},t)$ satisfies the usual advection-diffusion equation:

[TABLE]

The above equation differs from the form of Tiribocchi et al. [18] in that we have omitted the Ginzburg-Landau-like terms which they needed to study phase separation. In addition, we have neglected a non-linear term in the concentration current that leads to a term in Eq. (2.3) that is of higher order in the gradients than the convective non-linearity in that equation. Effectively in each equation we have retained the leading non-linearity in a gradient expansion. This will allow us to obtain a generalized Lorenz model at the end of this section.

The setup we envisage consists of two square parallel plates of side $L_{\perp}$, separated by a distance $L$ and placed parallel to the $X-Y$ plane. We work in the limit of $L_{\perp}\gg L$ (infinite aspect ratio).

The active matter gradient is maintained at a constant value $\frac{\Delta\phi_{0}}{L}$ across the plates. In this non equilibrium steady state (NESS), we have

[TABLE]

In this steady state there is a constant particle current $\mathbb{j}=-D\mathbb{\nabla}\phi=-D\frac{\Delta\phi_{0}}{L}\mathbb{\hat{z}}$ across the system. For the motionless state of the fluid, we require the pressure gradient to vanish.

The first question that we ask is whether this NESS is stable against small perturbations. If $\delta u_{i}(\mathbb{r},t)$ and $\delta\phi(\mathbb{r},t)$ are the perturbations around $u_{0}(\mathbb{r})$ and $\phi_{0}(\mathbb{r})$, then to linear order in $\delta u_{i}$ and $\delta\phi$, we have

[TABLE]

Imposing the incompressibility constraint $\partial_{i}\delta u_{i}=0$ from Eq. (2.1b) leads to

[TABLE]

Taking the Laplacian of Eq. (2.5) and using Eq. (2.6) leads to

[TABLE]

The concentration equation yields

[TABLE]

where $\delta w$ is the z-component of $\bm{\delta}\mathbb{u}$. Consequently, in Eq. (2.7) we use only the z-component which can be written as

[TABLE]

where $\nabla^{2}_{\perp}=\nabla^{2}-\frac{\partial^{2}}{\partial z^{2}}$. Eliminating $\delta\phi$ between Eqs. (2.8) and (2.9) we get

[TABLE]

The vertical velocity and its second derivative are taken to vanish at the two plates (stress-free boundary conditions) and because of the infinite extent in the $X-Y$ plane, we consider a periodic solution in the horizontal plane. The solution will have the form

[TABLE]

and inserting in Eq. (2.10), we have for each $n$,

[TABLE]

Instability would occur if for a given $n$ and $k_{\perp}$ it is possible for $\lambda_{n}(k_{\perp})$ to become zero or positive. It is immediately clear that for $\lambda_{n}$ to be zero or positive, we need $\zeta<0$. Hence for large negative values of the activity coefficient, we can have the state given by Eqs. (2.4a) and (2.4b) become unstable. Writing $-\zeta$ as $|\zeta|$, we find the condition of $\lambda_{n}$ becoming zero as

[TABLE]

Clearly, the lowest value of $n$ for which the instability can occur is $n=1$ and the critical value of $|\zeta|$ is

[TABLE]

We define $N\equiv|\zeta|\frac{(\Delta\phi_{0})^{2}}{D\nu}$ as the analogue of the Rayleigh number in fluids. This activity Rayleigh number $N$ has a minimum value when $k_{\perp}L=\pi$ and hence the critical value $N_{c}$ for the ”convective instability” is

[TABLE]

We have, thus, arrived at a situation which is completely analogous to hydrodynamic convection. The pure conduction state when a fluid is subjected to a constant temperature gradient loses stability above a critical Rayleigh number (the critical value is $27\frac{\pi^{4}}{4}$ for the boundary conditions that we have used) and a steady convective state becomes stable. The Rayleigh number is the dimensionless variable $Ra=\frac{\alpha g}{D\nu}\left(\frac{\Delta T}{L}\right)L^{4}$ with the $D$ standing for heat diffusivity (since the advected variable is heat), $\alpha$ is the expansion coefficient and $g$ is the acceleration due to gravity. The analogue of $\Delta\phi_{0}$ is $\Delta T$ and the analogue of $|\zeta|$ is $\alpha g$.

It should be noted that for physical reasons, only positive $\Delta T$ (heated from below) can cause convective instability while in the present scenario, any sign of the gradient will lead to an instability if $N$ exceeds the critical value.

In fluids, as the Rayleigh number is increased beyond a critical value, the steady convection becomes oscillatory as the roll-structure begins to undulate. With further increase of Rayleigh number, the oscillatory state makes a transition to a more complicated time dependence and eventually to a soft turbulence (chaos) followed by a hard turbulence characterized by definite scaling behavior. To see this transition in a tractable model, Lorenz carried out a Galerkin truncation of the hydrodynamic equations choosing the three most relevant modes. The resulting Lorenz model—a set of three ordinary non-linear differential equations has become one of the most studied examples in dynamical system theory for complex behavior. Motivated by Lorenz’s work, we construct a similar dynamical system model for this instability.

The simplest set of variables that we can take is those corresponding to a set of cylindrical rolls with the axis of the cylinder in the $y$-direction. This makes all the fields $\mathbb{u}(\mathbb{r},t)$ and $\phi(\mathbb{r},t)$ independent of the $y$-coordinate. The choice

[TABLE]

is the minimal nontrivial model. As shown in Ref. [20], this choice leads to the Lorenz model for driven active matter in the form

[TABLE]

In the above, $\sigma=\frac{\nu}{D}$ is the inverse Schmidt number and $r=\frac{N}{N_{c}}$ is the control parameter. It should be noted that the inverse Schmidt number in this model will be much greater than the Prandtl number in the usual Lorenz model since the particle diffusion coefficient $D$ is generally much smaller than the heat diffusion coefficient. The difference with the standard Lorenz model lies in the non-linear term in Eq. (2.17a). This term is the active matter contribution.

Using an ansatz, we arrive at the following Lyapunov function for the Lorenz model written down above

[TABLE]

with

[TABLE]

Clearly, $\dot{V}<0$ outside the ellipsoid $\sigma X^{2}+2\sigma rY^{2}+2\sigma r(Z-\frac{3}{2})^{2}=\frac{9}{2}\sigma r$ and hence all trajectories of the active matter Lorenz model remain bounded for all time. In the next section we will study the fixed points of this Lorenz model and their stability.

III Fixed Points and Bifurcations

The fixed points of Eqs. (2.17a) — (2.17c) are

The trivial fixed point

[TABLE] 2. 2.

Two non-trivial fixed points at

[TABLE]

The fixed point (3.1) corresponds to the state where there is no flow and a static concentration gradient exists between the plates. Linearizing about this fixed point, we have the perturbations $\delta X,\delta Y,\delta Z$ follow the dynamics

[TABLE]

One of the eigenvalues (growth rate of perturbation) is clearly $-2$ and the other two follow from

[TABLE]

which yields

[TABLE]

Clearly, this fixed point is stable for $r<1$ and unstable for $r>1$. The zero-velocity state consequently becomes unstable for $r>1$ and one has the steady convection roll state given by the fixed point (3.2). The values of $X_{0}$ and $Y_{0}$ correspond to the amplitudes of the roll-like solution while $Z_{0}$ is the Nusselt number for the problem. It describes the constant rate at which the active matter is transported from one plate to another.

Unlike the Rayleigh-Bernard case, the convective transport here is proportional to $(r-1)^{\frac{1}{2}}$ for $r>1$ which is different from the $r-1$ dependence of the usual Lorenz model. Similarly, for $r\gg 1$, the transport saturates in this case as opposed to the continuous increase in the standard situation.

We now look at the question of stability of fixed point (3.2). Writing $X=X_{0}+\delta X$, $Y=Y_{0}+\delta Y$ and $Z=Z_{0}+\delta Z$ and linearizing in $\delta X,\delta Y$ and $\delta Z$, we arrive at the following cubic for the eigenvalue $\lambda$

[TABLE]

For $r=1+\epsilon$ ($\epsilon\ll 1$ but positive) it is clear that all the three roots are real and negative. As one increases $r$ beyond the critical value of $1$, the stable fixed points change from nodes to spirals and eventually become unstable via a Hopf bifurcation at $r=r_{H}$. At the Hopf point, one root is negative and the other two roots become pure imaginary, i.e. $\pm i\omega_{H}$.

Before discussing what happens at $r=r_{H}$, we need to ask if a homoclinic bifurcation occurs for $1<r<r_{H}$, as in the standard Lorenz model [23]—[25]. To find the homoclinic bifurcation we recall that the origin is an unstable fixed point (for $r>1$) with a one-dimensional unstable manifold and a two-dimensional stable manifold (the eigenvalues are $\lambda_{1,2}$ from (3.5) and $\lambda_{3}=-2$). At the homoclinic point $r_{0}$, an orbit orbit leaves the origin along the unstable eigendirection, loops around the nearest stable fixed point and returns to the origin along the stable manifold. The homoclinic point $r_{0}$ as a function of $\sigma$ has been determined numerically (a typical orbit at $r=r_{0}=5.45$ and $\sigma=10$ is shown in Fig. 2) and have been plotted as the green curve in Fig. 3. The interesting thing to note is that $r_{0}<r_{H}$ for all values of $\sigma$. In fact for large $\sigma$, both $r_{H}$ and $r_{0}$ grow almost linearly with $\sigma$, with $r_{0}$ having a smaller slope as can be seen in Fig. 3. This means that a homoclinic bifurcation always occurs below the Hopf point and there could exist an unstable limit cycle in the region $r_{0}<r<r_{H}$, that is created from the homoclinic bifurcation at $r_{0}$ and vanishes at the Hopf bifurcation at $r_{H}$ (thus making it a sub-critical one). However this is not necessarily true for all $\sigma$, as we shall show soon.

Returning to the Hopf bifurcation point at $r=r_{H}$, we have

[TABLE]

We need to find $r_{H}$ and $\omega_{H}$ as functions of $\sigma$, using the $X_{0},Y_{0},Z_{0}$ given in Eqs. (3.2a)—(3.2c). This is difficult and so we use a large-$\sigma$ approximation since we expect $\sigma$ to be large for this problem. It is easy to check that for $\sigma\gg 1$, $r_{H}\approx\frac{\sigma}{4}$ and hence $r_{H}$ is large which in turn allows us to use asymptotic forms for $X_{0},Y_{0},Z_{0}$ which we write as follows (using $\epsilon=r-1$ as the natural variable)

[TABLE]

Using the above forms of $X_{0},Y_{0}$ and $Z_{0}$ in Eqs. (3.7a) and (3.7b), we obtain for $\sigma\gg 1$

[TABLE]

and

[TABLE]

The accuracy of the approximate $r_{H}$ of Eq. (3.9) can be assessed from Fig. 3 where we show the values of $r_{H}$ obtained from Eq. (3.9) and from a numerical solution of Eqs. (3.7a) and (3.7b) as a function of the inverse Schmidt number $\sigma$. It should be noted that for all $\sigma$, $r_{0}<r_{H}$.

The emergence of a limit cycle with frequency $\omega_{H}$ at $r=r_{H}$ implies that at $r=r_{H}$, the asymptotic solution for $X,Y,Z$ must have the form $(X,Y,Z)=(A,B,C)e^{\pm i\omega_{H}t}$ where $A,B$ and $C$ are constants. The fact that we have an instability at $r=r_{H}$ means that if we consider $r=r_{H}+\Delta r$ with $0<\Delta r\ll r_{H}$, then the amplitudes $A,B,C$ will evolve in time as $e^{\lambda t}$ with $\lambda>0$. We establish this behavior in Appendix A. To end this section, we check whether the growth saturates with increasing time (forward Hopf bifurcation). Mclaughlin and Martin [26] had used a solvability based approach for the original Lorenz model which did show that the bifurcation was backward for $\sigma=10$, but missed the fact that it could be forward for very high $\sigma$. This change with $\sigma$ was captured by the renormalization group treatment of Das et al. [27]. In this work, we show that a Krylov-Bogolyubov approach (generally reserved for two-dimensional dynamical systems) is the most direct route to the amplitude equation, the details of which are provided in Appendix A. Here we simply quote the final result (A.18) for the slow dynamics of the amplitude $A(t)$ of the limit cycle about the fixed point $(X_{0},Y_{0},Z_{0})$ as

[TABLE]

The bifurcation is forward (supercritical) for $\sigma\gg 1$ and Eq. (3.11) shows that for $0<\Delta r\ll r_{H}$, the radius of the limit cycle will grow as $\sqrt{\Delta r}$ independent of $\sigma$. The correction term in Eq. (3.11) has a positive sign which indicates that below a certain critical $\sigma$, the Hopf bifurcation will be backward (subcritical). The case $\sigma=10$ reported in Ref. [20] is consistent with a backward Hopf bifurcation. There, a limit cycle was observed for $r<r_{H}$ provided the initial conditions were sufficiently close to the center of the limit cycle. We note that transitions between different dynamical states have been studied recently in a variety of active matter systems, like fluid flocks with inertia [28], droplet growth in scalar active matter [29] and hydrodynamics of active defects [30]. We now present the detailed numerical results in Sec. IV.

IV Numerical Results and Discussion

Similar to the Lorenz model, a homoclinic bifurcation occurs at $r=r_{0}$ in our system. The bifurcation can be identified by the transition shown in Fig. 4. At the critical value $r=r_{0}$ there is a homoclinic orbit as shown in Fig. 2 where the trajectory comes back to the origin. Using this and elementary iterative techniques, we numerically calculate $r_{0}$ as a function of $\sigma$. The plot is shown in Fig. 5. As can be seen, for large $\sigma$, $r_{0}$ grows almost linearly. Plotting this curve in Fig. 3, we see that $r_{0}<r_{H}$ for all $\sigma$. Hence we always have a homoclinic bifurcation below the Hopf bifurcation.

For low values of $\sigma$, our system is expected to behave like the standard Lorenz model. For $r<1$ the origin is globally stable and the only fixed point. At $r=1$ a pitchfork bifurcation occurs and for $1<r<r_{H}$ the two non-trivial fixed points given by Eq. (3.2) are locally stable and trajectories starting sufficiently close to them spiral down to them as in Fig. 6. However the homoclinic bifurcation at $r_{0}$ leads to the formation of an unstable limit cycle as shown for $\sigma=10$ in Fig. 7 that exists for $r_{0}<r<r_{H}$. Integrating numerically in this region, (using a Runge-Kutta $(4,5)$ scheme with an absolute and relative error tolerance of $10^{-6}$) we get chaotic trajectories coexisting with the spirals and the unstable limit cycle as shown in Fig. 8. Since the limit cycle is unstable, even initial conditions very near it either settle on the fixed points or on the chaotic attractor. The limit cycle grows smaller on increasing $r$ from $r_{0}$ and vanish at $r=r_{H}$. For $r>r_{H}$, there are no periodic orbits or stable fixed points in the phase space. All trajectories end up on the Lorenz-like strange attractor shown in Fig. 9. The same behaviour is shown for $\sigma=12,15,20,$ etc.

For higher $\sigma$ the behaviour of the system is starkly different. From the theoretical work in Sec. III we expect the Hopf bifurcation to be forward. For $r<r_{H}$, the non-trivial fixed points given by Eq. (3.2) are stable. Although a homoclinic bifurcation still occurs in the $r_{0}<r<r_{H}$ region, there is no periodic orbit, unlike the low $\sigma$ case. Trajectories either spiral to the stable fixed point or end up on the Lorenz-like strange attractor as shown in Fig. 10. At $r=r_{H}$ a stable high-period limit cycle is born whose transition to chaos is quite unique and is shown in detail for $\sigma=30$ in Fig. 11. The limit cycle, shown in Fig. 11a, has a time period of $8.94$ which corresponds to $64\frac{2\pi}{\omega_{H}}$, where $\omega_{H}$ is supposed to be the Hopf frequency of Eq. (3.10). As $r$ is increased beyond $r_{H}$, the limit cycle undergoes a succession of period halving bifurcations, as can be seen in Fig. 11b, 11c and 11d, until we have a $\frac{2\pi}{\omega_{H}}$ period limit cycle at some $r$ as shown in Fig. 11e. As $r$ is further increased, this limit cycle starts to double its period as shown in Fig. 11f, 11g, 11h, 11i, 11j and ultimately becomes the Lorenz-like attractor as in Fig. 11k and 11l. Throughout this process, the limit cycle is globally stable and there exists no chaos in the region $r>r_{H}$ except after the Lorenz-like attractor is born. Beyond this point, the strange attractor becomes the only dominant feature of the system. This behaviour can be seen for $\sigma=22,23,25,\dots,30,40,50,\dots$ and even $\sigma=80$

Interestingly, as $\sigma$ is increased the zone of the reverse bifurcations narrows and the transition of the $\frac{2\pi}{\omega_{H}}$-period limit cycle into the strange attractor gets more delayed (occurs at even higher values of $r$) and in the $\sigma\to\infty$ limit, the system shows no chaos since (2.17a) gives

[TABLE]

and we are essentially left with a two-dimensional flow which can never be chaotic.

The critical sigma to separate the low and high cases is calculated numerically to be approximately $21$. For $\sigma>21$ the high-sigma route is followed whereas for $\sigma<21$ the low $\sigma$ behaviour is observed.

All these results are summed up in the illustrated plot in Fig. 1.

V Conclusion

We have considered a setup where a definite gradient of active matter is maintained in a given direction across the fluid medium in which it is suspended. As the gradient is increased the fluid medium develops “convective” instabilities and exhibits an infinite sequence of bifurcations. We have used a low-order Galerkin truncation to capture the dynamics of the system exactly in a manner analogous to that employed by Lorenz [4] for studying atmospheric flows. This new version of the Lorenz model, so obtained, has non-linear terms in all the three equations. We find that this Lorenz model can actually show one of the most desirable routes to turbulence—from a trivial steady-state to a non-trivial one, followed by a Hopf bifurcation to a periodic state and then a sequence of period doubling bifurcations until the flow loses all periodicity and shows the existence of a stable strange attractor. Very surprisingly, there is a small twist—the period doubling bifurcations leading to chaos as the control parameter is increased have a mirror image of bifurcations as the control parameter is decreased and this leads to the coexistence of a strange attractor and stable fixed points on the other side of the Hopf bifurcation point. For inverse Schmidt numbers below a critical value, the bifurcations are identical to what one sees in the traditional Lorenz model.

Appendix A Slow dynamics of limit cycle amplitude

We begin by rewriting our Lorenz equations as a single differential equation. Using the notation $D\equiv\frac{d}{dt}$, we define variables $u,v,w$ as

[TABLE]

and write Eqs. (2.17a)—(2.17c) as

[TABLE]

Combining into one equation, we have

[TABLE]

We also note

[TABLE]

We start exploring the region $r=r_{H}+\Delta r$ by ignoring the non-linear terms in Eq. (A.3) and splitting the linear part as a part at $r=r_{H}$ and a part of $O(\Delta r)$. We get

[TABLE]

Knowing that at $\Delta r=0$, the solution is $Ae^{i\omega_{H}t}$ where $A$ is a constant, we write the solution for $\Delta r\ll r_{H}$ as $u=A(t)e^{i\omega_{H}t}$ where $A(t)$ is a slowly varying function of time which implies $\frac{dA}{dt}$ is $O(\Delta r)$. The first term on the left-hand side of Eq. (A.5) factors as $(D^{2}+\omega_{H}^{2})(D+3+\sigma)$ [see Eqs. (3.6), (3.7a) and (3.7b)] and in the second term the $D$ operator can be replaced by $+i\omega_{H}$. Thus we find,

[TABLE]

In the first term, dropping the second and higher-order derivatives of $A(t)$, we have

[TABLE]

Now using the relations given in Eq. (3.8), we approximate the second term of Eq. (A.6) as

[TABLE]

to the leading order in $\sigma$ using $r_{H}\approx\frac{\sigma}{4}$. Thus, Eq. (A.6) becomes

[TABLE]

which clearly shows that A increases exponentially with time for $\Delta r>0$ and hence there is an instability at $r=r_{H}$.

To complete the story, we need to know the first non-linear contribution to Eq. (A.8) and determine if the growth saturates. To do this we need to go back to Eqs. (A.3) — (A.4b) and from the non-linear term on the right-hand side we extract the part with a time-dependence $e^{i\omega t}$. For this, we work with $u=Ae^{i\omega t}$ ignoring the time-dependence of $A$ in taking the derivatives on the right hand side of Eqs. (A.3) — (A.4b). This is allowed because the non-linear term itself will be $O(A^{3})$ which will make $A^{2}$ of $O(\Delta r)$ if the growth saturates. We first find the lowest order $v$ and $w$ ($v_{0}$ and $w_{0}$) in terms of $u$ from the linear terms in Eqs. (A.4a) and (A.4b). In the $\sigma\gg 1$ approximation, which keeps the algebra simple,

[TABLE]

The next step requires us to find the $O(A^{2})$ terms which we denote by $u_{1},v_{1},w_{1}$ from Eqs. (A.4a), (A.4b) and (A.2a). Since the non-linearity is quadratic, the expected structure is

[TABLE]

The constants $B_{2},B_{3}$ and $C_{2},C_{3}$ can be found from the non-linear part of Eqs. (A.4a) and (A.4b) and $B_{1},C_{1}$ from Eq. (A.2a). To keep the algebra simple, we show the results in the $\sigma\gg 1$ limit only. The derivation will be outlined for $v_{1}$ and the results for $w_{1}$ and $u_{1}$ written down. For this purpose, we note the following large $\sigma$ approximations

[TABLE]

Starting with Eq. (A.4a), we write

[TABLE]

We can now read off the coefficients $B_{2}$ and $C_{2}$ of Eq. (A.10b) from (A.12) as ($\sigma\gg 1$)

[TABLE]

Similarly, starting with Eq. (A.4b), we arrive at

[TABLE]

To find $u_{1}$ we turn to Eq. (A.2a) and write the second-order terms as

[TABLE]

From Eqs. (A.13) and (A.14), it is clear that among the three terms on the right-hand side of Eq. (A.15), the first one is $O(\sigma^{\frac{1}{2}})$ and the remaining ones are $O(\sigma^{-\frac{1}{2}})$. Hence for $\sigma\gg 1$

[TABLE]

and consequently,

[TABLE]

We return to Eq. (A.3) to extract the $e^{i\omega t}$ term from the right-hand side and note that the right-hand side of Eq. (A.3) which yields such a term is:

[TABLE]

The dominating term in the above expression for $\sigma\gg 1$ is the first term in which the coefficient of $e^{i\omega t}$ is $-\frac{1}{2}\sigma^{\frac{5}{2}}A^{3}\left[\left(-\frac{4}{\sigma^{\frac{3}{2}}}\right)+\left(-\frac{2}{\sigma^{\frac{1}{2}}}\right)\frac{2}{\sigma}\right]=4A^{3}\sigma$. Using this approximation for the right-hand side of Eq. (A.3), with the left-hand side obtained from Eq. (A.8), the final amplitude equation becomes

[TABLE]