Info

- The k mechanism. The first term in (15.70) shows the major influence of the opacity on pulsations. Positive values of KT and KP imply a higher opacity in the contracted stage of a pulsating star. More energy is retained and converted into mechanical energy. This produces pulsation driving at the expense of the radiative flux. In the expansion phase, the release of energy would be increased. This is just the k mechanism (Sect. 15.1).

- The partial ionization. In current stellar envelopes, one has Kramers' opacity, thus Kp = 1 and Kt = _4.5 with C = 1/Vad = 5/2. Thus, _ (C + kp) = 4/5 and this favors stability. However, if C becomes large as in a ionization zone (cf. Fig. 7.1 for Vad), the function _ (kt/C + kp) may be negative and has a destabilizing effect. Simultaneously, if C becomes very large, it also reduces the positive term 4/C and this also favors the instability. These destabilizing effects of a large C are sometimes called the y effect.

- Spherical geometry. The term _4/(3f[) is always destabilizing, it is absent in a plane parallel geometry. In view of the definition (7.57) of II, this term expresses the change of density during compression. In plane parallel geometry, the radiation would "see" the same column density of gas at all stages during a pulsation cycle. In spherical geometry, compression increases the column density radiation has to go through. Thus, more energy is retained at compression in spherical geometry and this is destabilizing.

- Radiative damping. The term 4/C is always stabilizing. It is due to the radiative damping which removes heat and acts like a viscous drag on the pulsations. One can realize the physical meaning of this term from the fact that 1/C = Vad expresses the change of T during compression. In an isothermal medium C ^ ^ this term would be zero as there is no radiative energy loss by pulsations. As seen above, a large C reduces this stabilizing effect.

In Fig. 15.4, the three terms of relation (15.70) are represented as well as their sum. The term _ (kt/C + kp) is stabilizing in the outer layers and destabilizing in the deeper envelope. If there would only be the dominant term kt , the destabilizing region would essentially cover the left zone of the k peak (cf. Fig. 15.1), where this term is positive. However, in regions deeper than the opacity peak (r3 _ 1) ^ 0, C ^ ™ according to (15.52) and the role of kt, which would favor stability there, vanishes. With kp still positive, this effect shifts the instability zone to deeper layers, as already seen in Fig. 15.1. The geometrical term _4/(3Ti) is always destabilizing, while the radiative damping is stabilizing as seen above.

15.4.1.1 Ionization Effects

The value of C (15.51) and (15.52), which varies a lot during ionization, plays a great role in the stellar stability. The value of C affects the coefficients A, B, D. Partial ionization by increasing C favors the instability, but it also has another effect: surprisingly, it limits the amplitude of the pulsations.

Let us consider a star in which the pulsations become large enough so that a significant mass shell is fully ionized in the contracted stage and recombined in the

4/C

" A í \ S\J 1

~ (KV'C+Jip)^,--'

i i i 'i

Fig. 15.4 The three terms of condition (15.70) for stability in a Cepheid model of 7 [email protected] and their sum. The stellar surface is on the left, the white area is the atmosphere with optical depth < 2/3, the gray area is the optically thick envelope. Adapted from Schaller [512]

Fig. 15.4 The three terms of condition (15.70) for stability in a Cepheid model of 7 [email protected] and their sum. The stellar surface is on the left, the white area is the atmosphere with optical depth < 2/3, the gray area is the optically thick envelope. Adapted from Schaller [512]

expanded stage. This would imply that C is small in these two leading phases, and we have seen that a small C favors stability. Thus, ionization, apart from having a destabilizing role, may also become an effect which may limit the growth of linear pulsations when the amplitudes become large. In addition, non-adiabatic effects, which lead to a heat leakage from pulsations, also limit the pulsation amplitudes.

15.4.1.2 A Straightforward Approach

The effects controlling the instability can also be understood in a rather direct approach. One has from the expression of i' given by (15.42), i' = 4r'- Kpp'- (kt -4)t'. (15.71)

Using the adiabatic approximations to express p' and t' as functions of q ', we have t' = (r3 - 1)q' and p' = r1Q', we get, ignoring the displacement term, i'^ (4 - kt)(r - 1)q' - KpT1Q' = [(4 - kt)(r - 1) - kq] q'. (15.72)

Let us consider a compression phase, thus q' is positive. For Kramers' opacity with kt = -3.5 and kq = 1, one has a positive i', meaning that more luminosity is leaving the shell than entering into it. Thus, no mechanical energy is available for pulsation driving and the star is stable, this is the general case. Thus, positive values of (-kt) and (-kq) contribute to the stability, in agreement with (15.70). Conversely, if (-kt) and (-kq) are negative, less luminosity is going out, more is retained in the shell and exerts pulsation driving. The destabilizing effect of partial ionization is also present through the term (r3 - 1), the y effect.

15.4.2 The Damping Timescale of Pulsations

We can gain further insight in expanding the solutions of (15.55) as functions of the coefficient K (15.44), which determines the non-adiabatic effects, s = so + s1K + s2K2 + (15.73)

We have the roots s0 = ±i\^Bo0 (15.60) for the adiabatic case with K = 0. Inserting (15.73) in the cubic, we get keeping only the first-order term in K,

3 s^ + c02Bs1 + o03D + OoAsl = 0 , s1 (—3o02B + oQB) = —o^D + c3AB , (15.74)

where the above expression of s20 is used. This yields s1 = — 2B (AB— D) . (15.75)

If AB — D > 0, which is our condition (15.63), we have s1 < 0 which implies a decreasing amplitude in time and thus stability with respect to pulsation. This justifies the above interpretation of this condition. AB — D > 0 implies the stability condition (15.70), provided the dynamic stability with B > 0 is also satisfied. If not, the amplitudes increase with time.

From the representation of the pulsation (15.54) and (15.73), it is clear that the damping time rd of the pulsation is of the order of (s1K), i.e.,

1 mpo mkto Eh

where we used (15.44) and the law of perfect gas. There, Eth is the thermal content of the layer of mass m. Thus, the damping time is the thermal diffusion timescale, i.e., the time needed for the luminosity to radiate all the heat content of the layer. For the whole star, this is the Kelvin-Helmholtz timescale (Sect. 3.2.4). A similar timescale intervenes for the growth of pulsations due to the non-adiabatic effects.

The second-order term s2 in (15.73), which we do not calculate here, would be imaginary. This indicates that there is only a second-order correction to be brought to the adiabatic frequency calculated above (15.60). Thus, as there is no first-order correction, the adiabatic frequency is a good approximation.

The timescale characterizing the non-adiabatic effects for the whole star is given by Td (15.76). The relative importance of non-adiabacy for a single layer during one pulsation cycle is given by the ratio (P/rd) of the pulsation period to the damping time. If the period is short, the non-adiabatic effects are small. Since rd « Eh/l0, this means that the non-adiabatic effects will be small when the heat content of the layer is high. Thus, the deep layers have in general little driving or damping effects for the K mechanism, so these effects are determined by the outer layers.

15.4.3 Secular Instability: Conditions on Opacities and Nuclear Reactions

For now, we have the adiabatic solution sad or so, the pulsational solutions with s1. There is also a so-called secular solution sS to (15.55). It is obtained if s is small, i.e., if the timescale of the perturbation is long. The cubic becomes in this case

and the solution is

If D > 0 for r1 > 4/3, the root is real and negative, which implies stability. From (15.78), the timescale of the instability is also

as for the damping time seen above (15.76), i.e., the timescale for the whole star is the Kelvin-Helmholtz time.

What is this secular instability ? It is a disequilibrium between the luminosity available from below at some level and the ability of the opacity to let this flux go out. If radiation cannot escape, there is an expansion at the timescale ts, or a gravitational contraction in the opposite case. The process is very slow with respect to current pulsations, while it is fast with respect to the nuclear timescale. This is why it is called a secular instability. It was already studied by Jeans [274] and by Ledoux [319]. The criterion for stability (for 71 > 4/3) is just D > 0, i.e.,

For a perfect gas with a = 5 = 1, the stability condition becomes kt + 4 kp > 0 . (15.81)

For the Kramers law with kt = -4.5 and kp = 1, the stability condition is not verified. This is normal, since we have not included the energy production. A radiative envelope deprived of energy source contracts at its Kelvin-Helmholtz time.

15.4.3.1 Why Are the Main Nuclear Burning Phases Stable ?

The best known secular instability at the Kelvin-Helmholtz timescale is the fast evolution toward the red giants at the end of the MS phase. It is interesting to see that it can find a consistent explanation in this context. The complete form of the criterion including the energy production rate £ in the equations lead to the following

15.5 Relations to Observations: Cepheids 391 stability condition [319] for a perfect gas, kt + 4KP > _3eQ _ £T , (15.82)

This criterion shows that for stability there must be some condition on the dependence on (q, T) of the opacity with respect to the (q, T) dependence of the production rate. For a given opacity, the positive sensitivity of the nuclear reactions to q and T must be strong enough to ensure secular stability. Conversely, for given nuclear reactions, the growth of the opacity with density and temperature must be as large as possible. For Kramers opacity and nuclear reactions with eQ = 1 and eT > 4, the above condition shows that the secular stability is largely ensured. This is why stars are secularly stable in current nuclear phases.

In a stage where gravitational contraction dominates, the rate of energy production eg is given by (20.6) for a mono-atomic perfect gas. Assuming homologous contraction with r/r = _(1/3)q/q, one gets

Here, eQ = _1 and eT = 1, one sees that the criterion (15.83) is not verified and the star is secularly unstable. In addition, if due to low T in the outer envelope, the gas is neutral, then kp « 0 according to Sect. 8.6.1 and this also favors the instability. This is quite consistent, when a star starts contracting it cannot get out of the instability until nuclear reactions take over with their strong dependence on T.

15.5 Relations to Observations: Cepheids

There are many kinds of variable stars: the main categories are shown in Fig. 15.5. Their properties have been reviewed for example by Gautschy and Saio [203, 204]. Some stars in this astrophysical zoo show radial oscillations like the Cepheids, other ones have non-radial oscillations like the slowly pulsating variables (SPB) and the solar-like variables. As an illustration, we discuss here some properties of the Cepheids. The non-radial oscillations of solar-like variables are extensively discussed in the context of Chap. 16.

15.5.1 The Period-Luminosity-Color Relations

There is a fundamental relation between the periods and average densities of pulsating stars, which leads to the period-luminosity-color (PLC) relations and in turn to

Wolf Rayet Star Diagram

Fig. 15.5 HR diagram with the principal categories of variable stars according to current literature superposed on evolutionary tracks of the Geneva group [513]. These are the Wolf Rayet stars (WR), luminous blue variables (LBV), red supergiants (RSG), P Cephei (P Ceph), Cepheids (Cepheids), slowly pulsating B stars (SPB), RV Tauri (RV Tau), Miras, irregular (Irr), blue variable subdwarfs (sdBV), RR Lyrae, 5 Scuti, roAp, y Doradus (y Dor) and the solar-like stars studied in asteroseismology. The variable white dwarfs are not represented

Fig. 15.5 HR diagram with the principal categories of variable stars according to current literature superposed on evolutionary tracks of the Geneva group [513]. These are the Wolf Rayet stars (WR), luminous blue variables (LBV), red supergiants (RSG), P Cephei (P Ceph), Cepheids (Cepheids), slowly pulsating B stars (SPB), RV Tauri (RV Tau), Miras, irregular (Irr), blue variable subdwarfs (sdBV), RR Lyrae, 5 Scuti, roAp, y Doradus (y Dor) and the solar-like stars studied in asteroseismology. The variable white dwarfs are not represented the famous period-luminosity (PL) relation for Cepheids, which is of major use for the distance calibration in the Universe.

Apart from second-order effects, the non-adiabatic periods are the same as the adiabatic values (Sect. 15.4.2), thus one can safely use the adiabatic relations. From (15.22) and (15.60), the fundamental period of pulsation P0 is

For the entire star of radius R, mass M and mean density q we get

This is a fundamental relation, which shows that the stellar pulsation periods (and the dynamical time 1.28) vary like the inverse of the mean density. For the Sun with a mean density of 1.4 g cm-3, the period is of the order of 55 minutes, for a red giant it may reach a hundred days and for a white dwarf it is about 1 s. A decrease of 71, as due for example to radiation pressure, produces a slight increase of the periods. More generally, it may be noted that a relation of the timescale with --1/2 exists for all kinds of configurations in gravitational equilibrium, from the Universe as a whole to a planet.

The relation between period and mean density is often written as

where Q is called the pulsation constant. Q evidently depends on the overtone considered: one writes Q0 for the fundamental mode, Q1 for the first overtone, Q2 for the second, and so on (when we write Q, this applies to any pulsation mode). Table 15.1 gives some typical values of the pulsation constant.

Table 15.1 Examples of pulsation constants Q

Type of stars

ß values (days)

polytrope n =

0

0.116

polytrope n =

1.5

0.071

polytrope n =

3

0.0383

Cepheids

0.035-0.050

ß Cepheids:

ßo

0.0375

ßi

0.027

ß2

0.022

Sources: Polytropes from Ledoux and Walraven [317]. Cepheids from Saio and Gautschy [504]. ß Cepheids from Lesh & Aizenman [324].

The typical period ratio of the first overtone to the fundamental mode is P1/P0 = 0.755, for the second overtone to the fundamental mode P2/P0 = 0.605 and for the third overtone to the fundamental mode P3/P0 = 0.506 for a typical MS model of 1.6M©. These fractions are evidently not 1/2, 1/3, 1/4 because the sound wave traveling through the star spends relatively more time between the surface and the first node, since the sound speed scales like VT (C.27).

The pulsation periods of any harmonics scale with mass and radius like

R3/2

One can eliminate the mass M with the mass-luminosity relation (3.30) and the radius R with (C.15),

eff which gives

Te4ff

Teff

In the mass domain of the Cepheids from about 3 to 12 M0, one typically has a = 33 and thus

log Po «log Qo - 3 log Teff + 0-60 log —- + const.

or log L = 1-67 log Po + 5 log Teff - 167 log Qo + const'• (15.92) L0

The periods are given in days, if the constant Q is expressed in days. This is the period-luminosity-color (PLC) relation characterizing radially pulsating stars. It depends on the value of Q for the corresponding mode of pulsation and on the stellar structure considered. Caution must be given to the fact that Q varies through the HR diagram, even for the Cepheids (see Fig. 15.7 right), the above relation is only a rough approximation. If we take a typical value of Q0 = 0-0383, appropriate to a polytrope of index n = 3 (which provides a simplified stellar model, cf. Sect. 24.5.1), we can draw iso-period lines over the whole HR diagram, as illustrated in Fig. 15.6 for a fundamental pulsation mode with the assumption of a constant Q0. The range of densities is such that the periods extend from about 1 h for the Sun to a few hundreds of days for the most luminous red supergiants. This figure is useful to anticipate the order of magnitude of the pulsation periods of stars to be observed.

15.5.2 Physics of the Instability Strip

The Cepheids occur in a so-called "instability strip", i.e., a relatively narrow band in the HR diagram crossed by the stars with a mass between about 3 and 12 M0 as illustrated in Fig. 15.5. Cepheids pulsate in the fundamental mode (Q0) with periods from about 2 to 50 days. Let us discuss the reasons why pulsating stars are present in a narrow range of Teff. One can schematically distinguish three zones in the stellar envelope:

- The outer zone near the stellar surface has a very low density. Even if it covers a significant radius, this zone has very small mass and heat content. The luminosity

Teff

Fig. 15.6 HR diagram with lines of constant period in the fundamental mode for a pulsation constant Qo = 0.0383. The periods are given in days (label "d") or in hours (label "h"). The spectral types are given along the Main Sequence and at the top of the figure for the giants and super-giants. The number along each track indicates the corresponding mass. The location of the Sun is indicated, as well as the location of Cepheids, 5 Scuti variables and P Cephei, see also Fig. 15.5. Courtesy from G. Burki [78]

Teff

Fig. 15.6 HR diagram with lines of constant period in the fundamental mode for a pulsation constant Qo = 0.0383. The periods are given in days (label "d") or in hours (label "h"). The spectral types are given along the Main Sequence and at the top of the figure for the giants and super-giants. The number along each track indicates the corresponding mass. The location of the Sun is indicated, as well as the location of Cepheids, 5 Scuti variables and P Cephei, see also Fig. 15.5. Courtesy from G. Burki [78]

perturbations £' only have negligible energy exchanges there and thus they remain essentially constant, despite this zone may be formally non-adiabatic.

- At moderate depths, there is an intermediate non-adiabatic zone with significant mass and heat contents, so that this zone may produce a large driving or damping of the pulsation.

- At large depths, the heat content is so large that the exchange of energy with the pulsation has no effect on its heat content and the pulsation may be regarded as adiabatic.

The instability strip is determined mainly by the location in depth of the intermediate zone. Let us consider, as an example, a 7 M0 star moving horizontally from the left to the right in the HR diagram due to its internal evolution. At Teff > 8000 K, the regions of partial ionization of H, He are in the outermost zone, the partial ionization of He+ ^ He++ (T « 4 x 104 K) being the leading effect. Elements hydrogen and helium are fully ionized in the intermediate zone, so that a Kramers-like opacity dominates. Coefficients kt and kp are negative and the star is stable.

As further evolution reduces Teff, the He+ region first and then the He and H partial ionization regions enter the intermediate non-adiabatic zone. The adiabatic exponents r decrease and the opacity regime may have positive kt and kp, two reasons making the star unstable.

A further decrease of Tef brings the zones of partial ionizations in the deep adia-batic zone, where the destabilizing factors have little influence. Thus we understand that it is only over a limited range of Teff values that the stars are unstable: this is the origin of the instability strip. This strip covers a broad range of magnitudes.

In practice, most Cepheids occur in the part of the blue loops where the evolution is the slowest, this is generally the blue extremity of the loops. Thus, most Cepheids are present at the luminosity where the extremity of the blue loops lie in the strip. As shown by Fig. 15.5, different kinds of variable stars, in addition to the Cepheids, occupy various parts of the instability strip, which extends down to the white dwarf sequence.

It is rather straightforward to derive the luminosity-Teff relation in the instability strip [147]. The non-adiabatic zone driving the instability lies at a depth Am in mass below the surface, where the ratio Eth/P of the thermal content of the upper layers to the period P is of the same order as the luminosity L (cf. Sect. 15.4.2),

The pressure at the base of the layer of thickness Am is a fraction Am/M of the central pressure (1.20). In addition, one easily finds the pressure (cf. 24.23) at some level within a hydrostatic envelope of perfect gas with a Kramers opacity of the form k = k0q T-3 5. Equaling these two estimates of the pressure at the depth Am gives, ignoring the numerical coefficients,

R4 V ^muK0/ V L J This provides the following scaling for Am,

Here, T is a constant corresponding to the ionization temperature considered, i.e., about 4 x 104 K for He+. For a constant Q, we get from the two expressions of Am the following relation for stars in the instability strip,

independently of M. Since by definition L = 4rcR2oTe4f, the scaling with R5/3 implies that Tef - R-1/3, so that

There is little change of Teff of Cepheids with their average luminosity; this explains why the instability strip is almost vertical in the HR diagram. The above simple developments give a relation for the instability strip, log Teff = -0.05 log L + const" . (15.99)

This slope is remarkably close to that obtained from numerical models, which give for the blue edge of the instability strip a relation log Teff = -0.036 log L + 3.925 [504]. The theoretical width of the strip is about 0.05 dex in Teff, the observational width being slightly larger, of the order of 0.08 dex in effective temperature [507].

15.5.3 The Period-Luminosity Relation

Since we have the PLC relation (15.92) and also a relation L — Teíí for the instability strip (15.99), one can eliminate Teff from the two and get a period- luminosity (PL) relation for the Cepheids:

log L = 1.34 log Po — 1.34 log Qo + const"'. (15.100)

The relation obtained from numerical models is [504]

for the blue edge, while for the red edge the slope is 1.244 and the constant 2.326. It is satisfactory that the simple analytical developments give a slope close to that from numerical models. These relations are illustrated in Fig. 15.7 left. At a given luminosity, the period at the red edge is about 25% longer than at the blue edge, because the average density is lower. The average observed PL relation for galactic Cepheids is [507]

My = —(3.087±0.085)log P — (0.914±0.098) . (15.102)

This corresponds to a slope of 1.235 in the logL vs. logP diagram, in relatively good agreement with the theoretical values. This relation plays a considerable role in the calibration of the extragalactic distance scale. The numerical models of Fig. 15.7 (left) present no metallicity effect in the PL relation over the range of the metallic-ities from the SMC to the Galaxy (Z = 0.004 to Z = 0.020). While several works claimed there is no observable Z effect in the PL relation, recent results by Sandage and Tammann [507] show significant differences with Z in the observed PL relations, so that a specific empirical PL relation has to be used for each metallicity. The extreme case is that of the Cepheids of Population II, the W Virginis stars, which

Fig. 15.7 Left: period-luminosity relation for the blue and red edges from numerical models of composition X = 0.70 with Z = 0.02-0.04. Right: the pulsation constant Qo for the Cepheids as a function of the periods, for the blue and red edge (same symbols as on the left and same composition). Adapted from Saio & Gautschy [504]

Fig. 15.7 Left: period-luminosity relation for the blue and red edges from numerical models of composition X = 0.70 with Z = 0.02-0.04. Right: the pulsation constant Qo for the Cepheids as a function of the periods, for the blue and red edge (same symbols as on the left and same composition). Adapted from Saio & Gautschy [504]

are old low-mass stars with M < 1 M0. A Cepheid and a W Virginis at the same location in the HR diagram have the same radius, however, their mean densities are very different. The much lower density of the W Virginis stars leads to a period about a factor 3.5 longer at a given luminosity.

We have considered the Q0 values as constant. However, the Q values in general depend on the stellar structure and the Cepheids of different luminosities have slightly different structures, thus the Q0 values also depend on the luminosity or period of the Cepheids. Figure 15.7 (right) shows these effects: Q0 increases with periods with a small difference between Cepheids at the blue or red edge. The maximum difference reaches about 50% which is far from negligible for the period-luminosity relation. The physical reason for the variation rests on the behavior of the ratio of the sound speed cS ~ VT in the outer layers to the average sound speed in the star: at higher luminosity (and longer periods), the Cepheids have a larger thin envelope with low temperature so that the above ratio is smaller. Therefore, the pulsation constants Q0 are larger for more luminous Cepheids.

15.5.4 Light Curves

There is a great variety of Cepheid light curves with different bumps, in particular there are systematic differences according to the luminosity. The amplitudes around the mean luminosity vary a lot, from a few percents to 1 magnitude. Also, at a given location in the HR diagram the amplitudes are not necessarily the same. Figure 15.8 shows the typical variations of magnitude, velocity, radius and Teff of a Cepheid. The velocity curve is the mirror image of the brightness curve with the same degree of skewness, the peak in luminosity corresponds to the largest expansion velocity. The variation of radius is such that its derivative gives the velocity curve. The maximum v -2.6 -2.4 -2.2

6000

5500

! _L 1-1-1-1-1-1—

--1-1 >*'nI

brightness

_ velocity ^^^^

radius

i i . 1 -r

Pilase

Pilase

Fig. 15.8 Schematic representation of the variations of V magnitude, radial velocity, radius with respect to the minimum radius and effective temperature of a classical Cepheid (5 Ceph) over one period. Adapted from J.P. Cox [146]

of brightness occurs slightly after the minimum radius, i.e., at the point with the average velocity on the steep descending part of the velocity curve. This is the so-called phase lag. It is the variation of Teff which essentially determines the brightness variations, because of the dependence in Tf of the luminosity. Teff is maximum, as a result of gas compression, when the radius is minimum. Thus, globally there is a correspondence between the minimum radius and the maximum luminosity. However, in details the star reaches its maximum luminosity when the opaque layers associated to the zone of partial ionization are the thinnest and the lag is due to the most external location of the H-ionization shell slightly after the minimum radius.

The non-linear effects which determine the amplitude of the pulsation are not yet fully explored. The same is true for the large differences of amplitudes between Cepheids in the Galaxy and the Magellanic Clouds and the amplitude differences from star to star. There is also a long-standing mass-discrepancy problem, in the sense that the masses derived from the pulsation theory are lower than those inferred from the evolutionary theory. Improved opacities have reduced the problem [203], however not totally.

Cepheids are slowly rotating stars, however, some of them may have had fast rotation during the MS phase. For such cases, rotation has two indirect effects on the properties of Cepheids. First, rotation during the MS phases increases the mass of the stellar cores and makes the stars more luminous when they are in the Cepheid instability strip. The radius is also increased, the mean density lower and the period longer. Thus, rotation produces a shift both in luminosity and in period. There is a second effect of rotation: rotational mixing produces a self-enrichment of the surface layers in helium. This may influence the driving and damping of the pulsation with possibly some consequences for the amplitudes of the pulsations. On the whole, we can say that the consequences for the Cepheids of the various rotational instabilities and mixing processes are far from being well explored.

Was this article helpful?

0 0

Post a comment