CHINESE ASTRONOMY AND ASTROPHYSICS
ELSEVIER
Chinese Astronomy and Astrophysics 35 (2011) 274–284
Unstable Non-radial Modes of Upper Main-sequence Stars† 1
2
XIONG Da-run1
DENG Li-cai2
Purple Mountain Observatory, Chinese Academy of Sciences, Nanjing 210008 National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100012
Abstract Using a non-local time-dependent theory of convection, the linear stability is analyzed for low-degree (l = 1 − 4) g39-p4 modes of the evolutionary sequences of stars with masses from 3M to 30M . The influence of convection on the pulsation stability of the non-radial modes is studied. The results show that this influence is not negligible. However, if the details of the single star and the single mode are not concerned, the location and scope of instability strip are not affected significantly by convection. Key words: convection—stars: oscillations—stars: variable
1. INTRODUCTION There are a series of pulsating variables from top to bottom along the main-sequence on the H-R diagram—from the OB-type star (ζ OPh) with line-profile variations, βCep-type stars and slowly pulsating B-type to δ Scuti stars. Until the appearance of OPAL and OP opacities the excitation mechanism of βCep was an unfathomable mystery for astronomers. Cox et al.[1] , Moskalik & Dziembowski[2] and Dziembowski & Pamyatnykh[3] calculated the non-adiabatic oscillations of OB-type stars by using the OPAL opacity. It was ultimately affirmed that the oscillations of βCep are excited by the opacity bump of the element Fe at temperature about 2× 105 K. It was also found that the main-sequence and early post mainsequence stars between βCep and δ Scuti stars are pulsationally unstable in g-modes[3,4] . The details of the investigation of oscillations for the upper main-sequence stars may be found in Pamyatnykh’s review paper and refrences therein[5] . †
Supported by National Natural Science Foundation Received 2010–06–02; revised version 2010–10–12 A translation of Acta Astron. Sin. Vol. 52, No. 2, pp. 115–125, 2011
[email protected]
0275-1062/11/$-see frontfront matter © 2011 B.V. All Science rights reserved. cElsevier 0275-1062/01/$-see matter 2011 Elsevier B. V. All rights reserved. doi:10.1016/j.chinastron.2011.07.005 PII:
XIONG Da-run and DENG Li-cai / Chinese Astronomy and Astrophysics 35 (2011) 274–284
275
The coupling between convection and oscillations was neglected in all the past theoretical calculations of non-adiabatic oscillations for upper main-sequence stars. This is mainly because that until the present we still lack a complete convection theory, which is not only simple enough for understanding but also convenient enough for application. This, surely, is a very difficult subject. Fortunately, the convection zone is very shallow and convection has no significant influence on the general qualitative properties of oscillations for these hot stars. Therefore it is possible to determine the approximate location of the pulsation instability strip on the H-R diagram, even the coupling between convection and oscillations is neglected. Although the convection zone of these hot stars is not well developed, the convective flux, as we shall see later, still occupies a few to 30 percent of total energy flux in their excitation zone of oscillations. Convection not only decreases the excitation strength of radiative κ-mechanism, but also affects the pulsation stability of stars by means of the convective energy transport (thermodynamic coupling between convection and oscillations) and turbulent pressure as well as turbulent viscosity (dynamical coupling between convection and oscillations). The goal of the present paper is to investigate the influence of convection on the pulsation stability for upper main-sequence stars. A brief description for treatment of non-local and time-dependent convection in non-radial oscillations of stars is given in Section 2, the results of numerical calculations are presented in Section 3, the influence of convection on the pulsation stability of stars is discussed in Section 4, the comparison with other theoretical works is given in Section 5 and the last section is a brief summary.
2. TREATMENT OF NON-LOCAL AND TIME-DEPENDENT CONVECTION FOR NON-RADIAL OSCILLATIONS OF STARS The goal of the present paper is to investigate the influence of convection on the pulsation stability of stellar non-radial oscillations. A non-local and time-dependent convection theory available for treatment of non-radial oscillations is the key to this subject. Differing from the radial oscillations of stars, the exchange of energy and momentum occurs not only in the radial direction, but also in the horizontal one. The phenomenological time-dependent mixing-length theory is helpless for this subject. Differing from the phenomenological mixing-length theory, our non-local and time-dependent convection theory is a dynamical theory of the correlation functions deduced from the hydrodynamic equations and turbulence theory. When convection sets in, any a physical quantity X in the hydrodynamical equations (i.e., conservations of mass, momentum and energy, and additionally an equation of radiation transfer in stars), such as the gas temperature, density, pressure, inter¯ and turbulent nal energy, entropy and velocity etc, is expressed as the sum of its average X fluctuation X , X = X + X .
(1)
Writing down all the quantities in the above form, putting them into the hydrodynamic equations, and making Taylor expansions for all fluctuations while keeping only the first order terms and ignoring all second and higher order terms, then averaging the whole set of equations, we can derive the equations of average motion. Taking the conservations of momentum and energy as examples and going through the described manipulation steps,
276
XIONG Da-run and DENG Li-cai / Chinese Astronomy and Astrophysics 35 (2011) 274–284
we have
1 Dui + ∇k (g ik P + ρui uk ) + g ik ∇k φ = 0 , Dt ρ ρ CP
DT DP −B = ρ(εN + ε1 ) − ∇k (Frk + ρ CP uk T ) , Dt Dt
(2) (3)
where ρ, T , P , N and ui denote, respectively, density, temperature, pressure, nuclear energy generation rate per unit mass and the i-th inverse variance component of velocity, φ is gravitational potential, CP is the specific heat at constant pressure, B = − (∂ ln ρ/∂ ln T )P is the gas expansion coefficient, and 1 = 1.56GMr ρ¯x3 /c1 r2 P¯ is the turbulent viscous dissipation rate per unit mass of fluid. When convection is absent, the molecular viscosity can be neglected. When convection turns up, however, turbulence gets energy from buoyant force and stellar oscillatory motion (happening mainly in the low wave number region of turbulent spectrum). Via the cascading process of turbulence, turbulent kinetic energy is transformed from low to high wave number regions in the spectrum, and is eventually dissipated into thermal energy through molecular viscosity (happening mainly in the highest wave number region of the turbulent spectrum). As a result of such a process, the setting in of convection may greatly enhance viscous dissipation. g ik is the metric tensor. The implicit summing rule is adopted here, i.e., a pair of identical sub- and super-scripts means summation going from 1 to 3. When convection happens, it follows from Eqs.(2) and (3) that the equation of average motion (2) has a Reynold’s stress term ρui uk , and the average energy equation (3) has a convective energy flux term F¯C = ρ¯C¯P uk T in addition to radiative energy flux. They are the second-order correlations of turbulent velocity and temperature fluctuations. Subtracting the corresponding average equations from the original hydrodynamic equations, we can get the dynamic equations of turbulent velocity (ui ) and temperature fluctuation (T /T ). Then we can derive the dynamic equations of auto- and cross-correlations, ui uk , 2 Z = T /T and V i = ui T /T . After simplifying the third-order correlations, all these dynamic equations of correlation functions, together with the average hydrodynamic equations, form a complete set of dynamic equations for stellar structure and oscillations. Our theory, differing from the phenomenological mixing-length theory, originates from the hydrodynamic equations and turbulent theory, therefore it has a more solid hydrodynamic base and can describe the dynamic behavior of turbulent convection more exactly than the phenomenological mixing-length theory. It can be seen from Eqs.(2) and (3) that the turbulent Reynold’s stress ρui uk and convective energy flux F¯C = ρ¯C¯P uk T are, respectively, a second-order tensor and three-dimensional vector. Therefore, our theory can be used to deal with the exchange of energy and momentum in both the radial and horizontal directions in the non-radial oscillations of stars. For details of description of our nonlocal and time-dependent convection theory, please refer to our previous works[6,7]. In the present paper, convection is treated by the nonlocal convection theory mentioned above in the calculations of both equilibrium models and non-adiabatic oscillations. In our theoretical scheme of non-local convection theory, the complete set of stellar equilibrium structure is composed of 10 equations, in which four are the traditional equations of stellar structure. However, it can be seen from Eqs.(2) and (3) that the turbulent pressure appears in the equation of static mechanical equilibrium (2) and the convective flux FC appears in the
XIONG Da-run and DENG Li-cai / Chinese Astronomy and Astrophysics 35 (2011) 274–284
277
conservation equation of energy (3). Hence three equations of auto- and cross- correlation functions of turbulent velocity and temperature have to be added to the equations of static structure of stars. All of them are second-order equations, so the total number of equations for stellar static structure is 10. For non-radial and non-adiabatic oscillations of stars, apart from above 10 equations, the Poisson equation of gravitational potential and the dynamical equation of turbulent velocity-temperature correlation in horizontal direction have to be supplemented. All of them are second-order equations, so the complete equations of non-radial and non-adiabatic oscillations of stars are in 14th-order. They become 6thorder when the coupling between convection and oscillations is neglected. Therefore, when nonlocal convection is taken into account, the calculations of both equilibrium models and non-adiabatic oscillations become very much complex.
3. RESULTS OF NUMERICAL CALCULATIONS Adopting the scheme described in the above section, we calculated the linear non-adiabatic oscillations of l = 1−4 and g39-p4 modes for the evolutionary sequences of stars with masses from 3M to 30M and the solar abundance (X=0.70, Z=0.02). The non-radial oscillations of spherical degree l = 1 − 4 for p4-g9 modes, and l = 1 − 2 for g10-g39 modes are calculated. The OPAL opacity8] and a MHD equation of state[9−11] modified by us are applied in the present work. Fig.1 illustrates the fractional radiation flux Lr /L as a function of depth for the equilibrium models of main-sequence and early post main-sequence stars with masses 3∼25 M . We choose only a small region of temperature, lgT from 5.1 to 5.5, because the convective flux is almost negligible in other regions of temperature for these hot stars, whose convection is less developmental. Figs.1(a) and (b), respectively, refer to nonlocal and local models. It is known that the curves of Lr /L for the nonlocal models become shallower and wider in comparison with the local ones. This is due to the effect of nonlocal diffusion. The change of profile of fractional convection flux Lc /L is not too large, while the turbulent velocities penetrate more deeply into the radiation zone than the convection flux does. It can be known from Eqs.(2) and (3) that convection affects the pulsation stability of stars by means of convective energy transport (the thermodynamic coupling) and turbulent pressure and viscosity (the dynamic coupling). The influences of change of the profile of turbulent velocities and convective flux in equilibrium models due to the nonlocal convective diffusion on the pulsation stability cannot be despised. In consideration of this reason, the nonlocal convection theory is used in calculations of both non-adiabatic oscillations and equilibrium models. In the evolution of stars from left to right on H-R diagram, the effective temperature of star decreases and convection intensifies. Correspondingly the curves become deeper and wider from upper to lower in Fig.1. The fractional convective flux Lc /L increases from a few to about 30 percent and this cannot be neglected. 3.1 Low-degree p- and g-modes (g4-p4 modes) We have found only the pulsationally unstable g1-, f- and p1-modes. There are no unstable g4-g2 modes and p-modes higher than first order in our calculations. Figs.2(a)-(d) show the distributions of pulsationally stable (small solid circles) and unstable g1- (open triangles), f- (open circles) and p1-modes (crosses) on H-R diagram respectively for l = 1−4,
278
XIONG Da-run and DENG Li-cai / Chinese Astronomy and Astrophysics 35 (2011) 274–284
that corresponds to the βCep instability strip. The p1-modes are located in the lower left area of the instability strip, g1-modes—in the upper right region of instability strip and f-modes lie between them. For l = 1 the extension of the instability strip is the largest, from M ≈ 7M to ≈ 30M in H-R diagram. With the increase of l it seems that the instability strip becomes short and narrow. It reduces to M from ≈ 8M to ≈ 20M for l = 4. In the meantime g1-modes move toward the high-temperature direction, from complete separation of f-mode to complete covering with f-modes when l increases from 1 to 4.
Fig. 1
Fractional radiation flux Lr /L vs. depth (lgT ) for the evolutionary models of main-sequence and
early post-main-sequence stars with masses M =5, 10, 15, 20 and 25 M . Note that the zero-lines in the upper panels are shifted by 0.5 relative to the lower one. (a) nonlocal convective models, (b) local convective models
3.2 Intermediate (g5-g19) and High-order (g20-g39) g-modes Fig.3 shows the distribution of pulsationally unstable intermediate- (panel (a)) and high-order (panel (b)) g-modes on the H-R diagram. The instability strip is a narrow and long band inclined slightly from upper left to lower right. It intersects the zero-age mainsequence at M ∼ (3 − 5)M for the intermediate-order g-modes or M ∼ (3 − 4)M for the high-order g-modes. Towards high luminosity, the instability strip moves into the evolved and post main-sequence. The width of instability strip has some uncertainty, it depends on the choice of overshooting parameter. The overshooting distance is taken to be 0.25 Hp in the present work. It can be seen from Figs.3(a) and (b) that with the increase of radial order the instability strip of g-modes tends to move towards low temperature. Fig.4 shows the period-luminosity relations of unstable modes for all the main-sequence and early post main-sequence stars. It can be clearly seen from Fig.4 that the unstable modes are divided into two completely separate groups: a group composed of unstable p1-, f- and g1-modes located in the low-right area of Fig.4, and the other group composed of unstable intermediate- and high-order g-modes located in the top of Fig.4. These two groups are separated by the stable g2-g4 modes. The first group is βCep stars, the pulsation periods of which are ∼ 0.1 − 1 days, and the second group is the slowly pulsating OB stars, the
XIONG Da-run and DENG Li-cai / Chinese Astronomy and Astrophysics 35 (2011) 274–284
279
pulsation periods of which are ∼ 0.5 − 20 days. On the H-R diagram, the βCep stars are a group of less evolved main-sequence stars located in the higher temperature region of instability strip than that of the slowly pulsating OB stars (refer to Figs.2 and 3). It can be seen from Fig.3 that the high-luminosity OB stars (lgL/L ≥ 3.5 or M > ∼ 8 M ) which are unstable in g-modes are high-evolved main-sequence or post main-sequence stars, and the < slightly low-luminosity stars (lgL/L < ∼ 3.5 or M ∼ 8 M ) which are unstable in g-modes are the slowly pulsating B stars and are located in main-sequence and have pulsation periods P ≈ 0.5 ∼ 7 days.
Fig. 2 Pulsationally stable (small solid circles) and unstable low-order modes on the H-R diagram. Open triangles, open circles and crosses are, respectively, the unstable g1-, f- and p1-modes. All the g2, g3, g4-modes and all the p-modes higher than first order are pulsationally stable. (a) l = 1, (b) l = 2, (c) l = 3, (d) l = 4
4. INFLUENCE OF CONVECTION ON THE PULSATION STABILITY OF UPPER MAIN-SEQUENCE STARS Fig.5 illustrates the accumulated work of g20-mode as a function of depth lgT for two 5M
280
XIONG Da-run and DENG Li-cai / Chinese Astronomy and Astrophysics 35 (2011) 274–284
Fig. 3 Pulsationally unstable intermediate-order (panel (a)) and higher-order g-modes (panel (b)) on H-R diagram for l = 1 and 2. The coupling between convection and oscillations is taken into account.
stars with different effective temperatures. The solid and dashed lines are, respectively, the accumulated work curves with and without coupling between convection and oscillations, and the fractional flux Lr /L (dotted lines) is also drawn in the figure. Our accumulated work is defined by integration from the stellar center to the surface, therefore the value of the accumulated work at the surface is exactly equal to the amplitude growth rate per period η = −2πωi /ωr , where ωi and ωr are, respectively, the imaginary and real parts of the complex circular frequency ω = iωi + ωr . Our calculations show that W (M0 ) and η agree with each other very well, the difference between them is, in general, within 1%. This is a straightforward and effective justification for the correctness of numerical calculations. It can be seen from Fig.5 that the influence of convection on the pulsation stability is not negligible, although the fractional convective flux Lc /L = 1 − Lr /L is not too large. Convection, in general, changes only the size of growth rate of pulsation amplitude (see Fig.5(a)). However, sometimes its sign is also changed, such as in Fig.5(b). By the way it can be seen from Fig.5 that the excitation region is located in the zone of Fe absorbing bump around T ∼ (2.0 − 2.5)105 K. Fig.6 shows the logarithm of absolute value of amplitude growth rate of g19-mode as a function of the effective temperature for two evolutionary sequences of stars with masses 4M and 8M . It can be seen from Fig.6 that convection changes the size of amplitude growth rate in certain degree, but does not change its sign for most stars. However, the sign of amplitude growth rate may be changed for the stars around the boundary of instability strip, while the coupling between convection and oscillations is or is not taken into account. This means that the pulsation stability changes in essential aspect: variation from pulsationally unstable to stable, or on the contrary, from pulsationally stable to unstable. Fig.7 shows the distribution of pulsationally unstable intermediate- (panel (a)) and
XIONG Da-run and DENG Li-cai / Chinese Astronomy and Astrophysics 35 (2011) 274–284
281
Fig. 4 Pulsationally unstable non-radial modes on the plane of lgP −lgL/L , where P is the pulsation period and L is the luminosity of stars. Open circles and crosses on the upper area of the figures are the unstable intermediate- and high-order g-modes, respectively, for l = 1 and 2, and inverse-triangles, circles and triangles on the low-right region of figures are, respectively, g1, f and p1-modes. Small solid and large open signs are, respectively, for l = 1 and 2. (a) The coupling between convection and oscillations (CCO) is taken into account, (b) CCO is neglected.
Fig. 5 Normalized accumulated work vs. depth for g20-mode of two evolutionary models of 5M star. The solid and dashed lines are, respectively, the cases with and without convective coupling. The dotted lines are the fractional radiation flux Lr /L. (a) lgL/L = 2.8356, lgTe = 4.2073, (b) lgL/L = 2.9926, lgTe = 4.1570
282
XIONG Da-run and DENG Li-cai / Chinese Astronomy and Astrophysics 35 (2011) 274–284
high-order g-modes (panel (b)) on the H-R diagram when the coupling between convection and oscillations is neglected. By comparison with Fig.3, it can be found that there are some differences between Fig.3 (with convective coupling) and Fig.7 (without convective coupling) in certain degree. The differences of intermediate-order g-modes (Fig.3(a) and Fig.7(a)) are larger than that of high-order modes (Fig.3(b) and Fig.7(b)). The differences are larger for low-order modes. However, the total number of the pulsationally unstable g-modes is too large, so if the details of single stars and single oscillation modes are not concerned, the locations and scopes of instability strip are, in general, close to each other in cases with or without convective coupling. This is true especially for high-order g-modes.
Fig. 6 Logarithm of the absolute value of amplitude growth rate vs. the effective temperature of stars for g19-mode of the evolutionary models of 4M (a) and 8M (b) stars. The circles and triangles are, respectively, for the cases with and without convective coupling. The small solid and large open signs are, respectively, the pulsutionally stable (η < 0) and unstable (η > 0) models.
5. COMPARISON WITH OTHER WORKS Our theoretical pulsation instability strip of βCep is roughly consistent with those of Dzimbowski & Pamyatnykh[3] as well as Dziembowski, Moskalik & Pamyatnykh[4]. However, there are obvious differences among them. 5.1 βCep Instability Strip (1) There is an upper boundary in our βCep instability strip, which extends to lgL/L ≈ 5.4 or M ≈ 30M for l = 1. The maximum luminosity of instability strip tends to decrease with increase of l. It decreases to lgL/L ≈ 5.0 or M ≈ 20M for l = 4. Dziembowski, Moskalik & Pamyatnykh’s instability strip[4] is more extended in low temperature than ours does. This may be due to the fact that the coupling between convection and oscillations was neglected in their work. (2) We do not find pulsationlly unstable g2-g4 modes and all the p-modes higher than the first order. In addition, the unstable g1-mode moves toward high temperature and
XIONG Da-run and DENG Li-cai / Chinese Astronomy and Astrophysics 35 (2011) 274–284
283
becomes close to the instability region of f-mode with increase of l (Figs.2(a)-(d)). However, there are pulsationally unstable p2, g2 and even g3 in Dzimbowski & Pamyatnykh’s paper[3] . We do not know the reasons causing these differences, because we do not know the details of their work. After comparison of theoretical results with observations, it is found that our theoretical βCep instability strip agrees with observations very well[12,13] , and it is not inferior to that of Dzimbowski & Pamyatnykh[3]. O-type βCep stars have not been found except HD-34656 in observations of variable stars. This seems to be favorable for our theoretical results.
Fig. 7 The pulsationally unstable intermediate-order (panel (a)) and higher-order g-modes (panel (b)) on the H-R diagram for l = 1 and 2, as in Figs.3(a) & (b), but the coupling between convection and oscillations is neglected.
The reason of existence of the βCep instability strip is the same as that of Cepheid instability strip. When the effective temperature of stars is too high, the region of Fe absorbing bump is too shallow and its mass is too small, then its excitation is not intense enough to compensate the radiation damping in the interior of star. Unless the effective temperature of star is low enough, the region of Fe absorbing bump is buried deeply enough, and the excitation becomes powerful enough to compensate the radiation damping in interior, then the star can be excited by itself. When the stellar temperature decreases further, the excitation due to Fe absorbing bump increases, and this is accompanied with the damping of the surface dissipation layer overlaying the excitation region of Fe absorbing bump increases quickly. The star becomes pulsationally stable again, while the damping of surface absorbing layer exceeds the excitation resulted from Fe absorbing bump. This is the reason why there is a red boundary in βCep instability strip. It is different from the reason of the existence of the red boundary of Cepheid instability strip. The latter is due to convection. It is also the reason why people cannot understand the existence of the red edge of βCep instability strip. 5.2 Intermediate- and High-order g-modes Dziembowski, Moskalik & Pamyatnykh[4] did not give the theoretical results in detail
284
XIONG Da-run and DENG Li-cai / Chinese Astronomy and Astrophysics 35 (2011) 274–284
and drew only a rough scope of instability strip for intermediate- and high-order g-modes. The red boundary of instability strip is defined by the limited width of main sequence. The accurate quantitative comparison is impossible. Comparing our Figs.3(a) & (b) with Fig.3 in Pamyatnykh’s paper[3] , it can be seen that they are consistent roughly. The instability strip of the slowly pulsating B-type stars covers a region from lgL/L ≈1.5 (M ≈ 3 M ) to lgL/L ≈4.0 (M ≈ 9 M ). When the luminosity is further higher, the instability strip will move to the outside of the main sequence and enter into the post main sequence. Please note that the convective overshooting has been considered in our evolutionary model of stars, so our main sequence is wider than theirs. 6. CONCLUSIONS Comparing the results of non-adiabatic oscillations with and without convective coupling, it is found that the effect of convection on pulsation stability is not negligible, although the surface convection zone is not well developed. The coupling between convection and oscillations, in general, changes only the amplitude of growth rate in most cases. Sometimes the sign of amplitude of growth rate is also changed. If the details of the single stars and the single oscillation modes are not concerned, the location and scope of the instability strip on H-R diagram do not change for both cases with and without the convective coupling for the upper main-sequence stars. We do not find pulsationally unstable g2-g4 modes and all the p-modes higher than the first order for upper main-sequence stars with the solar abundance. On the plane of lgL/L−lgP the pulsationally unstable modes can be evidently divided into two groups. One is the βCep stars unstable in p1-, f- and g1-modes, and another is the slowly pulsating OBtype stars pulsating in intermediate- and high-order g-modes. The low-luminosity members of the latter group (lgL/L <∼ 3.5 or M <∼ 8M ) are located in the main sequence and they are slowly pulsating B-type stars, while its high-luminosity members have moved to post main sequence. These results are, in general, consistent with those of Dziembowski, Moskalik &Pamyatnykh[4]. References 1
Cox A. N., Morgan S. M., Rogers F. J., Iglesias C. A., 1992, ApJ, 393, 272
2
Moskalik P., Dziembowski W. A., 1992, A&A, 256, L5
3
Dziembowski W. A., Pamyatnykh A. A., 1993, MNRAS, 262, 204
4
Dziembowski W. A., Moskalik P., Pamyatnykh A. A., 1993, MNRAS, 265, 588
5
Pamyatnykh A. A., 1999, AcA, 49, 119
6
Xiong D. R., 1989, A&A, 209, 126
7
Xiong D. R., Cheng Q. L., Deng L., 1997, ApJS, 108, 529
8
Rogers F., Iglesias C. A., 1992, ApJS, 79, 507
9
Hummer D. G., Mihalas D., 1988, ApJ, 331, 794
10
Mihalas D., D¨ appen, W., Hummer D. G., 1988, ApJ, 331, 815
11
D¨ appen W., Mihalas D., Hummer D. G., Mihalas B. W., 1988, ApJ, 332, 261
12
Deng L., Xiong D. R., 2001, MNRAS, 327, 881
13
Tian B., Me H., Deng L., Xiong D. R., Cao H. L., 2003, CJAA, 3, 125
14
Pigulski A., Kolaczkowski Z., 1998, MNRAS, 298, 753