## A note on Cepheid variables: reproducing the study of the period-luminosity laws

In the 13th year of our life, when the skies of our city were still tolerably dark, we observed two stars $\beta$ Persei (Algol) and $\delta$ Cephei from our balcony which faced north on every night they were visible. Thus, we reproduced the beginnings of one of the great, but less-appreciated, stories in science that took place in the late 1700s. The story of variable stars, i.e. stars whose brightness changes over time, goes back even earlier. Bullialdus in the late 1600s observed the star $o$ Ceti (Mira) for 6 years and established that it varies dramatically in its brightness (magnitude) over a period of 333 days. Presciently he declared that the variability was because of dark regions on the star that periodically presented themselves to the earth even as the sunspots. Newton in his Principia adopted the same explanation for the variability in Mira and subsequently after further observations Wilhelm Herschel too believed this to be the cause for stellar variability. Today we know that it is not the case for the Mira variables but then star spots do play a role for rotating variables like BY Draconis.

The first evidence for other mechanisms of variability came up over 100 years later when the astronomers Edward Pigott and young, deaf Goodricke studied $\beta$ Persei (Algol) and $\delta$ Cephei. They concluded that the variability of the former might be due to an eclipse by a darker body. Thus, the theory of eclipsing binaries was born and was subsequently confirmed for Algol-type variables. For $\delta$ Cephei Goodricke proposed the starspot mechanism. $\eta$ Aquilae and $\delta$ Cephei discovered by Pigott and Goodricke were the founder members of a whole class of regularly varying stars like it which were subsequently discovered: the Cepheids. The Cepheids, unlike the Miras tended to be of shorter period — several days rather than of an order closer to an year. Further, the Miras, unlike the Cepheids, while having a general periodicity did not strictly adhere to it both in terms of temporal period and actual magnitudes attained at maxima or minima. Like the Algol-type variables the Cepheids instead showed a strict periodicity in their light variation (the light curve) and also the maximum and minimum brightness they attained (amplitude). However, unlike the Algol-type variables but like the Miras they continuously varied in their brightness. This form of variability completely baffled the astronomers and physicists, who nevertheless came to believe in some kind of eclipse mechanism shortly after their discovery based on the strict periodicity they shared with the Algol and $\beta$ Lyrae type eclipsing variables. This idea remained in force for more than a century. In the late 1800s the German physicist August Ritter carried out one of the first detailed studies on gas spheres. In course of that he almost prophetically predicted an alternative mechanism for variability of stars i.e. pulsation; however, it was totally ignored by astronomers.

But at the turn of that century Schwarzschild’s study of the Cepheid $\eta$ Aquilae showed that its visual light variability was also accompanied by color (spectral type) variability. In 1908 Leavitt, while studying variables from the Small Magellanic cloud, discovered that for some variables with a regularly periodic light curves their period was directly correlated with their apparent luminosity (given that all stars in a Magellanic cloud were expected to be at roughly the same distance from our solar system). By 1912 she had established that this period-luminosity law specifically applied to Cepheids. Shortly thereafter Shapley marshalled various lines of evidence to suggest that the variability of Cepheids was likely not due to eclipses but due to pulsations. In the coming years Eddington developed the first physical theory for pulsating Cepheids by modeling them as thermodynamic engines that operate by pulsation, wherein the energy gain and dissipation compensate each other. Subsequent observations of the radial velocity of Cepheid derived from spectroscopy by Baade and others helped confirm the pulsation theory and also refined it further even as the physics of the stars was being better understood. In this context we should point out the key contribution of the sadly forgotten Hindu astronomer Hari Keshab Sen from Prayag, who was also a student of Vedic literature and studied the origin of the Hindu concepts of infinity in the upaniṣad-s. In 1948 while at the Harvard observatory, which was the great center of variable star research, he presented the first refined model for Cepheid variability: He showed that to explain the Cepheid variability one needed to develop a correct model for the internal structure of the star. He showed that for a star the ratio of the central gas density to its mean gas density was of key importance to explain variability. His work showed that this is particularly high Cepheids and using this in the model could explain the asymmetric amplitude of the Cepheid variability.

In years following the discovery of the period-luminosity law it was used by Hertzsprung to calculate the distance of the Magellanic clouds from the solar system. While this attempt had its errors it was the beginning of one of the most important developments in all of science: the use of Cepheids to measure the size of our universe, starting from the realization that the Magellanic clouds lay outside our own galaxy. The first step in correctly applying the period-luminosity law to measure the universe was the realization that the Cepheids were not a monolithic group but included a range of distinct families each obeying their own period-luminosity relationship:
1) W Virginis or type-II Cepheid stars: These are old stars with mass equal to or lower than the Sun but are distended having radii tens of times the Sun. They tend to be of the spectral types from mid F to early K and about -.5:-2.5 in absolute magnitude. They are mostly found in the central galactic spherical component and globular clusters.

2) RV Tauri variables: These tend to be more luminous giant versions of the W Virginis stars with an absolute magnitude in the -4:-5 range and are redder with more G to later K-type spectrum. Their light curves show alternating deep and shallow minima and concomitantly their period-luminosity laws are much more lax than the classic or type-II Cepheids.

3) The RR Lyrae variables: They are sub-giants usually of the A-F spectral class again with a low mass of about half the Sun. They too are old stars which have previously lost gas as red giants before climbing down to just above the main-sequence in the Hertzprung-Russell diagram. They too are present in the spherical center of the galaxy and globular clusters.

4) $\delta$ Scuti variables: These are whitish stars in the A to early F spectral class and are characterized by a low amplitude and short-period pulsation (roughly < 7 hours). They are close to the main-sequence.

5) Classic Cepheids (DCEP): These are also known as type-I Cepheids and are prototyped by the original members like $\eta$ Aquilae and $\delta$ Cephei. They are massive, giant stars: 4-20 Sun masses and with typical absolute magnitudes of -1.75:-6.4, i.e ~500:31,000 times the brightness of the Sun (the brightest being 200,000 times as bright as the Sun). They are young, metal-rich stars which in galaxies like ours are mostly in regular elliptical orbits around the galactic center. Within the classical Cepheids, we might further recognize the DCEPS group which includes representatives like the current Pole Star $\alpha$ Ursae Minoris which show an amplitude of visual magnitude pulsation less than .5 with generally symmetric light curves.

These Cepheid-like stars define a strip across of the Hertzsprung-Russell diagram known as the instability strip which mark the region in it where an evolving star develops period-luminosity law-type pulsation. Like any oscillator a pulsating variable star can pulsate in the fundamental mode of oscillation or as a higher overtone. The mode can be understood thus: if it is pulsating in the fundamental mode then all of the stellar sphere moves in and out at the same time -i.e. all of the star contracts and expands in unison. If it is pulsating the first overtone then we imagine a central shell, the nodal sphere, within the star where the gas is at rest. When the gas inside that shell is contracting then gas outside it expands and vice versa. In the second overtone pulsation we have two concentric nodal spheres with the gas inside the outer and innermost shells contracting when that inside the middle shell expands and vice versa. Both fundamental mode and overtone pulsation are seen across Cepheids types.

In the remainder of this note we shall explore the period-luminosity law for classical Cepheids using the wealth of modern data from ground- and satellite-based observations. To start we take a data-set of 674 Cepheids from the Milky Way for which we have reliable period, mean visual (V band, $m_V$) magnitude and parallax-derived distance information. Using the distance in parsecs $d$ we calculate the absolute magnitude ( $M_V$) of the stars using the formula:
$M_V=m_V+5-5\log_{10}(d)$

We then plot $M_V$ against the base 10 logarithm of the period (in days) of the star to see if we can recapitulate the period-luminosity law (Figure 1). In Figure 1 the regular classic Cepheids are in blue. The DCEPS group is further shown in red. Those DCEP that do have an overtone pulsation are colored violet and those DCEPS with overtone pulsation are colored green.

Figure 1

What we see in Figure 1 is that there is a general positive period-luminosity correlation. But by no means it can be termed a strong law
with a low correlation, $r^2=.22$. Nevertheless we can observe a few interesting things: the DCEPS are not found at the high end of the luminosity range but they are dominant at the lower end of the luminosity range. Similarly, the overtone pulsators are more at the lower side of the range. In any case with a correlation like this one can hardly expect to use Cepheids as an accurate measure of distance. So what is the issue here? The answer to this lies in the observation that was first made in the late 1800s that light from a star does not reach us unimpeded: it is scattered by dust in interstellar space. This dust selectively scatters the lower wavelength light towards the blue end of the spectrum more than that from red end. Details of this extinction of starlight were worked out only in the last century. It can be quantified by using the measure $E_{B-V}$ which is the observed excess of the difference between the B band and the V band filter magnitudes relative to the real $B-V$ color index.

We have these $E_{B-V}$ values for a high-quality data-set which is a subset of 383 brighter stars of our starting data-set. For this we additionally have the magnitudes measured in several photometric bands in addition to the V band. Hence, to start we make the same plot as above of $M_V$ computed from distance in parsecs (Figure 2).

Figure 2

We note that even though we have only half the data as in the first plot and we have lost some of the lower luminosity stars the results are similar. This suggests that the what we observed above is generally consistent across the data: a low correlation $r^2=.23$.

Now we make use of the $E_{B-V}$ values for these stars applying a multiplication factor of 3.23 and using the below formula:
$M_V=m_V-3.23E_{B-V}+5-5\log_{10}(d)$

Now we plot this extinction corrected $M_v$ against base 10 logarithm of the period (Figure 3).

Figure 3

Now with a correlation of $r^2=.92$ we get an unambiguous period-luminosity law for the classical Cepheids ranging from -1.75:-6.4 $M_v$. The tendency of the low-amplitude DCEPS to be excluded from the high end of the magnitude-period range is also clear.

In the above approach we solved the issue of extinction of starlight by interstellar dust on a star-by-star basis using $E_{B-V}$ for each star. To tackle this directly rather than for each case the astronomer Madore defined a new magnitude function called Wesenheit $W$ that is free of the extinction effect. Wesenheit can be calculated by measuring the magnitude in multiple photometric bands e.g. V and and the near-infrared I band (806 nm midpoint). For illustration with our current data-set we shall use the V and I band Wesenheit. For that we first calculate the near-infrared absolute magnitude from the apparent $m_I$ magnitude of the star which has been obtained using the I band filter:
$M_I=m_I+5-5\log_{10}(d)$

We then couple this with the uncorrected $M_V$ calculated from $m_V$ as above using the below Wesenheit formula:
$W=M_I-1.55(M_V-M_I)$

We then plot this against $\log_{10}(P)$ as before (Figure 4).

Figure 4

We see that we now get an even better period-Wesenheit law than with the case-by-case extinction correction: $r^2=.98$. Thus Wesenheit is a very effective means of capturing the Cepheid period-luminosity relationship in an extinction independent manner.

Figure 5

Finally, if we look at a histogram of extinction-corrected absolute magnitudes for these Milky Way Cepheids (Figure 5) we find an interesting distribution with a clear right skew and hint of bimodality or a thick right tail. The skew is likely in part an artefact of losing the lower luminosity overtone pulsators. There is a clear peak at $M_V=-3$, second peak or shoulder at around $M_V=-3.8$ that leads to a fat tail. This suggests that at least in this dataset classical Cepheids themselves might have at two sub-categories based on luminosity.

We next turn to classical Cepheids in the Large Magellanic cloud and the Tarantula nebula. For this we use the data collected by the ground-based Visible and Infrared Survey Telescope for Astronomy with a 4.1 meter mirror operating from Chile and released in 2012. This data-set contains 334 Cepheid stars from the LMC and the Tarantula Nebula with high quality period and magnitude data in the V, I and the $K_s$ (deeper near-infrared: 2190 nm midpoint) bands. It also has information regarding whether the star is pulsating in the fundamental mode or has overtone pulsation. However, there is no direct distance data: the LMC being extra-galactic is out range for direct parallax measurement. Hence, we shall use this data to perform a modern recapitulation of the efforts starting with Hertzsprung. However, in this case like Leavitt did we can assume that all the stars are at a generally similar distance from our solar system as they are from a single external galaxy. Thus, to investigate the period-luminosity law we may directly plot apparent magnitude $m_V$ against $\log_{10}$ period (Figure 6).

Figure 6

We see some interesting points right away: 1) We get a much better period-luminosity relationship even with no correction for starlight extinction: correlation $r^2=.58$. The LMC being an external galaxy would be subject to roughly similar extinction across all its stars. Hence, the deviations caused by extinction have less consequence for the basic period-luminosity relationship in this case. 2) As we saw in the case of the Milky Way, the overtone pulsators (plotted as red asterisks) occupy the lower-end of period-luminosity range while fundamental tone pulsators the higher end (blue dots). Notably they share a clear zone of overlap in the middle of the range. 3) Closer examination of the two pulsation modes hints that the two groups of stars might actually lie on separate lines.

We next utilize the $K_s$ magnitude to plot a period-luminosity for it (Figure 7). Being in the deeper near-infrared we expect it to suffer much less extinction than the visual band of the spectrum.

Figure 7

Sure enough we see that the $K_s$ magnitude shows much less dispersion than the V band magnitude. Strikingly, what we saw as less-clearly with the visual magnitude, comes out plainly in the $K_s$ band: not only are the overtone (red asterisks) and fundamental mode (blue) pulsators separated by the range occupied in the period-luminosity relationship, they also lie on two distinct, roughly parallel lines which appear shifted by an approximately constant factor.

We next utilize the V-I Wesenheit as per the above-stated formula to study the same thing (Figure 8).

Figure 8

Again, Wesenheit being free of the extinction factor gives a clean period-luminosity relationship and clearly shows the separation of the fundamental and overtone mode pulsators into two nearly parallel period-luminosity relationships.

Now that we have the period-luminosity laws for the LMC Cepheids in place we next use it calculate the distance of the LMC. For this we first need to calculate the absolute magnitude M for these stars from the period. For this we can try out any of a number of formulae:

$M_V=-2.43(\log_{10}(P)-1)-4.05$
This formula was based on relatively small sample of Cepheids studied using the Hubble space telescope.

$M_V= -2.76(\log_{10}(P)-1)-4.16$
This formula using the visual band magnitude was derived by Madore and Freedman in their seminal article on classical Cepheids.

$M_I= -3.06(\log_{10}(P)-1)-4.87$
This formula was provided for the infrared I band by Madore and Freedman in the same article.

$M_V=\dfrac{-\log_{10}(P)+.657}{.394}$
This formula was based on calibration of Milky Way Cepheids using the Hipparchos space observatory’s parallax measurements.

Once we get $M_V$ or $M_I$ with these formulae we next compute the distance modulus:
$\mu=m_V-M_V$
From the distance modulus we can easily get distance in parsecs as:
$d=10^{\tfrac{\mu+5}{5}}$

Since we saw considerable dispersion of the $m_V$ values (Figure 6) and given that the corresponding $M_V$ computation by these formulae will have the issues of extinction we shall instead resort other formulae that can correct for this using both $m_V$ and $m_I$ to compute distance modulus directly:
$\mu=m_V+3.746(\log_{10}(P)-1)-2.523(m_V-m_I)+5.959$
This is the formula of the team of the astronomer Sandage. It seems to underestimate $\mu$.

$\mu=m_V+3.255(\log_{10}(P)-1)-2.45(m_V-m_I)+5.899+.08$
This is the formula of the Freedman team which is what we shall use for our purpose. Having calculated $\mu$ we get distances for the stars from it and plot the distribution (Figure 9).

Figure 9

Thus, we get the median distance modulus $\mu=18.5$ and median distance for the LMC as 50,086.8 parsecs. However, we notice that there is a clear bimodal distribution of the computed distances with peaks around 54,000 and 42,000 parsecs. Why do we see this and which of these is correct? The answer lies in the earlier period-luminosity plots where we saw two distinct laws for fundamental mode and overtone pulsators. This was the key finding detailed by the astronomer Böhm-Vítense who suggested the overtone pulsators are shifted by $\Delta \log_{10}(P) \approx .15$ with respect to the fundamental mode Cepheids. Hence, to apply a uniform period-luminosity law we have to transform the observed period $P$ of the overtone pulsators to their fundamental mode period $P_0$ by the formula:
$P_0=\dfrac{P}{.716-0.027\log_{10}(P)}$

If we use this transformation in the above formula of the Freedman team we get the distance distribution plot seen in Figure 10.

Figure 10

We now get a sharp, single-peaked distribution with a mean distance modulus of $\mu=18.64$ and median distance of 53660 parsecs for the LMC. These values are within the range of values published in the last two decades. Most authors these days favor a distance modulus of $\mu=18.5$. Thus, even though the period-luminosity law looks like an obvious way to measure the distance of the universe and even estimate the expansion rate of the universe in the form of Hubble’s constant, actually doing so presents several challenges some of which we have just seen in the above experiments.