text:
# Synergetic System Analysis for the Delay-Induced Hopf Bifurcation in the Wright Equation

## I Introduction

Within the last decades synergetics has provided powerful concepts and methods to describe selforganization processes in various branches of science Haken1 ; Haken2 ; Haken3 ; Haken4 ; Haken5 . The spontaneous formation of spatial, temporal or functional patterns in complex systems has been successfully investigated by working out general principles and by mapping them onto universal mathematical structures. The important result is due to the fact that in the vicinity of a dynamical instability the high-dimensional set of nonlinear evolution equations modelling a complex system on a microscopic or a mesoscopic scale can be approximately reduced to a low-dimensional set of order parameter equations describing the evolving pattern formation on a macroscopic scale. To obtain such a simplified, reduced description of selforganization processes, the synergetic system analysis proceeds as follows. A linearization of the evolution equations around a stationary solution shows that a dynamical instability is always accompanied by a time-scale hierarchy between numerous fast modes $`s`$ and few slow modes $`u`$. A rigorous treatment of the full nonlinear evolution equations in the vicinity of the dynamical instability leads to a characteristic interdependence between both hierarchy levels which may be illustrated by a circular causality chain. On the one hand the slaving principle of synergetics states that the numerous fast modes $`s`$ quasi-instantaneously take values which are prescribed by the few slow modes $`u`$ according to $`s(t)=h(u(t))`$ with the center manifold $`h(u)`$. On the other hand an adiabatic elimination of the fast enslaved modes $`s`$ yields equations for the slow order parameters $`u`$ which depend, in general, on the center manifold $`h(u)`$ due to the nonlinear feedback.

In its original formulation, the synergetic system analysis was developed for complex systems which can be modeled by ordinary and partial differential equations as well as their stochastic generalizations. Some time ago, the general concepts and methods of synergetics have been extended to delay differential equations to deal with dynamical instabilities which are induced by the finite propagation time of signals in feedback loops Wolfgang . Taking into account the infinite-dimensional character of a delay system Hale ; Krasovskii , the adiabatic elimination of the stable modes leads to a low-dimensional set of order parameter equations which turn out to be of the form of ordinary differential equations, i.e. they no longer contain memory effects. The predictions of the synergetic system analysis have been quantitatively tested by investigating the delay-induced Hopf bifurcation of the electronic system of a first-order phase-locked loop (PLL) Wolfgang . The periodic orbit which results from the corresponding order parameter equation near the bifurcation point has been confirmed by both a multiple scaling procedure and numerical simulations Michael1 ; Michael2 . Although this application exemplarily proves the order parameter concept for delay systems, it does not allow to draw conclusions about the slaving principle. As the lowest nonlinear term in the scalar delay differential equation of the PLL is a cubic one, the center manifold does not influence the order parameter equation of the Hopf bifurcation in the lowest order. In order to check for delay systems both ingredients of the circular causality chain, i.e. the order parameter concept and the slaving principle, it is thus indispensible to study a scalar delay differential equation with a quadratic nonlinearity.

A candidate for such a study is provided by

$`{\displaystyle \frac{d}{dt}}z(t)=R\left[z(t)z(t\tau )^2\right].`$ (1)
With vanishing time delay $`\tau `$ it represents a system which is named after the Belgian mathematician P.F. Verhulst from the 19th century verhulst1 . It is used as a simplified model for the population dynamics of a species in an environment with limited food supply verhulst2 . The synergetic system analysis for the Verhulst system with time delay (1) has been already performed in Ref. Michael1 . There also the well-known equation of Wright Wright

$`{\displaystyle \frac{d}{dt}}z(t)`$ $`=`$ $`Rz(t\tau )\left[1+z(t)\right]`$ (2)
has been treated, where $`R`$ denotes a system parameter and $`\tau `$ a delay time. This delay differential equation is mentioned by Wright Wright  as arising in the application of probability methods to the theory of asymptotic prime number density. Cunningham Cunningham  depicts it as a “growth equation” representing a mathematical description of a fluctuating population of organisms under certain environmental conditions. In addition, it may describe the operation of a control system working with potentially explosive chemical reactions, and quite similar equations arise in economic studies of business cycles. Performing an appropriate scaling of time

$`t=\tau t^{},z^{}(t^{})=z(\tau t^{})`$ (3)
converts the Wright equation (2) to its standard form with the control parameter

$`R^{}=\tau R.`$ (4)
Thus varying the delay time $`\tau `$ corresponds to changing the control parameter $`R^{}`$. By omitting the prime ’ for the respective quantities, the standard form of the Wright equation reads

$`{\displaystyle \frac{d}{dt}}z(t)`$ $`=`$ $`Rz(t1)\left[1+z(t)\right].`$ (5)
In this paper we restrict ourselves to analyse this standard form of the Wright equation.

The Wright equation (5) shows a delay-induced instability, namely a Poincaré-Andronov-Hopf bifurcation at the critical value

$`R_c={\displaystyle \frac{\pi }{2}}`$ (6)
of the control parameter $`R`$. In Ref. Marsden  it is shown that the oscillatory solution in the vicinity of this instability, i.e. the emerging limit cycle, can be calculated approximately using the method of averaging. This approximation reads in the lowest order

$`z(t)`$ $`=`$ $`A\sqrt{R{\displaystyle \frac{\pi }{2}}}\mathrm{cos}\left({\displaystyle \frac{\pi }{2}}t\right)+𝒪\left(R{\displaystyle \frac{\pi }{2}}\right),`$ (7)
where the amplitude $`A`$ has the value

$`A=\sqrt{{\displaystyle \frac{40}{3\pi 2}}}.`$ (8)
In Section II we start with a linear stability analysis of the Wright equation (5) which confirms, of course, the delay-induced Poincaré-Andronov-Hopf bifurcation when the control parameter $`R`$ approaches the critical value (6). Near this instability we perform a nonlinear synergetic treatment in Section III and study in detail how the center manifold influences the order parameter equation. In Section IV the resulting order parameter equation is transformed to the normal form of a Hopf bifurcation, where the emerging periodic orbit is determined one order higher than the lowest-order result (7) and (8). The numerical investigations of Section V confirm the emerging periodic orbit, furthermore we discuss the global bifurcation scenario of the Wright equation (5).

## II Linear Stability Analysis

The solution of the delay differential equation (5) for times $`t0`$ depends on the initial values of the function $`z(t)`$ in the entire interval $`[1,0]`$. In order to properly define such an initial value problem, Hale Hale  and Krasovskii Krasovskii  proposed to transform the equation of motion (5) for a function $`z(t)`$ in the usual state space $`\mathrm{\Gamma }`$ to the extended state space $`𝒞`$ of continuous complex valued functions $`z_t(\mathrm{\Theta })`$ which are defined on the interval $`[1,0]`$:

$`{\displaystyle \frac{d}{dt}}z_t(\mathrm{\Theta })`$ $`=`$ $`\left(𝒢z_t\right)(\mathrm{\Theta })=\{\begin{array}{cc}{\displaystyle \frac{d}{d\mathrm{\Theta }}}z_t(\mathrm{\Theta }),\hfill & 1\mathrm{\Theta }<0\\ [z_t],\hfill & \mathrm{\Theta }=0\end{array}.`$ (11)
Following the notation of Ref. Wolfgang , we introduced not only the new function $`z_t𝒞`$, which is connected to the original function $`z(t)\mathrm{\Gamma }`$ through

$`z_t(\mathrm{\Theta })=z(t+\mathrm{\Theta }),1\mathrm{\Theta }0,`$ (12)
but also the nonlinear functional

$`\left[z_t\right]`$ $`=`$ $`{\displaystyle \underset{k=1}{\overset{2}{}}}{\displaystyle _1^0}𝑑\mathrm{\Theta }_1\mathrm{}{\displaystyle _1^0}𝑑\mathrm{\Theta }_k\omega ^{(k)}(\mathrm{\Theta }_1,\mathrm{},\mathrm{\Theta }_k){\displaystyle \underset{l=1}{\overset{k}{}}}z_t(\mathrm{\Theta }_l)`$ (13)
with the two scalar densities

$`\omega ^{(1)}(\mathrm{\Theta }_1)`$ $`=`$ $`R𝜹(\mathrm{\Theta }_1+1),`$ (14)
$`\omega ^{(2)}(\mathrm{\Theta }_1,\mathrm{\Theta }_2)`$ $`=`$ $`R𝜹(\mathrm{\Theta }_1+1)𝜹(\mathrm{\Theta }_2).`$ (15)
The stationary states of this system

$`z_{\mathrm{stat}}^\mathrm{I}=0,z_{\mathrm{stat}}^{\mathrm{II}}=1`$ (16)
are candidates for the reference state from which we start our further investigations. For the main body of the article we focus our attention on the stationary state $`z_{\mathrm{stat}}^\mathrm{I}`$ and choose it as the reference state. The other stationary state $`z_{\mathrm{stat}}^{\mathrm{II}}`$ will be discussed in Section V together with the global bifurcation scenario of the Wright equation (5).

Then we linearize the system (11) with respect to the stationary state $`z_{\mathrm{stat}}^\mathrm{I}=0`$ by using the decomposition

$`z_t(\mathrm{\Theta })=z_{\mathrm{stat}}^\mathrm{I}+\zeta _t(\mathrm{\Theta }),1\mathrm{\Theta }0.`$ (17)
This leads to the following linearized equation of motion for the deviation $`\zeta _t(\mathrm{\Theta })`$ from the stationary state $`z_{\mathrm{stat}}^\mathrm{I}=0`$:

$`{\displaystyle \frac{d}{dt}}\zeta _t(\mathrm{\Theta })`$ $`=`$ $`\left(𝒢_L\zeta _t\right)(\mathrm{\Theta })=\{\begin{array}{cc}{\displaystyle \frac{d}{d\mathrm{\Theta }}}\zeta _t(\mathrm{\Theta }),\hfill & 1\mathrm{\Theta }<0\\ [\zeta _t],\hfill & \mathrm{\Theta }=0\end{array},`$ (20)
where the linear functional is given by

$`\left[\zeta _t\right]`$ $`=`$ $`{\displaystyle _1^0}𝑑\mathrm{\Theta }\omega (\mathrm{\Theta })\zeta _t(\mathrm{\Theta })`$ (21)
with the scalar density

$`\omega (\mathrm{\Theta })`$ $`=`$ $`{\displaystyle \frac{\delta [z_t]}{\delta z_t(\mathrm{\Theta })}}|_{z_t(\mathrm{\Theta })=z_{\mathrm{stat}}^\mathrm{I}}=R𝜹(\mathrm{\Theta }+1).`$ (22)
Inserting the solution ansatz

$`\zeta _t(\mathrm{\Theta })=\varphi ^\lambda (\mathrm{\Theta })e^{\lambda t},1\mathrm{\Theta }0`$ (23)
into (20) leads to the eigenvalue problem of the infinitesimal generator $`𝒢_L`$:

$`\lambda \varphi ^\lambda (\mathrm{\Theta })=\left(𝒢_L\varphi ^\lambda \right)(\mathrm{\Theta }),1\mathrm{\Theta }0.`$ (24)
Taking into account the definition of $`𝒢_L`$ in (20), the eigenfunction $`\varphi ^\lambda (\mathrm{\Theta })`$ is determined to be

$`\varphi ^\lambda (\mathrm{\Theta })=N_\lambda e^{\lambda \mathrm{\Theta }},1\mathrm{\Theta }0,`$ (25)
and the eigenvalue $`\lambda `$ follows from

$`\lambda =L(\lambda ),`$ (26)
where $`L(\lambda )`$ is defined by

$`L(\lambda )={\displaystyle \underset{1}{\overset{0}{}}}𝑑\mathrm{\Theta }\omega (\mathrm{\Theta })e^{\lambda \mathrm{\Theta }}.`$ (27)
Using the scalar density (22), we obtain the following transcendental characteristic equation:

$`Re^\lambda \lambda `$ $`=`$ $`0.`$ (28)
Thus the spectrum of the linear operator $`𝒢_L`$ has the following properties Hale :

* It consists of a countable infinite number of eigenvalues which cumulate for $`\mathrm{}(\lambda )\mathrm{}`$.
* It is confined by an upper threshold for the real parts of the eigenvalues.
* At the bifurcation point, i.e. the instability, some of the eigenvalues reach the imaginary axes and thus become unstable.

Further properties of the eigenvalues of the characteristic equation (28) follow from the Hayes theorem which can be found in Ref. Hayes . It states that all solutions of the transcendental equation

$`p+qe^\lambda \lambda =0`$ (29)
possess a negative real part if and only if (a) $`p<1`$ and (b) $`p<q<\sqrt{a_1^2+p^2}`$. Here $`a_1`$ represents the solution of the transcendental equation $`a_1=p\mathrm{tan}(a_1)`$ which lies in the interval $`[0,\pi )`$. For the special case $`p=0`$ one can show that $`a_1`$ is equal to $`\pi /2`$. The shaded region in Fig. 1 represents that region of the parameter space $`p,q`$ where both conditions of the Hayes theorem are fulfilled. The upper boundary line stems from (a) $`p<1`$ and (b1) $`p<q`$, whereas the lower boundary line follows from (a) $`p<1`$ and (b2) $`q<\sqrt{a_1^2+p^2}`$.

Comparing (28) with (29) we obtain the identification $`p=0`$ and $`q=R`$. Changing the control parameter $`R`$ from $`0`$ to $`\pi /2`$, the corresponding point in the parameter space $`p,q`$ moves along the $`q`$-axis from the point $`q=0`$ to $`q=\pi /2`$. At this critical value it reaches the boundary of the shaded stability region, i.e. no longer all solutions of the characteristic equation (28) have a negative real part. Therefore an instability occurs at $`R_c=\pi /2`$.

Fig. 2 confirms this result by illustrating the movement of the ten solutions of the characteristic equation (28) with the largest real part when the control parameter is increased from $`0`$ to $`\pi /2`$. The eigenvalues were obtained with a Newton algorithm and the control parameter $`R`$ was increased in equidistant steps. For $`R=0`$ there exists only one real eigenvalue $`0`$ as the linearized delay differential equation (20) degenerates to an ordinary differential equation. For $`R>0`$ this eigenvalue remains real and becomes negative. Furthermore, a countable infinite number of conjugate complex eigenvalues and another real eigenvalue emerge from an infinite negative real part. At the value $`R=1/\text{e}`$ both real eigenvalues meet at the point $`(1/0)`$. They are converted to a pair of conjugated complex eigenvalues for $`1/\text{e}<R<\pi /2`$. These two complex conjugated eigenvalues have zero real part at the instability $`R_c=\pi /2`$, which thus represents a Hopf bifurcation. We can further analyze this instability by introducing the smallness parameter

$`\epsilon ={\displaystyle \frac{RR_c}{R_c}}R=R_c(1+\epsilon )`$ (30)
for the deviation from the critical control parameter $`R_c=\pi /2`$. In particular, we can determine both eigenvalues $`\lambda _u^\pm (\epsilon )`$ with nearly vanishing real part at $`\epsilon 0`$ from the characteristic equation (28):

$`\lambda _u^\pm (\epsilon )={\displaystyle \frac{R_c^2}{1+R_c^2}}\epsilon \pm iR_c\left(1+{\displaystyle \frac{1}{1+R_c^2}}\epsilon \right)+𝒪\left(\epsilon ^2\right).`$ (31)
In the vicinity of the instability $`\epsilon 0`$ we read off from Fig. 2 that only the two eigenvalues (31) have nearly vanishing real part, all other eigenvalues have a large negative real part:

$`\mathrm{}[\lambda _u^\pm (\epsilon 0)]0;\mathrm{}[\lambda _s^j(\epsilon 0)]<0,j=1,\mathrm{},\mathrm{}.`$ (32)
This characteristic property of the linearized system (20) leads to a the time-scale hierarchy

$`T_u^\pm ={\displaystyle \frac{1}{\mathrm{}[\lambda _u^\pm (\epsilon 0)]}}T_s^j={\displaystyle \frac{1}{\mathrm{}[\lambda _s^j(\epsilon 0)]}},j=1,\mathrm{},\mathrm{}.`$ (33)
Thus the infinite-dimensional extended state space $`𝒞`$ decomposes in an 2-dimensional subspace $`𝒰`$ of the linear unstable modes and a remaining infinite-dimensional subspace $`𝒮`$ of the linear stable modes Wolfgang . As a consequence, the extended state function $`z_t(\mathrm{\Theta })`$ can be decomposed near the instability according to

$`z_t(\mathrm{\Theta })=z_{\mathrm{stat}}^\mathrm{I}+u_t(\mathrm{\Theta })+s_t(\mathrm{\Theta })=u_t(\mathrm{\Theta })+s_t(\mathrm{\Theta }),1\mathrm{\Theta }0,`$ (34)
as we have $`z_{\mathrm{stat}}^\mathrm{I}=0`$. Here $`u_t(\mathrm{\Theta })`$ and $`s_t(\mathrm{\Theta })`$ denote the respective contributions of $`z_t(\mathrm{\Theta })`$ in the subspaces $`𝒰`$ and $`𝒮`$. In order to project into these subspaces we need the linear unstable modes

$`\varphi ^{\lambda _u^\pm }(\mathrm{\Theta })=N_{\lambda _u^\pm }e^{\lambda _u^\pm \mathrm{\Theta }},1\mathrm{\Theta }0`$ (35)
of the system (20) which have been already determined in (25). However, their knowledge is not sufficient as the infinitesimal generator $`𝒢_L`$ is not selfadjoint. Therefore we need also the linear unstable modes

$`\psi ^{\lambda _u^\pm }(s)=N_{\lambda _u^\pm }e^{\lambda _u^\pm s},0s1`$ (36)
of the adjoint system

$`{\displaystyle \frac{d}{dt}}\zeta _t^{}(s)`$ $`=`$ $`\left(𝒢_L^{}\zeta _t^{}\right)(s)=\{\begin{array}{cc}{\displaystyle \frac{d}{ds}}\zeta _t^{}(s),\hfill & 0<s1\\ ^{}[\zeta _t^{}],\hfill & s=0\end{array},`$ (39)
where the linear functional

$`^{}\left[\zeta _t^{}\right]`$ $`=`$ $`{\displaystyle _0^1}𝑑s\omega (s)\zeta _t^{}(s)`$ (40)
contains also the scalar density (22). Indeed, the solution ansatz

$`\zeta _t^{}(s)=\psi ^\lambda (s)e^{\lambda s},0s1`$ (41)
converts (39) to the eigenvalue problem

$`\lambda \psi ^{}{}_{}{}^{\lambda }(s)=\left(𝒢_L^{}\psi ^\lambda \right)(s),0s1,`$ (42)
which is adjoint to (24). Note that $`\psi ^\lambda (s)`$ and $`\zeta _t^{}(s)`$ are elements of the dual extended state space $`𝒞^{}`$ which consists of continuous complex valued functions on the interval $`[0,1]`$. The relation between both extended state spaces $`𝒞`$ and $`𝒞^{}`$ is defined by the bilinear form Wolfgang

$`\left(\psi ^{}|\varphi \right)=\psi ^{}(0)\varphi (0){\displaystyle \underset{1}{\overset{0}{}}}𝑑\mathrm{\Theta }{\displaystyle \underset{0}{\overset{\mathrm{\Theta }}{}}}𝑑s\psi ^{}(s\mathrm{\Theta })\omega (\mathrm{\Theta })\varphi (s).`$ (43)
Using this bilinear form one can show that the eigenfunctions (35) and (36) are biorthonormal:

$`(\psi ^{}{}_{}{}^{\lambda _i}|\varphi {}_{}{}^{\lambda _j})=\delta _{ij},i,j=\pm .`$ (44)
This determines the yet unknown normalization constants to be

$`N_{\lambda _u^\pm }`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{1+\lambda _u^\pm }}},`$ (45)
so they reduce near the instability because of (31) to

$`N_{\lambda _u^\pm }`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{1\pm iR_c}}}+𝒪\left(\epsilon \right).`$ (46)
Furthermore, the bilinear form (43) allows to define the projector into the 2-dimensional subspace $`𝒰`$ of the unstable modes:

$`(𝒫_u)(\mathrm{\Theta })`$ $`=`$ $`{\displaystyle \underset{i=\pm }{}}\varphi ^{\lambda _u^i}(\mathrm{\Theta })\left(\psi ^{\lambda _u^i}|\right).`$ (47)
Correspondingly, the projector into the remaining infinite-dimensional subspace $`𝒮`$ of the stable modes reads

$`𝒫_s`$ $`=`$ $`(𝒫_u).`$ (48)
Applying the projector $`𝒫_u`$ to $`z_t𝒞`$ leads to $`u_t𝒰`$ according to

$`u_t(\mathrm{\Theta })=\left(𝒫_uz_t\right)(\mathrm{\Theta })={\displaystyle \underset{i=\pm }{}}u^i(t)\varphi ^{\lambda _u^i}(\mathrm{\Theta }),`$ (49)
where the amplitudes of the linear unstable modes $`\varphi ^{\lambda _u^\pm }(\mathrm{\Theta })`$ are defined by

$`u^\pm (t)=\left(\psi ^{\lambda _u^\pm }|z_t\right).`$ (50)
Later on, these amplitudes represent the order parameters which indicate the emergence of an instability. Analogously, the projector (48) leads to the stable modes

$`s_t(\mathrm{\Theta })=\left(𝒫_sz_t\right)(\mathrm{\Theta }),1\mathrm{\Theta }0.`$ (51)
## III Nonlinear Synergetic Analysis

After having performed a linear stability analysis around the reference state $`z_{\mathrm{stat}}^\mathrm{I}=0`$ in the vicinity of the instability $`R_c=\pi /2`$, we now return to our original nonlinear evolution equation (11) in the extended state space $`𝒞`$. We proceed by decomposing the generator $`𝒢`$ into its linear part $`𝒢_L`$ and a remaining effective nonlinear part:

$`{\displaystyle \frac{d}{dt}}z_t(\mathrm{\Theta })`$ $`=`$ $`\left(𝒢_Lz_t\right)(\mathrm{\Theta })+X_0(\mathrm{\Theta })^{\text{eff}}[z_t],1\mathrm{\Theta }0.`$ (52)
Here we introduced the scalar function

$`X_0(\mathrm{\Theta })`$ $`=`$ $`\{\begin{array}{cc}0,\hfill & 1\mathrm{\Theta }<0\\ 1,\hfill & \mathrm{\Theta }=0\end{array}`$ (55)
and the effective nonlinear functional

$`^{\text{eff}}[z_t]`$ $`=`$ $`{\displaystyle _1^0}𝑑\mathrm{\Theta }_1{\displaystyle _1^0}𝑑\mathrm{\Theta }_2\omega ^{(2)}(\mathrm{\Theta }_1,\mathrm{\Theta }_2)z_t(\mathrm{\Theta }_1)z_t(\mathrm{\Theta }_2)`$ (56)
with the scalar density (15). Using the projectors (47) and (48) and their properties (49)–(51), we can investigate the respective contributions of the order parameters $`u^\pm (t)`$ and the linear stable modes $`s_t𝒮`$ to the nonlinear dynamics (52). Thus we obtain the following system of coupled nonlinear mode equations:

$`{\displaystyle \frac{d}{dt}}u^\pm (t)`$ $`=`$ $`\lambda _u^\pm u^\pm (t)+\psi ^{}{}_{}{}^{\lambda _u^\pm }(0)^{\text{eff}}[{\displaystyle \underset{j=\pm }{}}\varphi {}_{}{}^{\lambda _u^j}u_{}^{j}(t)+s_t],`$ (57)
$`{\displaystyle \frac{d}{dt}}s_t(\mathrm{\Theta })`$ $`=`$ $`\left(𝒢_Ls_t\right)(\mathrm{\Theta })+\left(\left(𝒫_u\right)X_0\right)(\mathrm{\Theta })^{\text{eff}}\left[{\displaystyle \underset{j=\pm }{}}\varphi {}_{}{}^{\lambda _u^j}u_{}^{j}(t)+s_t\right].`$ (58)
It is still exact and describes completely the nonlinear dynamics. However, a solution to these equations can be only found by means of an approximation method. Such a well-established approximative solution is provided by the slaving principle of synergetics Haken1 ; Haken2 ; Haken3 ; Haken4 ; Haken5 .

To this end we start with the time-scale hierarchy (33) near the instability which leads to the fact that the dynamics of the stable modes $`s_t𝒮`$ evolves much faster than the order parameters $`u^\pm (t)`$. In Ref. Wolfgang  it has been shown for a quite general class of delay differential equations that such a time-scale hierarchy leads to a slaving of the stable modes, i.e. the numerous fast modes $`s_t𝒮`$ quasi-instantaneously take values which are prescribed by the few slow order parameters $`u^\pm (t)`$. In our context, the slaving principle states mathematically that the dynamics of the stables modes $`s_t𝒮`$ is determined by the center manifold $`h(\mathrm{\Theta },u^+,u^{})`$ according to

$`s_t(\mathrm{\Theta })=h(\mathrm{\Theta },u^+(t),u^{}(t)).`$ (59)
Inserting this ansatz in (58) leads to an implicit equation for the center manifold $`h(\mathrm{\Theta },u^+,u^{})`$:

$`{\displaystyle \underset{i=\pm }{}}{\displaystyle \frac{h(\mathrm{\Theta },u^+(t),u^{}(t))}{u^i(t)}}(\lambda _u^iu^i(t)+\psi ^{}{}_{}{}^{\lambda _u^i}(0)^{\text{eff}}[{\displaystyle \underset{j=\pm }{}}\varphi {}_{}{}^{\lambda _u^j}u_{}^{j}(t)+h])`$
$`=\left(𝒢_Lh\right)(\mathrm{\Theta })+\left(\left(𝒫_u\right)X_0\right)(\mathrm{\Theta })^{\text{eff}}\left[{\displaystyle \underset{j=\pm }{}}\varphi {}_{}{}^{\lambda _u^j}u_{}^{j}(t)+h\right].`$ (60)
It can be approximately solved in the vicinity of the instability as follows. We assume that the order parameters $`u^\pm (t)`$ possess a certain dependence on the smallness parameter (30) which is typical for a Hopf bifurctation:

$`u^\pm (t)=𝒪(\epsilon ^{1/2}).`$ (61)
Furthermore, we perform for the center manifold $`h(\mathrm{\Theta },u^+,u^{})`$ the lowest-order ansatz

$`h(\mathrm{\Theta },u^+,u^{})={\displaystyle \underset{j_1=\pm }{}}{\displaystyle \underset{j_2=\pm }{}}H_{j_1j_2}(\mathrm{\Theta })u^{j_1}(t)u^{j_2}(t),`$ (62)
as $`r=2`$ is the order of the effective nonlinear functional (56). From (61) and (62) follows then in lowest order of $`\epsilon `$ that the effective nonlinear functional $`^{\text{eff}}`$ in (60) can be approximated by

$`^{\text{eff}}\left[{\displaystyle \underset{j=\pm }{}}\varphi {}_{}{}^{\lambda _u^j}u_{}^{j}(t)+h\right]{\displaystyle \underset{j_1=\pm }{}}{\displaystyle \underset{j_2=\pm }{}}F_{j_1j_2}^{\mathrm{eff}}u^{j_1}(t)u^{j_2}(t),`$ (63)
where the coefficients $`F_{j_1j_2}^{\mathrm{eff}}`$ read

$`F_{j_1j_2}^{\text{eff}}={\displaystyle _1^0}𝑑\mathrm{\Theta }_1{\displaystyle _1^0}𝑑\mathrm{\Theta }_2\omega ^{(2)}(\mathrm{\Theta }_1,\mathrm{\Theta }_2)\varphi ^{\lambda _u^{j_1}}(\mathrm{\Theta }_1)\varphi ^{\lambda _u^{j_2}}(\mathrm{\Theta }_2).`$ (64)
Taking into account (15) and (35), these coefficients turn out to be

$`F_{++}^{\text{eff}}=F_{}^{\text{eff}}{}_{}{}^{}=RN_{\lambda _u^+}^2e^{\lambda _u^+},F_+^{\text{eff}}=F_+^{\text{eff}}=RN_{\lambda _u^+}N_{\lambda _u^{}}e^{\lambda _u^+}.`$ (65)
As a consequence, we conclude from (60) in lowest order of $`\epsilon `$ that the coefficients $`H_{j_1j_2}(\mathrm{\Theta })`$ of the center manifold (62) are given by

$`H_{j_1j_2}(\mathrm{\Theta })`$ $`=`$ $`F_{j_1j_2}^{\text{eff}}K_{j_1j_2}(\mathrm{\Theta }),`$ (66)
where the coefficients $`K_{j_1j_2}(\mathrm{\Theta })`$ follow from

$`K_{j_1j_2}(\mathrm{\Theta })=\left(\left[𝒢_L\text{λ}\right]^1\left(𝒫_uX_0X_0\right)\right)(\mathrm{\Theta })`$ (67)
with the abbreviation

$`\text{λ}={\displaystyle \underset{k=1}{\overset{2}{}}}\lambda _u^{j_k}.`$ (68)
In Ref. Wolfgang  it is shown that the operator $`\left[𝒢_L\text{λ}\right]^1`$ has the explicit representation

$`\left(\left[𝒢_L\text{λ}\right]^1\chi \right)(\mathrm{\Theta })={\displaystyle \underset{0}{\overset{\mathrm{\Theta }}{}}}𝑑se^{\text{λ}(\mathrm{\Theta }s)}\chi (s)+[L(\text{λ})\text{λ}]^1\left(\chi (0){\displaystyle \underset{1}{\overset{0}{}}}𝑑\mathrm{\Theta }{\displaystyle \underset{0}{\overset{\mathrm{\Theta }}{}}}𝑑se^{\text{λ}(\mathrm{\Theta }s)}\omega (\mathrm{\Theta })\chi (s)\right)e^{\text{λ}\mathrm{\Theta }},`$ (69)
where $`L(\lambda )`$ is already defined in (27). After some calculation, which also involves (22), (35), (36), (47), and (55), it thus follows that the coefficients (67) are given by

$`K_{j_1j_2}(\mathrm{\Theta })={\displaystyle \underset{j=\pm }{}}{\displaystyle \frac{N_{\lambda _u^j}^2}{\lambda _u^j\text{λ}}}e^{\lambda _u^j\mathrm{\Theta }}{\displaystyle \frac{e^{\text{λ}\mathrm{\Theta }}}{L(\text{λ})\text{λ}}}.`$ (70)
Thus we obtain together with (22), (27), and (68):

$`K_{++}(\mathrm{\Theta })=K_{}^{}(\mathrm{\Theta })`$ $`=`$ $`{\displaystyle \frac{N_{\lambda _u^+}^2e^{\lambda _u^+\mathrm{\Theta }}}{\lambda _u^+}}+{\displaystyle \frac{N_{\lambda _u^{}}^2e^{\lambda _u^{}\mathrm{\Theta }}}{\lambda _u^{}2\lambda _u^+}}+{\displaystyle \frac{e^{2\lambda _u^+\mathrm{\Theta }}}{Re^{2\lambda _u^+}+2\lambda _u^+}},`$ (71)
$`K_+(\mathrm{\Theta })=K_+(\mathrm{\Theta })`$ $`=`$ $`{\displaystyle \frac{N_{\lambda _u^+}^2e^{\lambda _u^+\mathrm{\Theta }}}{\lambda _u^{}}}{\displaystyle \frac{N_{\lambda _u^{}}^2e^{\lambda _u^{}\mathrm{\Theta }}}{\lambda _u^+}}+{\displaystyle \frac{e^{(\lambda _u^++\lambda _u^{})\mathrm{\Theta }}}{Re^{(\lambda _u^++\lambda _u^{})}+\lambda _u^++\lambda _u^{}}}.`$ (72)
This completes the lowest-order result for the center manifold $`h(\mathrm{\Theta },u^+,u^{})`$ which is given by (62), (65), (66), (71), and (72).

Thus we can now consider the order parameter equation (57). In lowest order in $`\epsilon `$ we take into account (15), (39), (56), (59), and (62) so that it reduces to

$`{\displaystyle \frac{d}{dt}}u^\pm (t)=\lambda _u^\pm u^\pm (t)RN_{\lambda _u^\pm }{\displaystyle \underset{l=1}{\overset{2}{}}}\left[{\displaystyle \underset{j=\pm }{}}\varphi ^{\lambda _u^j}(\vartheta _l)u^j(t)+{\displaystyle \underset{j_1=\pm }{}}{\displaystyle \underset{j_2=\pm }{}}H_{j_1j_2}(\vartheta _l)u^{j_1}(t)u^{j_2}(t)\right],`$ (73)
where we set

$`\vartheta _l=\{\begin{array}{cc}1,\hfill & l=1\\ 0,\hfill & l=2\end{array}.`$ (76)
Note that the order parameter equation (73) turns out to be an ordinary differential equation, i.e. it no longer contains memory effects. Furthermore, we observe that the center manifold explicitly enters the order parameter equation (73) as a direct consequence of the quadratic nonlinearity of the Wright equation (5). We remark that this effect, which is essential for the present synergetic analysis, was neglected in the neurophysiological study in Ref. Peter . In the subsequent section we show how the order parameter equation (73) is converted to the normal form of a Hopf bifurcation.

## IV Normal Form

Now we perform a nonlinear transformation of the order parameters which eliminates those terms which are irrelevant for the normal form of a Hopf bifurcation. As far as the so-called near identity transformation and the theory of normal forms in general is concerned, we refer to the Refs. Marsden ; NF1 ; NF2 . The terms in (73) which are relevant for the normal form of a Hopf bifurcation read

$`{\displaystyle \frac{d}{dt}}u^\pm (t)`$ $`=`$ $`\lambda _u^\pm u^\pm (t)+q_0^\pm u^\pm (t)^2+q_1^\pm u^\pm (t)u^{}(t)+q_2^\pm u^{}(t)^2+k_1^\pm u^\pm (t)^2u^{}(t),`$ (77)
as we can neglect quartic terms and nonresonant cubic terms due to the rotating wave approximation. The respective coefficients in Eq. (77) are given by

$`q_0^\pm `$ $`=`$ $`RN_{\lambda _u^\pm }\varphi ^{\lambda _u^\pm }(1)\varphi ^{\lambda _u^\pm }(0)=RN_{\lambda _u^\pm }^{}{}_{}{}^{3}e^{\lambda _u^\pm },`$ (78)
$`q_1^\pm `$ $`=`$ $`RN_{\lambda _u^\pm }\left[\varphi ^{\lambda _u^\pm }(1)\varphi ^{\lambda _u^{}}(0)+\varphi ^{\lambda _u^\pm }(0)\varphi ^{\lambda _u^{}}(1)\right]=RN_{\lambda _u^\pm }^{}{}_{}{}^{2}N_{\lambda _u^{}}\left(e^{\lambda _u^\pm }+e^{\lambda _u^{}}\right),`$ (79)
$`q_2^\pm `$ $`=`$ $`RN_{\lambda _u^\pm }\varphi ^{\lambda _u^{}}(1)\varphi ^{\lambda _u^{}}(0)=RN_{\lambda _u^\pm }N_{\lambda _u^{}}^{}{}_{}{}^{2}e^{\lambda _u^{}},`$ (80)
$`k_1^\pm `$ $`=`$ $`RN_{\lambda _u^\pm }\{\varphi ^{\lambda _u^\pm }(1)[H_+(0)+H_+(0)]+\varphi ^{\lambda _u^{}}(1)H_{++}(0)+`$ (81)
$`\varphi ^{\lambda _u^\pm }(0)[H_+(1)+H_+(1)]+\varphi ^{\lambda _u^{}}(0)H_{++}(1)\},`$
where we did not write down the explicit form of $`k_1^\pm `$ for simplicity. Then the previous order parameters $`u^\pm (t)`$ are transformed to new order parameters $`v^\pm (t)`$ by the near identity transformation

$`u^\pm (t)`$ $`=`$ $`v^\pm (t)+\alpha _0^\pm v^\pm (t)^2+\alpha _1^\pm v^\pm (t)v^{}(t)+\alpha _2^\pm v^{}(t)^2,`$ (82)
with the yet-unknown coefficients $`\alpha _0^\pm `$, $`\alpha _1^\pm `$, and $`\alpha _2^\pm `$. As the $`u^\pm (t)`$ are small quantities in the vicinity of the instability, the same holds for the $`v^\pm (t)`$. Inserting (82) in (77), we obtain a system of ordinary differential equations of the form

$`M(t){\displaystyle \frac{d}{dt}}\left(\begin{array}{c}v^+(t)\\ v^{}(t)\end{array}\right)=\left(\begin{array}{c}w^+(t)\\ w^{}(t)\end{array}\right),`$ (87)
where the matrix $`M(t)`$ is defined by

$`M(t)`$ $`=`$ $`\left(\begin{array}{cc}1+2\alpha _0^+v^+(t)+\alpha _1^+v^{}(t)& \alpha _1^+v^+(t)+2\alpha _2^+v^{}(t)\\ \alpha _1^{}v^{}(t)+2\alpha _2^{}v^+(t)& 1+2\alpha _0^{}v^{}(t)+\alpha _1^{}v^+(t)\end{array}\right).`$ (90)
For simplicity, we don’t write down the explicit form of $`w^+(t)`$ and $`w^{}(t)`$, but we note that they contain $`v^+(t)`$ and $`v^{}(t)`$ at least in first order. Thus we obtain from (87)

$`{\displaystyle \frac{d}{dt}}\left(\begin{array}{c}v^+(t)\\ v^{}(t)\end{array}\right)=M^1(t)\left(\begin{array}{c}w^+(t)\\ w^{}(t)\end{array}\right),`$ (95)
with the inverse matrix

$`M(t)^1`$ $`=`$ $`{\displaystyle \frac{1}{\text{Det}M(t)}}\left(\begin{array}{cc}M_{22}(t)& M_{12}(t)\\ M_{21}(t)& M_{11}(t)\end{array}\right),`$ (98)
where the determinant has the form

$`\text{Det}M(t)=1+v^+(t)\left(2\alpha _0^++\alpha _1^{}\right)+2v^+(t)v^{}(t)\left(\alpha _0^+\alpha _0^{}\alpha _2^+\alpha _2^{}\right)+2v^+(t)^2\left(\alpha _0^+\alpha _1^{}\alpha _1^+\alpha _2^{}\right)+c.c..`$ (99)
Expanding the right-hand side of (95) in powers of $`v^+(t)`$ and $`v^{}(t)`$ up to the third order, we yield

$`{\displaystyle \frac{d}{dt}}v^\pm (t)`$ $`=`$ $`\lambda ^\pm v^\pm (t)+\left(q_1^\pm \alpha _0^\pm \lambda ^\pm \right)v^\pm (t)^2+\left(q_0^\pm \alpha _1^\pm \lambda ^{}\right)v^+(t)v^{}(t)+\left[q_2^\pm +\alpha _2^\pm (\lambda ^\pm 2\lambda ^{})\right]v^{}(t)^2`$ (100)
$`+[k_1^\pm +q_0^\pm (\alpha _1^{}\alpha _0^\pm )q_0^{}\alpha _1^\pm +q_1^\pm \alpha _1^\pm +2q_2^\pm \alpha _2^{}2q_2^{}\alpha _2^\pm +\alpha _1^\pm \alpha _1^{}\lambda ^\pm +2\alpha _2^\pm \alpha _2^{}(2\lambda ^\pm \lambda ^{})`$
$`+\alpha _0^\pm \alpha _1^\pm (\lambda ^\pm +2\lambda ^{})]v^\pm (t)^2v^{}(t).`$
Now we can fix the yet-unknown coefficients $`\alpha _0^\pm `$, $`\alpha _1^\pm `$, and $`\alpha _2^\pm `$ of the near identity transformation (82) in such a way that all quadratic terms vanish. This leads to the conditions

$`\alpha _0^\pm ={\displaystyle \frac{q_0^\pm }{\lambda _u^\pm }},\alpha _1^\pm ={\displaystyle \frac{q_1^\pm }{\lambda _u^{}}},\alpha _2^\pm ={\displaystyle \frac{q_2^\pm }{2\lambda _u^{}\lambda _u^\pm }}.`$ (101)
Thus Eq. (100) reduces to the normal form of a Hopf bifurcation

$`{\displaystyle \frac{d}{dt}}v^\pm (t)=\lambda _u^\pm v^\pm (t)+b^\pm v^\pm (t)^2v^{}(t),`$ (102)
where the Hopf parameter $`b^\pm `$ is given by

$`b^\pm =k_1^\pm +{\displaystyle \frac{q_0^\pm q_1^\pm (4\lambda _{u}^{\pm }{}_{}{}^{2}\lambda _{u}^{}{}_{}{}^{2})+q_1^\pm q_1^{}(2\lambda _u^\pm \lambda _u^{}\lambda _{u}^{}{}_{}{}^{2})+2q_2^\pm q_2^{}\lambda _u^\pm \lambda _u^{}}{\lambda _u^\pm \lambda _u^{}(2\lambda _u^\pm \lambda _u^{})}}.`$ (103)
Taking into account (30), (31), (35), and (46) as well as (65), (66), (71), and (72) together with (78)–(81), this Hopf parameter $`b^\pm `$ reads in the vicinity of the instability:

$`b^\pm ={\displaystyle \frac{R_c}{5(1+R_c^2)^{\frac{3}{2}}}}\left[(3R_c1)\pm i(R_c+3)\right]+𝒪\left(\epsilon \right).`$ (104)
Performing the ansatz

$`v^\pm (t)`$ $`=`$ $`r(t)e^{\pm i\phi (t)},`$ (105)
the normal form (102) is transformed to polar coordinates

$`{\displaystyle \frac{d}{dt}}r(t)`$ $`=`$ $`r(t)\left[\mathrm{}\left(\lambda _u^\pm \right)+\mathrm{}\left(b^\pm \right)r(t)^2\right],`$ (106)
$`{\displaystyle \frac{d}{dt}}\phi (t)`$ $`=`$ $`\pm \left[\mathrm{}\left(\lambda _u^\pm \right)+\mathrm{}\left(b^\pm \right)r(t)^2\right].`$ (107)
Thus taking into account (31) and (104) near the instability, we result in the oscillatory solution

$`r_{\mathrm{stat}}`$ $`=`$ $`\sqrt{{\displaystyle \frac{\mathrm{}\left(\lambda _u^\pm \right)}{\mathrm{}\left(b^\pm \right)}}}=\sqrt{{\displaystyle \frac{5R_c}{3R_c1}}}\sqrt[4]{1+R_c^2}\sqrt{\epsilon }+𝒪\left(\epsilon \right),`$ (108)
$`{\displaystyle \frac{d}{dt}}\phi (t)`$ $`=`$ $`\pm \left[\mathrm{}\left(\lambda _u^\pm \right)+\mathrm{}\left(b^\pm \right)r_{\mathrm{stat}}^2\right]=R_c{\displaystyle \frac{R_c}{3R_c1}}\epsilon +𝒪\left(\epsilon ^2\right).`$ (109)
In order to compare this result with numerical simulations, we have to convert this oscillatory solution back to the original state space $`\mathrm{\Gamma }`$. At first we observe that we obtain for $`z_t𝒞`$ from (34), (49), (59), and (62) near the instability:

$`z_t(\mathrm{\Theta })={\displaystyle \underset{j=\pm }{}}\varphi ^{\lambda _u^j}(\mathrm{\Theta })u^j(t)+{\displaystyle \underset{j_1=\pm }{}}{\displaystyle \underset{j_2=\pm }{}}H_{j_1j_2}(\mathrm{\Theta })u^{j_1}(t)u^{j_2}(t).`$ (110)
Taking into account the near identity transformation (82) together with (35), this yields up to the first order in $`\epsilon `$:

$`z_t(\mathrm{\Theta })=N_{\lambda _u^+}e^{\lambda _u^+\mathrm{\Theta }}v^+(t)+N_{\lambda _u^{}}e^{\lambda _u^{}\mathrm{\Theta }}v^{}(t)+a_0(\mathrm{\Theta })v^+(t)^2+a_1(\mathrm{\Theta })v^+(t)v^{}(t)+a_2(\mathrm{\Theta })v^{}(t)^2,`$ (111)
where the coefficients $`a_0(\mathrm{\Theta })`$, $`a_1(\mathrm{\Theta })`$, $`a_2(\mathrm{\Theta })`$ read

$`a_0(\mathrm{\Theta })`$ $`=`$ $`\left[N_{\lambda _u^+}e^{\lambda _u^+\mathrm{\Theta }}\alpha _0^++N_{\lambda _u^{}}e^{\lambda _u^{}\mathrm{\Theta }}\alpha _2^{}+H_{++}(\mathrm{\Theta })\right],`$ (112)
$`a_1(\mathrm{\Theta })`$ $`=`$ $`\left[N_{\lambda _u^+}e^{\lambda _u^+\mathrm{\Theta }}\alpha _1^++N_{\lambda _u^{}}e^{\lambda _u^{}\mathrm{\Theta }}\alpha _1^{}+H_+(\mathrm{\Theta })+H_+(\mathrm{\Theta })\right],`$ (113)
$`a_2(\mathrm{\Theta })`$ $`=`$ $`\left[N_{\lambda _u^+}e^{\lambda _u^+\mathrm{\Theta }}\alpha _2^++N_{\lambda _u^{}}e^{\lambda _u^{}\mathrm{\Theta }}\alpha _0^{}+H_{}(\mathrm{\Theta })\right].`$ (114)
Due to the relation (12) between $`z(t)\mathrm{\Gamma }`$ and $`z_t𝒞`$, we conclude from (111):

$`z(t)=N_{\lambda _u^+}v^+(t)+N_{\lambda _u^{}}v^{}(t)+a_0(0)v^+(t)^2+a_1(0)v^+(t)v^{}(t)+a_2(0)v^{}(t)^2.`$ (115)
Near the instability we obtain from (46)

$`N_{\lambda _u^\pm }={\displaystyle \frac{1}{\sqrt[4]{1+R_c^2}}}e^{\pm i\psi _1}+𝒪(\epsilon )`$ (116)
with some phase $`\psi _1`$, whereas (105)–(107) leads to

$`v^\pm (t)=r_{\mathrm{stat}}e^{\pm i\phi (t)}`$ (117)
with the radius (108) and the phase

$`\phi (t)=\mathrm{\Omega }(\epsilon )t+\phi _0.`$ (118)
Here the frequency turns out to be

$`\mathrm{\Omega }(\epsilon )=R_c{\displaystyle \frac{R_c}{3R_c1}}\epsilon .`$ (119)
Furthermore, we yield from (112)–(114) by taking into account (65), (66), (71), (72), and (116) in the lowest order of $`\epsilon `$

$`a_0(0)={\displaystyle \frac{1}{\sqrt{5(1+R_c^2)}}}e^{i\psi _2},a_1(0)=0,a_2(0)={\displaystyle \frac{1}{\sqrt{5(1+R_c^2)}}}e^{i\psi _2},`$ (120)
where $`\psi _2`$ denotes some phase. Thus we obtain the following result for $`z(t)\mathrm{\Gamma }`$ near the instability

$`z(t)=c_0(\epsilon )+c_1(\epsilon )\mathrm{cos}\left[\phi (t)+\psi _1\right]+c_2(\epsilon )\mathrm{cos}\left[2\phi (t)+\psi _2\right]+𝒪\left(\epsilon ^{\frac{3}{2}}\right),`$ (121)
where the respective coeffcients read

$`c_0(\epsilon )=0,c_1(\epsilon )=2\sqrt{{\displaystyle \frac{5R_c}{3R_c1}}}\sqrt{\epsilon },c_2(\epsilon )=2{\displaystyle \frac{\sqrt{5}R_c}{3R_c1}}\epsilon .`$ (122)
Now we compare the oscillatory solution (7), (8), which was obtained from using the method of averaging, with ours (118), (119), (121), (122) by taking into account the critical value (6) of the control parameter. We conclude that both results coincide in the lowest-order $`\epsilon ^{1/2}`$, but our result is even correct up to the order $`\epsilon `$.

From the near identity transformation (82) as well as from (108) and (117) we conclude that the order parameters $`u^\pm (t)`$ turn out to be of the order $`\epsilon ^{1/2}`$. This result is consistent with our original assumption (61) which was the basis of our approximate solution of the implicit equation for the center manifold (60) in the vicinity of the instability. Thus our synergetic system analysis is justified a posteriori by self consistency.

Note that the same perturbative result (118), (119), (121), (122) for the oscillatory solution above the Hopf bifurcation can be derived with the multiple scaling method Michael1 . It represents a technical procedure to deduce the normal form, once the bifurcation type is known, by using the knowledge how the respective quantities depend on the smallness parameter $`\epsilon =(RR_c)/R_c`$. Although the multiple scaling method has been originally developed for ordinary differential equations ms1 ; ms2 ; ms3 , it can be also applied to delay differential equations (see, for instance, the treatment in Ref. elena ).

## V Numerical Investigation

In order to numerically verify our analytical result, we integrated the underlying delay differential equation of Wright (5). By doing so, we varied the control parameter $`R`$ in the vicinity of the instability $`R_c=\pi /2`$ in such a way that the smallness parameter $`\epsilon =(RR_c)/R_c`$ took 200 equidistant values between $`10^5`$ and $`10^1`$. We used a Runge-Kutta-Verner method of the IMSL library as an integration routine and performed a linear interpolation between the respective values. In particular in the immediate vicinity of the instability the phenomenon of critical slowing down led to a transient behavior. To exclude this, we iterated the discretized delay differential equation for each value of the control parameter at least $`10^6`$ times. Afterwards we calculated the power spectrum with a complex FFT so that the basic frequency $`\mathrm{\Omega }`$ of the oscillatory solution could be determined with high resolution. Then we performed a real FFT with the period $`T=2\pi /\mathrm{\Omega }`$ of the simulated periodic signal $`z(t)=z(t+T)`$:

$`z(t)={\displaystyle \frac{a_0}{2}}+{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}\left[a_k\mathrm{cos}\left(k\mathrm{\Omega }t\right)+b_k\mathrm{sin}\left(k\mathrm{\Omega }t\right)\right].`$ (123)
The Fourier coefficients follow from integrations with respect to one period $`T=2\pi /\mathrm{\Omega }`$:

$`a_k={\displaystyle \frac{2}{T}}{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}𝑑tf(t)\mathrm{cos}\left(k\mathrm{\Omega }t\right),k=0,1,\mathrm{},\mathrm{};b_k={\displaystyle \frac{2}{T}}{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}𝑑tf(t)\mathrm{sin}\left(k\mathrm{\Omega }t\right),k=1,\mathrm{},\mathrm{}.`$ (124)
From (123) follows then the spectral representation

$`z(t)=c_0+{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}c_k\mathrm{cos}\left(k\mathrm{\Omega }t+\varphi _k\right)`$ (125)
with the quantities

$`c_0={\displaystyle \frac{a_0}{2}},c_k=\sqrt{a_k^2+b_k^2},\varphi _k=\text{artan}{\displaystyle \frac{b_k}{a_k}},k=1,\mathrm{},\mathrm{}.`$ (126)
Thus our analytical result (118), (121) can be interpreted as the first terms within a spectral representation (125), where the frequency $`\mathrm{\Omega }=2\pi /T`$ and the Fourier coefficients $`c_0`$, $`c_1`$, $`c_2`$ are given by (119) and (122). Numerically analysing the Hopf bifurcation with the FFT, the results for $`\mathrm{\Omega }`$, $`c_0`$, $`c_1`$, $`c_2`$ are plotted in Fig. 3 versus the smallness parameter $`\epsilon `$. Comparing the respective numerical and analytical results, we observe some deviations for small and for large values of the smallness parameter $`\epsilon `$. The former are due to the phenomenon of critical slowing down, i.e. the system stays longer in the transient state when the instability is approached, and the latter arise from the neglected higher order corrections in the analytical approach. Therefore we restricted our numerical analysis to the intermediate interval $`[10^5,10^1]`$ of the smallness parameter $`\epsilon `$. In Tab. 1 we see that the analytical and numerical determined quantities agree quantitatively very well. Thus our synergetic system analysis for the delay-induced Hopf bifurcation in the Wright equation is numerically verified.

For the sake of completeness we have also investigated oscillatory solutions for values of the control parameter $`R`$ which are larger than the critical one $`R_c=\pi /2`$. Figure 4 shows that all these periodic solutions oscillate around the stationary state $`z_{\mathrm{stat}}^\mathrm{I}=0`$ which becomes unstable at $`R_c=\pi /2`$. It turns out that a global bifurcation occurs for $`R_c^g=3.247`$ as then the oscillatory solution comes close to the other stationary state $`z_{\mathrm{stat}}^{\mathrm{II}}=1`$ which turns out to be linear unstable for all values of the control parameter $`R>0`$. Indeed, performing a linear stability analysis according to Section II around the stationary state $`z_{\mathrm{stat}}^{\mathrm{II}}=1`$ leads to the characteristic equation

$`R\lambda =0,`$ (127)
so we have from Eq. (28) the identification $`p=R`$ and $`q=0`$ (compare the shaded stability region in Fig. 1).

## VI Summary and Outlook

In this article a linear stability analysis of the Wright equation (5) around the stationary state $`z_{\mathrm{stat}}^\mathrm{I}=0`$ showed that a delay-induced Hopf bifurcation occurs at the critical value $`R_c=\pi /2`$ of the control parameter $`R`$. Within a subsequent nonlinear synergetic analysis we adiabatically eliminated the stable modes and derived the normal form of this Hopf bifurcation. It is explicitly influenced by the center manifold in the lowest order, as the Wright equation (5) has a quadratic nonlinearity. Solving the normal form we obtained a periodic solution above the Hopf bifurcation which was numerically verified.

In contrast to the corresponding analysis of the electronic system of a first-order phase-locked loop with time delay Michael2 , this paper confirms not only the order parameter concept for delay systems. It also represents a successful test for the slaving principle of synergetics, i.e. for the influence of the center manifold on the order parameter equations. Thus the validity of the circular causality chain of synergetics (see Fig. 5) has been demonstrated for the Hopf bifurcation of a delay differential equation.

It remains to investigate the circular causality chain also for other bifurcations. The Floquet theory for delay differential equations and thus the linear stability analysis for a periodic reference state was already established in the Refs. Christian1 ; Christian2 . However, a corresponding synergetic system analysis is still missing which derives the order parameter equations and the normal forms for bifurcations of oscillatory solutions Michael3 .

## Acknowledgements

We are thankful to Michael Bestehorn and Rudolf Friedrich for contributing various useful comments at an initial stage of this work. Furthermore we thank Hermann Haken and Arne Wunderlin for teaching us synergetics for many years. Finally, A.P. is grateful for the hospitality of Günter Wunner at the I. Institute for Theoretical Physics of the University of Stuttgart as this article was finished there.


format: text

text:
# Initial conditions for quintessence after inflation

## I Introduction

In order to fit the impressive array of observational data now available, it is widely believed that a model of the Universe should feature a period of inflation in its early history, leading to the generation of density perturbations via quantum fluctuations, and that the present expansion rate should be accelerating. While the latter phenomenon is often considered due to a cosmological constant, it is at least as attractive to presume that the acceleration is driven by the same mechanism usually exploited to give early Universe inflation, namely domination by the potential energy of a scalar field. Such models are known as quintessence models.

An important class of quintessence models are known as tracking models track ; RP ; nucl2 ; scal ; equip , where the late-time evolution of the field has an attractor behaviour rendering it fairly independent of initial conditions. However, despite the existence of tracking behaviour the details of the initial conditions may yet be important. For example, in several models of quintessence the tracking solution cannot be achieved too early; big bang nucleosynthesis is spoiled if the quintessence field has too large a density at that time nucl2 ; nucleo . More obviously, a quintessence scenario will not work unless the initial energy density is at least as large as is required by the present. It is therefore useful to have further guidance as to the likely initial conditions for quintessence.

The simplest assumption concerning quintessence is that it is a fundamental scalar field (rather than a low-energy composite field), and as such was already present during the inflationary epoch in the early Universe. If the quintessence field is sufficiently weakly coupled that it is not affected by the inflaton decays ending inflation, its possible initial conditions are restricted by its dynamics during the inflationary period. Note that we are not considering the situation where the quintessence field and the inflaton are the same field quintinf .

In this paper we carry out a comprehensive study of the influence of inflationary dynamics on the quintessence field. We choose parameters so that the quintessence sector matches the observed acceleration of the Universe, and the inflaton sector generates suitable perturbations to initiate structure formation. We investigate both classical and quantum dynamics of the quintessence. We assume a flat Universe throughout.

## II Models and observational constraints

### II.1 Quintessence

The quintessence model is defined by a scalar field $`Q`$ evolving in a potential $`V(Q)`$. Constraints on such models from observation have recently been considered by various authors obscon ; BHM02 ; the precise details are not important for our considerations but for definiteness we use the results of Bean et al. BHM02 , who recently combined constraints on such models from type Ia supernovae, CMB peak positions and large-scale structure surveys. They found one-sigma constraints on the quintessence density parameter and its pressure–density ratio of

$$0.57\mathrm{\Omega }_Q0.72\mathrm{and}1w_Q0.85.$$
(1)
With this, and assuming the Hubble constant to be given by $`h=0.72\pm 0.08`$ HKP , it is possible to find constraints on the parameters of the quintessence models. To locate viable regions of parameter space, we approximate these three constraints as Gaussian distributed (in the case of $`w_Q`$ a half-Gaussian centred at $`w_Q=1`$) and independent. This defines a confidence region in the space of parameters, and using numerical simulations we can translate the confidence region into other variables such as the density of non-relativistic matter today and the parameters of the quintessence potential.

Our main discussion will focus on quintessence models with inverse power-law potentials, as originally introduced by Ratra and Peebles RP

$$V(Q)=V_0\left(\frac{Q}{m_{\mathrm{Pl}}}\right)^\beta .$$
(2)
In order to simplify some later expressions we define the notation $`\stackrel{~}{X}=X/m_{\mathrm{Pl}}`$, where $`X`$ is any field or variable. Using the procedure described above and projecting the result onto the $`\beta `$$`V_0`$ plane gives the allowed region shown in Fig. 1. We see that $`\beta 1`$ and $`V_05\times 10^{124}m_{\mathrm{Pl}}^4`$ at 68% confidence, with a broader region at 95% confidence. In defining our models, we allow $`\beta `$ to vary and for each choice fix $`V_0`$ at its best-fit value.

### II.2 Inflation

We take inflation to be driven by a scalar field $`\varphi `$ (the inflaton) in a potential $`W(\varphi )`$. Because of the small value of $`V_0`$, the inflaton dominates the evolution and therefore we can use the standard constraints on its parameters. The inflation model is to be chosen so as to reproduce the amplitude of perturbations needed to generate the observed structures, which requires LL

$$\frac{512\pi }{75}\frac{W^3}{m_{\mathrm{Pl}}^6\left|\mathrm{d}W/\mathrm{d}\varphi \right|^2}4\times 10^{10}.$$
(3)
This is to be evaluated at $`\varphi _{\mathrm{obs}}`$, the value of the field corresponding to the epoch when the scale of the observable Universe crossed the horizon.

In this paper we will consider potentials where inflation ends by violation of the slow-roll conditions, so the slow-roll approximation

$`\epsilon (\varphi )`$ $`={\displaystyle \frac{m_{\mathrm{Pl}}^2}{16\pi }}\left({\displaystyle \frac{\mathrm{d}W/\mathrm{d}\varphi }{W}}\right)^21`$ (4)
$`\eta (\varphi )`$ $`={\displaystyle \frac{m_{\mathrm{Pl}}^2}{8\pi }}{\displaystyle \frac{\mathrm{d}^2W/\mathrm{d}\varphi ^2}{W}}1`$ (5)
gives the end of inflation and the number of $`e`$-foldings

$$N\frac{8\pi }{m_{\mathrm{Pl}}^2}_{\varphi _{\mathrm{end}}}^{\varphi _{\mathrm{obs}}}\frac{W(\varphi )}{\mathrm{d}W/\mathrm{d}\varphi }d\varphi 50$$
(6)
gives us $`\varphi _{\mathrm{obs}}`$.

Fifty or so $`e`$-foldings is the minimum amount of inflation capable of giving our observed Universe, but crucial for the considerations of this paper is that typically one expects enormously more inflation. If one only imposes the constraint that the initial energy density be below the Planck energy, most models allow very large amounts of inflation. However the early stages of such inflation may be dominated by quantum dynamics (the so-called stochastic inflation regime SI ). The maximum value of the inflaton field we consider, $`\varphi _{\mathrm{max}}`$, is given by the quantum-to-classical evolution transition

$$\frac{\mathrm{\Delta }_{\mathrm{qu}}}{\mathrm{\Delta }_{\mathrm{cl}}}=\frac{H}{2\pi }\left|\frac{H}{\dot{\varphi }}\right|1,$$
(7)
where the $`\mathrm{\Delta }`$ represent the typical quantum and classical evolution per Hubble time.

For simplicity we focus on a inflaton model with a power-law potential

$$W(\varphi )=W_0\left(\frac{\varphi }{m_{\mathrm{Pl}}}\right)^\alpha .$$
(8)
The former constraints lead to

$`\stackrel{~}{\varphi }_{\mathrm{end}}`$ $``$ $`\mathrm{Max}[{\displaystyle \frac{\alpha }{4\sqrt{\pi }}},\sqrt{{\displaystyle \frac{\alpha (\alpha 1)}{8\pi }}}],`$ (9)
$`\stackrel{~}{\varphi }_{\mathrm{obs}}`$ $``$ $`\sqrt{{\displaystyle \frac{50\alpha }{4\pi }}+\stackrel{~}{\varphi }_{\mathrm{end}}^2},`$ (10)
$`\stackrel{~}{\varphi }_{\mathrm{max}}`$ $``$ $`\left({\displaystyle \frac{3\alpha ^2m_{\mathrm{Pl}}^4}{128\pi W_0}}\right)^{1/(\alpha +2)},`$ (11)
$`W_0`$ $``$ $`{\displaystyle \frac{3\times 10^8\alpha ^2m_{\mathrm{Pl}}^4}{512\pi }}\stackrel{~}{\varphi }_{\mathrm{obs}}^{(\alpha +2)}.`$ (12)
## III Quintessence field evolution during inflation

The models are now fully defined and we can study the dynamics of each field, beginning with the classical evolution.

### III.1 Classical evolution

We assume that the quintessence field is completely uncoupled from everything else, including the inflaton. The equations of motion are the Euler–Lagrange equations for the inflaton and quintessence fields

$$\ddot{\varphi }=3H\dot{\varphi }\frac{\mathrm{d}W}{\mathrm{d}\varphi };\ddot{Q}=3H\dot{Q}\frac{\mathrm{d}V}{\mathrm{d}Q}$$
(13)
and the Einstein equation

$$H^2=\frac{8\pi }{3m_{\mathrm{Pl}}^2}\left[\frac{1}{2}\dot{\varphi }^2+W(\varphi )+\frac{1}{2}\dot{Q}^2+V(Q)\right].$$
(14)
Because of the very small value of $`V_0`$, the evolution of the Universe is determined by the inflaton. The quintessence field may dominate if its value is very small, but, as we will see this situation would end very quickly. Therefore, using the slow-roll approximation we have the solution

$$\stackrel{~}{\varphi }(N)\sqrt{\stackrel{~}{\varphi }_{\mathrm{ini}}^2\frac{\alpha N}{4\pi }},$$
(15)
where $`N=\mathrm{log}(a/a_{\mathrm{ini}})`$. The subscript “ini” always means the value at some time $`t_{\mathrm{ini}}`$. Knowing the dynamics of the Universe we can now study the quintessence field.

The classical behaviour of the quintessence field is rather straightforward. If the field is in a region where $`\left|\mathrm{d}V/\mathrm{d}Q\right|`$ is small, its classical evolution will be highly suppressed by the friction arising from the inflationary expansion and to a good approximation the field will retain its initial value. Conversely, if $`\left|\mathrm{d}V/\mathrm{d}Q\right|`$ is large, the field will quickly roll down until the potential becomes too flat. More precisely, using the slow-roll approximation for the quintessence field we have

$$\dot{\stackrel{~}{Q}}\frac{\beta V_0}{3H(t)m_{\mathrm{Pl}}^2}\stackrel{~}{Q}^{\beta 1}.$$
(16)
Assuming $`H(t)`$ constant (roughly correct for a few Hubble times) the solution is

$$\stackrel{~}{Q}(t)\left[\frac{\beta (\beta +2)V_0}{3Hm_{\mathrm{Pl}}^2}(tt_{\mathrm{ini}})+\stackrel{~}{Q}_{\mathrm{ini}}^{\beta +2}\right]^{1/(\beta +2)}.$$
(17)
We define

$$\stackrel{~}{Q}_{\mathrm{min}}=\left[\frac{\beta (\beta +2)V_0}{3H^2m_{\mathrm{Pl}}^2}\right]^{1/(\beta +2)}$$
(18)
and see that the evolution during a Hubble time has the following approximate description:

$$\begin{array}{ccc}\hfill \stackrel{~}{Q}(t_{\mathrm{ini}}+H^1)& \stackrel{~}{Q}_{\mathrm{ini}}\hfill & \mathrm{if}\stackrel{~}{Q}_{\mathrm{ini}}\stackrel{~}{Q}_{\mathrm{min}}\hfill \\ \hfill \stackrel{~}{Q}(t_{\mathrm{ini}}+H^1)& \stackrel{~}{Q}_{\mathrm{min}}\hfill & \mathrm{if}\stackrel{~}{Q}_{\mathrm{ini}}\stackrel{~}{Q}_{\mathrm{min}}\hfill \end{array}.$$
(19)
Therefore, after a few Hubble times the classical evolution of the quintessence field is to remain constant.

### III.2 Quantum fluctuations

Although the large friction from the Hubble expansion renders the classical evolution of the quintessence field negligible, the same is not necessarily true of the quantum fluctuations. Indeed, as we shall see, the effect of these in the quintessence potential is to drive typical regions of the Universe to quite large values of the quintessence field, corresponding to low energy densities. This possibility was first mentioned in Ref. SS .

Let us now recall briefly the way we treat these fluctuations. Following Ref. LL , we choose the spatially-flat slicing and we split the quintessence field into an unperturbed part and a perturbation

$$Q(𝐱,t)=Q(t)+\delta Q(𝐱,t).$$
(20)
We quantize the perturbation and expand its Fourier components

$$\delta \widehat{Q}_𝐤(t)=w_k(t)\widehat{a}_𝐤+w_k^{}(t)\widehat{a}_𝐤^{},$$
(21)
where $`\widehat{a}_𝐤`$ is the annihilation operator. In the linear approximation we have to solve

$$\ddot{w}_k+3H\dot{w}_k+\left[(k/a)^2+m_{\mathrm{eff}}^2\right]w_k=0,$$
(22)
where $`m_{\mathrm{eff}}^2=\mathrm{d}^2V/\mathrm{d}Q^2`$. If the field is effectively massless and $`H`$ is approximately constant, the power spectrum of the fluctuations is given by

$$𝒫_k(t)=\frac{L^3k^3}{2\pi ^2}\left|w_k(t)\right|^2\left[1+\left(\frac{k}{aH}\right)^2\right]\left(\frac{H}{2\pi }\right)^2.$$
(23)
This means that before horizon-crossing $`𝒫_k(t)1/a^2(t)`$ and afterwards it becomes constant. Using

$$\delta \widehat{Q}^2(𝐱,t)=_0^{\mathrm{}}𝒫_k(t)\frac{\mathrm{d}k}{k},$$
(24)
and the fact that the fluctuations become classical after they cross the horizon, we find that the quintessence field receives quantum kicks whose size is estimated as $`H/2\pi `$ per Hubble time.

The effect of these quantum fluctuations can be studied using the Fokker–Planck formalism SI2 , which allows one to follow the probability distribution $`f(Q,t)`$ of the quintessence field during inflation. A simple derivation is as follows. We take discrete steps $`\mathrm{\Delta }t=H^1\epsilon `$ during which there are a random jump of magnitude $`\mathrm{\Delta }Q_{\mathrm{qu}}=\sqrt{\mathrm{\Delta }t/H^1}H/2\pi =\sqrt{\mathrm{\Delta }\epsilon }H/2\pi `$ due to quantum fluctuations and a classical step $`\mathrm{\Delta }Q_{\mathrm{cl}}=\dot{Q}\mathrm{\Delta }t=\epsilon V^{}/3H^2`$, where the prime denotes a derivative with respect to $`Q`$. During this time an interval $`\mathrm{d}Q`$ shrinks by a factor $`1+\epsilon V^{\prime \prime }/3H^2`$. We find the following equation:

$$\begin{array}{cc}\hfill f(Q,t+\mathrm{\Delta }t)\mathrm{d}Q=& \frac{1}{2}[f(Q\mathrm{\Delta }Q_{\mathrm{cl}}+\mathrm{\Delta }Q_{\mathrm{qu}},t)\hfill \\ & +f(Q\mathrm{\Delta }Q_{\mathrm{cl}}\mathrm{\Delta }Q_{\mathrm{qu}},t)]\hfill \\ & \times (1+\epsilon V^{\prime \prime }/3H^2)\mathrm{d}Q.\hfill \end{array}$$
(25)
Subtracting $`f(Q,t)\mathrm{d}Q`$, dividing by $`\epsilon `$ and taking the limit $`\epsilon 0`$ leads to

$$\begin{array}{cc}\hfill \frac{f(Q,t)}{t}=& \frac{H^3(t)}{8\pi ^2}\frac{^2f(Q,t)}{Q^2}\hfill \\ & +\frac{1}{3H(t)}\frac{}{Q}\left(\frac{V(Q)}{Q}f(Q,t)\right).\hfill \end{array}$$
(26)
On the right-hand side, the first term produces diffusion and the second one produces a drift. Note that unlike the usual stochastic inflation situation, here $`H(t)`$ is determined externally by the inflaton field evolution, and the fluctuations in $`Q`$ do not back-react on the expansion rate.

To be complete we need the boundary condition

$$\underset{Q0}{lim}\left[\frac{H^3(t)}{8\pi ^2}\frac{f(Q,t)}{Q}+\frac{V(Q)/Q}{3H(t)}f(Q,t)\right]=0,$$
(27)
the requirement that $`f(Q,t)0`$, and normalization $`_0^{\mathrm{}}f(Q,t)dQ=1`$ which implies $`lim_Q\mathrm{}f(Q,t)=0`$.

From now on we focus on the inverse power-law potential for the quintessence. First of all, one may wonder whether the massless condition we used for the typical quantum jump is justified. Indeed, for values smaller than $`Q_{\mathrm{lim}}`$, where

$$\stackrel{~}{Q}_{\mathrm{lim}}=\left[\frac{\beta (\beta +1)V_0}{H^2m_{\mathrm{Pl}}^2}\right]^{1/(\beta +2)},$$
(28)
the quintessence field has an effective mass $`m_{\mathrm{eff}}`$ which is no longer negligible compared to the Hubble rate.111For power-law potentials, this field value is almost the same as $`Q_{\mathrm{min}}`$, the value below which the classical evolution is important. For all wavelengths $`a/k`$ bigger than $`1/m_{\mathrm{eff}}`$, Eq. (22) implies

$$w_k(t)a^{3/2}(t)$$
(29)
and therefore the power spectrum of $`P_k(t)`$ decreases as $`1/a^3(t)`$ leading to a smaller value of the quantum jump after horizon-crossing. Nevertheless, because the typical quantum jump $`H/2\pi `$ (typically $`10^5m_{\mathrm{Pl}}`$) is dramatically bigger than $`Q_{\mathrm{lim}}`$ (with $`\beta =1`$, typically $`10^{39}m_{\mathrm{Pl}}`$), the massless condition is broken only during very short periods of time (much less than a Hubble time) and so the massless approximation is an excellent one.

The Fokker–Planck equation becomes

$$\begin{array}{cc}\hfill \frac{f(Q,t)}{t}=& \frac{H^3(t)}{8\pi ^2}\frac{^2f(Q,t)}{Q^2}\hfill \\ & \frac{}{Q}\left(\frac{\beta V_0m_{\mathrm{Pl}}^\beta }{3H(t)}Q^{\beta 1}f(Q,t)\right).\hfill \end{array}$$
(30)
For $`Q>Q_{\mathrm{min}}`$ we can drop the terms coming from the classical evolution since we have seen it is negligible. For $`Q<Q_{\mathrm{min}}`$ the effect of the potential becomes important and it acts as a wall preventing the quantum fluctuations driving the quintessence field to the origin. Therefore, to a good approximation we have to solve

$$\frac{f(Q,t)}{t}=\frac{H^3(t)}{8\pi ^2}\frac{^2f(Q,t)}{Q^2},$$
(31)
with the boundary condition $`\mathrm{d}f(Q0,t)/\mathrm{d}Q=0`$ (no flux through the origin) and maintaining $`f(Q,t)0`$ and $`_0^{\mathrm{}}f(Q,t)dQ=1`$. This is equivalent to a random walk with a wall. Starting with a (half) Gaussian distribution centred on $`Q=0`$ and with variance $`\sigma ^2(t)`$, and assuming $`H(t)`$ is constant during $`\mathrm{\Delta }t`$, we have

$$\sigma ^2(t_{\mathrm{ini}}+\mathrm{\Delta }t)=\sigma _{\mathrm{ini}}^2+\frac{H^3}{4\pi ^2}\mathrm{\Delta }t.$$
(32)
More precisely, the equation to solve is

$$\frac{\mathrm{d}\sigma ^2}{\mathrm{d}t}=\frac{H^3(t)}{4\pi ^2},$$
(33)
which has solution

$$\sigma ^2(\stackrel{~}{\varphi })\sigma _{\mathrm{ini}}^2+\frac{16W_0}{3\alpha (\alpha +2)m_{\mathrm{Pl}}^2}\left(\stackrel{~}{\varphi }_{\mathrm{ini}}^{\alpha +2}\stackrel{~}{\varphi }^{\alpha +2}\right).$$
(34)
After enough time the evolution is roughly independent of the initial condition $`\sigma _{\mathrm{ini}}`$ and we can set it to zero.

We can now ask what one expects the mean value of the quintessence field to be in our region of the Universe. Because we are presently interested only in the mean value, we should only consider perturbations on scales larger than our present horizon, which are generated between the initial inflaton value $`\varphi _{\mathrm{ini}}`$ and $`\varphi _{\mathrm{obs}}`$. Because of the presence of the wall, the net effect of the quantum fluctuations is to diffuse the quintessence field to larger values. The extent of this diffusion depends on how much inflation occurs before the last fifty $`e`$-foldings; if more inflation occurs then there are more ‘steps’ in the diffusion and also the early steps are larger. We have also computed the distribution obtained by averaging over all possible initial inflaton values between $`\varphi _{\mathrm{obs}}`$ and $`\varphi _{\mathrm{max}}`$ assuming a flat probability.

Actually we are more interested in the probability distribution $`f_q`$ of $`q=\mathrm{log}_{10}(Q/m_{\mathrm{Pl}})`$. We have $`f_q(q,t)=f(Q(q),t)\mathrm{d}Q/\mathrm{d}q`$ which has a peak at $`q_{\mathrm{peak}}=\mathrm{log}_{10}(\sigma (t)/m_{\mathrm{Pl}})`$. In Fig. 2 we show some distributions $`f_q(q,t_{\mathrm{obs}})`$ in the particular case of $`\alpha =2`$; we show the probability distributions for two possible initial conditions for the inflaton, and then the probability distribution averaged over a uniform initial distribution for $`\varphi _{\mathrm{ini}}`$ between $`\varphi _{\mathrm{obs}}`$ and $`\varphi _{\mathrm{max}}`$. Fig. 3 shows the position of the peak in the distributions as a function of $`\alpha `$; as we can see, the result is roughly independent of $`\alpha `$, and as expected the more inflation there is the further the distribution diffuses to large $`Q`$.

We wish to know the quintessence value after inflation to assess when the solution begins tracking. In the slow-roll approximation, the solution for the tracker is

$$\stackrel{~}{Q}(a)\left(\frac{\beta (\beta +2)V_0}{24\pi \rho _{\mathrm{f0}}(1+w_\mathrm{f})}\right)^{1/(\beta +2)}\left(\frac{a}{a_0}\right)^{3(1+w_\mathrm{f})/(\beta +2)},$$
(35)
where $`\rho _{\mathrm{f0}}`$ is the present value of the energy density and $`w_\mathrm{f}`$ the pressure–density ratio of the dominant fluid (either radiation or non-relativistic matter). Some examples of trackers are plotted in Fig. 4. After inflation the quintessence field remains constant as long as its energy density is lower than the tracker’s. We can therefore find the probability distribution $`f_z`$ for $`\mathrm{log}_{10}(z_{\mathrm{tr}}+1)`$, where $`z_{\mathrm{tr}}`$ is the redshift at which the quintessence reaches the tracking behaviour, by comparing the value of $`Q`$ at the end of inflation with the tracking solution. Using Fig. 3 and Fig. 4 we can estimate at which redshift the tracker is reached for some values of $`\alpha `$ and $`\beta `$. We notice that the smaller $`\beta `$ is, the later the tracker is reached for a given initial value of $`Q`$. In Fig. 5 we show some examples of distributions for the redshift at which tracking begins. The discontinuity in the distribution comes from the fact that at about radiation–matter equality there is a transition between two trackers.

The striking feature of this figure is how low the expected redshift of tracking is. While some previous papers have advocated equipartition of the quintessence energy density as an initial condition equip , leading to prompt tracking, we find that tracking is postponed until the late stages, and in particular well after nucleosynthesis. Indeed, the bulk of the probability is not only after nucleosynthesis but after decoupling too; however there is only a small probability that tracking has not begun by the present, which would not lead to acceptable quintessence.

Finally we must say that our results are valid for small initial values of the quintessence field before inflation (in fact formally zero). If this were not the case, the probability distribution for the mean value of the quintessence field at the end of inflation will broaden to larger values and the tracker will be reached even later or not at all.

## IV Conclusions

We have analyzed the quantum dynamics of the quintessence field during a period of early Universe inflation. Due to quantum fluctuations, even if the initial value of $`Q`$ in a certain region of the Universe is small, it is rapidly diffused to large field values and hence low energy density. We have found that typically tracking behaviour begins only at quite a late stage of evolution, well after nucleosynthesis and quite likely after decoupling too.

Although we have discussed specific models of both inflation and quintessence, we expect our results to be quite general; as far as inflation is concerned we need only the assumption that there are significantly more than fifty $`e`$-foldings in total and a standard value of $`H`$ during inflation, while for other models of quintessence where the potential diverges at the origin the result should also be the same as the precise form of the potential is dynamically irrelevant. In particular, while we have only considered models where the tracking density during nucleosynthesis is negligible, these considerations may reinstate models whose tracking density would be unacceptable during nucleosynthesis.

###### Acknowledgements.

M.M. was supported by the Fondation Barbour, the Fondation Wilsdorf and the Janggen-Pöhn-Stiftung, and A.R.L. in part by the Leverhulme Trust. We thank Pier-Stefano Corasaniti, Ruth Durrer, Anne Green, Lev Kofman, Dmitri Pogosyan and Emmanuel Zabey for useful discussions.