Characterisation of chaotic dynamics in pulsating stars

Vistas in Astronomy, Vol. 33, pp. 399-411, 1990 Printed in Great Britain. All rights reserved.

CHARACTERISATION

0083--6656/90 $0.00 + .50 © 199(} Pergamon Press plc.

OF

IN

CHAOTIC

PULSATING

M. Auvergne, DASGAL,

DYNAMICS

STARS

M.J. Goupil, A. Baglin & T. Serre

Observatoire F-92195

Paris,

de

Meudon

URA

Cedex,

CNRS

335

France

1. I n t r o d u c t i o n Most of the variable stars have been thought more or less periodic, but with the advent of increasingly accurate observations and extended time series, many stars have been found to be multiperiodic or irregular. For a review on irregular stellar variability, see Perdang (1985). Two examples of irregular pulsating star light curves, a cepheid HR 7308 and a DA white dwarf G191-16 are shown in Fig. 1 and 2. Observational and theoretical arguments indicate that nonlinear effects are at work in pulsating stars and in some cases dominate the dynamics. Chaos has been suggested to be responsible for some of the observed irregularities and its existence has been shown possible in stellar conditions (Dappen and Perdang 1984, Kovacs and Buchler 1989, Aikawa 1987, Goupil et al. 1988). Therefore the question is raised: how to recognize a chaotic dynamics from observational data of variable stars, i.e. from noisy and time-limited series of luminosity or velocity measurement s.

Jo

-~o -20

-~

. :..-.-"

{!

~8 3~0

'"

. ,:~

.

"':-

~0

"o

""

~0

'

H.~-2.~ .00

""" ~

~

i

.Ik-SEP

l A06 ~ 3

' '

Fig. 1. Velocity curve of HR 7308 from Burki et al. 1982

In the theory of dissipative dynamical systems, the chaotic ones are characterized by two major features: the exponential divergence of trajectories in phase space and/or the fractal structure of the attractor. They are associated with two quantities, the Lyapounov exponents which measure the divergence and the dimension of the attractor. A third characteristic is the route to chaos, the way by which the system evolves from a stable state to a chaotic one. Up to now, three routes leading to chaos have been identified, which are: the Pommeau-Manneville scenario through intermittency, the Ruelle-Takens one through quasiperiodicity and phase locking, and the period doubling cascade of Feigenbaum. The last scenario is the easiest to identify through the existence of subharmonics in the Fourier spectrum. 399

400

M. Auvergne et al.

In the first section of this paper we introduce the fractal structure of a chaotic attractor by means of a simple nonlinear model of a pulsating star. In the second part we present some of the dimension and Lyapounov exponent computation techniques and we discuss the effect of noise. In the fourth section some results on white dwarfs stars are presented and finally a method of ODE reconstruction is presented. 0.4

.

.

.

.

I

"

'

'

'

I

.

.

.

.

I

Am

a

0.2

t i m e (s) -0.2

,

,

,

I

,

,

.

t

5000

0

I

.

.

i

l

104

,

I

1.5 104

Fig. 2. Light curve of G191-16 from Vauclair et al. 1989

2. A s i m p l e n o n l i n e a r m o d e l for s t e l l a r p u l s a t i o n The four basic equations of motion, continuity, energy and transfer governing the motion of a radiative envelope can be written as: OP O---M-

I (GM 4~rr~ ~ Or OM

1

47rv2p

OL OT OM - - C P Ot aT OM

a2r~ + Ot 2 /

10P p Ot

3~L 64v2acr4T 3

where p, P, T, Cp, ~, L, M have their usual meaning. Their linearization has been extensively used to study the instability of stellar models. In the classical one zone formalism introduced by Baker (1966), all gradients are neglected except the luminosity gradient, and the P D E are approximated by a single ODE. We obtain a third order differential equation: da x d2 z ~dz dt ~ + A K ~ o - d t 2 + (3F1 - 4 ) W o - ~ - + B K w ~ z = 0

with z = "6,~ ,

e P ) ,lad . ' F I = (\ ~dlogp

_q~ra, A = (F3 - 1)(s + 4), B = (Fa - 1 ) ( a n - 4 s ) , where a and n are the opacity gradients, assuming a power law dependence of the O 3 0 -----

(i)

Chaotic Dynamics in Pulsating Stars

401

absorption coefficient: ~ = x o p n T - * . K L is the ratio of the luminosity AmC.Two to the internal energy per period. It measures tlie departure from adiabaticity. K = 0 gives the adiabatic approximation. The criterion for instability of this model depends on A, B and K as related to the energy terms, the second and fourth in Eq. (i). Opacity and thermodynamic exponents enter those quantities and the instability mechanism has been named after them '~ and ~' mechanisms. - -

F1

4/3 ~

..... :-

0

1 Y

Fig. 3. Variation of r~ as a function of the ionisation degree

When Eq. (a) is used to model a layer in a partial ionization zone, the physical behavior can be complemented by taking into account the variation of r~ with the ionisation degree zt (Fig. 3). We also have rz < ~ for y = 0.5 (Cox and Giuli, 1968). In a crude way, but preserving the trend of rz, it can be written as (Auvergne and Baglin 1984): 3F1 - 4 = 71(1 + 72)z ~ 3 v~214 -y~ = 3rl(u = 0.5) - 4 and 72 = ~3+,~,/4),. u is the abundance by number of the ionizing element and ,7 = ~ (x is the ionization potential).

We obtain the equation:

d3x

d2z

dr---Y + ~ + ~ ( z + ~

2)

+~=0

This equation is a peculiar case of the equation studied by Moore and Spiegel (1966) and Baker et al. (1971) in a simplified model of radiating convective elements. With reasonable values of the parameters the adiabatic version of Eq. (2) has three fixed points, z = 0 and • = ±=s. Two real eigenvalues are associated to the origin one positive and one negative. The full system (3) exhibits chaotic oscillations. A Poincard section of the trajectories by the plane = = 0 shows the fractal structure of the attractor (Fig. 4). Though very crude, this model favors the idea of the presence of chaotic dynamics in variable stars, directly linked to the nonfinear behavior of the matter of stellar interiors, as confirmed by more sophisticated simulations.

(2)

402

M. A u v e r g n e et al.

!

,

!

'

I

25

20

/

/

./

S

15

oj

10

~,..o

5

i

~0

! ,

I

I

I

I

-200

I

,

0

2O0

Fig. 4. Poincarg map for solution of Eq. (2) 3. C h a r a c t e r i z a t i o n o f c h a o t i c t i m e series We present here several methods developed in the framework of the theory of dynamical systems. When applied to variable stars we have to take into account some peculiarities of astrophysical data: - the dynamics is known through the luminosity (L ~ R2T~41:) or velocity variations. - the data are noisy ( S / N < 50o), the time series short. In addition we have no control on the object. For instance in experimental studies of turbulence, parameters as Reynolds or Rayleigh numbers allow one to control the state of the dynamics. Such a parameter is inaccessible in astrophysical conditions.

2.1 P h a s e s p a c e r e c o n s t r u c t i o n Generally a dynamical system is only known through a single observable s(t), and to calculate, for instance, a Poincaxd map we must be able to construct a set in a n dimensional space equivalent to the original one. A theorem due to Takens (1980) states that a set equivalent to the original one is obtained through the following algorithm known as the time delay method. We construct the set of vectors z(t~) = {s(t), s(t + r), s(t + 2r) ..... s(t + (n - 1)r}, where n is the dimension of the Euclidian space used for the reconstruction. Takens shows that this process is an embedding i.e. that the original set and the embedded set are related through a diffeomorphism. We have the condition n > 2D + 1 where D is the dimension of the original set. This result is independent of the delay T ¢= 0. However if r is too small we have s(t) ~ s(t + r) . . . . . s(t + (n - 1)r) and all points are close to the first bisector. Experience shows that a 'good' value is of the order of 3o% of the characteristic time scale.

403

C h a o t i c D y n a m i c s in P u l s a t i n g S t a r s

2.2 D i m e n s i o n computations Roughly speaking, the dimension of a set is related to the amount of information needed to locate a point on it. Several 'fractal dimensions' have been defined corresponding to the nature of the information one wants to obtain. We present here only three of them, those which lead to algorithmic implementation. For a review on dimensions of chaotic attractors see Farmer et al. (1983). 2.2.1. T h e b o x counting algorithm. The most simple, but tedious, algorithm is the box counting algorithm. Suppose a set in R" defined by points as produced by numerical integrations of ODE or by an observational or experimental apparatus. Dividing R" with hypercubes of size l, it is easy to see that the average number of points inside a box varies as i n where D is the dimension of the set. Suppose a line uniformly covered by points, the number of points in a segment of length I is proportional to I. In the case of a surface it is proportional to l 2. So the number of boxes Nbo=(l) o, l - D and the dimension is given by: D = lim ~-o

log Nbox (l) log 1/I

This box counting algorithm is computationally inefficient because time consuming. On the other hand D does not take into account the density of points on the attractor.

2.2.2 Correlation integral In this second method, the computed quantity is: 1

C(1) = 1v-¢¢ lim g ( Y -

N

1) ~

0 ( 1 - I l z l - zil[)

i,]=l

where 0 is the Heaviside function (Grassberger and Procaccia 1983). c(1) counts the number of couples of points such that l < Iz~ - zjl) (Fig. 5). For small values of l one can show that: c ( o ~ t°

Since the dimension of the attractor is not known, the dimension of the embedding space is also unknown. Therefore c ( 0 is computed for increasing values of , , When the slope of c(l) becomes independent of n one has reached the dimension of the set. C(l) is independent of the particuliar choice of the norm. The max norm or the norm d(z, y) = ~ ]z~ - Y~l give fastest numerical computations.

2.2.3 Effect of noise on C(/) A one percent white noise, mimicking the atmospheric fiuctuations, added to the solution of Eq. (2) swamps the fractal structure (Fig. 6), and we observe that the attractor becomes thick. This corresponds to a diffusive effect in the embedding space.

(~)

M. Auvergne et al.

404

CO

.

.

.

.

I

.

.

.

.

i

.

.

.

.

r

•

""

"

"

I

.

.

.

.

Dtmoneion

.

•

. +

÷

. ÷ w

•

•

~'+

..~

÷

* * * + + + + 4 - +

+ + + 4 - + +

* + +

+

÷

I

,

,

,

,

J

,

~

i

f

I

,

a

,

i

,

]

. I°w2t

Fig. 5. Slope of log(C(0) computed for a solution of Eq. (2) with 1% noise (stars) and without noise (crosses). The dimension of the embedding space is ,, = 5. The dimension of the attractor is 2.08. I

I

I

'~

25 ,s

I

g. 8•

o ".1 r

20

15

10

.'fIe=**

5

-200

0

•

200

Fig. 6. Same attractor as Fig. 4. One percent noise has been added. So we expect that for I smaller than a value related to the variance of the noise the slope of c(l) will be equal to the dimension embedding space. But there are exceptions. Osborne and Provenzale (1989) have shown that, applying the correlation integral method to colored noise (with a spectral power law (p(,~) ~

Chaotic Dynamics in Pulsating Stars

405

w - " ) the resulting dimension is decreasing with increasing a. For a N dimensional noisy curve, with a between 1 and 2.5, a correlation dimension between Dc > 10 to Dc = 1.4 is found. Theoretical values axe between 10 and 1.6. Colored noises of this kind are found, for instance in photometric observations of variable stars, due to the atmospheric scintillation and transparency fluctuations, with an exponent a ~ 1. (A simple and efficient algorithm to model colored noise is the one proposed by Baxnsley et al. in the book Fractal Images.) 2.2.4 H a u s d o r f f d i m e n s i o n

To define the Hausdorff dimension of a set in R", consider a covering with n dimensional cubes of variable edge length e~. We define the quantity HD(~ ) by: HD(,) = i n f ~

D with the condition el < ' i

Let HD = lira HD(e). Hausdorff showed that a non zero limit only exists for D = e~o

DH. Maxtinez et al. (1990 and this meeting) propose an implementation of this definition based on a minimum spanning tree. A set of points being given, we want to connect them in such a way that we have the smallest sum of link lengths. This is the minimun spanning tree (Prim, 1957, Whitney, 1972). If one identifies the length of the largest links with e, the minimun spanning tree is related to HD(~). Increasing v the density of points on the set, ~ decreases in average and the authors expect that the total length of the tree lim l(v) behaves as lira HD(e). With this method the authors have computed the dimension of clusters of galaxies. They found that Dc < D~ and suggest that the concept of multifractals would be more appropriate to the description of the matter distribution in the universe. 2.3 L y a p o u n o v e x p o n e n t s

The Lyapounov exponents describe how a volume of phase space is shrunk by the flow. If for a dissipative system the volume always decreases, some directions expand when others contract (Fig. 7). Two methods have been proposed to calculate Lyapounov exponents on experimental data (Wolf et al, 1985). We present here the method proposed by Eckmann and Ruelle (1982) and Eckmann et al. (1986). 2.3.1 C h a r a c t e r i s t i c e x p o n e n t s : a d e f i n i t i o n a n d a s c h e m a t i c a l g o r i t h m

We suppose that a dynamical system is described by an ordinary differential equation: = F(z) The linearized equation associated with (5) is: d --u(t) = (D~,(t)F)u(t) dt and a solution of (5) is z(t) = I'(z(o)). A solution of the linear equation is u(t) = (D~(o)/t)u(0). We define the matrix T~ (tangent map) by: Of~ or ~d T'x(o)-- (D~F)T~(o) T~ = OzQj

(4)

406

M. A u v e r g n e et al.

p(%.~

Fig. 7. Contraction of a ball of radius r by a flow f . and compute lim (Tt'"rt~ 112t =. A= x-x -zJ t~oo

The eigenvalues of At are equal to exp(A~) where A~ are the exponents.If the phase space is R 3, the rate of growth of a volume element gives the sum of the three exponents A1 + A2 + As, a surface element gives A1 + ~2 and a vector gives the largest exponent A1.

Fig. 8.

407

Chaotic Dynamics in Pulsating Stars

Concerning real data the difficulty is the estimation of the tangent map. Let 16' be the flow; we have ¢6t(z(0 ) = z(t + at). We choose the points z'(t) inside a ball of radius r such that (Fig. 8): d(~(~), ~'(t)) _ ~ ~nd d(~(, + at), ~'(t + at)) _<

where d is a distance. The linear map T~' fit such that:

= D=/6'

(5)

is determined by a least square

T~'(~(t) - ~'(,)) ~ ~(, + at) - ~'(t + at)

Successive matrices T2 6' are computed. In the last step, one successively determines orthogonal matrices Q~ and upper triangular matrices R1 such that QO _ I and T~Q °

=

Q1R1

TjQj-1

=

Q~R'

The Lyapounov exponents are given by ~m=o K-~ log RTk = ,% up to a normalization factor. Remarks: - The length of the original time series must be large enough so that there is a fair number of points in the ball of radius r, satisfying (6). - This method sometimes gives only the positive exponents. This is true for periodic orbits. If the system is perturbed, transients also give negative exponents. In the case of a chaotic attractor some negative exponents can be obtained as for the Rossler attractor. This means also that we have more information on a physical system if it is chaotic than if it is periodic. This is also the case if the system is stochastically perturbed (we will come back to this point later) A Lyapounov dimension can be defined. The exponents being sorted in decreasing order: c(k) DA=k+--

IAk+ll

with c(/) = ~ Ai, and/~ given by c(/¢) > 0 and c(k + 1) < o. It has been conjectured by Kaplan and Yorke that we have in general: Du = D^. The equality holds in many cases but counter examples are known. 2.3.2 E f f e c t o f noise o n L y a p o u n o v e x p o n e n t s It has been reported by Conte and Dubois (1987) that noise reduces the absolute value of Lyapounov exponents. This effect is particularly important for positive exponents because they generally are small. The noise introducing a diffusion of the trajectories in the phase space, the exponential divergence needs larger time interval to be detected and the estimation of the tangent map becomes imprecise. It has been also stressed by Badii et al. (1988) that an extra positive Lyapounov exponent is introduced by filtering data, if the cutoff frequency is of the same order as a negative exponent. Consider a dynamical system defined by Eq. (4) and model a low pass filter by an extra equation:

(6)

M. Auvergne et al.

408

where u(t) is one component of ,. y is a source term in Eq. (7) and all time scales of y shorter than are smoothed out in the new observable quantity z. An extra Lyapounov -,1 is added. Suppose that our system has three exponents {.~1,0, As}. If -,7 < ~ then Dl=2+-and if,7 < IA~I and ,1 < A1, then

D~:2+-this implies:

~ <~

This is also true when the calculated dimension is the information dimension. This result says that we have to be cautious with filtering methods. The variation of the dimension with ,1 gives also an estimation of the Lyapounov exponents.

!

o

7 0.005

Fig. 9. Power spectrum of PG 1351

3. A p p l i c a t i o n to variable w h i t e dwarfs Some white dwarf stars are variable, and their pulsations are triggered by the ,~ mechanism. They are highly stratified objects. A dense degenerate core of heavy elements (carbon, oxygen, nitrogen) is surrounded by layers of pure and partially ionized helium and hydrogen. Despite the fact that white dwarfs are faint objects, we have chosen them because they have the shortest periods. Their values are between two hundred and eight hundred seconds and one can obtain, on one night of observations, about forty periods. They are known also to have complicated light curves. Photometric observations are too noisy and the time series too short (one to three hundred periods are necessary), to obtain firm conclusions when the dimension is computed. The period doubling cascade being often found in physicochemical systems (laser, weak turbulence, chemical reaction, biology) we have searched for subharmonics

Chaotic Dynamics in Pulsating Stars

409

in the power spectra of light curves of white dwarfs. We have found such subharmonics in two large amplitude stars, PG 1351 and G191-16. The power spectrum of PG 1351 is shown in Fig. 9 and the autocorrelation of the power spectrum of G191-16 (Fig. 10) shows the existence of subharmonics at frequencies v/2 and v/4 (Goupil et al 1988 and Vauclair et al. 1989).

0.5 c(,)

0.4 0.3

I

0.2

o.1

0

0

/ 2

~

V0 --1.12 10-a 4

6

Fig. 10. Autocorrelation of the power spectrum of G191-16

Those subharmonics are not always visible, indicating a more complicated dynamics. In the case of G191-16 they are detected when the amplitude is large confirming their link with nonlinear effects. 4. T r a n s i e n t s in t h e v i c i n i t y o f a n a t t r a c t o r a n d i n f o r m a t i o n p r o d u c t i o n We have already said that when dimension computations are applied to astrophysical data, a firm conclusion is difficult to obtain, because of the high noise level and the short length of the time series. In that case, the power spectrum is the most tractable indicator, (see also the wavelet transform, Goupil et al. 1990). But from a theoretical point of view, it is not a good indicator for the specific analysis of chaotic motions. In a quasiperiodic motion, the number of frequencies is equal to the dimension of the attractor. This is no longer true for a chaotic one. A chaotic dynamics produces information because two indistinguishable initial conditions evolve into distinguishable states after a finite time. The Lyapounov exponents measure this effect. In some cases the volume of the phase space which is explored by the system is larger than for a quasiperiodic attractor involving the same number of independent modes. This is why the dimension of a chaotic attractor is also larger. Another way to explore larger and larger volumes of phase space in order to obtain more and more physical information about a system, is to introduce a perturbation. To illustrate this point one considers the two equations:

M. Auvergne et al.

410 - the adiabatic oscillator, +w2z

=

O,

Zo

=

(8)

0.5

- a Van der Pol type oscillator, +

w 2 z _ r n ( A 2 w ~ _ ~2 _ w 2 z ~ ) ~

w i t h rn

=

0.8 a n d A

=

0.5

(9)

They have the same solution for t ~ oo, z(t) = A s i n ( w t ) . If we want to be able to distinguish between the two system, when z(t) is only known, we must obtain a knowledge of the trajectories far from the periodic orbit. A perturbation at time ~ is relaxed by Eq. (9). On the other hand Eq. (8) keeps the memory of the perturbation over an infinite time. Now, suppose that from the signal s(t), we want to determine an equation which has a(t) as a solution. We write:

~,(~) = ~ ai¢(~) J where ~i(t) = ~(t + it) and the ~(z) axe nonlinear functions of the *k as monomials or Legendre polynoms. Then the coefficients a; are determined by a least square fit. The construction of a vector field (Eq. (10)) from one single observable usually is not unique. For instance the knowledge that a periodic orbit is stable is not contained in the signal a(*). To construct a vector field which takes this additional information into account, one must use the additional knowledge provided by transient trajectories in the vicinity of the periodic attractor. For astrophysical applications a stochastically perturbed system would be interesting. This work is in progress. References

Aikawa, T. 1987: Astrophys. Sp. Sc. 139, 281 Auvergne, M. ~: Baglin, A.: 1985 Astron. Astrophys. 142,388 Auvergne, M. ~: Baglin, A.: 1986 Astron. Astrophys. 168, 118 Badii, R., Broggi, G., Derighetti, B., Ravani, M., Ciliberto, S., Politi, A. & Rubio, M.A.: 1988, Phys. Rev. Let. 60, 979 Baker, N.H.: 1966, in Stellar Evolution, eds. R.F. Stein t¢ A.G.W. Cameron, Plenum Press Baker, N.H., Moore, D.W. & Spiegel, E.A.: 1971, Quart. J. Mech. Appl. Math. X X I V , 391 Barnsley, M.F., Devaney, R.L., Mandelbrot, B.B., Peitgen, H.O., Saupe, D. & Voos, R.F.: 1988, The Science of Fraetal Images, eds. tt.O. Peitgen ~z D. Saupe, Springer-Verlag Burki, G., Mayor, M. &: Benz, W.: 1982, Astron. Astrophys. 109,258 Grassberger, P. & Procaccia, I.: 1983, Phys. Rev. A 28, 2591 Goupil, M.J., Auvergne, M. & Baglin, A.: 1988, Astron. Astrophys. 196, L13

(10)

Chaotic Dynamics in Pulsatin9 Stars

411

Conte R. & Dubois M.: 1987, in Proceedings of the IVth Workshop on Nonlinear Evolution Equations and Dynamical Systems ed. J.P. Leon, World Scientific Cox, J.P. & Giuli, R.N.: 1968, in Principle of Stellar Structures, Gordon and Breach Dappen, W. & Perdang, J.: 1984, Mere. Soc. Astron. Italiana. 55,299 Eckmann, J.P. & Ruelle, D.: 1982, Rev. Mod. Phys. 57, 617 Eckmann, J.P., Oliffson Kamphorst, S., Ruelle, D. & Ciliberto, S.: 1986, Phys. Rev. A 34, 4971 Farmer, J.D., Ott, E. & Yorke, J.A.: 1983 Physica D 7, 153 Goupil, M.J., Auvergne, M. & Baglin, A.: 1990, in Proceedings of the Colloquium Confrontation between Evolution and Pulsation, Bologna Kovacs, G. & Buchler, J.R.: 1989, Astrophys. J. 346, 1026 Moore, D.W. & Spiegel, E.A.: 1966, Astron. Astrophys. 116, 341 Martinez, V.J., Jones, B.J.T., Dominguez-Tenreiro, R. &: Van de Weygaert, R.: 1990, preprint Osborne, A.R. & Provenzale, A.: 1989, Physica D 35,357 Perdang, J.: 1985, in Chaos in Astrophysics, eds. J.R. Buchler et al., D. R.eidel Publ. Co. Prim, R.C.: 1957, Bell System Tech. J. 36, 1389 Takens, F.: 1980, in Dynamical System and Turbulence Warwick Lecture Notes Math. 898, eds. D.A. Rand & L.S. Young Vauclair, G., Goupil, M.J., Baglin, A., Auvergne, M. & Chevreton, M.: 1989, Astron. Astrophys. 215, L17 Whitney, V.K.M.: 1972, Comm. ACM 15,273 Wolf, A., Swift, J.B., Swinney, II.L. & Vastano, J.A.: 1985, Physica D 16, 285