Dynamics of a delay turbidostat system with contois growth rate

1.   Introduction

The chemostat is a laboratory bio-reactor used for the continuous culture of microorganism (Figure 1). The apparatus is of both ecological and mathematical interest since it can be used to represent a lot of microorganism systems, such as the simple lake and the waste-water treatment, and plays an important role in micrology and population ^[8,9,30].

Recently, a large of research devoted to modifying the chemostat for higher economical values by controlling the dilution rate of the chemostat. The chemostat with the feedback control of the dilution rate is established, which is referred to as turbidostat by biologists ^[14] (Figure 1). In the turbidostat, the optical sensor measures the concentration of the microorganism and the signal feedback control the dilution rate. For example, Flegr ^[5] analyzed the coexistence of two species in the turbidostat by numerical analysis, and later Leenheer and Smith ^[13] also demonstrated results of Flegr by theoretical analysis. They showed that a turbidostat with monotone response functions permits a unique coexistence equilibrium and if it is locally asymptotically stable then it is globally attracting. Furthermore, they also obtained that coexistence is not achievable in the turbidostat with more than two competing organisms. Li ^[14] established a mathematical model of competition in a turbidostat for an inhibitory nutrient and achieved the stability of the equilibrium and the existence of limit cycles. Li and Chen ^[16] studied a mathematical model of the turbidostat with impulsive state feedback control and obtained the existence and the stability of periodic solutions of order one.

The uptake function in the model of microorganism continuous culture is the growth rate of the microorganism including many forms, such as, Monod ^[19], Moser ^[20] and Contois ^[4]. The most classical and common uptake function is Monod one, which is the function of concentration of microorganism and the constant half saturation term. However, Contois ^[4] presented experimental results to show that it is not necessarily a constant. Thus, the uptake function of Contois was introduced, which depends on the ratio substrate to microorganism. Jost ^[24] explored the relation of the Contois growth rate in microbiology and the ratio-dependent functional response in population ecology introduced by Arditi and Ginzburg ^[1]. As such, many researchers also call it ratio-dependent growth rate just as in population ecology ^{[10,11,17,34]}. Especially, Hu et al. ^[11] investigated the existence and the stability of the periodic solution of order one for the turbidostat system with ratio-dependent growth rate and impulsive state feedback control. The Contois model gave predictions that were in excellent agreement with experiment measurement. Especially, the Contois model was found to accurately describe the processing of industrial wastewaters. Nelson investigated a series of chemostat models with Contois growth rate. For example, Nelson and Sidhu ^[21] investigated the Hopf bifurcation and degenerate Hopf bifurcation of the chemostat model with Contois growth rate and the variable yield. Alqahtani et al. ^[22] studied the chemostat model with variable yield and contois growth kinetics with substrate inhibition and showed the Hopf bifurcation, degenerate Hopf bifurcation and Bogdanov-Takens bifurcation.

Researchers recognized that time delays are natural in the ecological systems ^[12,27,29]. Smith and Waltman ^[30] showed that there are two obvious sources of delays in the cultivation of microorganism, one of which is the possibility that the microorganism stores the nutrient. Caperon ^[3] first introduced the delay into the chemostat model by some experiments. Bush and Cook ^[2] corrected Caperon's model and established a chemostat model with delay in the growth rate of microorganism. The time delays yield some complex impacts to the chemostat model and biologists gave the reasonable explanations for the observations and experimental datas by investigating the chemostat model with delays ^{[18,15,23,26,31]}. For example, Li et al. ^[15] considered the global dynamics of the chemostat model with two populations of microorganisms competing for two perfectly complementary nutrients when distributed delays are involved. Ruan and Wolkowicz ^[26] studied the existence and the stability of the Hopf bifurcation of the chemostat model with a distributed delay and found that the periodic solutions become unstable if the dilution rate is increased. Wang and Wolkowicz ^[31] analyzed the local stability of the equilibria and the globally asymptotical stability of the single species survival equilibrium when $n$ species compete in the chemostat with time delay for a single resource, hence the competitive exclusion principle holds. Yao et al. ^[32] discussed the delayed turbidostat model with Monod growth rate and obtained the phenomenon of oscillation. In the turbidostat, additional factor causing delay is due to the process that the sensor controls the dilution, for example, Tagashira and Hara ^[28] and Yuan et al. ^[33] considered the turbidostat models of the dilution rate with delayed feedback control.

Li and Xu ^[17] investigated the following chemostat system with Contois growth rate and time delay

where $S(t)$ and $x(t)$ present the concentration of substrate and the concentration of microorganism at time $t$ , respectively, $S_0>0$ stands for the input concentration of the nutrient, $\alpha+\beta S(t)$ $(\alpha>0, \beta>0)$ is variable yield, $m>0$ and $a>0$ are growth parameters of the microorganism, $\tau>0$ is the time delay of the growth response of the microorganism, and $d>0$ presents the flow volume. They considered the local and global stability of the equilibria and Hopf bifurcation.

In this paper, motivated by the works of Bush and Cook (^[2]), Contois (^[4]) and Li and Xu (^[17]), one focuses on the dynamics of a turbidostat model with Contois growth rate and time delay which follows as

where $S(t)$ , $x(t)$ and parameters $S_0$ , $m$ , $a$ , $\tau$ are defined as in (1.1), $\gamma>0$ is yield constant and $d+kx(t)$ $(d>0, k>0)$ presents the dilution rate of the turbidostat.

The aim of this paper is to investigate the qualitative properties of system (1.2) including the stability of the equilibria and the bifurcations. The outline of this paper is as follows. In Section 2, I analyze the existence and local stability of the equilibria, and the transcritical bifurcation at the boundary equilibrium and the Hopf bifurcation at the positive equilibrium. In Section 3 mostly focuses on the stability and type of the bifurcating periodic solutions induced by Hopf bifurcation to system (1.2). I in Section 4 further illustrate our main results by numerical simulations.

2.   Stability of the equilibria and bifurcations

In this section, we investigate the existence and stability of the equilibria and the bifurcations at equilibria including the transcritical bifurcation and Hopf bifurcation of system (1.2).

It is convenient to introduce dimensionless variables. In particular, we set

System (1.2) becomes

where we still denote $\tilde{S}, \tilde{x}, \tilde{k}, \tilde{a}, \tilde{m}, \tilde{t}$ with $S, x, k, a, m, t$ , respectively. First of all, we discuss the existence of the equilibria of system (2.1) and have the following result.

Lemma 2.1. System (2.1) always has a boundary equilibrium $E_0: (1, 0)$ ; If $m>1$ , then system (2.1) has a unique positive equilibrium $E_1: (S^*, x^*)$ , where $S^* = \frac{1+k}{m+k}$ when $a = 1$ ,

when $a\neq1$ and $x^* = 1-S^*$ .

Proof. From the right part of system (2.1), we can obtain the boundary equilibrium easily, which we denote as $E_0: (1, 0)$ . In order to get the positive equilibrium, we need to discuss the roots of the equation $\frac{mS}{a+(1-a)S}-1-k(1-S) = 0$ on the interval $(0, 1)$ . We define

Therefore, we just need to discuss the zeros of the function $f(S)$ on the interval $(0, 1)$ , which need to be discussed in the following three cases.

(ⅰ). If $a = 1$ , then we have $f(S) = (k+m)S-1-K$ and $f(S)$ has a zero $S = \frac{1+k}{m+k}$ on the interval $(0, 1)$ when $m>1$ , and $f(S)$ has no zero on the interval $(0, 1)$ when $m\leq1$ . Thus, we obtain $S^* = \frac{1+k}{m+k}$ when $a = 1$ and $m>1$ .

(ⅱ). If $0 < a < 1$ , then $f(S)$ is a quadratic function of $S$ . Since the coefficient of quadratic term $k(1-a)>0$ and $f(0) = -a(1+k) < 0$ , $f(S) = 0$ has two roots

and it is obvious that $S_1 < 0, S_2>0$ . By simplifying the inequality $S_2 < 1$ , we can get that $S_2 < 1$ if and only if $m>1$ . Thus, we have $S^* = S_2$ when $0 < a < 1$ and $m>1$ .

(ⅲ). If $a>1$ , then $f(S)$ still is a quadratic function of $S$ with the coefficient of quadratic term $k(1-a) < 0$ and $f(0) = -a(1+k) < 0$ . The discriminant of $f(S)$ is

We can take $\Delta$ as the function of $a$ . Therefore, the discriminant of $\Delta$ is $\tilde{\Delta}: = km(m-1)(1+k)$ . For this case we need further discussion in the following three cases.

(a). If $m = 1$ , then $\tilde{\Delta} = 0$ and the double root of $\Delta = 0$ is $-k$ . Therefore, $\Delta>0$ holds when $a>1$ . We have $f(S) = 0$ has two roots $S_1 = \frac{a(1+k)}{k(a-1)}>1$ and $S_2 = 1$ . Thus, $f(S)$ has no zero on the interval $(0, 1)$ when $a>1$ and $m = 1$ .

(b). If $m < 1$ , then $\tilde{\Delta} < 0$ , which implies $\Delta>0$ . We have $f(S)$ has two roots $S_1$ and $S_2$ presented in (2.3) and $S_2 < S_1$ . By simplifying $S_2>1$ , we can get that $S_2>1$ if and only if $m < 1$ . Thus, $f(S)$ has no zero on the interval $(0, 1)$ when $a>1$ and $m < 1$ .

(c). If $m>1$ , then $\tilde{\Delta}>0$ . Since $4km-2k+2m-2>0$ and $(k-m+1)^2\geq0$ , the two roots of $\Delta = 0$ denoted by $a_1$ and $a_2$ satisfy $a_1 < a_2\leq0$ . Therefore, $\Delta>0$ holds when $a>1$ and $m>1$ . We obtain $f(S)$ has two roots presented in (2.3) and $S_1>S_2$ . Further, we can get $0 < S_2 < 1$ if and only if $m>1$ and $a>1$ . Suppose that $S_1 < 1$ , then we get $m>1$ , which is a contradiction. Thus, we have $S^* = S_2$ when $a>1$ and $m>1$ .

This completes the proof.

From the above analysis we have the existence of equilibria of system (2.1). In the following, we consider qualitative properties of system (2.1) including stability of the equilibria and the bifurcations. We discuss system (2.1) in two cases: (ⅰ) $\tau = 0$ and (ⅱ) $\tau>0$ .

Case (ⅰ): $\tau = 0$ , i.e. there is no time delay in system (2.1), then system (2.1) becomes

We first investigate the stability of the equilibria $E_0$ and $E_1$ and have the following result.

Theorem 2.2. If $m < 1$ , then $E_0$ is a stable node; If $m>1$ , then $E_0$ is a saddle and $E_1$ is a stable node; If $m = 1$ , then $E_0$ is a saddle-node.

Proof. The Jacobian matrix at $E_0$ is

The determinant, trace and discriminant are respectively

If $m < 1$ , then $D_0>0$ and $T_0 < 0$ . Thus, $E_0$ is a stable node. If $m>1$ , then $D_0 < 0$ . Thus, $E_0$ is a saddle. If $m = 1$ , then $D_0 = 0$ and $T_0 = -1$ . This is the degenerate case, which needs further discussion.

Translating the equilibrium $E_0$ to the origin by $\tilde{S} = S-1$ , system (2.4) becomes the following system

where we still use $S$ to present $\tilde{S}$ . Using the linear transformation $u = x$ , $v = S+x$ and time-rescaling $t_1 = -t$ to normalize the linear part of system (2.5), we can change system (2.5) into the following

where we still denote $t_1$ as $t$ . By the implicit function theorem, there is a unique function $v = \phi(u)$ such that $\Psi(u, v) = 0$ . We can obtain $v = \phi(u) = 0$ by solving $\Psi(u, v) = 0$ . Substituting $v = 0$ into $\Phi(u, v) = 0$ , we get

Thus $E_0$ is a saddle-node when $m = 1$ by Theorem 7.1 in Zhang et al. (^[35]). Moreover, the parabolic sector of the saddle-node lies on the right-hand plane in the $(u, v)$ -coordinates and the two hyperbolic sectors lie on the left-hand plane since $a+k>0$ .

We continue to prove the properties of $E_1$ .

The Jacobian matrix at $E_1$ is

The determinant, trace and discriminant of $J_1$ are respectively

Thus, $E_1$ is a stable node if it exists.

We complete the proof.

It is indicated in Lemma 2.1 that system (2.4) has either exact one equilibrium $E_0$ when $m\leq1$ or exact two equilibria $E_0$ and $E_1$ when $m>1$ . In the following we reduce the system to a 1-dimensional system on a center manifold and display the mechanism for $E_1$ to arise.

Theorem 2.3. System (2.4) experiences a transcritical bifurcation at $E_0$ when $m = 1$ .

Proof. Let $\varepsilon = m-1$ and translate the equilibrium $E_0$ to the origin by $\tilde{S} = S-1$ , system (2.4) becomes the following system

where we still denote $\tilde{S}$ as $S$ .

Applying the transformation $S = -u+v$ and $x = u$ to diagonalize the linear part of system (2.7), we can change system (2.7) into the suspended system

By the center manifold theory, system (2.8}) has a two-dimensional center manifold $W^c: v = h(u, \varepsilon)$ near $O$ , which is tangent to the plane $v = 0$ at $O$ in the $(u, v, \varepsilon)$ -space. In order to obtain the second-order approximation of function $h$ , we set

Substituting (2.9) into the equality $\dot{v} = h_u\dot{u}$ and comparing the coefficient of $u^2$ , $\varepsilon^2$ and $u\varepsilon$ , we have $a_1 = b_1 = c_1 = 0$ . Hence the center manifold is $v = o(|u, \varepsilon|^2)$ and the restricted system of (2.8}) on the center manifold (2.9) is

The expression (2.10) shows that a transcritical bifurcation occurs at $E_0$ as $\varepsilon$ varies through the bifurcation value $\varepsilon = 0$ (^[6]). More concretely, when $\varepsilon < 0$ , $E_0$ is stable and the other equilibrium appears on the negative $u$ -axis; when $\varepsilon = 0$ , the two equilibria coincide at $E_0$ , which is a saddle-node; when $\varepsilon>0$ , $E_0$ remains an equilibrium but is unstable while a stable equilibrium $E_1$ arises.

We complete the proof.

Case (ⅱ): $\tau>0$ .

In this case we will investigate the effect of time delay on the system. The time delay can cause the loss of stability of $E_1$ and can induce periodic solutions. Note that if $m>1$ , then system (2.1) has a unique positive equilibrium $E_1$ . To further consider the local stability of $E_1$ and the Hopf bifurcations induced by the delay, we set $\tilde{x} (t) = x(t)-x^*, \tilde{S} (t) = S(t)-S^*$ and still denote $\tilde{x}, \tilde{S}$ as $x, S$ . Then system (2.1) can be changed into the following system

where

The following characteristic equation can be achieved from the linear system of system (2.11)

To investigate the stability of the equilibrium $E_1$ and Hopf bifurcation of (2.1), we must study the distribution of the roots of (2.12).

Lemma 2.4. If $a_{11}b_{12}+a_{12}b_{11}-a_{11}b_{13} < 0$ , then $\pm iw_0(w_0>0)$ are the eigenvalues of (2.12) when $\tau = \tau_j$ . The values of $w_0$ and $\tau_j$ can be presented as follows

Proof. If $\lambda = iw$ $(w>0)$ is a root of (2.12), then we have

Separating the real and imaginary parts, we obtain

which implies the following equation

Denote

Since

and the condition $a_{11}b_{12}+a_{12}b_{11}-a_{11}b_{13} < 0$ , we have the constant term of $h(w)$ is negative. In addition, we can obtain the coefficient of quadratic term of $h(w)$ is positive. In fact,

Thus, $h(w)$ has a unique positive real root

Substituting $w_0$ into (2.13), we conclude that

Hence, (2.12) has a pair of purely imaginary roots $\pm iw_0$ as $\tau = \tau_j,$ $j = 0, 1, 2\cdots$ .

This completes the proof.

Let $\lambda(\tau) = \alpha(\tau)+i\beta(\tau)$ denote the root of (2.12) near $\tau = \tau_j$ satisfying $\alpha(\tau_j) = 0$ and $\beta(\tau_j) = w_0$ . Then we have the following transversality condition.

Lemma 2.5. The following transversality condition holds

Proof. Differentiating (2.12) with respect to $\tau$ yields that

Thus,

Since $h(w)$ is a quartic function of $w$ and the leading coefficient of which is positive, furthermore, from the analysis of Lemma 2.4, we obtain that the four roots of equation $h(w) = 0$ are a pair of conjugate complex roots, a negative real root and a positive real root $w_0$ , respectively. Hence $h(w)$ is monotonically increasing at $w_0$ , i.e., $h'(w_0)>0$ , which implies $h'(w_0)/\{2w_0^3b_{13}^2+2w_0(a_{11}b_{13}-a_{12}b_{11})^2\}>0$ . Thus, we have

The proof is completed.

Now we have the following result about the distribution of the roots of the exponential polynomial (2.12) by Corollary 2.4 in ^[25].

Lemma 2.6. Suppose that $w_0$ and $\tau_j$ $(j = 0, 1, 2\cdots)$ are defined by (2.14) and (2.15), respectively. We have the following results

($i$ ) If $a_{11}b_{12}+a_{12}b_{11}-a_{11}b_{13}\geq0$ , then all the roots of (2.12) have negative real parts for all $\tau>0$ .

($ii$ ) If $a_{11}b_{12}+a_{12}b_{11}-a_{11}b_{13} < 0$ and $\tau = \tau_j$ , then (2.12) has a pair of simple imaginary roots $\pm iw_0$ . Furthermore, if $\tau\in [0, \tau_0)$ , then all the roots of (2.12) have negative real parts; if $\tau\in (\tau_j, \tau_{j+1}), j = 0, 1, 2\cdots$ , then (2.12) has at least one root with positive real parts.

From Lemma 2.4, Lemma 2.5, Lemma 2.6 and the Hopf bifurcation theorem, we have the following result of $E_1$ .

Theorem 2.7. For system (2.1), we have

($i$ ) if $a_{11}b_{12}+a_{12}b_{11}-a_{11}b_{13}\geq0$ , then the unique positive equilibrium $E_1$ is asymptotically stable for all $\tau\geq 0$ .

($ii$ ) if $a_{11}b_{12}+a_{12}b_{11}-a_{11}b_{13} < 0$ , then the equilibrium $E_1$ is asymptotically stable for $\tau\in [0, \tau_0)$ and unstable for $\tau>\tau_0$ . Hopf bifurcation occurs when $\tau = \tau_j, j = 0, 1, 2\cdots$ .

3.   Direction and stability of the hopf bifurcation

In this section, we consider the direction and stability of the bifurcating periodic solutions of system (2.1) induced by the Hopf bifurcation by using the normal form theory and the center manifold theorem by Hassard et al. (^[7]). We compute (see Appendix for details of the computation)

where the first four coefficients $g_{20}$ , $g_{11}$ , $g_{02}$ and $g_{21}$ that we need for determining the properties of the Hopf bifurcation are of the following forms

in which the terms $W_{11}^{(1)}(0)$ , $W_{11}^{(2)}(0)$ , $W_{11}^{(1)}(-1)$ , $W_{11}^{(2)}(-1)$ , $W_{20}^{(1)}(0)$ , $W_{20}^{(2)}(0)$ , $W_{20}^{(1)}(-1)$ , $W_{20}^{(2)}(-1)$ , $q_1$ , $q_1^*$ and $\bar{D}$ are calculated in Appendix.

Now using the four coefficients we obtain the values of the parameters $\mu_2$ , $\beta_2$ and $T_2$

Thus, using the quantities above, the properties of the Hopf bifurcation are determined by the following theorem.

Theorem 3.1. The properties of the Hopf bifurcation are determined by the parameters $\mu_2$ , $\beta_2$ and $T_2$ , where $\mu_2$ determines the direction of the Hopf bifurcation: if $\mu_2>0$ $(<0)$ , then the Hopf bifurcation is supercritical (subcritical); $\beta_2$ determines the stability of the bifurcating periodic solutions: if $\beta_2 < 0$ $(>0)$ , the bifurcating periodic solutions are stable (unstable); and $T_2$ determines the period of the bifurcating periodic solutions: if $T_2>0$ $(<0)$ , the period increases (decreases).

4.   Numerical simulation and discussion

We discussed a turbidostat system with Contois growth rate and time delay in this paper. We investigated the qualitative properties of system (1.2) including the existence and the stability of the equilibria and the bifurcations. In the case of no time delay, the boundary equilibrium $E_0$ is globally asymptotically stable if it is the unique equilibrium, and the positive equilibrium $E_1$ arising from a transcritical bifurcation is globally asymptotically stable if it exists. In the case of system with delay, the stability of $E_1$ is changed and the Hopf bifurcation occurs by choosing the delay as the bifurcation parameter. Furthermore, the stability and direction of the bifurcating periodic solutions are discussed by using the normal form and center manifold theorem. More concretely, when the delay is greater than the critical value $\tau_0$ , the positive equilibrium $E_1$ loses its stability and the stable periodic solution appears. The qualitative analysis and theoretical results show that the delay can change the topological structures of the system and produces more complicated dynamic behaviors.

We next offer an example to illustrate the feasibility of our results. When system (2.1) without delay, i.e. $\tau = 0$ , we set $a = 1$ and $k = 0.6$ , the system has the unique equilibrium $E_0$ , which is a stable node as $m = 0.5$ ; the system has a saddle $E_0$ and a stable node $E_1 = (0.76, 0.24)$ as $m = 1.5$ ; the system has the unique equilibrium $E_0$ , which is a saddle-node as $m = 1$ (Figure 2). However, when $m = 1.5$ , $a = 1$ and $k = 0.6$ , by choosing the time delay $\tau$ as the bifurcation parameter, we obtained the critical value of bifurcation $\tau_0\approx 7$ . Thus, the positive equilibrium $E_1$ is asymptotically stable when $\tau = 6 < \tau_0$ , which is supported by Figure 3. The positive equilibrium $E_1$ is unstable and a stable bifurcating periodic solution occurs from $E_1$ when $\tau = 8 >\tau_0$ , which can be seen clearly in Figure 4. According to (4.10) we can compute

Further, we can get

Then, in accordance with (4.22), we can obtain

Therefore, when $\tau = 8$ , $\mu_2>0$ and $\beta_2 < 0$ , then the Hopf bifurcation for system (2.1) is supercritical and the stable bifurcating periodic solutions occur from the positive equilibrium $E_1$ . From Figure 3 and Figure 4, we can find that the dynamics of system (2.1) change when $\tau$ locates near $\tau_0$ . The bifurcation diagram of $x-\tau$ is presented in Figure 5, from which we find that even for parameter values not chosen, the stable periodic solutions occur in a large region of time delay.

Remark 1. This is a remark about Figure 2. In fact, considering the biological background, the set $\Omega = \{(S, x): 0\leq S\leq1, x\geq0\}$ is positively invariant with respect to system (2.4). Further, $E_0$ is globally asymptotically stable with respect to $\Omega$ if $m < 1$ ; $E_0$ is globally attractive with respect to $\Omega$ if $m = 1$ ; and $E_1$ is globally asymptotically stable with respect to $\Omega$ if $m>1$ . These results can be proved easily by Liapunov-LaSalle invariant principle and Poincare-Bendixson Theorem respectively.

Acknowledgments

The author is very grateful to the anonymous referees for their careful reading and valuable suggestions, which have notably improved the quality of this paper.

Conflict of interest

All authors declare no conflicts of interest in this paper.

Appendix

We obtained in Section 2 that system (2.1) undergoes the Hopf bifurcation at the positive equilibrium $E_1$ when $\tau = \tau_j$ . In the following the properties of the Hopf bifurcation are determined. We first consider system (2.1) by the transformation $\tilde{y}_1(t) = S(\tau t),$ $\tilde{y}_2(t) = x(\tau t),$ $\tau = \tau_j+\mu$ and still denote $\tilde{y}_1(t), \tilde{y}_2(t)$ as $y_1(t), y_2(t)$ . Then system (2.1) is equivalent to a functional differential equation defined in $\mathcal{C} = C([-1, 0], R^2)$

where $y(t) = (y_1(t), y_2(t))^T \in R^2$ , and $L_\mu : \mathcal{C}\rightarrow R^2$ and $h : R\times \mathcal{C}\rightarrow R^2$ are given respectively by

and

where

and $\varphi = (\varphi_1, \varphi_2)^T\in \mathcal{C}.$

By the Riesz representation theorem, there exists a matrix whose components are bounded variation functions $\eta(\theta, \mu)$ for $\theta\in [-1, 0]$ such that

Indeed, we may choose

where $\delta$ is the Dirac delta function.

For $\varphi\in C^1 ([-1, 0], R^2)$ , we further define the operators $A$ and $B$ as

Then system (4.1) is equivalent to

where $y_t(\theta) = y(t+\theta)$ for $\theta\in[-1, 0]$ .

For $\psi \in C^1([0, 1], (R^2)^*)$ , we define $A^*$ of $A$ as

In order to normalize the eigenvectors of $A$ and $A^*$ , we define a bilinear inner product

where $\eta(\theta) = \eta(\theta, 0)$ .

From the above discussion we assume $q(\theta)$ and $q^*(s)$ are eigenvectors of $A$ and $A^*$ corresponding to $iw_0\tau_j$ and $-iw_0\tau_j$ . Suppose that $q(\theta) = (1, q_1)^Te^{iw_0\tau_j \theta}$ is the eigenvector of $A(0)$ corresponding to $iw_0\tau_j$ , then $Aq(0) = iw_0\tau_jq(0)$ . On the basis of the definitions of $A(0)$ , $L_\mu \varphi$ and $\eta(\theta, \mu)$ , the following conclusion can be achieved

We can get

Similarly, by the definition of $A^*$ , we have

and

We can then obtain the value of $D$ by $\langle q^*, q\rangle = 1$ . Further it holds from (4.3) that

Hence

In the following part, we compute the coordinates describing center manifold $\mathcal{C}_0$ at $\mu = 0$ by using the theory in ^[7]. Define

Then on the center manifold $\mathcal{C}_0$ , we have

where $z$ and $\bar{z}$ are local coordinates for center manifold $\mathcal{C}_0$ in the direction of $q^*$ and $\bar{q}^*$ . It is easy to see that $W$ is real if $y_t$ is real. Thereby, we next only consider real solutions of (4.2). If $y_t\in \mathcal{C}_0$ is the solution of (4.2), since $\mu = 0$ , which implies that

where

In the following we need to compute the coefficients $g_{20}$ , $g_{11}$ , $g_{02}$ and $g_{21}$ . Note that $y_t(\theta) = (y_{1t}(\theta), y_{2t}(\theta))^T = W(t, \theta)+zq(\theta)+\overline{z q(\theta)}$ and $q(\theta) = (1, q_1)^T e^{iw_0\tau_j \theta}$ . It then follows that

It follows together with the definition $h(\mu, \varphi)$ that

By substituting (4.9) into the above equation and comparing the coefficients with (4.8) we have

In order to determine $g_{21}$ we need to compute $W_{20}(\theta)$ and $W_{11}(\theta)$ . In the light of (4.2) and (4.7), we obtain

where

Expanding the previous series and comparing the coefficients, we obtain that

Case 1: We first consider the case $\theta\in[-1, 0)$ , it follows from (4.11) that

Comparing the coefficients with (4.8), we obtain

It follows from (4.13), (4.15) and the definition of $A$ that

Note that $q(\theta) = (1, q_{1})^T e^{iw_0\tau_j \theta}$ , then we have

where $E = (E^{(1)}, E^{(2)})^T\in R^2$ and $F = (F^{(1)}, F^{(2)})^T\in R^2$ are all constant vectors. In what follows, we will seek appropriate $E$ and $F$ .

Case 2: We now consider the case $\theta = 0$ . From (4.11), we have

It then follows from $H(z, \bar{z}, \theta)$ and (4.12) that

where

From (4.13) and the definition of $A$ , we obtain

where $\eta(\theta) = \eta(0, \theta)$ .

Substituting (4.17), (4.18) into (4.19) and noting that

and

We conclude that

that is

Thus

Similarly, we have

which enables us to assert that

Still now we can determine $W_{20}(\theta)$ and $W_{11}(\theta)$ from (4.17). Therefore, all $g_{ij}$ in (4.10) can be determined. Furthermore, we can compute the following values