Fairy circles and temporal periodic patterns in the delayed plant-sulfide feedback model

1.   Introduction

Spatial patterns are widespread in ecological and chemical systems, such as salt marshes ^[1,2,3], predator-prey systems ^[4,5,6], the Brusselator model ^[7,8,9], Sel'kov model ^[10,11], Lengyel-Epstein system ^[12,13] and Degn-Harrison system ^[14,15]. It was proposed that spatial patterns may unveil underlying mechanisms that drive ecological resilience, and may serve as signals for environmental changes in many ecological systems (see ^[16,17,18]). It is assumed that self-organization theory is of great significance in helping us to understand spatial patterns, which was first proposed by Alan Turing in his seminal work ^[19]. Mathematically, asymptotic dynamics, which corresponds to persistent patterns, has been intensively studied, and the underlying mechanism was proposed as scale-dependent alternation between facilitation and inhibitory interactions, also called scale-dependent feedback ^[20]. There are many interesting persistent patterns which have been discovered, such as stable fairy circles and ring type patterns (^{[21,22,23,24,25,26]}). There also is an increasing recognition that dynamics on ecological time scales, called transient patterns, may be of some significance. The impermanence of transient patterns means that an ecological system in a transient state may change abruptly, even without any underlying change in environmental conditions, while the possibility of transient patterns implies that an ecological system may remain far from its asymptotic patterns for some period of time ^[27]. In a recent study conducted by Zhao et al. ^[28], transient fairy circles of S. alterniflora were observed in salt marsh pioneer zones in the Yangtze estuary north branch in eastern China coasts. The authors proposed that hydrogen sulfide ($\text{H}_2\text{S}$) may serve as a feedback regulator to drive such transient fairy circles, and proposed the following mathematical model to demonstrate the sulfide feedback mechanism:

The plant biomass concentration at the location $x$ and time $t$ is denoted by $P(x, t)$, where $x$ is a point on the plane. Assume the plant population growth follows the Verhulst-Pearl logistic pattern with the intrinsic growth rate $r$ and the carrying capacity $K$. The sulfide concentration is denoted by $S(x, t)$. Organic matters including plant biomass in intertidal salt marshes can produce hydrogen sulfide (H$_{2}$S). Assume the plant population produces hydrogen sulfide proportionally with the effective production rate $\epsilon$. Hydrogen sulfide is toxic to plants and can lead plants to die off. It was assumed that the loss of plants is to increase with an increase of the sulfide concentration and plant biomass, represented by $cPS$. As field studies about salt marsh plant species S. mariqueter and S. alterniflora show that plant lateral expansion through vegetative growth can be described by random walk, the diffusion process without drifts was employed to model plant dispersal with diffusion coefficient $D_P$. $D_S$ is the planar dispersion rate of the sulfide concentration. In addition, plants produce dissolved organic carbon which promote bioactivities of sulfate-reducing bacteria which, in turn, promote sulfide enrichment. This effect is represented by the term $\frac{k_s}{P(x, t)+k_s}$. The parameter $d$ is the maximum escape rate of sulfide through the mud-air interface. The parameter $\xi$ is dimensionless, which controls the time scale between the plant biomass and sulfide concentration. For a more detailed description of the model (1.1), we refer the reader to ^[28,29,30]. The conclusion of the study ^[28] is that transient fairy circle patterns in intertidal salt marshes can both infer the underlying ecological mechanisms and provide a measure of ecosystem resilience.

In this study, we would like to investigate how the self-correcting mechanism of the plant population influence transient patterns of fairy circles in salt marshes. For a single population growth, Hutchinson considered that for $\frac{K-P}{K}$ formally describing a self-regulatory mechanism, the only formal conditions that must be imposed on the biologically possible mechanism is that they operate so rapidly that the lag, $\tau$, is negligible between $t$ when any given value of $P$ is reached, and the establishment of the appropriately corrected value of the effective reproductive rate $r\frac{K-P}{K}$ ^[31]. Now, we would like to consider the "negligible time lag" for the plant population. It should be noticed that the time lag may not be negligible since different populations may have different growth properties. Particularly, for plant populations within constricted environments, the self-correcting adjustment with time lag may be necessary for understanding the underlying ecological mechanisms. Thus, the time delay will be incorporated in the self-regulatory mechanism. That is, the logistic growth term $rP(x, t)\left(1-\frac{P(x, t)}{K}\right)$ for the plant population will be replaced by $rP(x, t)(1-\frac{P(x, t-\tau)}{K})$. Thus, we will have a delayed plant-sulfide feedback model as follows:

For simplicity, we introduce the following non-dimensional variables

Dropping the hats of $\hat{t}$ and $\hat{\tau}$, we get the system with Neumann boundary conditions and initial conditions as follows

where $\nu$ is the outward unit normal vector on $\partial\Omega$.

The aim of this article is to conduct a detailed analysis about the effect of the time delay on the dynamics of the system (1.3). Our analysis shows that there is a critical value of the time delay. When the time delay is greater than the critical value, the system will have temporal periodic solutions. Since there are no analytical methods to study transient patterns yet, we will use numerical simulations to study transient patterns. We show that there are transient fairy circles for any time delay. However, there are different types of fairy circles and rings occurring in this system.

The rest of this paper is organized as follows. In Section 2, we mainly study the stability of the positive steady state of the system (1.3) by discussing the distribution of the eigenvalues, and give the sufficient conditions for the occurrence of Hopf bifurcations induced by the time delay. In Section 3, by using the center manifold theory and normal form theory for partial differential equations, we analyze the properties of Hopf bifurcations, and obtain the formulas determining the direction of Hopf bifurcation and stability of bifurcating periodic solutions. In Section 4, we conduct numerical studies to demonstrate transient patterns, and give some simulations to illustrate our theoretical results. The paper is closed with a brief discussion.

2.   Stability and Hopf bifurcation

Clearly, the system (1.3) has two nonnegative constant steady states: $E_0 = (0, 0)$ and $E^* = (u^*, v^*)$, where

We claim that $0 < u^* < k$. In fact,

Notice that $a, \; b, \; k > 0$, and we can obtain that $E^* = (u^*, v^*)$ is the unique positive steady state of the system (1.3).

Let $\Omega\subset\mathbb{R}^2$ be a bounded domain with smooth boundary $\partial\Omega$, and denote

In the phase space $\texttt{C} = C([-\tau, \; 0], \; X)$, we can rewrite the system (1.3) as

where $D = \text{diag}(d_1, d_2)$, $U_t(\cdot) = U(t+\cdot)$, and $F:\texttt{C}\to X$ is defined by

By calculation, we have the linearization of system (2.1) at $E^*$, which can be written as

where $L: \texttt{C}\to X$ is given by

and

Let

be all the eigenvalues of the operator $-\Delta$ on $\Omega$ with the Neumann boundary conditions. Using the techniques in ^[32], we obtain that the characteristic equation of the system (2.2) can be expressed as

where $I_2$ is the $2\times2$ identity matrix and $Q_n = \mu_n\text{diag}(d_1, d_2)$. That is, each eigenvalue $\lambda$ should satisfy the equation as follows

where

When $\tau = 0$, the characteristic equation (2.4) becomes

We give the following conclusion about the stability of the positive steady state $E^*$ as $\tau = 0$.

Theorem 2.1. When $\tau = 0$, all of the roots of the Eq (2.4) have negative real parts. That is, the positive steady state $E^*$ of the system (1.3) without delay is locally asymptotically stable.

Proof. If $T_n+B > 0, \; D_n+M_n > 0$ for all $n\in\mathbb{N}_0$, then all roots of the Eq (2.5) have negative real parts. By calculation, we have that $0 < u^* < k$ and $v^* > 0$. Recall that $k > 0, \; a > 0$, and $b > 0$, we can get

for any $n\in\mathbb{N}$. So, we complete the proof.

In the following, we will study the effect of delay on the stability of $E^*$. Recall that $D_n+M_n > 0$ for any $n\in\mathbb{N}_0$, we can see that $0$ is not the root of (2.4). Now, we will check whether there exist the critical values of $\tau$ such that (2.4) has a pair of simple purely imaginary eigenvalues for some $n\in\mathbb{N}_0$. Let $\pm \text{i}\omega(\omega > 0)$ be the solutions of the $(n+1)$-th Eq (2.4), then we get

Separating the real and imaginary parts, and we obtain that $\omega$ and $\tau$ should satisfy

It yields

where

Denote $Z = \omega^2$, and then (2.7) becomes

(2.9) has two roots:

From the proof of Theorem 2.1, we get $D_n+M_n > 0$ for any $n\in\mathbb{N}$. So, the sign of $D_n^2-M_n^2$ is the same to that of

To study the existence of positive roots of (2.9), we only need to discuss the sign of $D_n-M_n$, $T_n^2-B^2-2D_n$, and $(T_n^2-B^2)(T_n^2-B^2-4D_n)+4M_n^2$. It is not difficult to get the following Lemma.

Lemma 2.2. For the Eq (2.9), the following conclusions hold.

1) If $T_n^2-B^2-2D_n > 0$ and $D_n-M_n > 0$ or $(T_n^2-B^2)(T_n^2-B^2-4D_n)+4M_n^2 < 0$ for any $n\in\mathbb{N}_0$, then the Eq (2.9) has no positive root.

2) If there exists some $n\in\mathbb{N}_0$ such that $D_n-M_n < 0$, then the Eq (2.9) has a positive root $Z_n^+$.

3) If there exists some $n\in\mathbb{N}_0$ such that $T_n^2-B^2-2D_n < 0$, $D_n-M_n > 0$, and $(T_n^2-B^2)(T_n^2-B^2-4D_n)+4M_n^2\geq0$, then the Eq (2.9) has two positive roots $Z_n^\pm$.

Next, we discuss some sufficient conditions for Lemma 2.2 to hold. Using the first equation of (2.8), we have

It is easy to see that $T_0^2-B^2-2D_0 > 0$ when

holds. Similarly, we can get from (2.10) that $D_0-M_0 > 0$ provided that

To illustrate the existence of positive roots of the Eq (2.9), we mainly study three cases as follows.

Case Ⅰ: $b > \underline{b}$. If $b > \underline{b}$, we can obtain that $D_0-M_0 < 0$. It is easy to see that $D_n-M_n \rightarrow \infty$ as $n\rightarrow \infty$. So there must exist a $N\in\mathbb{N}$ such that $D_n-M_n < 0$ for $n < N$ and $D_n-M_n\geq0$ for $n\geq N$. From Lemma 2.2, we know that the Eq (2.9) has positive roots $Z_n^+$ for $n < N$.

Case Ⅱ: $\frac{u^*(1+u^*)d_2}{kd_1} < b < \min\{\underline{b}, \overline{b}\}$. If $b < \underline{b}$, then $D_0^2-M_0^2 > 0$. Furthermore, we can obtain from (2.10) that $D_n^2-M_n^2 > 0$ for any $n\in\mathbb{N}_0$ when $b > \frac{u^*(1+u^*)d_2}{kd_1}$. Notice that $T_0^2-B^2-2D_0^2 < 0$ when $b < \overline{b}$, and then from (2.8), we know that there exists some $N_1\in\mathbb{N}$ such that $T_n^2-B^2-2D_n^2 < 0$ for $n < N_1$, and $T_n^2-B^2-2D_n^2\geq0$ for $n\geq N_1$. Through calculation, we can obtain

where

Notice that $\lim_{n\rightarrow \infty} (T_n^2-B^2)(T_n^2-B^2-4D_n)+4M_n^2 = \infty$, and we can find a $N_2\in\mathbb{N}$ such that $(T_n^2-B^2)(T_n^2-B^2-4D_n)+4M_n^2\geq 0$ for $n\geq N_2$. Therefore, when $N_2 < N_1$ and $\frac{u^*(1+u^*)d_2}{kd_1} < b < \min\{\underline{b}, \overline{b}\}$ are satisfied, the Eq (2.9) has two positive roots $Z_n^\pm$ for $N_2\leq n < N_1$. Otherwise, the Eq (2.9) does not have a positive root for any $n\in\mathbb{N}_0$.

Case Ⅲ: $\max\{\overline{b}, \frac{u^*d_2(1+u^*)}{kd_1}\} < b < \underline{b}$. We have from Case Ⅰ and Case Ⅱ that $D_n^2-M_n^2 > 0$ and $T_n^2-B^2-2D_n^2 > 0$ for any $n\in\mathbb{N}_0$ when $\max\{\overline{b}, \frac{u^*d_2(1+u^*)}{kd_1}\} < b < \underline{b}$. Therefore, the Eq (2.9) has no positive root as $\max\{\overline{b}, \frac{u^*d_2(1+u^*)}{kd_1}\} < b < \underline{b}$.

Combining Lemma 2.2 with the above analysis, we have the following result.

Corollary 2.3. Denote $\mathcal{D}_1 = \{n\in\mathbb{N}_0|T_n^2-B^2-2D_n^2 < 0\}$, $\mathcal{D}_2 = \{n\in\mathbb{N}_0|(T_n^2-B^2)(T_n^2-B^2-4D_n)+4M_n^2\geq 0\}$, $\omega_n^\pm = \sqrt{Z_n^\pm}$.

1) If $b > \underline{b}$, then there exists some $N\in\mathbb{N}_0$ such that the Eq (2.7) has a positive root $\omega_n^+$ for $n < N$.

2) If $\frac{u^*(1+u^*)d_2}{kd_1} < b < \min\{\underline{b}, \overline{b}\}$ and $\mathcal{D}_1\bigcap\mathcal{D}_2\neq\emptyset$, then the Eq (2.7) has two positive roots $\omega_n^{\pm}$ for $n\in\mathcal{D}_1\bigcap\mathcal{D}_2$.

3)If $\max\{\overline{b}, \frac{u^*(1+u^*)d_2}{kd_1}\} < b < \underline{b}$, or $\frac{u^*(1+u^*)d_2}{kd_1} < b < \min\{\underline{b}, \overline{b}\}$ and $\mathcal{D}_1\bigcap\mathcal{D}_2 = \emptyset$, then the Eq (2.7) has no positive roots for any $n\in\mathbb{N}_0$.

For simplicity, we define the following set

It is easy to see that Eq (2.7) has a pair of positive roots $\omega_n^\pm$ for $n\in\Gamma\subset\mathbb{N}_0$. Then the Eq (2.4) has a pair of purely imaginary roots $\pm \text{i}\omega_n^\pm$ when $\tau$ takes the critical values $\tau_{n, j}^\pm$, which can be determined from (2.6), given by

for $j\in \mathbb{N}_0$. Let $\lambda(\tau) = \alpha(\tau)+\text{i}\omega(\tau)$ be the root of (2.4) satisfying $\text{Re}\lambda(\tau_{n, j}^\pm) = 0$ and $\text{Im}\lambda(\tau_{n, j}^\pm) = \omega_n^\pm$.

Lemma 2.4. Assume that the condition 2 or 3 of Lemma 2.2 holds, then

1) $\mathit{\text{Re}}\lambda'(\tau_{n, j}^\pm) = 0$, when $(T_n^2-B^2)(T_n^2-B^2-4D_n)+4M_n^2 = 0$.

2) $\mathit{\text{Re}}\lambda'(\tau_{n, j}^+) > 0$, $\mathit{\text{Re}}\lambda'(\tau_{n, j}^-) < 0$, when $(T_n^2-B^2)(T_n^2-B^2-4D_n)+4M_n^2 > 0$.

Proof. Differentiating the two sides of the Eq (2.4) with respect to $\tau$, we have

Thus,

Following the techniques in Cooke and Grossman ^[33], and using the Eqs (2.4) and (2.6), we obtain that

Note that

and we can complete the proof.

Combining Theorem 2.1 and Lemma 2.4, we have that $\tau_{n, 0}^+ < \tau_{n, 0}^-$ holds true for $n\in\Gamma$. Then, we define the smallest critical value such that the stability of $E^*$ will change, which can be given by

Combined the above analysis with Corollary 2.4 in Ruan & Wei ^[34], we know that all of the roots of (2.4) have negative real parts when $\tau\in [0, \tau^*)$, and the $(n+1)$-th equation of (2.4) has a pair of simply purely imaginary roots when $\tau = \tau_{n, j}^\pm$. Moreover, we see that (2.4) has at least one pair of conjugate complex roots with positive real parts when $\tau > \tau^*$. Based on the above discussion, we can obtain the following conclusion about the stability of $E^*$.

Theorem 2.5. For $\tau^*$ defined in (2.15), the following statements about system (1.3) hold true.

1) If $T_n^2-B^2-2D_n > 0$ and $D_n-M_n > 0$ or $(T_n^2-B^2)(T_n^2-B^2-4D_n)+4M_n^2 < 0$ hold for any $n\in\mathbb{N}_0$, then the positive steady state $E^*$ is locally asymptotically stable for any $\tau\geq0$.

2) If $D_n-M_n < 0$ holds for some $n\in\mathbb{N}_0$, then

(a) the positive steady state $E^*$ is locally asymptotically stable for $\tau \in [0, \tau^*)$, and unstable for $\tau > \tau^*$.

(b) the system (1.3) undergoes a Hopf bifurcation at $E^*$ when $\tau = \tau_{n, j}^+$ for $j\in\mathbb{N}_0$.

3) If $T_n^2-B^2-2D_n < 0$, $D_n-M_n > 0$, and $(T_n^2-B^2)(T_n^2-B^2-4D_n)+4M_n^2 > 0$ hold for some $n\in\mathbb{N}_0$, then the positive steady state $E^*$ is locally asymptotically stable for $\tau \in [0, \tau^*)$. Moreover, the system (1.3) undergoes a Hopf bifurcation at $E^*$ when $\tau = \tau_{n, j}^\pm$ for $j\in\mathbb{N}_0$, where $\tau_{n, j}^\pm$ is defined in (2.14).

Remark 2.6. From Lemma 2.4 and Theorem 2.5, we have that

provided that $T_n^2-B^2-2D_n < 0$, $D_n-M_n > 0$, and $(T_n^2-B^2)(T_n^2-B^2-4D_n)+4M_n^2 > 0$ are satisfied. In this case, the stability switch may exist.

3.   Stability and direction of the Hopf bifurcations

In this section, we study the stability and direction of the Hopf bifurcations by applying the center manifold theorem and the normal formal theory of partial functional differential equations ^[32,35]. First, the system (1.3) can be represented as an abstract ODE system. Second, on the center manifold of the ODE system corresponding to $E^*$, the normal form or Taylor expansion of the ODE system will be computed. Then, using the techniques in ^[36], the coefficients of the first 4 terms of the normal form will reveal all the properties of the periodical solutions. These analytical results will also be used in numerical studies.

In this section, we choose $\Omega = [0, l\pi]\times[0, l\pi]$. To write the system (1.3) as an ODE system, we define a function space

where $u_1(\cdot, t) = u(\cdot, \tau t)-u^*, \; u_2(\cdot, t) = v(\cdot, \tau t)-v^*$, and $U(t) = (u_1(\cdot, t), u_2(\cdot, t))^T$. Then the system (1.3) can be written as

in the function space $\texttt{C} = C([-1, 0], X)$, where $D = \text{diag}(d_1, d_2)$, $L(\tau)(\cdot):\texttt{C}\to X$ and $f:\texttt{C}\times\mathbb{R}\texttt{C}\to X$ are given, respectively, by

with

$f_1(\varphi, \tau) = a_1\varphi_1(0)\varphi_2(0)+a_2\varphi_1(0)\varphi_1(-1),$

$f_2(\varphi, \tau) = a_3\varphi_1^2(0)+a_4\varphi_1(0)\varphi_2(0)+a_5\varphi_1^3(0)+a_6\varphi_1^2(0)\varphi_2(0)+\mathcal{O}(4),$

for $\varphi = (\varphi_1, \varphi_2)^T\in\texttt{C}$,

where

Let $\tau = \tau^*+\mu$, and then (3.1) can be rewritten as

where

for $\varphi\in\texttt{C}$.

From the analysis in Section 2, we know that system (3.2) undergoes Hopf bifurcation at the equilibrium $(0, 0)$ when $\mu = 0$ (i.e., $\tau = \tau^*$). We assume that when $\mu_{n_0} = \frac{j_0^2+k_0^2}{l^2}$, (2.4) has roots $\pm\text{i} \omega^*$ as $\tau = \tau^*$. Moreover, we also have that $\pm \text{i} \omega^*\tau^*$ are simply purely imaginary eigenvalues of the linearized system of (3.2) at the origin:

as $\mu = 0$ and all other eigenvalues of (3.3) at $\mu = 0$ have negative real parts.

The eigenvalues of $\tau D\Delta$ on $X$ are $-\tau d_1\frac{j^2+k^2}{l^2}$ and $-\tau d_2\frac{j^2+k^2}{l^2}, \; j, k\in\mathbb{N}_0$, with corresponding eigenfunctions $\beta_{j, k}^1(x) = (\gamma_{j, k}(x), 0)^T$ and $\beta_{j, k}^2(x) = (0, \gamma_{j, k}(x))^T$, where $x = (x_1, x_2)$, $\gamma_{j, k}(x) = \frac{\cos\frac{jx_1}{l}\cos\frac{kx_2}{l}}{\sqrt{\int_0^{l\pi}\cos^2\frac{jx_1}{l}\text{d}x_1\int_0^{l\pi}\cos^2\frac{kx_2}{l}\text{d}x_2}}.$

We define a space as $M_{j, k} = \text{span}\{\langle\varphi, \beta_{j, k}^i\rangle\beta_{j, k}^i:\varphi\in\texttt{C}, i = 1, 2\}, \; j, k\in\mathbb{N}_0$, and the inner product $\langle\cdot, \cdot\rangle$ is defined by

Then, on $M_{j, k}$, the Eq (3.3) is equivalent to the ODE on $\mathbb{R}^2$:

Now, we compute the normal form in the center manifold. There are several steps. We first compute eigenvectors of the infinitesimal generator of the semigroup defined by the linearized system at $\tau = \tau^*$. From the Riesz representation theorem, there exists a bounded variation function $\eta_{j, k}(\mu, \theta)$ for $\theta\in[-1, 0]$, such that

for $\varphi\in C([-1, 0], \mathbb{R}^2)$. In fact, we can choose

Let $A$ denote the infinitesimal generator of the semigroup defined by (3.4) with $\mu = 0, \; j = j_0, \; k = k_0$ and $A^*$ denote the formal adjoint of $A$ under the bilinear form

for $\phi\in C([-1, 0], \mathbb{R}^2)$ and $\psi\in C([0, 1], \mathbb{R}^{2^T})$. Then, we know that $\pm \text{i}\omega^*\tau^*$ are simply purely imaginary eigenvalues of $A$, and they are also eigenvalues of $A^*$. By direct computations, we get $q(\theta) = q(0)\text{e}^{\text{i}\omega^*\tau^*\theta} = (1, q_1)^T\text{e}^{\text{i}\omega^*\tau^*\theta}(\theta\in[-1, 0])$ is eigenvector of $A$ corresponding to $\text{i}\omega^*\tau^*$, where

Similarly, we have $q^*(s) = \text{e}^{-\text{i}\omega^*\tau^* s}(1, q_2)(s\in[0, 1])$ is eigenvector of $A^*$ corresponding to $\text{i}\omega^*\tau^*$, where

Let $\Phi = (\Phi_1, \Phi_2) = (\text{Re}q, \text{Im}q)$ and $\Psi^* = (\Psi_1^*, \Psi_2^*)^T = (\text{Re}q^*, \text{Im}q^*)^T$. Denote

where

Let $\Psi = (\Psi_1, \Psi_2)^T = (\Psi^*, \Phi)_{j_0, k_0}^{-1}\Psi^*$, $(\Psi, \Phi)_{j_0, k_0} = I_2$, and $I_2$ is a $2\times2$ identity matrix.

We now write the reduced equation on the center manifold. The center subspace of linear equation (3.3) with $\mu = 0$ is given by $P_{CN}\texttt{C}$, where

with $\beta_{j_0, k_0} = (\beta_{j_0, k_0}^1, \beta_{j_0, k_0}^2)$ and $c\cdot\beta_{j_0, k_0} = c_1\beta_{j_0, k_0}^1+c_2\beta_{j_0, k_0}^2$ for $c = (c_1, c_2)^T\in\texttt{C}$. Let $P_S\texttt{C}$ denote the stable subspace of linear equation (3.3) with $\mu = 0$, and then $\texttt{C} = P_{CN}\texttt{C}\bigoplus P_S\texttt{C}$.

Using the decomposition $\texttt{C} = P_{CN}\texttt{C}\bigoplus P_S\texttt{C}$ and following ^[32], the flow of (3.2) with $\mu = 0$ in the center manifold is given by the following formulae:

with $h(0, 0, 0) = 0$ and $Dh(0, 0, 0) = 0$.

Let us write the reduced equation in complex form. Set $z = y_1-\text{i}y_2$ and $\Psi(0) = (\Psi_1(0), \Psi_2(0))^T$, and then $q = \Phi_1+\text{i}\Phi_2$ and $\Phi(y_1(t), y_2(t))^T\cdot \beta_{j_0, k_0} = (qz+\overline q\overline z)\cdot \beta_{j_0, k_0}/2$. Thus, (3.7) can be written as

where

From (3.8) and (3.9), we obtain that $z$ satisfies

where

Now, let us compute $g(z, \overline z)$. Set

Let $(\psi_1, \psi_2) = \Psi_1(0)-\text{i}\Psi_2(0)$. From (3.7), (3.9) and (3.10), we can get the following quantities:

and

To obtain $g_{21}$, we need to compute $w_{11}$ and $w_{20}$. The calculation of $w_{11}$ and $w_{20}$ is somewhat tedious. Let $A_U$ denote the generator of the semigroup generated by the linear system (3.3) with $\mu = 0$. From (3.9) and (3.10), we have

where

Let

Furthermore, by comparing the coefficients, we obtain that

By using the chain rule,

and we obtain, from (3.11) and (3.12), that

As $2\text{i}\omega^*\tau^*$ and 0 are not characteristic values of (3.3), (3.14) has unique solutions $w_{20}$ and $w_{11}$ in $P_S\texttt{C}$, given by

Using the definition of $A_U$, we get, from the first equation (3.13) and (3.14), that for $\theta\in[-1, 0]$,

Therefore

where $E$ is a 2-dimensional vector in $X$. According to the definition of $\beta_{j, k}^i(i = 1, 2)$ and $q(\theta)(\theta\in[-1, 0])$, we have

From (3.14), we get that

Representing $E$ and $f_{z^2}$ by series:

We get from (3.16) that

So, $E_{j, k}$ could be calculated by

where

with

Similarly, we get

where

with

Then, the coefficient $g_{21}$ is completely determined.

Let $\lambda(\tau) = \alpha(\tau)+\text{i}\omega(\tau)$ denote the eigenvalues of (3.3). Thus we can compute the following quantities:

According to the Hopf bifurcation theory (see ^[36]), we write our derivation as a theorem which is well-known.

Theorem 3.1. The quantity $\mu_2$ determines the direction of the Hopf bifurcation (forward if $\mu_2 > 0$, backward if $\mu_2 < 0$). The quantity $\beta_2$ determines the stability of the bifurcating periodic solutions (stable if $\beta_2 < 0$, unstable if $\beta_2 > 0$). The quantity $T_2$ determines the period of the bifurcating periodic solutions (the period increases if $T_2 > 0$, decreases if $T_2 < 0$).

4.   Numerical studies

In this section, we perform two sets of numerical studies. According to theoretical derivations in the previous sections, we choose parameter values to produce numerical simulations. We observe several different shapes of transient fairy circles, and stable temporal periodic solutions.

Example 4.1. We choose $x = (x_1, \; x_2)\in\Omega = [0, 2\pi]\times[0, 2\pi], \; d_1 = 1, \; d_2 = 3, \; a = 0.1, \; b = 0.4, \; k = 1$. This set of parameters satisfies Case Ⅰ of Section 2. Through calculation, we get $E^*(0.7016, 0.2984)$, and only for

(2.4) has pure roots $\pm\text{i}\omega_n^+, \; n = 0, 1, 2$, where

and the corresponding bifurcation values $\tau_{k, j}^+$ are

respectively. When $\tau < \tau^* = 1.9885$, the positive equilibrium $E^*(0.7016, 0.2984)$ of (1.3) is asymptotically stable (see Figure 1). Before the system reaches this stable equilibrium solution, there are transient fairy circles. The plant population and sulfide concentration have similar-shaped fairy circles.

By the formulas derived in the previous section, we get $c_1(0)\approx -0.5506 + 1.3320\text{i}$. Because $\text{Re}c_1(0) < 0$, we know that when $\tau > \tau^* = 1.9885$, there exist orbitally stable periodic solutions (see Figure 2). We observe that there are transient fairy circles before the system tends to stable periodic solutions. However, the shapes of the fairy circles for the plant population and sulfide concentration are different. In the numerical simulations for Figures 1 and 2, the initial conditions are

for $t\in[-\tau, 0]$.

Example 4.2. We choose $x = (x_1, \; x_2)\in\Omega = [0, 2\pi]\times[0, 2\pi], \; d_1 = 0.5, \; d_2 = 0.1, \; a = 0.2, \; b = 0.2, \; k = 1$. This set of parameters satisfies Case Ⅱ of Section 2. Through calculation, we get $E^*(0.4142, 0.5858)$, and only for

(2.4) has pure roots $\pm\text{i}\omega_n^{\pm}, \; n = 0, 1, 2$, where

and the corresponding bifurcation values $\tau_{k, j}^+$ are

respectively. When $\tau < \tau^* = 2.8577$, the positive equilibrium $E^*(0.7016, 0.2984)$ of (1.3) is asymptotically stable (see Figure 3). We observe that there are transient fairy circles before the system reaches the stable equilibrium point. The shapes of the fairy circles for the plant population and sulfide concentration are different.

By the formulas derived in the previous section, we get $c_1(0)\approx -1.3262 + 1.2114\text{i}$. Because $\text{Re}c_1(0) < 0$, we know that when $\tau > \tau^* = 2.8577$, there exist orbitally stable periodic solutions (see Figure 4). We observe that there are transient fairy circles. The shapes of the fairy circles for the plant population and that for the sulfide concentration are the same as that in the equilibrium point case, respectively, although the shapes of the fairy circles for the plant population and sulfide concentration are different. In the numerical simulations for Figures 3 and 4, the initial conditions are

for $t\in[-\tau, 0]$.

5.   Discussion

Transient spatial patterns in ecosystems have gained more attention in recent research. Several mathematical models were proposed to understand fairy circles and rings observed in salt marshes. One conclusion drawn from this research was that transient fairy circle patterns can infer the underlying ecological mechanisms and provide a measure of resilience of salt marsh ecosystems. The underlying mechanism proposed was plant-sulfide feedbacks instead of scale-dependent feedbacks. It is assumed that any population of a species in a location has a self-regulatory mechanism, and the time of mechanism action has some time lag. In this research, we consider how the time delay in the self-regulatory mechanism influences fairy circle formations in plant-sulfide feedbacks. Based on a mathematical model ^[28], we proposed a delay plant-sulfide feedback system. We performed a detailed investigation of the plant-sulfide feedback system subject to Neumann boundary conditions, and identified the parameter ranges of stability of the positive equilibrium and the existence of Hopf bifurcation. There is a critical value of the time delay. When the time delay is smaller than the critical value, the system will approach the spatial homogeneous equilibrium state asymptotically. At the critical value, Hopf bifurcations occur. When the time delay is greater than the critical value, the system will have temporal periodic solutions. Since we do not have any method to analyze transient patterns yet, we numerically demonstrated that there always are transient fairy circles for any time delay. We also found that there are different shapes of fairy circles or rings. This confirms that transient fairy circle patterns in intertidal salt marshes can infer the underlying ecological mechanisms and provide a measure for ecosystem resilience.

In ^[28], simulations were conducted with several points as initial values, and several circles progressed to one circle and to a uniform distribution. It is easy to see that the number and shape of fairy circles may depend on initial values. In our numerical study, we observed several different patterns of transient fairy circles. The number and shape of these fairy circles seem also to depend on initial distributions. This may be reasonable in reality since a natural species population usually starts with various situations and then progresses to its asymptotical patterns.

It is conventional that a natural population will obey the self-regulatory mechanism which is different from the scale-dependent feedback mechanism. The scale-dependent feedback mechanism assumes that scale-dependent feedbacks between localized facilitation and large-scale inhibition induce spatial self-organization. The self-regulatory mechanism assumes that the population growth rate is limited by its total population. The former may be considered as a special case of the latter. However, when mathematical models are constructed for those mechanisms, different terms should be taken in the equations. This is the reason we incorporate the self-regulatory mechanism with time delay to the plant-sulfide feedback mechanism.

In ^[28], two additional models, the nutrient depletion model and scale-dependent model, were also proposed to explain transient fairy circle patterns. We may consider to study time delayed versions of those models in order to compare how different mechanisms with natural time delay influence transient spatial patterns in the future.

In general, it is difficult to characterize transient patterns in dynamical systems. One direction to attempt may be time-transformation. Given that we want to know the dynamics of a system during a finite period of time, we make a time-transformation that changes this finite period of time to infinity, and then study the transformed system. However, it is required that the time period is given. It seems that we need to develop new analytical tools for transient patterns. This is an interesting mathematical question in dynamical systems for the future.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgment

This research is supported by the National Natural Science Foundation of China (Nos.11901172, 12271144) and Fundamental Research Fund for Heilongjiang Provincial Colleges and Universities (Nos. 2021-KYYWF-0017, 2022-KYYWF-1043).

Conflict of interest

The authors declare that there is no conflict of interest.