A&A 490, 743-752 (2008)
DOI: 10.1051/0004-6361:200809891

Direct numerical simulations of the ${\rm\kappa}$-mechanism

II. Nonlinear saturation and the Hertzsprung progression

T. Gastine - B. Dintrans

Laboratoire d'Astrophysique de Toulouse-Tarbes, Université de Toulouse, CNRS, 14 avenue Edouard Belin, 31400 Toulouse, France

Received 2 April 2008 / Accepted 24 July 2008

Abstract
Context. We study the $\kappa$-mechanism that excites radial oscillations in Cepheid variables.
Aims. We address the mode couplings that manage the nonlinear saturation of the instability in direct numerical simulations (DNS).
Methods. We project the DNS fields onto an acoustic subspace built from the regular and adjoint eigenvectors that are solutions to the linear oscillation equations.
Results. We determine the time evolution of both the amplitude and kinetic energy of each mode that propagates in the DNS. More than 98% of the total kinetic energy is contained in two modes that correspond to the linearly-unstable fundamental mode and the linearly-stable second overtone. Because the eigenperiod ratio is close to 1/2, we discover that the nonlinear saturation is due to a 2:1 resonance between these two modes. An interesting application of this result concerns the reproduction of Hertzsprung's progression observed in Bump Cepheids.

Key words: hydrodynamics - instabilities - waves - stars: oscillations - stars: variables: Cepheids - methods: numerical

1 Introduction

In Gastine & Dintrans (2008) (hereafter Paper I), we modelled the $\kappa$-mechanism in classical Cepheids using a simplified approach, that is, the propagation of acoustic waves in a layer of gas where the opacity bump is represented by a hollow in radiative conductivity. To maintain the mode excitation, the two following conditions apply:

• a sufficient width and amplitude are required for the conductivity hollow;

• this hollow must be located in a precise zone called the ``transition region''.

Our parametric approach enabled us to check the reality of these two conditions and to obtain the instability strips from a linear-stability analysis. By starting from the most favourable situations (i.e. the most unstable linear modes), we performed the corresponding 1D and 2D nonlinear Direct Numerical Simulations (hereafter DNS). These DNS confirmed with a noteworthy agreement both the growth rates and spatial structures of linearly-unstable modes. The nonlinear saturation was reached and a quantitative study of the involved mode couplings remained to be done.

This is the principle aim of this paper. To study the nonlinear interactions between modes, we adapt a powerful method used in, e.g., aeroacoustics (Salwen & Grosch 1981; Wang et al. 2006) or numerical astrophysics (Bogdan et al. 1993, hereafter BCM93). It is based on a projection of DNS fields onto a basis shaped from the regular and adjoint eigenvectors that are solutions to the linear oscillation equations. This method has already been used in a simplified version by one of us to study internal waves in DNS (Dintrans & Brandenburg 2004; Dintrans et al. 2005). In the quoted investigation the eigenvalue problem was adiabatic, and therefore Hermitian, and the authors only accounted for projections onto regular eigenfunctions. This reduction is not applicable in our $\kappa$-mechanism simulations since the transition region requires low densities close to the surface and then high diffusivities (Paper I). Eigenmodes are highly non-adiabatic and imply non-Hermitian oscillation operators with non-orthogonal eigenfunctions. Solving the adjoint problem is therefore mandatory to determine both mode amplitudes and energy.

By using projections onto the two respective sets of eigenvectors, the regular and adjoint ones, the time evolution of each acoustic mode propagating in DNS is completed. The kinetic energy in each acoustic mode is available that highlights the energy transfer between modes. One of the main results of this work is that more than $98\%$ of the total kinetic energy is contained in both the fundamental mode and the second overtone. Because this second overtone is linearly stable, we show that its large amplitude results from a 2:1 resonance with the fundamental mode. In our numerical experiments, this resonance is efficient because the ratio of the two involved periods, that is P₂ / P₀, is close to 1/2.

This nonlinear saturation based on a 2:1 resonance is an interesting result because this mechanism has been proposed to explain the secondary bump in light curves of some classical Cepheids named ``Bump Cepheids''. The bump position is correlated with the oscillation period that leads to the well-known Hertzsprung progression (Payne-Gaposchkin 1954; Hertzsprung 1926). The bump first appears on the descending branch of the light curve of Population I Cepheids with periods of about 6-7 days and then travels up this curve to reach its maximum for 10-11 day Cepheids. For longer periods, it moves down in the ascending branch and disappears for periods longer than 20 days (Bono et al. 2000; Tsvetkov 1990; Whitney 1983). Two main theories have been developed to explain this phenomenon:

• Whitney (1956) suggested an ``echo mechanism'' that Christy (1968) developed in more detail. In this model, two radial-pressure waves (called Christy's waves) moving in opposite directions are generated during compression phases in the HeII ionisation region. The pulse running downward reflects on the stellar core and reaches the surface at the next period during which the bump appears;

• the second explanation, known as the ``resonance mechanism'', has been proposed by Simon & Schmidt (1976, hereafter SS76). The bump results from a resonance between the fundamental mode and the second overtone that is possible for a period ratio P₂ / P₀ close to 1/2.

Three papers of Whitney (1983) and Aikawa & Whitney (1984, 1985) compared the echo and resonance theories. The basic idea was to consider both approaches as two complementary sides of the same physical process. Despite promising results on polytropes (Whitney 1983), the dynamical approach showed that the acoustic-ray formalism did not reproduce the ``running waves'' in Christy's diagram. This implies that most of the energy is contained in standing waves rather than progressive waves (Aikawa & Whitney 1984,1985).

The resonance mechanism has also been investigated using the amplitude-equations formalism by Buchler & Goupil (1984) and Klapp et al. (1985). They established the dominant role of a 2:1 resonance because phases and amplitudes of modes show characteristic variations with the period-ratio P₂/P₀ (see also Buchler et al. 1990; Kovacs & Buchler 1989; Simon & Lee 1981). Despite these results favoring the resonance phenomenon, the echo-resonance controversy remains topical (Bono et al. 2000,2002; Fadeyev & Muthsam 1992).

It is interesting to know whether our model is able to reproduce this Hertzsprung progression. Up to 400 DNS are completed to cover a significant range in the P₂/P₀ ratio, while studying the bump position in luminosity curves. This allows us to accurately reproduce the expected bump-progression with the ratio value: (i) if P₂ /P₀ > 1/2, the bump is located in the descending branch; (ii) if P₂ /P₀ < 1/2, the bump is located in the ascending branch.

In Sects. 2 and 3, we introduce the general-oscillations equations and the associated adjoint problem, respectively. We develop our projection method and provide results in Sect. 4. Section 5 concerns the Hertzsprung progression and we outline our conclusions in Sect. 6.

   
2 The pulsation model and corresponding DNS

We focus on radial modes propagating in Cepheids and restrict our study to the 1D case. Our system, which represents a local zoom around an ionisation region, is composed of a layer of width d filled with a monatomic and perfect gas with $\gamma=c_p/c_v=5/3$ (c_p and c_v are specific heats). Gravity $\vec{g}=-g \vec{e}_z$ and kinematic viscosity $\nu$ are assumed to be constant. Following Paper I, the ionisation region is described by a parametric hollow in radiative conductivity that corresponds to a bump in opacity. We recall below the temperature-dependent profile that is adopted for the radiative conductivity

$$% K_0(T)=K_{{\rm max}}\left[1+{\cal A}\frac{-\pi/2+\arctan(\sigma T^+T^-)}{\pi/2+\arctan(\sigma e^2)}\right],$$ (1)

with

$$% {\cal A}=\frac{K_{{\rm max}}-K_{{\rm min}}}{K_{{\rm max}}},\quad T^{\pm}=T-T_{{\rm bump}}\pm e,$$ (2)

where $T_{{\rm bump}}$ is the position of the hollow in temperature and $\sigma$, e and ${\cal A}$ denote its slope, width and relative amplitude, respectively.

   
2.1 Instability strips from the linear-stability analysis

We are interested in small perturbations about the hydrostatic and radiative equilibria. The layer is fully radiative and the diffusion approximation implies that

$$% \vec{F}'=-K_0New AT'-K' New AT_0$$ (3)

for the radiative flux perturbation (hereafter the ``0'' subscripts mean equilibrium quantities, while primes denote Eulerian ones).

Following Paper I, the depth d of the layer is selected as the length scale, and the top density  $\rho_{{\rm top}}$ and temperature  $T_{{\rm top}}$ as density and temperature scales, respectively. The velocity scale is thus $(c_pT_{{\rm top}})^{1/2}$. Finally, gravity and fluxes are provided in units of $c_p T_{{\rm top}}/d$ and $\rho_{{\rm top}}(c_p T_{{\rm top}})^{3/2}$, respectively. The linearised perturbations obey the following temperature, momentum and continuity dimensionless equations:

$$% \left\lbrace \begin{array}{rcl} \lambda T' &=& \displaysty... ...}- \frac{{\rm d} \ln \rho_0}{{\rm d}z} u, \end{array} \right.$$ (4)

where $R\equiv \rho'/\rho_0$ denotes the density perturbation, u the velocity and $F_{{\rm bot}}$ the imposed bottom flux. Here we seek normal modes using a time-dependence of the form  $\exp(\lambda t)$, with $\lambda = \tau + \hbox{i}\omega$ ($\tau$ is the damping or growth rate of the mode and $\omega$ its frequency). Finally, by assuming rigid walls at both limits of the domain for the velocity, a perfect conductor at the bottom and a perfect insulator at the top for the temperature, we obtain the following boundary conditions

$$% \left\lbrace\begin{array}{l} u=0 ~{\rm for}~ z=(0,1), \\ \\... ...\rm for}~ z=0, \\ \\ T'=0 ~{\rm for}~ z=1. \end{array}\right.$$ (5)

Equations (4), (5) define an eigenvalue problem

$$% A \vec{\Psi}_n = \lambda_n \vec{\Psi}_n,$$ (6)

where $\lambda_n$ is the eigenvalue associated with the (complex) eigenvector $\vec{\Psi}_n=(T',\ u,\ R)^T$. Here n defines the mode order, that is, the number of nodes of the eigenfunction u while the matrix A is a differential operator provided in Eq. (7), where the differential notation $D\equiv {\rm d}/{\rm d}z$ is used for clarity

 

$$% \tiny A~=~ \left( \begin{array}{cccc} {\displaystyle\frac{\gamm... ...\frac{\gamma-1}{\gamma}} T_0 D\\ \\ 0 & -D-D\ln\rho_0 & 0 \end{array}\right).$$ (7)

	Figure 1: Instability strip for the fundamental mode in the plane  $(T_{{\rm bump}},\ {\cal A})$ for $e=0.4,\ \sigma =7$ and $\nu =1$ $\times$ 10-4. The top (at z=1) density and temperature are $\rho =2.5$ $\times$ 10 -3 and T=1, respectively, while gravity is g=7 everywhere. The white cross denotes a particular simulation further studied for $T_{{\rm bump}}=2.1$ and ${\cal A}=70\%$.
Open with DEXTER

Figure 1 displays the instability strip for the fundamental mode. In this parametric survey, we fix the slope $\sigma$ and width e of the conductivity hollow, whereas its amplitude ${\cal A}$ and temperature  $T_{{\rm bump}}$ vary (Eq. (1)). For every couple  $(T_{{\rm bump}},\ {\cal A})$, both equilibrium fields and solutions to the eigenvalue problem (6) are completed using the LSB code (Linear Solver Builder, Valdettaro et al. 2007). Unstable modes with $\tau > 0$ are identified among all eigenvalues in each spectrum and only the fundamental mode is excited by the $\kappa$-mechanism in these simulations.

2.2 The associated nonlinear DNS

To confirm the reality of the instability strip and study the mode saturation, we perform a DNS of the nonlinear problem. We start from a linearly-unstable setup discovered in the parametric survey (the white cross in Fig. 1 of which the oscillation spectrum is provided in Table 1) and advance hydrodynamic equations in time using the Pencil-Code. Because we investigate the 1D case in this work, we cannot apply the 2D Alternate Direction Implicit (ADI) scheme developed in Paper I. We therefore code a 1D version of this scheme that consists of a semi-implicit Crank-Nicholson algorithm where nonlinearities are dealt with using a Jacobian factorisation (see e.g. Press et al. 1992, Sect. 16.6).

Table 1: First linear eigenvalues of the unstable setup used to start the DNS (this setup is denoted by the white cross in the instability strip in Fig. 1). Note that only the fundamental n=0 mode is linearly unstable.

	Figure 2: a) Temporal power spectrum for the momentum in the $(z,\ \omega )$-plane. b) The resulting spectrum after integrating over depth. c) Comparison between normalised momentum profiles for $n=(0,\ 2)$ modes according to the DNS power spectrum (solid black line) and to the linear-stability analysis (dotted blue line).
Open with DEXTER

To determine which modes are present in the DNS once it has reached the nonlinear-saturation regime, we perform in Fig. 2a a temporal Fourier transform of the momentum field  $\rho u(z,t)$ and plot the power spectrum in the $(z,\ \omega )$-plane. With this method, acoustic modes are extracted because they emerge as ``shark fin profiles'' about definite eigenfrequencies (see Dintrans & Brandenburg 2004). In Fig. 2b, we integrate $\widehat{\rho u}(z,\omega)$ over depth to obtain the mean spectrum. Several discrete peaks corresponding to normal modes appear but the fundamental mode close to $\omega_0=5.439$ clearly dominates.

The linear eigenfunctions can be compared to the mean profiles computed from a zoom about given eigenfrequencies in the DNS power spectrum shown in Fig. 2. Figure 2c displays such profiles about $\omega_0=5.439$ and $\omega_2=11.06$(solid black lines), whereas associated eigenfunctions are overplotted in dotted blue lines. The agreement between the linear-stability analysis and the DNS is remarkable.

In summary, Fig. 2 shows that several overtones are present in this DNS, even for long times. Because these overtones are linearly stable (see Table 1), some underlying energy transfers occur between modes. To investigate these intricate couplings in more detail, we need to follow the evolution of each mode separately using the projection method developed in the next section.

   
3 The adjoint problem

As said in the Introduction, the radiative diffusivity in our model is large close to the surface due to the instability criterion (the so-called ``transition region'' criterion). Eigenvectors strongly differ from adiabatic ones, and therefore the quasi-adiabatic approximation fails in this case. The dissipative effects must be fully taken into account and the oscillations operator is non-Hermitian. In other words, the matrix A provided in Eq. (7) is not self-adjoint and its eigenvectors  $\vec{\Psi}_j$ are not mutually orthogonal. This implies that they cannot be used to arrange an acoustic subspace on which the physical fields of the DNS are projected. We thus consider the adjoint problem that has the same spectrum, while its eigenvectors are orthogonal to A-ones because of the biorthogonality property (see Appendix A.1).

3.1 The adjoint matrix A${^\dag }$

We start from the following inner-product

$$% \left\langle \vec{X}, \vec{Y} \right\rangle=\int_0^1 \vec{X}^\dag\vec{Y} {\rm d}z,$$ (8)

where $\dag$ denotes the Hermitian conjugate. The adjoint of the operator A is defined by (e.g. Löwdin 1983)

$$% \left\langle \vec{X}, A\vec{Y} \right\rangle =\left\langle A^\dag\vec{X}, \vec{Y} \right\rangle.$$ (9)

The eigenvectors $\vec{\Phi}_i$ of $A^\dag$ are normalised to verify (see Appendix A.1)

$$% \left\langle\vec{\Phi}_i, \vec{\Psi}_j \right\rangle = \delta_{ij},$$ (10)

where $\vec{\Phi}_i$ are the eigenvectors of $A^\dag$ such that

$$% \left\lbrace \begin{array}{rcl} A \vec{\Psi}_j & = & \lam... ...c{\Phi}_i & = & \lambda_i^* \vec{\Phi}_i. \end{array} \right.$$ (11)

The adjoint-matrix calculation is detailed in Appendix A.2 and the resulting $A^\dag$ is provided in Eq. (16). This corresponds to the following oscillation equations for the adjoint problem:

$$% \left\lbrace \begin{array}{rcl} \lambda T' & = &\displays... ...m d}z}-\frac{F_{{\rm bot}}}{K_0}u\right). \end{array} \right.$$ (12)

3.2 The adjoint boundary conditions

These boundary conditions are detailed in Appendix A.3. While integrating Eq. (A.6) by parts, we obtain surface contributions that must vanish to fulfill the adjoint definition (9) (Bohlius et al. 2007). This implies the following set of adjoint boundary conditions

$$% \left\lbrace\begin{array}{l} u=0 ~{\rm for}~ z=(0,\ 1), \\ ... ...\rm for}~ z=0, \\ \\ T'=0 ~{\rm for}~ z=1. \end{array}\right.$$ (13)

3.3 Results

We solve the adjoint eigenvalue problem defined by Eqs. (12), (13) again using the LSB code, that is, we compute both the whole spectrum of eigenvalues  $\lambda^*_n$and their associated eigenvectors  $\vec{\Phi}_n$. An example is given in Fig. 3 where we plot the real part of both the regular (solid black line) and adjoint (dashed blue line) eigenfunctions of the fundamental n=0 mode.

	Figure 3: Real part of normalised eigenfunctions  $(u,\ R,\ T')$ for the fundamental mode (solid black line) and its adjoint (dotted blue line), for the equilibrium setup denoted by a white cross in Fig. 1.
Open with DEXTER

   
4 The projection method

Once the two sets of (regular) $\vec{\Psi}_n$ and (adjoint)  $\vec{\Phi}_n$ eigenvector are known, the DNS fields are projected at each timestep using the following procedure:

1.

We first arrange the physical fields as

$$% \vec{\Xi}_{{{\rm DNS}}}(z,t) = \left( \begin{array}{l} T'... ...\phantom{'} \\ R\phantom{'} \end{array}\right)_{{\rm DNS}},$$ (14)

where u(z,t) is the velocity and $T'(z,t),\ R(z,t)$ the temperature and density perturbations provided by

$$% T' (z,t)= T(z,t)-T_0(z),\quad R(z,t)={\rho'\over \rho_0}={\rho(z,t)\over \rho_0(z)}-1.$$ (15)

$$% A^\dag = \left( \begin{array}{cccc} \displaystyle \frac{... ...0 D -\frac{F_{{\rm bot}}}{K_0}\right) & 0 \end{array} \right)$$ (16)

2.

We decompose these fields onto the regular-linear eigenvectors

$$% \vec{\Xi}_{{{\rm DNS}}}(z,t) = \Re\left\lbrace\sum_{n=0}^\infty c_n(t) \vec{\Psi}_n (z) \right\rbrace,$$ (17)

where the c_n(t) coefficients are completed using adjoint eigenvectors

$$% c_n(t)=\left\langle \vec{\Phi}_n,\vec{\Xi}_{{{\rm DNS}}} \r... ... \int^1_0 \vec{\Phi}_n^\dag\ \vec{\Xi}_{{{\rm DNS}}} {\rm d}z.$$ (18)

This complex function c_n(t) defines the mode amplitude because its time evolution is related to the eigenfrequency ${\omega}_n$ by

$$% c_n(t)=\vert c_n(t)\vert\ {\rm e}^{\hbox{i} \phi_n(t)} \hbox{~with~} \phi_n(t)\propto \omega_n t.$$ (19)

3.

We perform Fourier transforms or phase diagrams of these c_n(t) coefficients to investigate the evolution of modes.

4.1 Time evolution of the complex-mode amplitudes c_n(t)

	Figure 4: Left panel: time evolution of the real part of the complex-mode amplitudes cn(t) for $n\in \llbracket 0,\ 3\rrbracket$ modes. The theoretical growth or damping curves are overplotted in dashed blue lines. Right panel: their respective Fourier transforms. The modes eigenfrequencies are overplotted in dashed blue lines.
Open with DEXTER

The left panels of Fig. 4 show the time evolution of the real part of the projection coefficients c_n(t) for the first four modes $n\in \llbracket 0,\ 3\rrbracket$. We overplot the curves $\propto$ $\exp(\tau_n\ t)$ derived from the linear-stability analysis (see Table 1). As expected, only the fundamental mode amplitude c₀(t) grows at each timestep. For other linearly-stable modes, amplitudes decay except for the n=2 one. The amplitude of this mode shows a transient decrease until time $t\simeq 10$ and then begins to increase. This behaviour is an estimate of the nonlinear coupling that occurs between this overtone and (at least) the fundamental mode.

We next perform Fourier transforms of these mode amplitudes (right panel of Fig. 4). Following Eq. (19), the c_n(t) amplitudes behave as $c_n(t)\propto \exp(\hbox{i} \omega_n t)$ because theoretical eigenfrequencies and peaks appearing in power spectra overlap one another.

4.2 Phase diagrams

	Figure 5: Phase diagrams for $n\in \llbracket 0,\ 3\rrbracket$ modes. The abscissa (ordinate) corresponds to the real (imaginary) part of the complex mode amplitude cn(t). The increasing orbits are plotted in blue, whereas the decreasing ones are in black.
Open with DEXTER

In each panel of Fig. 5, the imaginary part of the complex mode amplitude c_n(t) is plotted as a function of the real part for the four modes shown in Fig. 4. To avoid unreadable plots, we stop the time evolution at t=50. The decreasing orbits are plotted as black lines, whereas the increasing ones are in blue. In these diagrams, the radius of each trajectory is related to the modulus of the complex mode amplitude, which implies that a decreasing (increasing) radius is associated with a damped (excited) mode.

Once again, the fundamental mode appears to be continuously excited. We note that the excitation phenomenon is coherent in the sense that no discontinuity is observed in the mode orbits, that is, the spiral develops continuously. For the second overtone, after the linearly-transient decrease, the radius begins to increase significantly and this is still the signature of a nonlinear coupling with the fundamental mode. For the linearly-stable n=1 and n=3 modes, things are more intricate because a marginal increase is shown during the last orbits (also seen in Fig. 4). Unfortunately, these phase diagrams are not adapted to investigate this long-time dynamics because orbits will overlap one another.

4.3 The mode kinetic energy content

To address the mode couplings responsible for the nonlinear limit-cycle stability observed at late times, we compute the energy content of each mode separately.

In BCM93, the total amount of sound generated by the turbulent convection is weak because they find that only $0.22\%$ of the kinetic energy is stored in acoustic standing waves. Our problem however consists of an initially static radiative zone without convection and the velocity field that develops is only due to acoustic modes, that is,

$$% E^{\rm tot}_{\rm kin} = E_{{\rm waves}}=\sum_{n=0}^\infty E_n,$$ (20)

where E_n is the kinetic energy contained in the n-acoustic mode. One advantage of our projection method lies in its ability to compute this kinetic energy content E_n because Eq. (17) implies that

$$% u_{{{\rm DNS}}} = \Re \left\lbrace\sum_{n=0}^\infty c_n(t) u_n \right\rbrace,$$ (21)

where $u_{{\rm DNS}}$ and u_n are the velocity of the DNS and the regular eigenvector, respectively. By multiplying by $1/2\rho_0 u$, the kinetic energy density  $E_{{\rm kin}}(z)$ is obtained as

$$% E_{{\rm kin}}(z) \displaystyle\frac{1}{2}\rho_0\Re \left\lbrace\sum_{n=0}^\infty c_n(t) u_n \right\rbrace \Re \left\lbrace\sum_{p=0}^\infty c_p(t) u_p \right\rbrace$$ =

$$\displaystyle\sum_{n=0}^\infty \sum_{p=0}^\infty \frac{1}{2}\rho_0 \Re[c_n(t) u_n]\Re[c_p(t) u_p].$$ = (22)

After integrating over the domain, the total kinetic energy reads

$$% E^{\rm tot}_{{\rm kin}}=\sum_{n=0}^\infty \left(\sum_{p=0}^... ...int_0^1\rho_0 \Re[c_n(t) u_n]\Re[c_p(t) u_p] {\rm d}z\right).$$ (23)

It is possible to determine the contribution of each acoustic mode to the whole amount of kinetic energy because Eq. (20) allows us to define E_n as

$$% E_n=\sum_{p=0}^\infty \frac{1}{2}\int_0^1\rho_0 \Re[c_n(t) u_n]\Re[c_p(t) u_p] {\rm d}z.$$ (24)

Figure 6a displays the long-time evolution of the kinetic energy ratio $E_n /E^{\rm tot}_{{\rm kin}}$ for $n \in \llbracket 0, \ 6\rrbracket$. After the transient linear growth of the fundamental mode, a given fraction of energy is progressively transferred to upper overtones and the nonlinear saturation is achieved above $t\simeq 150$. These nonlinear couplings mainly involve the n=0 and n=2 modes because their energy ratios are dominant. However, the other overtones (i.e. $n\in\llbracket 1,3{-}6\rrbracket$) also participate in this nonlinear saturation process but their energy ratios remain below $1\%$. This is compatible with the weak peaks around the corresponding frequencies ${\omega}_n$ shown in the power spectrum in Fig. 2b or with the marginal increase of the radius orbits in Fig. 5. It is at first glance not surprising because a high-order mode is associated with smaller scales and therefore with a high damping rate (Table 1). This implies that the excitation of these overtones from the fundamental mode would require efficient underlying couplings that are not observed in this simulation.

	Figure 6: a) Kinetic energy ratios for $n \in \llbracket 0, \ 6\rrbracket$ in a logarithmic y-scale. b) Zoom in the n=0 (solid black line) and n=2 (dashed green line) mode energy contents.
Open with DEXTER

The fact that the n=2 mode is more excited than the n=1 one nevertheless contradicts this explanation because its damping rate is more than twice the n=1 one. In Fig. 6b, we expand Fig. 6a for the fundamental mode and this second overtone only. The major energy transfer well occurs between these two modes because they contain about $98\%$ of the total kinetic energy of this DNS. The reason for this favored coupling lies in the period ratio existing between these two modes (see Table 1): the fundamental period is $P_0=2\pi/\omega_0\simeq 1.155$, while the n=2 one is $P_2\simeq 0.568$ such that the corresponding period ratio is close to one half ( $P_2/P_0\simeq 0.491$). The n=2 mode therefore receives in a preferential way some energy from the n=0 mode every two periods and that corresponds, from a dynamical point of view, to a classical 2:1 resonance. Such a resonance is usual in celestial mechanics with, e.g., Jupiter's moons Io (P=1.769 d), Europa (P=3.551 d) and Ganymede (P=7.154 d) and it is well known that it helps to stabilise orbits. In our case, this stabilisation takes the form of a nonlinear saturation: the linear growth of the fundamental mode is balanced by the pumping of energy from the linearly-stable second overtone behaving as an energy sink. The final amplitudes give about $87\%$ of the kinetic energy in the fundamental mode and $11\%$ in the second overtone, the remaining 2 percent being in higher overtones as displayed in Fig. 6a.

   
5 The Hertzsprung progression

The nonlinear saturation in our excitation model therefore results from a 2:1 resonance between the fundamental mode and the second overtone. As shown in the Introduction, this mechanism has also been proposed to explain the secondary bump observed in the luminosity variations of Bump Cepheids (SS76). An interesting correlation between the value of the ratio P₂/P₀ and the bump position was also emphasised leading to the so-called ``Hertzsprung progression'' (Buchler et al. 1990; Kovacs & Buchler 1989; Simon & Lee 1981). Nevertheless this correlation with P₂/P₀ is not observed as the period P₂ deduced from luminosity curves is necessarily equal to P₀/2. Indeed, the second overtone period changes due to nonlinear couplings and finally reaches the value P₂= P₀/2 once the nonlinear saturation is achieved (i.e. the 2:1 resonance). As a consequence, only linear eigenperiods are relevant when studying the influence of the P₂/P₀ ratio on the bump position. Furthermore, Buchler et al. (1990) have shown that this correlation between the phase and the P₂/P₀ ratio is not reproduced in their models when plotting the phase as a function of P₀only. As our hydrodynamical model deals with dimensionless quantities, studying the phase variations according to P₂/P₀ is also more consistent.

	Figure 7: P2/P0 level lines overplotted on the instability strip. Dashed white lines denote ratio P2/P0 <1/2, while solid white lines correspond to P2/P0 >1/2. The critical level P2/P0=1/2 is emphasised as a red line, and the dashed box delimits the region where a large number of DNS are performed to study Hertzsprung's progression (see Figs. 8, 9).
Open with DEXTER

To determine whether our model is able to reproduce this Hertzsprung progression, we first overplot in Fig. 7 the P₂/P₀ level lines with the instability strip discovered in the linear-stability analysis (Sect. 2.1). The solid white level lines correspond to P₂/P₀ > 1/2, the dashed lines to P₂/P₀ < 1/2 and the red line is for the critical level P₂/P₀=1/2. The instability strip is clearly split in two separate parts and the most unstable modes belong to the P₂/P₀<1/2 region. If we now perform the DNS corresponding to those different setups, we can check if the luminosity bump position is really correlated with the period ratio value. Because we are dealing with the 1D case, the luminosity reduces to the radiative flux at the top of the domain, that is,

$$% L(t)\equiv F_{{\rm top}}(t),$$ (25)

and the luminosity curves are transformed into emerging radiative flux curves. Because the bottom flux and the top temperature are imposed in the DNS, these flux variations  $F_{{\rm top}}(t)$ happen about  $F_{{\rm bot}}$.

We then run a large number of DNS (up to 400) in the instability strip to cover a significant range in the P₂/P₀ ratio and cross the 1/2-line; all these DNS belong to the dashed box shown in Fig. 7. In each case, the simulation is performed above the nonlinear saturation ( $t_{{\rm end}}=500$). Because we know that the fundamental mode and the second overtone enter in the time evolution of the emerging flux  $F_{{\rm top}}(t)$, we fit the obtained flux variations using the following decomposition:

$$% f(t)=F_{{\rm bot}}+ A_0\cos(\omega_0 t + \phi_0) + A_2\cos(2\omega_0 t + \phi_2),$$ (26)

where the free parameters $[A_0,\ A_2,\ \phi_0,\ \phi_2]$ are computed using the Levenberg-Marquardt algorithm that is a nonlinear least-square method (e.g. Press et al. 1992, Sect. 15.5). The phase shift between the two signals is defined by

$$% \Delta \phi=\phi_2 -2 \phi_0,$$ (27)

and allows first to phase the flux curves and second, to know where the secondary bump is. If $\Delta \phi$ is less (greater) than $\pi$, this secondary bump will be located in the descending (ascending) branch after (before) the flux maximum.

The whole DNS are summarised in Fig. 8 where isocontours of $\Delta \phi$ are plotted in the same $(T_{{\rm bump}},\ {\cal A})$ instability strip plane used in Fig. 7. The red regions are associated with $\Delta \phi > \pi$, the blue ones with $\Delta \phi < \pi$ and the white ones with $\Delta \phi =\pi$. By overplotting the level lines of the period ratio P₂/P₀, the following correspondence appears:

 

$$\Delta \phi < \pi \hbox{~(bump \textit{after} the max.)}\Leftrightarrow P_2/P_0 > 1/2,$$

$$\Delta \phi > \pi \hbox{~(bump \textit{before} the max.)}\Leftrightarrow P_2/P_0 < 1/2,$$ (28)

which corresponds to the observed correlation (SS76).

	Figure 8: Isocontours of $\Delta \phi$ in the $(T_{{\rm bump}},\ {\cal A})$ plane (note that we rather plot $\Delta \phi - \pi$ to highlight the critical line $\Delta \phi =\pi$). This figure corresponds to the dashed box in Fig. 7 and the black circles display eleven selected simulations showed in Fig. 9.
Open with DEXTER

	Figure 9: Rephased time evolutions of the surface radiative flux about time t=500 for the eleven selected simulations plotted as black circles in Fig. 8. The ratio of eigenperiods P2/P0 is indicated in each case.
Open with DEXTER

To highlight this result, examples of light curves shaping an Hertzsprung progression are provided in Fig. 9. These eleven curves correspond to the black circles distributed across the P₂/P₀=1/2 level line in Fig. 8. The correspondences (28) appear clearly: the bump is located in the descending branch of the light curve if P₂/P₀ > 1/2, while values P₂/P₀ <1/2 lead to a bump in the ascending branch.

Observations also show that the amplitude of the surface velocity presents a double-peaked behavior along the Hertzsprung progression as its center corresponds to a minimum value for the velocity (Bono et al. 2000). In our simulations, no clear correlation between the amplitude ratio A₂/A₀ and the bump position appears and we guess that it may be related to the two following restrictions:

• rigid walls for the velocity have been assumed. As we are interested in the physics involved near the top of the box, the ratio A₂/A₀ is not really relevant there compared to, say, the phase difference between the two modes;

• the sound speed contrast across the box is relatively weak. As the frequency step is given by $\Delta\nu \sim \left(2 \int_L {\rm d}r/c_{\rm s}\right)^{-1}$, there are not enough variations to obtain significant changes in this ratio.

However our model is able to qualitatively reproduce the Hertzsprung progression, as the position of the secondary bump with regard to the main bump evolves in the expected way with the P₂/P₀ ratio. We therefore conclude that the 2:1 resonance supplying the nonlinear saturation in our model results in an ``observational'' bump in the light curve of which the position is linked to the period ratio P₂/P₀.

   
6 Conclusion

In this paper, we investigate the nonlinear saturation of the acoustic modes excited by the $\kappa$-mechanism in direct numerical simulations (DNS). This study follows on from our former work where we began our modelling of the $\kappa$-mechanism in Cepheids by looking the propagation of acoustic modes in a layer of gas in which the ionisation is modelled by a hollow in the radiative conductivity profile (Gastine & Dintrans 2008, Paper I).

To catch the evolution of linearly-unstable modes, we apply an original method consisting of the projection of the DNS fields onto the eigenmodes that are solutions to the linear oscillation equations (Dintrans & Brandenburg 2004). Because the oscillation operator is non-Hermitian, the regular eigenvectors are not mutually orthogonal and we consider the adjoint problem that has the same spectrum while its eigenvectors are orthogonal to the regular ones. Thanks to these two sets of eigenvectors, we compute the projection coefficients for each DNS snapshot and follow the temporal evolution of acoustic modes separately. Corresponding phase diagrams show that the orbit radius of the linearly-unstable fundamental mode is continuously growing with time. On the contrary, orbits associated with the linearly-stable overtones exhibit decreasing radii, except the n=2 one that is increasing at late times.

We derive the kinetic energy contributions of each mode propagating in the DNS. More than $98\%$ of the total kinetic energy is contained in two modes corresponding to the fundamental mode and the second overtone. This last mode, which represents about $10\%$ of the total kinetic energy, is involved in the nonlinear saturation of the instability through a 2:1 resonance with the fundamental mode. It means that the unstable fundamental mode gives energy to the stable second overtone that leads to the limit-cycle stability, as displayed in Fig. 6.

This 2:1 resonance is striking because the same mechanism is probably responsible for the observed bump in the luminosity curves of Bump Cepheids, leading to the well-known ``Hertzsprung's progression''. This progression links the bump position to the period ratio P₂/P₀ existing between the fundamental mode and the second overtone (SS76). We perform a large number of DNS covering a significant range in this ratio and extract the resulting luminosity curves. The 2:1 resonance does lead to an Hertzsprung progression as the bump arises in the ascending branch for P₂/P₀ < 1/2, while it is located in the descending branch for P₂ /P₀ > 1/2.

The study of the $\kappa$-mechanism in a purely-radiative layer is now completed with this work. The physical criteria for the instability as well as the nature of its nonlinear saturation are addressed. Future works will be devoted to the influence of convection on the mode stability. Indeed, a coupling between the convection and pulsations is suspected to play a major role in the disappearance of unstable modes close to the red edge of Cepheids' instability strip. Time-dependent theories of convection have difficulty in explaining the red-edge location since they rely on many unconstrained parameters (e.g. Wuchterl & Feuchtinger 1998; Yecko et al. 1998). Because DNS fully account for crucial nonlinearities, they are appropriate to properly investigate this dynamical coupling between the acoustic modes and convection. We thus plan to perform 2D DNS of the $\kappa$-mechanism in the presence of convection. The interesting cases will of course correspond to a (local) convective-turnover timescale of the same order as the period of the fundamental mode.

Acknowledgements

We thank the referee (G. Bono) for his fruitful comments. Calculations were carried out on the CalMip machine of the ``Centre Interuniversitaire de Calcul de Toulouse'' (http://www.calmip.cict.fr/) and Grid'5000, which is an initiative from the French Ministry of Research through the ACI GRID incentive action, INRIA, CNRS and RENATER and other contributing partners (see https://www.grid5000.fr).

Appendix A: The adjoint formalism

   
A.1 The biorthogonality relation

We consider a continuous linear operator H with eigenvalues $\lambda_i$ and eigenvectors  $\vec{\Psi}_i$ linked by

$$% H\vec{\Psi}_i = \lambda_i \vec{\Psi}_i.$$ (A.1)

The operator $H^\dag$ is its adjoint with eigenvalues $\mu_i$ and eigenvectors $\vec{\Phi}_i$ as

$$% H^\dag\vec{\Phi}_i = \mu_i \vec{\Phi}_i.$$ (A.2)

Using the selected inner product (8), we have

$$% \langle \vec{\Psi}_j, H^\dag\vec{\Phi}_i \rangle = \mu_i \langle \vec{\Psi}_j, \vec{\Phi}_i \rangle,$$ (A.3)

that implies

$$% \langle H \vec{\Psi}_j, \vec{\Phi}_i \rangle = \mu_i \langle \vec{\Psi}_j, \vec{\Phi}_i \rangle,$$ (A.4)

by using the adjoint definition (9). We thus obtain

$$% (\lambda^*_j-\mu_i)\langle \vec{\Psi}_j, \vec{\Phi}_i \rang... ...quad \langle \vec{\Psi}_j, \vec{\Phi}_i \rangle = \delta_{ij}.$$ (A.5)

This property is the biorthogonality satisfied by both the regular and adjoint eigenfunctions (e.g. Löwdin 1983).

   
A.2 The derivation of the adjoint matrix A${^\dag }$

The relation (9) that links an operator with its adjoint takes the following form in our problem

$$% \displaystyle \int_0^1 \left( \begin{array}{ccc} T'_{2} ... ...array}{c} T'_2 \\ u_2 \\ R_2 \end{array} \right)^* {\rm d}z,$$ (A.6)

where $( T'_2 \ u_2 \ R_2)^T$ is an adjoint eigenvector and $( T'_1 \ u_1 \ R_1)^T$ is a regular one. To determine the adjoint operator $A^\dag$, we calculate each element successively from the regular matrix in Eq. (7). We provide below an example of such a calculation for the first column of the adjoint matrix given in Eq. (16).

$\bullet$ Calculation of $A^{\dag }_{11}$

$A^{\dag }_{11}$ comes from the following inner-product

$$% I \!=\! \langle T'_{2}, A_{11} T'_1 \rangle \! =\! \displa... ..._0}\left[K_0 D^2\!+\!2DK_0 D \!+\!D^2 K_0\right]T'_1 {\rm d}z.$$ (A.7)

This integral is split in 3 parts as

$$\gamma\left[\underbrace{\int_0^1 \frac{T'^*_2 K_0}{\rho_0}\frac{{... ...\frac{{\rm d} K_0}{{\rm d}z}\frac{{\rm d} T'_1}{{\rm d}z}{\rm d}z}_{(2)}\right.$$                          I =

$$\left.+ \underbrace{\int_0^1\frac{T'^*_2}{\rho_0}\frac{{\rm d}^2 K_0}{{\rm d}z^2}T'_1 {\rm d}z}_{(a)}\right].$$ (A.8)

(1) and (2) are integrated by parts, leading to

$$% \begin{array}{rcl} (1) & = & \displaystyle\left[\frac{T'^*... ...{\rm d} K_0}{{\rm d}z}\right)T'_1 {\rm d}z}_{(c)}. \end{array}$$ (A.9)

While letting aside in a first step the integrated terms, the three terms (a)-(c) imply that

$$% \begin{array}{rcl} (a)+(b)+(c) & = &\displaystyle\int_0^1 T... ...\ \\ & = & \langle H^\dag T'_{2}, T'_1 \rangle . \end{array}$$ (A.10)

Applying Eq. (A.10), we finally get

$$% A^\dag _{11} = \frac{\gamma K_0}{\rho_0}\left\lbrace D^2-2D\ln\rho_0 D -D^2\ln\rho_0+(D\ln\rho_0)^2\right\rbrace.$$ (A.11)

$\bullet$ Calculation of $A^{\dag }_{21}$

$A^{\dag }_{21}$ comes from the following inner-product

$$% \langle T'_{2}, A_{12} u_1 \rangle = \int_0^1 T'^*_{2}\left[\frac{F_{{\rm bot}}}{K_0}-(\gamma-1)T_0 D\right]u_1 {\rm d}z.$$ (A.12)

Integrating by parts implies that

$$% \begin{array}{rcl} \langle T'_{2}, A_{12} u_1 \rangle &=& ... ...e{(\gamma-1)\left[ T'^*_2 T_0 u_1\right]_0^1}_{=0}. \end{array}$$ (A.13)

Using Eq. (5), the integrated part cancels and we obtain

$$% A^\dag _{21}=(\gamma-1)T_0 D +(2-\gamma)\frac{F_{{\rm bot}}}{K_0}\cdot$$ (A.14)

$\bullet$ Calculation of $A^{\dag }_{31}$

Because A₁₃=0 in Eq. (7), we have

$$% \langle T'_2, A_{13} R_1 \rangle = 0,$$ (A.15)

and therefore

$$% A^\dag _{31}= 0.$$ (A.16)

   
A.3 The boundary conditions satisfied by A${^\dag }$

To determine the boundary conditions for the adjoint problem, we follow a method used in Bohlius et al. (2007). Equation (A.9) contains the boundary conditions on temperature because the definition of the adjoint operator (9) imposes that every integrated term vanishes as

$$% \begin{array}{rcl} (1)+(2) &= & \displaystyle \left[\frac{... ...rac{T'^*_2 K_0}{\rho_0}\right)T'_1\right ] _0^1 =0. \end{array}$$ (A.17)

Applying Eq. (5), this equation implies the two following boundary conditions on temperature:

• At the top z=1:
  

  $$% \frac{T'^*_2 K_0}{\rho_0}\frac{{\rm d} T'_1}{{\rm d}z}=0 \Longleftrightarrow T'_2=0.$$ (A.18)

  
  

• At the bottom z=0:
  

  $$% \left[2\frac{T'^*_2}{\rho_0}\frac{{\rm d} K_0}{{\rm d}z}-\f... ...rm d}z}\left( \frac{T'^*_2 K_0}{\rho_0} \right)\right]T'_1=0,$$ (A.19)

  
  
  which simplifies to
  

  $$% K_0\frac{{\rm d} T'_2}{{\rm d}z}-\left(\frac{{\rm d} K_0}{{\rm d}z}+K_0\frac{{\rm d} \ln\rho_0}{{\rm d}z}\right)T'_2=0.$$ (A.20)

  
  

References

• Aikawa, T., & Whitney, C. A. 1984, ApJ, 282, 527 [NASA ADS] [CrossRef] (In the text)
• Aikawa, T., & Whitney, C. A. 1985, ApJ, 296, 165 [NASA ADS] [CrossRef]
• Bogdan, T. J., Cattaneo, F., & Malagoli, A. 1993, ApJ, 407, 316 [NASA ADS] [CrossRef] (BCM93) (In the text)
• Bohlius, S., Pleiner, H., & Brand, H. R. 2007, Phys. Fluids, 19, 094103 [NASA ADS] [CrossRef] (In the text)
• Bono, G., Castellani, V., & Marconi, M. 2002, ApJ, 565, L83 [NASA ADS] [CrossRef]
• Bono, G., Marconi, M., & Stellingwerf, R. F. 2000, A&A, 360, 245 [NASA ADS]
• Brandenburg, A., & Dobler, W. 2002, CoPhC, 147, 471 [NASA ADS]
• Buchler, J. R., & Goupil, M.-J. 1984, ApJ, 279, 394 [NASA ADS] [CrossRef] (In the text)
• Buchler, J. R., Moskalik, P., & Kovacs, G. 1990, ApJ, 351, 617 [NASA ADS] [CrossRef]
• Christy, R. F. 1968, QJRAS, 9, 13 [NASA ADS] (In the text)
• Dintrans, B., & Brandenburg, A. 2004, A&A, 421, 775 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Dintrans, B., Brandenburg, A., Nordlund, Å., & Stein, R. F. 2005, A&A, 438, 365 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Fadeyev, Y. A., & Muthsam, H. 1992, A&A, 260, 195 [NASA ADS]
• Gastine, T., & Dintrans, B. 2008, A&A, 484, 29 [NASA ADS] [CrossRef] [EDP Sciences] (Paper I) (In the text)
• Hertzsprung, E. 1926, Bull. Astron. Inst. Netherlands, 3, 115 [NASA ADS]
• Jones, E., Oliphant, T., Peterson, P., et al. 2001-2008, SciPy: open source scientific tools for Python, http://www.scipy.org/
• Klapp, J., Goupil, M. J., & Buchler, J. R. 1985, ApJ, 296, 514 [NASA ADS] [CrossRef] (In the text)
• Kovacs, G., & Buchler, J. R. 1989, ApJ, 346, 898 [NASA ADS] [CrossRef]
• Löwdin, P.-O. 1983, J. Math. Phys., 24, 70 [NASA ADS] [CrossRef] (In the text)
• Payne-Gaposchkin, C. 1954, Annals of Harvard College Observatory, 113, 151 [NASA ADS]
• Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical recipes in FORTRAN, The art of scientific computing (Cambridge: University Press), 2nd edn. (In the text)
• Salwen, H., & Grosch, H. E. 1981, J. Fluid Mech., 104, 445 [NASA ADS] [CrossRef]
• Simon, N. R., & Lee, A. S. 1981, ApJ, 248, 291 [NASA ADS] [CrossRef]
• Simon, N. R., & Schmidt, E. G. 1976, ApJ, 205, 162 [NASA ADS] [CrossRef], SS76 (In the text)
• Tsvetkov, T. G. 1990, Ap&SS, 173, 31 [NASA ADS] [CrossRef]
• Valdettaro, L., Rieutord, M., Braconnier, T., & Fraysse, V. 2007, JCoAM, 205, 382 [NASA ADS] (In the text)
• Wang, M., Freund, J. B., & Lele, S. K. 2006, Ann. Rev. Fluid Mech., 38, 483 [NASA ADS] [CrossRef]
• Whitney, C. A. 1956, AJ, 61, 192 [NASA ADS] [CrossRef] (In the text)
• Whitney, C. A. 1983, ApJ, 274, 830 [NASA ADS] [CrossRef]
• Wuchterl, G., & Feuchtinger, M. U. 1998, A&A, 340, 419 [NASA ADS]
• Yecko, P. A., Kolláth, Z., & Buchler, J. R. 1998, A&A, 336, 553 [NASA ADS]