© D. Gillet et al. 2019

1 Introduction

The kinematics and dynamics of the outer layers of pulsating stars, a fortiori those undergoing Blazhko effect (Blazhko 1907), is usually not known in detail. Among variable stars, RR Lyr is the brightest RR Lyrae star of the sky with 7 ≲ V ≲ 8. During its pulsation cycle, its photospheric radius relative variation ΔR_⋆∕R_⋆ = 17% (Fokin & Gillet 1997), while is spectral type evolves from A7III to F8III (Gillet & Crowe 1988). RR Lyr variability was discovered by the Scottish astronomer Williamina Fleming at Harvard in 1901 (Pickering et al. 1901). The pulsation cycle occurs on approximately 0.5668 d or, equivalently, 13.6 h. Furthermore, the light curve also presents light modulations with a variable period around 39 d. These amplitude and phase modulations are known as the Blazhko effect. Its physical origin has remained a mystery up to the present day, but recently several interesting ideas were proposed to explain the Blazhko phenomenon (Smolec 2016; Kovács 2016). All these phenomena are the rendering of atmospheric movements that are not simple outward and inward movements that may show light curves (see RR41 run of Fokin & Gillet 1997).

In order to study the atmospheric dynamics, the whole pulsation cycle has to be recorded in detail. The first detailed spectroscopic study of RR Lyr was done by Preston et al. (1965) during a Blazhko cycle, but this investigation was limited to the pulsation phases of rising light within the 13.6 h cycle. Chadid & Gillet (1996) carried on spectroscopic observation and showed metallic line doubling in RR Lyr. They suggested conducting some extra observations with a larger telescope to get a higher signal-to-noise ratio (S/N) and suitable time resolution. Preston (2011) and Chadid & Preston (2013) performed extra spectroscopic observations. They studied dynamics of various RR Lyrae stars, but did not focus on RR Lyr itself. Fossati et al. (2014) carried out a spectroscopic abundance analysis of RR Lyr over a pulsation cycle. This analysis was centered on the so-called “quiet phase” around the maximum radius. In particular, they show that the microturbulent velocity V_mic close to the photospheric level is not significantly affected by the pulsation and remains almost constant during the pulsation cycle while, there is an increase in V_mic around the bump phase and again on the rising branch. During these two phases, strong shocks propagate rapidly through the atmosphere. Fossati et al. (2014) confirm the presence of a strong peak of the turbulent velocity around the pulsation phase 0.9. Nevertheless, none of these historical surveys gave a description of the whole atmospheric dynamic structure during a pulsation cycle. A recent paper (Gillet et al. 2017) presents an additional spectroscopic survey dedicated to the “third apparition.” Since they saw an interesting behavior in the sodium doublet region, we decided to conduct a more intensive run in order to get a more detailed view of RR Lyr dynamics.

In order to find out the origins of the atmospheric layer dynamics, it is necessary to compare observations with relevant atmosphericmodels, but this is only possible if an observational survey with a high temporal resolution during the whole pulsation cycle is available. There are many models that describe photosphere evolution. Unfortunately, only a small number of them integrate a description of the atmosphere and the photosphere together with an adapted treatment of shock waves. Feuchtinger (1999) developed a nonlinear convective model of pulsating stars with shock waves based on the Vienna model. Computed light curve morphology is in good agreement with observed RR Lyrae light curves,including the well known RR Lyrae phase discrepancy problem (Simon 1985). Recently, Geroux & Deupree (2013) implemented a two-dimensional (2D) radiation hydrodynamic code that simulates the interaction of radial pulsation and convection. The code is able to match the observed light curve shape near the red edge of the RR Lyrae instability strip the in HR diagram. However, neither of these models is able to compute a spectral line profile.

While convection is a tricky phenomenon to model, only a few pulsation models are able to calculate line profiles. Hill (1972) built a nonlinear hydrodynamic model of atmosphere with shock waves which leads to a successful comparison of light curves, absorption line velocities, emission strengths, and shock wave dissipation for the RR Lyrae star X Ari. Stellingwerf (1975) described a nonlinear hydrodynamic model atmosphere with shock waves taking into account dissipation. It is able to compute light curves,growth rates, layer temperatures, and velocities during the pulsation cycle. More recently, Stellingwerf (2013) developed an adiabatic three-dimensional (3D) pulsation model for RR Lyrae stars that he applied to RR Lyr itself. It considers the outer boundary condition and the interaction between pulsation and turbulence, but this model does not include radiative effects such as radiation and ionization. In these conditions, turbulent mixing and the dynamics of an extended expanding atmosphere produce shocks with extreme intensities. However, such large intensities have not been observed yet.

A work dedicated to RR Lyr was done by Fokin & Gillet (1997). They built a series of nonlinear, nonadiabatic models including atmosphere and shock wave propagation. In particular, they used a purely radiative approach because convection is thought to strongly affect the helium ionization zone of RR Lyr according to its effective temperature (Xiong et al. 1998). The resulting Fe II and Ba II line profiles have been confirmed by high-resolution observations such as line doubling phenomena according to the Schwarzschild (1952) mechanism. Lines profiles are computed at any pulsation phase according to predicted physical parameters such as temperature, density, or turbulent velocity occurring within the line-forming region. This seems to be correct when lines are formed over large areas of the atmosphere (absorption lines). On the other hand, when lines are produced in the radiative wake of the shock, such as Hα emission, line profiles are only approximate because non-LTE effects are not fully taken into account. In addition, new high-quality spectral observations around the relevant pulsation phases may involve improvements in this kind of model.

There are currently too few pulsation atmospheric models that deal with both convection and shock waves that are able to compute synthetic lines profiles. For the coming years, the development of such models is strongly needed.

We present here simultaneous observations of profile variations of sodium and Hα lines in RR Lyr itself. In Sect. 2 we describe observation data and their reduction processes. Observations analyses including the sodium line velocity measurement are presented in Sect. 3. Atmospheric implications of sodium observations are discussed in Sect. 4 and the dynamical structure of the atmosphere is described in Sect. 5. Section 6 is dedicated to shock wave dynamics. Finally, we draw our conclusions in Sect. 7.

2 Observations and data analysis

2.1 Data acquisition

In order to describe the atmosphere and shock wave dynamics, data were assembled from different sources. On the one hand, to achieve a complete movie of atmospheric layers over the whole pulsation cycle, we used time series spectra focused on Na D lines. On the other hand, for the shock wave study, since a 300 s time resolution was required to complete velocity measurements, we reexamined previous published spectra. Data used in this paper were obtained with three different spectrographs:

• ESHEL I spectrograph: we used an echelle spectrograph called ESHEL attached to an automated 35 cm telescope at the Observatoire de Chelles (Lemoult et al., in prep.). The fiber-fed ESHEL spectrograph is described in Thizy & Cochard (2011) and is built by Shelyak Instruments¹. The spectral domain is spread over 4300−7100Å for a resolving power of R = 10 500 at λ5896.92. The detector used is an Atik 460EX camera (Sony ICX694 CCD) on which the resolution element represents about 3.2 pixels. The exposure time was 900 s providing a mean S/N of 75. A total of 47 h of exposure time were gathered between April 7 and 18, 2017. In order to cover a full cycle, we selected79 spectra for a total of 20 h of exposures. These 79 spectra achieve a complete and continuous RR Lyrae cycle with a maximum of 900 s, i.e., a 2% cycle step. Observations were started at the Blazhko phase ψ = 0.77 to finish at maximum Blazhko ψ = 0.00 in a single Blazhko cycle.

• ESHEL II spectrograph: in August 2014, we observed with a similar echelle spectrograph ESHEL attached to a 35 cm telescope at the Observatoire de Foncaude, France. The detector was a CCD sensor Kodak KAF-3200ME. The resolution element represents about 4.3 pixels. The aim was to collect mid-resolution spectra of RR Lyr over 45 min around the outward atmosphere acceleration at φ = 0.91. Eight 300 s spectra were obtained over the spectral visible domain (4300−7100 Å) and with a resolution power of about 10 500 and a S/N of 50.

• AURELIE spectrograph: on October 12, 2013, we used the  AURELIE (Gillet et al. 1994) spectrograph (1200 t mm⁻¹, Blaze 5000 Å with OG590 order filter) recorded with a EEV 42-20 CCD on T152 at the Observatoire de Haute-Provence (OHP), France. In order to collect high-resolution spectra of RR Lyr, 600 s spectra were acquired over the visible domain (6513–6713 Å) and with a resolution power of about 22 700 and a S/N of 86. The resolution element represents about 2.8 pixels.

• ELODIE spectrograph: attached to the 193 cm telescope at the OHP, it is described in Baranne et al. (1996). Observations were performed in 1994−1997 in the context of a survey led by D. Gillet and published in Chadid et al. (1999). We used some spectra from these observations to verify the Na D red component at different pulsation phases. Typical exposure times were between 5−8 min, leading to a S/N of about 57 per pixel for a resolving power of R = 44 000 on the whole visible domain (3900−6800 Å).

The observations are summarized in Table 1.

2.2 Data reduction

AURELIE data were reduced with the SpcAudace package (Mauclaire 2017), part of Audela² software (Klotz et al. 2012). ESHEL spectra processing was done using a dedicated echelle sub-package of the Audela software too. ELODIEspectrograph uses its own MIDAS pipeline (Baranne et al. 1996).

These packages perform classical operations through pipelines, such as preprocessing (bias and dark subtraction, flat-fielding), masking of bad pixels, spectrum extraction, and wavelength calibration. Line profiles were corrected from instrumental response and heliocentric velocity. All spectra were normalized to the local continuum around the Hα line.

Wavelengths are presented in the stellar rest frame of RR Lyr using its actual radial velocity V_* = −73.5 km s⁻¹ with respect to the solar system barycenter (Chadid & Gillet 1996). V_* is also called the γ-velocity. All wavelengths (λ) are given in angström (Å).

2.3 Ephemeris computation

In this paper, the pulsation phase is denoted φ, where φ = 0.0 corresponds to the luminosity maximum. The Blazhko phase is denoted ψ and has its maximum at ψ = 0.0 when the highest luminosity amplitude is observed during a Blazhko cycle.

The period and reference date of pulsation and Blazhko ephemerides evolve with time. Up-to-date parameters are thus required. Le Borgne et al. (2014) showed that the pulsation period varies alternately between two states, defined as a pulsation over long (longer than 0.56684 d) and short (shorter than 0.56682 d) periods. These states alternate with intervals of 13–16 yr. After a sudden change in period value, this interval remains constant for a few years. The physical origin of these two states is not known yet.

We thus used the ephemerides provided by Le Borgne et al. (2014) for the AURELIE and ESHEL II observations made in 2013 and 2014. For the 2017 data, a late epoch for the ephemeris provided by Le Borgne et al. (2014), the pulsation period P_P was computed using the difference between the two closest O–C minima from the GEOS RR Lyr web database³ and close to our 2017 observations. The reference HJD maximum light (HJD0_P) was also provided by this database. A summary of the ephemerides used for the different data sets is provided in Table 2.

The determination of the Blazhko phase is more problematic since the Blazhko period P_B does not follow the variation in the pulsation period, but presents rather erratic changes (Le Borgne et al. 2014). Moreover, during 2014, the photometric O−C were historically low (close to nil), preventing an easy Blazhko period computation.

Therefore, in order to find maximum light amplitude dates HJD0_B, we developed an innovative method, which will be described in a forthcoming paper (Gillet et al., in prep.), based on the equivalent width (EW) and main shock velocities driven by our many years of experience in performing spectral observations. For each photometric maximum (φ = 0) corresponds to a spectroscopic maximum (φ = 0.91) where Hα presents an emission line. Consequently, |EW| and |V_shock| (defined in Eq. (3)) reach high values. They reach a maximum when the light amplitude O-C correspond to their extrema, thus pointing out Blazhko maxima. We expect an uncertainty of ± 2 days on HJD0_B and of ± 0.2 d on P_B. Ephemerides were computed involving the period P_B = 39.0 d (Le Borgne et al. 2014). A summary of the Blazhko ephemeris used for the different data set is provided in Table 2.

3 Observations of Hα and sodium in RR Lyr

We analyzed profile variations for three spectral lines: Hα, He, and Na D1 during a complete pulsation period. This last line is a resonance line that, following Fokin & Gillet (1997), is formed in the deep layers of the photosphere, while the core of the Hα and helium lines are formed well above the photospheric layers, hereafter called the upper atmosphere. Thus, these three lines, by their different formation heights, are able to probe and to follow the propagation of the waves in the atmosphere along the pulsation cycle.

Because one observing night is not long enough to obtain spectra during the whole pulsation period (P_P ≃ 13 h 36 min), we had to recompose the data of 11 consecutive nights. It was shown by Gillet & Crowe (1988) that lines could vary depending on the Blazhko phase with a period of P_B ≃ 39 d. Thus, as 55 spectra out of 79 were obtained over three days, on this timescale the Blazhko effect must be secondary in the framwork of our data. Consequently, the use of data obtained in different periods will not change appreciably the different phases of the atmospheric dynamics and the effects of the shock wave propagation on these different phenomena, including their impact on the turbulence amplification.

3.1 Hα line evolution

The three successive appearances of hydrogen emission (Preston 2011) have been observed in detail. Between φ = 0.892−0.929, just before the luminosity maximum, an intense blueshifted emission is present, the well-known “first apparition.” A smaller red bump is detectable from φ = 0.637 to φ = 0.736. Panel a of Fig. 1 zooms in on this “second apparition,” which is the weak hump at the top of the blue line wing around φ ≈0.73. Finally, another small red bump corresponding to the “third apparition” is detected for φ = 0.181−0.387. These three apparitions are compared in Fig. 1. While only one line profile is plotted for each apparition in this figure, they are observed in several consecutive spectra (see Figs. 2 and 3).

The blueshifted emission intensity corresponding to the first apparition (Fig. 1b) represents up to 30% of the continuum, depending on the Blazhko phase. The duration of the emission of the Hα line for the first apparition is between 15 and 40 min, i.e., less than 5% of the pulsation cycle. This phenomenon is rather short compared to the second and third apparitions (Figs. 1a and c, respectively). Their durations are respectively about 1.5 and 3 h, i.e., 11 and 20% of the pulsation period.

The second apparition was first observed in RR Lyr by Gillet & Crowe (1988). They measured a 30 min duration in nonconsecutive observations. This emission may depend on the pulsation cycle and requires observations over several cycles to be clearly defined.

On the other hand, the third apparition was first reported in RR Lyr itself by Gillet et al. (2017), where the emission intensity and duration are consistent with our observations. However, these two apparitions intensities represent only a few percents of the continuum.

Figures 2 and 3 show the evolution of the Hα line profile over the full pulsation cycle. They start at the bottom left and go toward the top right. The Hα line becomes redshifted at φ = 0.618. It is followed by a line doubling starting at φ = 0.943 and ending at φ = 1.027, i.e., ~ 8% of the pulsation cycle. Then the Hα line remains blueshifted until φ = 1.461, where it reaches its laboratory wavelength, and then becomes redshifted again. The shock intensity depends on the Blazhko phase: the more important it is, the more visible the He I emission lines become.

3.2 Helium line evolution

Figure 4 shows the He I λ5876 line in emission at φ = 0.911, i.e., the same phase of the first apparition of the Hα line. This emission was first reported in RR Lyr by Gillet et al. (2013). Its intensity is 8% of the continuum and the line velocity reaches − 6.5 ± 0.8 km s⁻¹, which confirm the Gillet et al. (2013) measurements. While this emission is a short phenomenon of less than 30 min, i.e., ~ 3% of the pulsation period, it is followed by an absorption shape between φ = 0.925−0.954. We note that the low value S/N_{He line}≅3 prevents us from performing a more detailed analysis of such a weak feature. No emission of the He I λ6678 line is observed at the same phase. This is consistent with the Blazhko phase ψ = 0.77 for this observation on the night of April 8–9, 2017. Indeed, this Blazhko phase corresponds to a strong shock, but not an exceptional one. Stronger shocks are expected around ψ = 0 (Gillet 2013).

Fig. 1 Comparison of the three apparitions of Hα line profile of RR Lyr. Panel a: second apparition at φ = 0.74 represented with a spectrum dated on April 10, 2017. Panel b: first apparition at φ = 0.91 represented with a spectrum dated on April 09, 2017. Panel c: third apparition at φ = 1.25 represented with a spectrum dated on April 08, 2017. Vertical line indicates the Hα line laboratory wavelength.

3.3 Sodium line evolution

Figures 5 and 6 show the evolution of the D1 sodium line profile (λ5895.92) over the fullpulsation cycle. Two strong telluric lines are blanketing the center of the D2 sodium component at λ5889.95 line (Hobbs 1978). Consequently, we focused our analysis on the D1 sodium component at λ5895.92 line. This line is denoted Na D in the paper. Figures 5 and 6 show the evolution of the D1 line profile over the full pulsation cycle.

A line doubling starts somewhere for φ ∈ [0.892;0.925], i.e., very close to the pulsation phase of the first apparition of the Hα line. It remains visible along φ ∈ [0.855;1.600], i.e., over 75% of the pulsation cycle. The origin and evolution of these two components of the D1 line profile are analyzed in detail in Sects. 4 and 5.

3.4 Line velocity evolution

From the whole time series composed of 79 spectra, we selected 19 spectra that led us to describe every step of the layer motion evolution during one pulsation cycle. In order to probe atmospheric dynamics we extracted the radial velocities of the sodium and the Hα absorption lines. Measurements were done with respect to the stellar rest frame. Line center measurements were obtained by Gaussian fitting or first moment of the line profile, depending on the line shape. The mean uncertainty is ± 0.6 km s⁻¹. These Doppler velocities are listed in Table 3 and presented in Fig. 7 where curves are computed with a cubic spline.

The sodium blueshifted component has a very perturbed behavior, alternating receding and ascending movements. However, the redshifted velocity component remains constant around 50 km s⁻¹. A detailed analysis is presented in Sect. 4.

Like all spectral lines, the Hα absorption line evolution alternates ascending and receding movements over the whole pulsation cycle. Near φ = 0.892, the first emission arises. This is followed by a brief line doubling episode (φ = 0.943−1.027), during which the extremal blue (B) and red (R) radial velocities of these two components are V (Hα_B) = −63.2 km s⁻¹ and V (Hα_R) = +49.8 km s⁻¹. At φ = 1.037, the Hα line starts a progressive redshift from its maximum blueshifted velocity value − 54.6 km s⁻¹ to finish at + 46.9 km s⁻¹ near φ = 1.874. The Hα line dynamics is studied in detail in Sect. 4.

4 The very redshifted sodium component

Our observations (Figs. 5 and 6) show that a very redshifted component is present within the Na D line profiles near λ5897. It is clearly visible from φ = 0.911 to φ = 0.302, i.e., when it is not blended with the variable pulsational sodium component. Between 1994 and 2017, its velocity remains constant near + 50.3 ± 1.2 km s⁻¹ (see Table 4). The shape of this component also seems constant, i.e., it is not dependent on the pulsation phase. However, the line profiles of other metallic lines such as Fe II λ4923.921 and Mg I λ5183.604 (Chadid & Gillet 1996) as well as Balmer lines (Preston et al. 1965) observed at a resolving power of up to R = 44 000 do not present such a very redshifted component as observed for Na D lines. This points toward a nonphotospheric origin of this very redshifted Na D line component.

The maser emissions SiO or OH are considered as proof of the presence of a circumstellar envelope. Nevertheless, until now this type of emission has never been detected for RR Lyr. In addition, microwave CO lines, which could be the fossil remnants of an ancient shock ejection, have not been detected. Moreover, RR Lyr is not known to have strong mass loss due to spectacular shock ejections like those occurring in unstable red giant envelopes (Tuchman et al. 1978, 1979). As a result, it is unlikely that RR Lyr has a circumstellar envelope. The most likely explanation of this redshifted component is an interstellar medium (ISM) origin. It would be due to the absorption by interstellar dust and gas on the line of sight of RR Lyr.

The sodium splitting is seen in only 75% of the pulsation cycle, but not in the other 25%. The stellar component merges with the static ISM component during φ interval 0.563−0.874 since the Na stellar component is redshifted at φ = 0.338−0.892.

The corresponding large observed heliocentric velocity (V_helio = + 50.3 − 73.5 = −23.2 km s⁻¹, where V_* = −73.5 km s⁻¹) of the redshifted component is consistent with an interstellar origin because, as RR Lyr has a low galactic latitude (b = 12°), the radial velocity of the expected interstellar material could be large due to small projection effect. The observed FWHM of the red component remains approximately constant between 1994 and 2017 around 0.182 ± 0.02Å, i.e., 9.3 ± 0.9 km s⁻¹ (see Table 4). This line width corresponds to the sum of a thermal part and a turbulent part occurring within the ISM gas. It is consistent with the FWHM observed on the high-resolution spectral profiles of the interstellar Na I D line toward 80 early-type stars located in the local interstellar medium (Welsh et al. 1994).

According to the Gaia database, RR Lyr has a parallax of 3.64 ± 0.23 mas, i.e., a distance of 895 light-years (ly; Gaia Collaboration 2016). The Sun is located in a low-density region of the interstellar medium (about 0.05 atoms per cubic centimeter of neutral hydrogen, a near-perfect vacuum) called the Local Bubble (Frisch et al. 2011). This gigantic asymmetric cavity of 330–490 ly in diameter has a wall of denser matter at its surrounding. It is believed that the Local Bubble has been carved by fast stellar winds and supernova explosions within the last 10 Ma. Consequently, this cavity is partially filled with hot gas. RR Lyr is located in another expanding cavity called Loop III (Kun 2016). Between them, at about 200 ly from RR Lyr, Loop III and the Local Bubble are separated by a wall of about 150 ly in size, of colder and denser neutral gas. Because the gas density within the cavities is at least 10 times lower than the average ISM density of the rest of the Milky Way’s spiraldisk, we can expect that the very redshifted sodium component is formed within the Local Bubble wall.

Gillet et al. (2017) have performed a single observation (2017-03-26) of the sodium line at φ = 0.227. The profile shows a line doubling. Because the two components are similar (same FWHM and depth) it was thought that they were formed in a same atmosphere. Consequently, Gillet et al. (2017) suggested that the sodium doubling at φ = 0.227 is due to the presence of a shock within the RR Lyr atmosphere. In the present paper we present a series of 79 spectra sampled approximately every 15 min over the whole pulsation cycle, so we are able to show that the radial velocity of this component, when it is clearly visible, remains stable during the whole pulsation cycle for all line profile parameters. Consequently, we suggest that the redshifted component of the sodium profile has an interstellar origin. As a result, there would be no large velocity (+ 50.3 ± 1.2 km s⁻¹) sodium layer falling on the atmosphere around the pulsation phase φ ~ 0.30, which would enhance the atmospheric compression causing the appearance of the third Hα emission.

Fig. 2 Evolution of the Hα line profile based on RR Lyr observations obtained in April 2017. The pulsation phase varies from φ = 0.618 (bottom left) to φ = 0.083 (top right), while the Blazhko phase ψ overlays the interval 0.77–0.03. The exposure time is 900 s, but the time step is shorter. Weak absorption features are telluric lines. The vertical line indicates the Hα line laboratory wavelength.

Fig. 3 Evolution of the Hα line profile based on RR Lyr observations obtained in April 2017. The pulsation phase varies from φ = 0.093 (bottom left) to φ = 0.600 (top right), while the Blazhko phase ψ overlays the interval 0.77–0.03. The exposure time is 900 s, but the time step is shorter. Weak absorption features are telluric lines. The vertical line indicates the Hα line laboratory wavelength.

Fig. 4 Apparition of the He I λ5876 emission lineat φ = 0.911, i.e., when the Hα line is also in emission. The vertical line indicates the He I line laboratory wavelength.

Fig. 5 Evolution of the D1 sodium line profile (λ5895.92) based on RR Lyr observations obtained in April 2017. The pulsation phase varies from φ = 0.618 (bottom left) to φ = 0.083 (top right), while the Blazhko phase ψ overlays theinterval 0.77–0.03. The exposure time is 900 s, but the time step is shorter. There is a line doubling that stays over 75% of the cycle. The vertical line indicates the sodium line laboratory wavelength.

Fig. 6 Evolution of the D1 sodium line profile (λ5895.92) based on RR Lyr observations obtained in April 2017. The pulsation phase varies from φ = 0.093 (bottom left) to φ = 0.600 (top right), while the Blazhko phase ψ overlays theinterval 0.77–0.03. The exposure time is 900 s, but the time step is shorter. There is a line doubling which remains over 75% of the cycle. The vertical line indicates the sodium line laboratory wavelength.

Fig. 7 Curves of velocity with respect to the stellar rest frame of the sodium and Hα over the pulsation phase φ = 0.8 to φ = 2.2. Shown are the blue components, V (NaB) (blue filled line) and V (HαB) (blue dot-dashed line), and their respective red components, V (NaR) (green dotted line) and V (HαR) (red dashed line). The mean uncertainty is ± 0.6 km s−1.

5 Dynamical structure of the atmosphere

The pulsation of RR Lyr concerns essentially the photosphere and the atmospheric layers, which represent only a small percent of the mass of the star.

Here we consider the main dynamic phenomena occurring during a pulsation period of RR Lyr of average amplitude. It is assumed to take place about halfway between the minimum and the maximum Blazhko, i.e., around the Blazhko phase ψ = 0.75. The pulsation amplitude is minimum at Blazhko minimum and maximum at Blazhko maximum. The amplitude gradually increases between Blazhko minimum and Blazhko maximum. However, it is known that appreciable fluctuations are observed from one cycle to another. The Kepler-amplitude K_p of the light curve of RR Lyr (see Figs. 2 and 3 in Gillet 2013), which is correlated to the intensity of the main shock, strongly depends on the previous pulsation cycle. As a result, this random behavior induces appreciable effects on the atmospheric dynamics, such as changes in observed line profiles. For example, helium lines can be in emission or not. Apart from these fluctuations, observations and models show that the dynamics of the pulsation follows a general reproducible scenario from one cycle to another. In this section we decompose these atmospheric phenomena into six main steps that characterize the different dynamic processes occurring during a pulsation cycle (see Table 5).

5.1 Phase 1 (0.874–0.892): emergence of the main shock

Near the pulsation phase φ = 0.89, a blueshifted Hα emission component suddenly appears within the Hα profile (Fig. 2). When the Hα emission reaches its maximum intensity (φ = 0.911), He I and H II lines can also be in emission if the light curve amplitude, hence the main shock, is large enough. Helium emission lines are not easily observed because their intensity is low (Fig. 4). Consequently, their exhaustive detection requires spectra with a high S/N, i.e., the use of large-diameter telescopes.

As noted in the introduction, RR Lyrae models that consider an atmosphere and shocks are rare. Fokin & Gillet (1997) presented in detail two models of RR Lyr that give the same fine atmospheric structure (see Appendix A). At the end of the global compression of the star near the phase φ = 0.80, the κ − γ mechanism associated with the H-ionization zone provokes a local overpressure, followed by the generation of a wave w1 that quickly breaks into the shock wave s1. At approximately the same phase the second helium ionization zone, which is located far below the photosphere, also generates by a similar mechanism, a compression wave w2 which quickly turns into the shock wave s2.

As shown by Fokin & Gillet (1997), these two waves propagate first toward the interior of subphotospheric layers. In approximately the phase interval 0.75–0.90, the whole atmosphere falls on the star with an infall velocity exceeding that of the waves. Thus, all waves (or shocks) are first receding for the observer (Eulerian coordinates), but their motion is expanding outward in the Lagrangian frame. After phase 0.90, the Mach number of shocks becomes as high as 15–25 and with the deceleration of the atmosphere during its downward motion, the shock starts to propagate outward in the Eulerian frame.

From the sodium and Hα line profiles reported in this paper (see Figs. 3, 2, 5, and 6), it is not possible to distinguish independently the shock waves s1 and s2. Since the latter are formed well below the photosphere (T ~ 40 000 K), the optical depth is certainly too large for the Hα profile to be affected.

5.2 Phase 2 (0.892–0.929): the shock is radiative

For the observations presented in this paper, an Hα emission occurs at the pulsation phase φ = 0.892. The emission is visible during a short time interval from φ = 0.892 to φ = 0.929, i.e., during 5% of the pulsation cycle. The intensity of the shock is thus important enough to produce emission lines: it is a radiative shock.

The Hα emission is formed just behind the shock front within the radiative shock wake. This region is very narrow, typically a few kilometers depending on the shock velocity and physical condition within the atmosphere (Fadeyev & Gillet 2004). At what velocity (or Mach number) does the shock become radiative?

At φ = 0.892, the emission is blueshifted (and defined by V_e1 = −44.4 km s⁻¹, Fig. 2). It is possible to obtain an estimate of the shock front velocity from the velocity of the blueshifted Hα emission. The gas layers emitting the Balmer line radiations are located at the rear of the shock front in the hydrogen recombination zone. We use shock models of Fadeyev & Gillet (2004), which consider the structure of steady-state plane-parallel radiative shock waves propagating through a partially ionized hydrogen gas of pre-shock temperature T₁ = 3 000 K and pre-shock densities between ρ₁ = 10⁻⁹ and 10⁻¹² g cm⁻³. Thus, in the frame of the observer, the gas flow velocity in the hydrogen recombination zone is roughly one-half of the shock wave velocity: (1)

From these models it is also possible to estimate the ratio of V_e1 (denoted δV in Fadeyev & Gillet 2004) of the hydrogen first emission component Doppler shift to the gas flow velocity V_recomb in the hydrogen recombination zone: (2)

Combining these two equations, we obtain (3)

where λ_e1 is the wavelength of the maximum intensity (φ ≃ 0.90) of the Hα emission and λ₀ its laboratory wavelength. Since during the first apparition (pulsation phase φ ~ 0.91) the emission line is blueshifted, then V_shock < 0 and the shock wave is coming toward the observer. However, for practical purpose we use the absolute value of V_shock in this paper. The uncertainties of Eqs. (1)–(3) were estimated from the values given in Table 1 of Fadeyev & Gillet (2004) for the different shock velocities and atmospheric densities. For the considered shock wave model, the validity of Eq. (3) is verified for shock velocities from 40 to 90 km⁻¹, i.e., for Mach numbers ranging from 6.2 to 14.

If we apply this equation to the observed Hα line profile in thepulsation phase interval 0.89−0.93 (Fig. 2), we obtain an estimate of the shock front velocity. We derive the Hα emission position by Gaussian fitting or first moment method depending on the line shape. Velocity uncertainties, which are much larger than measurements uncertainties, are computed with Eq. (3) estimated from the Fadeyev & Gillet (2004) models. Since the emission line is blueshifted for the first apparition, the shock intensity, i.e., the absolute value of V_shock, is considered for the clarity of the analysis. Most of measurements come from the night of April 09, 2017. In order to improve the time resolution, we added an observation from the night of April 12, 2017, which corresponds to φ = 0.925, assuming that shock velocity does not vary too much between consecutive pulsation cycles. This assumption is verified because the measured shock velocity at this phase (+ 111 km s⁻¹) is an intermediate value between + 108 and + 158 km s⁻¹. These values are listed in Table 6.

With an assumed sound velocity of 10 km s⁻¹, the Mach number is between 11 and 13. Thus, with physical conditions present in the atmosphere of RR Lyr, the shock becomes radiative for M ≳ 10.

The average shock front velocity is around 120 km s⁻¹ (see Table 6). This value is much greater than the expected maximum velocity of + 40 km s⁻¹ of the infalling atmosphere (see Fig. 7 near φ = 0.82). Our estimate of the shock front velocity is thus consistent with the hypothesis that the intensity of the shock wave must be intense enough to reverse the supersonic infall motion of the sodium layer.

Our observations of the hydrogen emissions were made on the night of April 8, 2017. The corresponding Blazhko phase is ψ = 0.77, i.e., almost halfway between the minimum and the maximum Blazhko. This must explain why the observed shock velocity (from 133 to 108 km s⁻¹; see Table 6) is relatively moderate in contrast to those which can occur at maximum Blazhko, up to ≃ 170 km s⁻¹ (Gillet et al., in prep.).

The question is whether the Eq. (3) used in this paper can be applied to velocities as high as 110 km s⁻¹ or even higher. When the shock velocity increases, the hydrogen gas behind the shock front is rapidly fully ionized. Moreover, the optical thickness at the Hα central wavelength decreases when the shock velocity increases. Thus, the radiative cooling of the gas increases. In the limit case of very efficient cooling, the isothermal approximation can be made: the cooling is assumed to be instantaneous. Thus, the gas returns to its original pre-shock temperature (T₁) just behind the shock front: T₂ = T₁. This means that the cooling time is very short compared to the dynamic time of the shocked atmosphere.

However, as shown by Fadeyev & Gillet (2004), for photospheric densities and with a shock front velocity of 90 km s⁻¹, the compression rate within the hydrogen recombination zone only reaches a small percent of the isothermal case. Due to the lack of shock models with front velocities up to 200 km s⁻¹, we can only assume that the approximate Eq. (3) remains valid in the radiation flux-dominated regime.

Finally, we find that in the phase interval φ = 0.892−0.929 the shock velocity decreases all the time. Consequently, the observations do not show any acceleration phase of the shock, but only its deceleration phase in the atmosphere. The acceleration phase of the shock necessarily occurs between the moment when the compression waves w1 and w2 are created by the κ − γ mechanism and the beginning of the deceleration phase of the shock.

5.3 Phase 3 (0.320): maximum expansion of the sodium layer

From Table 3, the photospheric sodium component (blue component) reaches a velocity maximum of − 31.3 km s⁻¹ at φ = 0.991 after a continuous velocity increase from φ = 0.911. After φ = 0.911, the velocity continually decreases. This means that, contrary to the Hα layer where the Hα blueshifted absorption is produced, the sodium layer first accelerates before decelerating from φ = 0.911.

At φ = 0.320, the sodium layer stops its ascending motion and starts its slowly infalling motion. The sodium layer thus reaches its maximum expansion in the atmosphere at φ = 0.320, while the Hα layer, where the Hα absorption is formed, is still ascending. The sodium and Hα layers have now an opposite motion and consequently, from this time, an increasing rarefaction zone appears in the atmosphere.

The generation of rarefaction and compression waves during the development of the pulsation of the atmosphere have already been theoretically demonstrated by Hill (1972). Thus, it appears that the motion of the photospheric layers is rather complicated when a layer reverses its movement in relation to another.

5.4 Phase 4 (0.455): maximum expansion of the Hα layer

After the end of the Hα line doubling (φ = 0.027), the Hα profile is in absorption and single. Its radial velocity begins to rapidly decrease from φ = 0.074 (− 52.3 km s⁻¹) to φ = 0.455 where its value reaches ~ 0 km s⁻¹ (Table 3 and Fig. 3). Consequently, at φ = 0.455, the Hα line-forming region stops its expansion. At φ = 0.473, the Hα absorption line becomes redshifted (Fig. 3). Thus, this part of the atmosphere begins to reverse its movement toward the center of the star.

During the ascension end of the sodium layer, which was driven by the shock now located just above it, the radial velocity of the Hα blueshifted component continues to decrease slowly (from −21.1 to −18.8 km s⁻¹) as described in Table 3. This means that the expansion of the Hα layer located just behind the shock front is not yet complete. The Hα profile is now in absorption and single. Its radial velocity continues to decrease to ~ 0 km s⁻¹ at φ = 0.455 (Table 3 and Fig. 3). At φ = 0.473, the Hα absorption line becomes redshifted (Fig. 3). Thus, this part of the atmosphere begins to reverse its movement toward the center of the star.

5.5 Phase 5 (0.320–0.455): a two-step infalling motion

During the phase interval 0.320 ≲ φ ≲ 0.455 the sodium line is redshifted. As a result, the lower atmosphere, where the sodium component is produced, begins its infalling movement well before (φ = 0.320) that of the atmosphere in which the Hα line is formed (φ = 0.473). It is necessary to wait for 15% of the pulsation period for the hydrogen layer to begin to fall on the star.

Our observations show that the atmospheric layers do not follow a synchronized ballistic movement even in the lower part of the atmosphere. It is the propagation of the radiative (hypersonic) shock wave that causes this strong desynchronization of the pulsation movement. If we limit our observational analysis to the movement of the layers of sodium and hydrogen, we have an infalling movement in two successive steps. This means that during the phase interval 0.320 ≲ φ ≲ 0.455, the sodium and Hα layers, which are localized behind the shock, have opposite movements within the atmosphere.

Finally, it is expected that when a shock of large amplitude (hypersonic) occurs within a pulsating atmosphere, a strong desynchronization of the layer motions takes place in the atmosphere. It should be easy to validate this statement using high-resolution spectroscopic observations.

5.6 Phase 6 (~0.36): a strong photospheric compression

From our observations presented in a previous paper (Gillet et al. 2017), the Hα third emission line is observed approximately within the phase interval 0.20 ≲ φ ≲ 0.40. Due to the variation between pulsation cycles, the duration range of visibility of the third emission may undergo some variations on the order of 5–10% in the pulsation phase. In all cases, the presence of a third emission line indicates that a strong warming occurs in the atmosphere due to the intense compression induced by the infalling layers.

Due to the new pulsation cycle, as shown by the model, at φ = 0.20 the whole atmosphere located behind the shock wave is expanding. Only the highest layers of the atmosphere that have not yet been traversed by the shock are still falling back on the star. The ballistic motion of these layers was induced by the previous pulsation cycle.

From our observations we find that the Hα line-forming region stops its expansion atφ = 0.455. Because the third Hα emission disappears near φ = 0.40, this means that the ascending Hα region contributes all the time to the compression, inducing the Hα third emission line. On the other hand, this is not always the case of the sodium layer because it begins its infalling movement at φ = 0.320.

The maximum intensity of the Hα third emission should correspond to the maximum of the compression. For our observations made during the night of September 14, 2014 (Gillet et al. 2017), the maximum compression probably occurs around the pulsation phase φ = 0.357. Thus, the most intense part of the compression occurs after the beginning of the infalling motion of the sodium layer.

The atmospheric phenomena and the shock wave behavior are summarized in Table 7.

6 Shock wave dynamics

6.1 Hα line components and shock velocity

At φ = 0.943, the Hα line doubling phenomenon is already present (Fig. 2). It is observed up to φ = 0.027 when the redshifted component disappears. The velocity of the latter remains virtually constant: +50 km s⁻¹ from φ = 0.855 to φ = 0.027, hence just before the emission, during its presence, after its disappearance, and until the end of the line doubling phenomenon. In addition, the velocity of the blueshifted Hα component decreases all the time from − 58.5 km s⁻¹ (φ = 0.943) to − 52.3 km s⁻¹ (φ = 0.074). This means that the motion of the Hα layer located just behind the shock front undergoes a continuous but modest deceleration (the velocity decreases only by 11% over 13% of the pulsation cycle).

The estimation of the velocity of the shock from the radial velocity of the blueshifted Hα absorption component of the line doubling profile is not obvious. This component is not formed in the radiative wake of the shock wave, i.e., just behind the shock front; on the contrary, it is formed far beyond the shock, within atmospheric layers that undergo complex movements due to the possible presence of rarefaction waves. Under these conditions, only a complete model of the dynamics of the atmospheric layers, taking precisely into account the radiative dissipation of the shock and the dynamical interaction of layers, should allow us to obtain a pertinent estimate of the shock front velocity. Unfortunately, today such a model is not available.

6.2 Shock wave velocity evolution

All the observations presented in this paper were made between April 8 and 18, 2017, as described in Sect. 2. The Blazhko phase ψ was approximately 0.77 at the beginning of the observations and 0.03 at the end, i.e., at a Blazhko maximum. The Hα emission is clearly visible at the pulsation phases φ = 0.892 and φ = 0.911. These two spectra were obtained during the same pulsation cycle, which occurred on the night of April 8, 2017. Using Eq. (3), we obtain a shock velocity of 133–118 km s⁻¹, respectively, for these two phases (see Table 6).

The maximum velocity of the shock is variable from one Blazhko maximum to another. For example, by August 13, 2014, the Hα profile was observed (Gillet et al. 2017) very close to a Blazhko maximum (ψ = 0.04). It shows that the shock velocity decreases rapidly from 156 km s⁻¹ at φ = 0.903–101 km s⁻¹ at φ = 0.941. From Hα profiles observed near ψ≅0.10 by Gillet et al. (2013), Gillet & Fokin (2014) estimated the shock velocity from 132 km s⁻¹ at φ = 0.902 to 114 km s⁻¹ at φ = 0.927.

Finally, these estimates clearly show that the values of the shock velocity as well as the rate of deceleration are highly variable along the Blazhko cycle. Thus, it would be interesting to carry out a high-quality observational study to improve our knowledge on the variation in the shock velocities occurring in RR Lyr during a Blazhko cycle and in different Blazhko cycles.

In Fig. 3, the Hα third emission is clearly visible from phase φ = 0.202 to φ = 0.382. As discussed above, during this phase interval, the shock wave velocity is probably less than 50 km s⁻¹. Below this velocity, the shock is no longer energetic enough to produce an emission component in its wake in the Hα profile: the shock is no longer radiative.

For the 2017 run, the time resolution did not allow us to precisely see how |V_shock| evolves during a pulsation cycle. We thus used a previous set of spectra obtained during August 13, 2014 and October 12, 2013 with exposure times of 300 s to enable such a study. The |V_shock| values are given in Table 8. Hence, from the 2014–2013 data set, it was possible to compute a cubic spline interpolation within a 0.05 phase resolution. All measurements are listed in Table 8. Phased kinetic and radiative flux measurements are plotted with trend curves in Fig. 8 to show the global evolution.

In order to get an idea of the rate of decay of the shock velocity during its elevation in the atmosphere, we have the shock velocity estimate at ψ = 0.04 (August 13, 2014) from φ = 0.903 to φ = 0.941 in addition to those at φ = 1.04 (October 12, 2013, ψ = 0.13). These values are listed in Table 8 and displayed in Fig. 8a. These shock velocities come from two different Blazhko cycles, but they are all close to a Blazhko maximum. They cover the phase interval from 0.90 to 1.04, i.e., 14% of the pulsation period. They correspond to the appearance in the atmosphere of the main shock (φ = 0.90) until it disappears (φ > 1.04) as a radiative shock, i.e., producing the Hα and He I emissions within the shock wake. Although the shock velocities come from two different Blazhko cycles, we clearly identify a decrease in the shock velocity of a factor of about three. Thus, during this strong deceleration phase, the shock propagates over a distance of about two solar radii, i.e., over 30% of the stellar radius R_⋆. When the P Cygni profile of the D3 helium line is observed at phase φ = 1.04, the atmospheric extension of 16% is not very high, but high enough to already observe the characteristic P Cygni shape. Thus, we confirm that the shock is well detached from the photospheric disk at the pulsation phase φ = 1.04.

Uncertainties are computed using Eq. (3). They represent the dispersion of the resulting values of the models calculated by Fadeyev & Gillet (2004) for different shock velocities and unperturbed densities used as parameters in their Table 1.

As discussed in Sect. 5.1, from our observations of the sodium and Hα line profiles, it isnot possible to distinguish independently the shock waves s1 and s2. We only observe the resulting shock s4+s3+s3′′+s1+s2, called the main shock, obtained from the successive merging of s1, s2 and the secondary shocks s4, s3, and s3′. In their model RR41, Fokin & Gillet (1997) show the fusion of the two shocks s1 and s2. Moreover, they show that the amplitude of the shock s4+s3+s3′+s1 decreases rapidly first after merging s1 and again after merging s2 (see their Fig. 6). This weakening of the shock amplitude is consistent with our observations of the decreasing shock velocity when it rises in the atmosphere, as shown in Fig. 8. Thus, in the context of the physical and dynamic conditions taking place in the RR Lyr atmosphere, the shocks do not seem to present an accelerating phase above the photosphere.

6.3 Shock weakening regimes

During its propagation in the atmosphere, the shock is continuously decelerated (Fig. 8). It loses energy in two forms. One form is density dilution (1∕r²) as its radius r increases. The other involves radiative losses as soon as it reaches a hypersonic speed, i.e., when its Mach number is greater than approximately 5. Using the Rankine–Hugoniot equation for energy and assuming that the shock is strong (until the isothermal limit), the maximum value of the radiative flux F_r produced by the shock is (see Appendix) (4)

where ρ₁ is the density of the unperturbed gas in front of the shock and V_shock the shock front velocity. Thus, the flow of radiative losses increases rapidly with the shock front speed (to the third power). In the context of a realistic model, Fadeyev & Gillet (2001) estimated the irreversible loss occurring in the shock wave energy. Comparing the radiation flux emerging from the boundary of the shock layer with the total energy flux, they showed that the radiative losses increase very rapidly with increasing shock front velocity. They showed that 70% or more of the shock energy is lost in radiative form as soon as the shock front velocity is greater than 50 km s⁻¹ for physical conditions occurring in the atmosphere of pulsating stars with unperturbed gas of temperature between T₁ = 3 000 K and T₁ = 6 000 K with a density ρ₁ = 10⁻¹⁰ g cm⁻³.

Consequently, applying the equation above, the shock would lose, in first approximation, up to 18 times more radiation energy at φ = 0.903, 15 times at φ = 0.909, and only 6 times at φ = 0.928 with respect to φ = 1.040 (see Fig. 8b). Accordingly, the shock wave will probably lose most of its energy before the maximum brightness of the star. In 2017, when the maximum velocity of the shock was smaller (133 km s⁻¹ instead of 156 km s⁻¹), the radiative losses only represent 11 times those of the adiabatic regime. For the 2017 data we computed the radiative flux ratio with the 2014 value of F_r(1.04), i.e., V_shock(1.04), assuming that V_shock converges asymptotically to the same value as 2014, as shown in Fig. 8.

In 2014, the deceleration is not uniform between the pulsation phases 0.90 and 1.04 (Fig. 8a). A two-step deceleration phase clearly appears. First, during a very short interval of the pulsation cycle (5% from φ = 0.90 to 0.95), the shock velocity decreases by one-third of its value. This rate of fast decrease in the shock velocity occurs when the shock has a very large radiative loss, close to that of an isothermal shock. Then, up to φ = 1.04, but slower, the shock velocity decreases byabout a factor of 40% over 10% of the period. During this second phase, the radiative losses of the shock become secondary with respect to those induced by the geometric dilution due to the shock propagation. These two phases are clearly visible in Fig. 8.

In 2017, the maximum shock velocity did not reach a value (133 km s⁻¹) as large as in 2014 (156 km s⁻¹), i.e., a Mach number of 13 instead of 16. This difference seems sufficient to induce the transition between the two hydrodynamic regimes: theisothermal one, in which the radiative losses dominate, and the adiabatic one, in which the losses due to the geometric dilution are dominant.

Finally, when the intensity of the main shock wave is strong enough (isothermal phase), most of the shock wave energy is dissipated by radiative processes at the photospheric level, i.e., in the deep atmosphere. The remainder of the shock wave energy is dissipatedby dilution in the high atmosphere when the radiative nature of the shock is no longer dominant (adiabatic phase).

Fig. 8 Phased kinetic and radiative flux measurements with trend curve. Measurements points are the most probable values of |Vshock| and Fr (φ)∕Fr(1.04). Associated uncertainties are computed in Tables 8 and 6. Red dots are measured values from August 13, 2014 while the red triangle represents the night of October 12, 2013. Velocities from April 2017 are black squares with the corresponding trend curve (grey dotted line). Panel a: evolution of |Vshock| plotted with a cubic spline interpolation curve. In the phase interval φ = 0.90−0.92, Vshock decreases much more than in the phase interval φ = 0.92−1.04. Thus, for 40% of the period, the shock front velocity decreases by a factor of more than three. Panel b: evolution of maximum values of Fr(φ)∕Fr(1.04) ratio plotted with a cubic spline interpolation curve. In the phase interval φ = 0.90−0.92, the Fr ratio decreases much more than in the phase interval φ = 0.92−1.04. The flow of radiative losses increases rapidly with the speed of the shock front.

6.4 Turbulent state of the atmosphere

Gillet et al. (1998) discussed the turbulence amplification in the atmosphere of a radially pulsating star, showing that it is due to the global compression of the atmosphere during the pulsation, and that it is also affected by strong shock waves propagating from the deep atmosphere. As already noticed by Kolenberg et al. (2010) and Fossati et al. (2014), among others, by definition the microturbulence velocity V_mic characterizes a difference between observed and theoretical atmospheric parameters affecting the FWHM of an absorption spectral line: (5)

Therefore, it represents only the upper limit of atmospheric turbulence. Indeed, all line broadening effects such as velocity gradients or non-LTE effects are not always taken into account in atmosphere models.

A time variation in the microturbulence velocity V_mic in the atmosphere of RR Lyrae has been estimated from a comparison of the observed and theoretical full width at half maximum (FWHM) curves of the Fe II 4923.921Å absorption line by Fokin et al. (1999). This velocity V_mic reflects motions on scales smaller than the line-forming region. Fokin et al. (1999) found that turbulence in this region varies from 2 to 7 km s⁻¹, over a large phase interval (0.5 ≲ φ ≲ 1.0). The maximum (7 km s⁻¹) occurs at the pulsation phase φ = 0.7 (during the light bump) and a secondary peak (3.5 km s⁻¹) occurs at φ = 0.2. This prominent peak of V_mic starts to appear when the two secondary shocks s4 and s3, as noted by Fokin & Gillet (1997), merge at φ = 0.52 when V_mic is at its lowest level (just after the phase of the maximum radius). The physical origin of the shock s4 is due to the accumulation of several weak compression waves at the sonic point during the beginning of the atmospheric compression, while s3 comes from the stop of the hydrogen recombination front near the phase of maximum expansion. In the case of RR Lyr, the peak maximum of the microturbulence velocity at 7 km s⁻¹ takes place when s4+s3 merge with the shock s3′ assumed to be caused by the 1H-mode (Fokin & Gillet 1997). At this stage, the overall velocity of the shock s4+s3+s3′ is at a maximum, but still modest (55 km s⁻¹), although supersonic (Mach number M ≳ 5). Then the velocity of this shock slowlydecreases almost to 40 km s⁻¹ (Fokin & Gillet 1997, their Fig. 6). As a consequence the microturbulence velocity decreases until its rapid and abrupt rise induced this time by the arrival of the strong shocks s2 and s1 near φ = 0.91, which have their origin in the κ-γ mechanism associated with the H- and He-ionization zones. However, the effect of these two last shocks on the microturbulence velocity remains modest because they merge with s4+s3+s3′ in the upper part of the line-forming region of the Fe II line.

Recently, as part of an in-depth spectroscopic analysis of the Blazhko star RR Lyrae, Kolenberg et al. (2010) deduced a microturbulence velocity curve from the FWHM of the Fe II line at 4508.288Å. In a complementary work, Fossati et al. (2014) showed this curve (green solid line in their Fig. 6); it presents an asymmetrical peak starting abruptly at φ = 0.9, with a maximum near φ = 0.95, and then a slow decrease until φ = 0.2. It is centered just before the middle of the rising branch of the V-light curve. The minimum microturbulent velocity occurs at the maximum radius during the most stable pulsation phase, called the quiet phase, near φ = 0.33. The value of the microturbulence velocity obtained by these authors is much higher (between 15 and 35 km s⁻¹) compared to the values between 2 and 7 km s⁻¹ found by Fokin et al. (1999). This may originate, at least in part, from the use of a static atmosphere model that is not able to take into account the strong velocity gradients present in the star. Indeed, such a velocity gradient is an important source of line broadening within a pulsating atmosphere (see Fokin et al. 1996). From the range of line formation depths covered by the measured lines, Fossati et al. (2014) estimate the V_mic profiles in the corresponding optical depth range logτ_Ross = − 4 to logτ_Ross = 0, i.e., within the middle part of the atmosphere. The most striking fact is that the turbulent velocity continually decreases while the depth in the atmosphere increases. This is true for all pulsation phases – especially during the quiet phase and during the maximum compression induced by the global atmospheric pulsation – suggesting that the large-scale motions have no effect on the microturbulent velocity, a disturbing result. Between the bottom and the top of the line-forming region, the microturbulence velocity varies from 3 to 20 km s⁻¹ except during the peak where the microturbulent velocity may be as high as 50 km s⁻¹, i.e., largely supersonic. In addition, Fossati et al. (2014) also point out that the microturbulence velocity of the deepest atmospheric layer increases, while that of the layer just above remains unchanged, also a surprising behavior.

As shown by Gillet et al. (1998), the turbulence amplification during a compression of the atmosphere without the presence of waves results from the variation of the density within the line-forming region. Because in radially pulsating stars, such as RR Lyrae stars, the atmospheric curvature is very small, it is possible to assume that the atmospheric distortion is not locally spherical but almost radial. Thus, compression must induce a preferential amplification direction. Gillet et al. (1998) estimated that the global atmospheric compression can be almost considered as an adiabatic homogeneous axial compression. Consequently, in the framework of a rapid distortion, the solenodial RDT model of Jacquin et al. (1993) appears to give the better quantitative prediction. As an application, for the radially pulsating star δ Cephei, the predicted turbulence amplification induced by the global atmospheric compression is consistent with the solenodial RDT (Gillet et al. 1998).

In the case of the presence of strong shocks as in the case of RR Lyrae stars, i.e., when the compression rate becomes higher than 2 (Mach number ≳2), radiative effects take place and adiabatic turbulence amplification theories break down. In this case, these theories give a considerable overestimation of the amplification (Gillet et al. 1998). Consequently, in the limiting case of very strong shocks (isothermal shocks), the compression ratio η caused by the shock is η ≃ M². Basically, effects induced by radiative terms in conservation equations are certainly at the origin of the observed limit of the turbulence amplification. Nevertheless, the relevant theory does not yet exist and a new theoretical approach is required to take into accountthe radiative field within hypersonic shocks.

In this paper, we have established that the movements in the atmosphere are complex because the motion of atmospheric layers are not synchronous and are disturbed by the passage of strong shock waves. It is not possible to observe the steps of formation and acceleration of the main shock. Only its long deceleration phase is observable. During slightly less than 15% of the pulsation cycle (0.90 ≤ φ ≤ 1.04), the main shock velocity decreases from 150 to 60 km s⁻¹, i.e., from a Mach number of 15 to 6. This strong deceleration occurs in only 2 h, a rather short duration compared to the 13.6 h period. The radiative flux emitted by the shock then drops by a factor of 18 (Fig. 8b), demonstrating that its energy is essentially dissipated in radiative form. All atmospheric layers do not reach their maximum expansion at the same time (a 2-h difference between sodium and Hα as shown inFig. 7), generating rarefaction waves during the pulsation from φ = 0.30 to 0.50 (see Table 5). Then a violent atmospheric compression occurs in the lower part of the atmosphere at φ = 0.36 due to the rapid fall of the upper layers of the atmosphere.

Therefore, in the context of the extremely complex atmospheric pulsation of RR Lyr, it is not obvious to validate the use of successive static atmosphere models, assuming LTE and plane-parallel geometry, to interpret the different phenomena involved in this highly dynamic atmosphere. As a result, it seems difficult to get a good determination of the effective temperature through synthetic spectra fitting of the observed Hγ line, and of the surface gravity, which is obtained by imposing the ionization equilibrium for iron. Moreover, while imposing the equilibrium between the iron line abundance and equivalent widths to determine the microturbulent velocity, Fossati et al. (2014) introduced the need of a depth-dependent V_mic. Thus, in the context of a very disturbed atmosphere revealed by the observations reported in this article, it would be interesting to confirm the result of Fossati et al. (2014) regarding the depth-dependent of V_mic.

For δ Cephei, the maximum compression rate within the Fe II line-forming region induced by the global compression of the atmosphere near φ = 0.85 is of the same order of magnitude as that caused by the main shock (Gillet et al. 1998). Because the maximum amplitude of the latter in RR Lyr is about 5 times larger than that of δ Cephei, we must expect that the dominant amplification phenomenon could be the main shock and not the global atmospheric compression. Moreover, as shown at the beginning of this section, the intensity of the main shock of RR Lyr is extremely large (M ≳ 15) only in the lowest part of the atmosphere because the shock velocity decreases very rapidly when it rises in the atmosphere, due to the intense radiative loss occurring in the shock wake. Consequently, the greatest microturbulence velocities should be observed essentially in the lower part of the atmosphere. Furthermore, as shown by Gillet et al. (1998), the presence of several secondary shocks in the atmosphere, as expected by pulsation models (Fokin & Gillet 1997), can also be strong enough in the Fe II line-forming region to contribute to an amplification of the microturbulence velocity on the order of that induced by the global compression of the atmosphere. This is consistent with the result obtained by Fokin et al. (1999) because the wide maximum of microturbulence velocity spanning half a pulsation period occurs when secondary shocks s3 and s3’ are crossing the line-forming region. Also, the third Hα emission observed near the pulsation phase 0.3 (Gillet et al. 2017) is interpreted as the consequence of a large atmospheric compression which provokes a strong and rapid increase of density and thus an appreciable amplification of the microturbulence velocity. This atmospheric compression could explain the secondary peak (3.5 km s⁻¹) of microturbulence velocity over the phase interval 0.10 ≤ φ ≤ 0.35.

Finally, in order to confirm that the microturbulence velocity increases with height in the atmosphere, as found by Fossati et al. (2014), it would be necessary to perform for RR Lyr a computation using a convective pulsation model including an extended atmosphere and allowing the generation of shock waves of large amplitude.

7 Conclusion

The objective of this paper was to determine from spectral observations the main dynamic phenomena occurring in the atmosphere of RR Lyr during a typical pulsation cycle. As the intensity of the pulsation cycles varies greatly over time, in particular because of the Blazhko effect, we considered a mean amplitude cycle, i.e., observations with a Blazhko phase close to 0.75. The atmosphericphenomena that can occur during the cycles of extreme amplitude will be examined in a second paper.

This observational study was based on a series of 79 spectra sampled approximately every 15 min on the pulsation cycle. Thus, we are close to a temporal resolution of two hundredths (period ~13 h 36 min). This resolution seems quite sufficient to follow all the developments occurring in spectral line profiles such as Hα and Na D. Moreover, the resolving power used is relatively modest (R = 11 000), all the potential spectral information that could exist in the profiles seems to be highlighted. All the spectra of this series were obtained between April 7 and 18, 2017; therefore, during the same Blazhko cycle. Observations were started at the Blazhko phase ψ = 0.75 and were stopped at maximum Blazhko (ψ = 0.03). Since the most intense hydrogen emission occurs at the beginning of the observation series, the intensity of the main shock corresponds to the average pulsation cycle (ψ ~ 0.77).

It should be noted that these observations were made with a modest telescope (35 cm). As RR Lyr is the brightest RR Lyrae star (7 < V < 8), an effective exposure time of 900 s allows a S/N of around 100 per pixel depending on weather conditions. As the feasibility conditions of this spectroscopic study were satisfying for this star, its extension to other RR Lyrae stars requires at least a one-meter class telescope.

The observed sodium and Hα profiles presented in this paper did not allow us to resolve the shocks s1 and s2 individually. Only the main shock, the sum of these two shocks induced by the κ-mechanisms in the He and H-ionization subphotospheric layers, is confirmed. It is possible that observations with higher time resolution (exposure time less than 5 min), high S∕N ≳ 100 and high spectral resolution (R ≳ 50 000) may allow the detection of these two shocks separately. Indeed, the observation of the shock s2, which is initiated deep below the photosphere, is probably impossible due to the excessive opacity of the gas. To observe s2, it would have to have existed above the photosphere before it merged with s1, but this does not seem to be the case from our observations.

One of the most striking results of this observational study is that we do not observe the acceleration phase of the main shock as soon as the Hα emission produced in the shock wake becomes visible. Consequently, we directly observe the deceleration of the wave. The deceleration is faster when the shock is radiative than after when dilution phenomena become dominant. In the deep atmosphere (for 0.903 < φ < 0.941), the shock is radiative when its Mach number is larger than 10, whereas when in the high atmosphere (φ = 1.040) a Mach number around 6 is enough.

The sodium layer reaches its maximum expansion at φ = 0.320, while for the layer in which the Hα line core is formed, it is much later, at φ = 0.455. Thus, a rarefaction wave appears in the atmosphere. This phenomenon is clearly predicted by the models of Fokin & Gillet (1997).

The so-called third Hα emission is observed around φ ~ 0.36. It occurs when the ascending Hα layer and the highest infalling atmospheric layers of the previous pulsating cycle compress the gas located between them. At this pulsation phase (φ ~ 0.36) the sodium layer, located in the deep atmosphere, already has an infalling movement, which means that the sodium does not participate in the compression inducing the third Hα emission.

As previously shown by Gillet et al. (1998), the turbulence amplification of the atmospheric gas is mainly due to both the global compression of the atmosphere during its pulsation and to strong shock waves propagating through the atmosphere. In the case of RR Lyr, because amplitudes of shocks occurring within the atmosphere are large, the turbulence amplification is mostly caused by shocks, while that provoked by the global compression of the atmosphere seems secondary. This is the reason why the main peak of microturbulence velocity is very wide for RR Lyr conversely to the case of δ Cephei where it is quite narrow. Furthermore, because the main shock is always observed with a rapidly decreasing velocity during its propagation, the consecutive induced turbulence amplification should be highest in the lower part of the line-forming region of the Fe II line. It would be interesting to establish whether the microturbulence velocity increases or decreases with height in the atmosphere;this might be accomplished with a pulsation model including convection and an extended atmosphere and allowing the generation of shock waves of large amplitude.

In either case, in the future it would be constructive to know the variations in the shock intensity present in the atmosphere of RR Lyr, as well as those occurring in the various line profiles induced by atmospheric dynamics. However, this study requires a large observational survey because it is essential to integrate the effects due to the Blazhko phenomenon on atmospheric dynamics. This will be the goal of the second article in this series.

Finally, the data and work presented in this paper demonstrate further the increasing role of the amateur spectroscopy community in stellar surveys.

Acknowledgements

We thank Lux Stellarum and the French OHP-CNRS/ PYTHEAS for their support. We also thank D. Boureille for his help on smart small tools. The present study has used the SIMBAD database operated at the Centre de Données Astronomiques (Strasbourg, France) and the GEOS RR Lyr database hosted by IRAP (OMP-UPS, Toulouse, France), created by J.F. Le Borgne. Spectroscopic data have been analyzed with the SpcAudace software written by B.M. (ARAS group, France). We especially thank Christian Féghali for his very careful reading of the manuscript and for his pertinent remarks. We appreciate the helpful comments from the reviewer and editor. We gratefully acknowledge Helenka Kinnan for her very careful reading of the final version of this paper.

Appendix A The shock model of Fokin & Gillet (1997)

Fokin & Gillet (1997) computed more than 20 nonlinear, nonadiabatic, purely radiative pulsation models of the RR Lyrae prototype, RR Lyr. However, RR Lyr is located very close to the red edge of the instability strip, convection probably slightly affects the ionization zones of hydrogen and helium at the origin of the pulsation. These models go through the subphotospheric layers up to an extended atmosphere.

The shock waves are implicitly calculated with the Von Neumann-Richtmyer artificial viscosity that mimics the viscous dissipation of the kinetic energy of the shock. The consequence of this artificial viscosity is to widen the shock wake. Thus, it is not possible with this kind of numerical approach to resolve the fine structure of the shock because the spatial resolution is too weak (about 10³ km), or even to resolve the largest zone of the wake where photorecombinations occur (typically over a few km). As part of their numerical approach, Fokin & Gillet (1997) note that the differentiation between a wave and a shock wave is not always obvious through the numerical model. They used the word “shock” when the density jump exceeds 5 over four consecutive mass zones.

The effect of excessive damping of the pulsation and waves due to the artificial viscosity have been tested by Fokin & Gillet (1997). It appears that their models are not very sensitive to the damping and cutoff parameters used. Consequently, the characteristics of the shocks and their associated hydrodynamics phenomena are stable and are present from one cycle to another.

During each pulsation cycle, four main shocks, with a well-identified physical origin, are produced in the outer part of the stellar envelope. In the case of RR Lyr, a fifth shock (shock s3’) is also highlighted; however, its physical origin has not yet been clearly identified. Contrary to the static regime, the extension of the atmosphere increases by several times during the pulsation. It is interesting to note that the four main shocks are stronger than in classical Cepheids, but their general features, especially their physical origin, are qualitatively similar.

The two strongest shocks are the consequence of the recombination of the ionization zones of hydrogen and helium in the subphotospheric region. Fokin & Gillet (1997) call theses shocks s1 and s2, respectively.

The scheme of phenomena predicted by Fokin & Gillet (1997) model are summarized in Table A.1. Their models (including RR41) limit the description of the atmosphere to an external density to 10⁻¹⁴ g cm⁻³. Thus, the maximum atmospheric height above the photosphere only reaches 0.22R_ph near the phase of the maximum expansion (φ = 0.55).

Appendix B Maximum radiative flux produced by the shock

This section describes which hypotheses are applied on the Rankine–Hugoniot equation system to obtain an upper limit on the radiation flux relation (Eq. (4)) used in this paper.

Considering the total energy Rankine–Hugoniot relation, the radiation occurs through three different terms: the radiative pressure p_r, the specific radiative energy e_r, and the radiative flux F_r: (B.1)

In general, the exact value of the three radiative terms is difficult to determine, especially when the gas opacity is strongly frequency dependent. The first three angular moments of the specific intensity I_μν of the radiation field are (Mihalas & Mihalas 1984, Sects. 6.4–6.6) (B.2) (B.3)

and (B.4)

where μ ≡ cosθ is the cosine of the polar angle θ between the direction of the considered radiation beam and the normal to the emitting surface. Nevertheless, when the fluid is close to the thermodynamical equilibrium, i.e., when it is optically thick, then the evaluation of the angular moments I_μν is straightforward. In this case (see, e.g., Sampson 1965, Sect. 5.1 or Mihalas & Mihalas 1984, Sects. 6.8 and 6.9) the specific radiative energy is (B.5)

and (B.6)

where a_r ≡ 4σ_r∕c = 7.56 10⁻¹⁵ erg cm⁻³ K⁻⁴ is the radiation constant and σ_r = 5.67 10⁻⁵ erg s⁻¹ cm⁻² K⁻⁴ is the Stefan–Boltzmann constant. At the strictly thermodynamic equilibrium limit, the radiation field is isotropic. Consequently, since F_r represents the flux, i.e., the sum of photon contributions in all directions, this means that the F_r value is zero. Nevertheless, in the shock wake region located in front of the region where recombinations take place, the radiative flux cannot be in equilibrium because the gas is partially optically thin. This region includes the precursor, the shock front, and the first part of the wake where the thermalization occurs between heavy particles and electrons. Ina first approximation, we can consider the wake recombination region as a blackbody at temperature T and producing an equilibrium radiative flux, (B.7)

in the shock front direction. Because all atoms do not recombine at the same temperature and electronic density, T must be considered as an effective temperature integrating large gradients of the physical variables present in the recombination region. Thus, from these thermodynamical equilibrium relations, we have (B.8)

It should be noted that if the radiation field is larger in one direction (that of the outward flow), then p_r will exceed its isotropic limit of 1∕3e_r (see, e.g., Mihalas & Mihalas 1984, Sect. 66). In the extreme limit case occurring in plane waves, and more approximately in the outermost layers of the very extended atmosphere, we have p_r = e_r (see Mihalas & Mihalas 1984, Sect. 66). Thus, in any case, p_r remains close to e_r. Finally, at thermodynamical equilibrium, Eq. (B.8) shows that the radiative flux is always the main radiative contributionto the energy Rankine–Hugoniot relation.

This is explained by the fact that the speed of light is usually much higher than the relative post-shock velocity u ≡ V_shock −v. This is almost always the case in stellar atmospheres, where u rarely exceeds 300 km s⁻¹ (0.1% of c). We note that for extremely fast shocks, such as supernovae explosions, where the initial shock front velocity is around 5% of the speed of light, the radiative pressure and energy are nearly equal to the gas pressure and energy. Thus, it is only for very fast flows of this kind that the three radiative terms must be considered.

Consequently, hereafter the only radiative term that we consider in the total energy Rankine–Hugoniot relation is the radiative flux F_r: (B.9)

The gas energy e_g is the sum of the energy e_t involved into the translational degrees of freedom of atoms and the energy e_i stocked into their internal degrees of freedom (excitation, dissociation, ionization): (B.10)

In an ideal gas, the energy e_t in translational degrees of freedom is given by (B.11)

where γ is the ratio of specific heat.

If E, D, and I respectively denote the excitation, dissociation, and ionization specific energies stored in the degrees of freedom of atoms, the internal gas energy e_i is (B.12)

If indexes 1 and 2 denote the gas just before and just after the shock front, and if we neglect the gas energy in the unperturbed gas, then we have (B.13)

If we consider the medium 2, the shock cools very quickly, and then we have T₂ = T₁. In this case (isothermal shock), all the energy stored in the internal degrees of freedom of atoms is totally dissipated after the shock front. In addition, if the ballistic movement of the atmosphere has already been completed before the arrival of the next shock, then the velocity of the gas falling down on the star at the level of the shock is zero (u₁ = V_shock). In these conditions, the above Rankine–Hugogniot relation reduces to (B.14)

This means that the shock cannot radiate more radiative energy flux than (B.15)

References

All Tables

All Figures

Fig. 1 Comparison of the three apparitions of Hα line profile of RR Lyr. Panel a: second apparition at φ = 0.74 represented with a spectrum dated on April 10, 2017. Panel b: first apparition at φ = 0.91 represented with a spectrum dated on April 09, 2017. Panel c: third apparition at φ = 1.25 represented with a spectrum dated on April 08, 2017. Vertical line indicates the Hα line laboratory wavelength.

In the text

	Fig. 2 Evolution of the Hα line profile based on RR Lyr observations obtained in April 2017. The pulsation phase varies from φ = 0.618 (bottom left) to φ = 0.083 (top right), while the Blazhko phase ψ overlays the interval 0.77–0.03. The exposure time is 900 s, but the time step is shorter. Weak absorption features are telluric lines. The vertical line indicates the Hα line laboratory wavelength.
In the text

	Fig. 3 Evolution of the Hα line profile based on RR Lyr observations obtained in April 2017. The pulsation phase varies from φ = 0.093 (bottom left) to φ = 0.600 (top right), while the Blazhko phase ψ overlays the interval 0.77–0.03. The exposure time is 900 s, but the time step is shorter. Weak absorption features are telluric lines. The vertical line indicates the Hα line laboratory wavelength.
In the text

	Fig. 4 Apparition of the He I λ5876 emission lineat φ = 0.911, i.e., when the Hα line is also in emission. The vertical line indicates the He I line laboratory wavelength.
In the text

Fig. 5 Evolution of the D1 sodium line profile (λ5895.92) based on RR Lyr observations obtained in April 2017. The pulsation phase varies from φ = 0.618 (bottom left) to φ = 0.083 (top right), while the Blazhko phase ψ overlays theinterval 0.77–0.03. The exposure time is 900 s, but the time step is shorter. There is a line doubling that stays over 75% of the cycle. The vertical line indicates the sodium line laboratory wavelength.

In the text

Fig. 6 Evolution of the D1 sodium line profile (λ5895.92) based on RR Lyr observations obtained in April 2017. The pulsation phase varies from φ = 0.093 (bottom left) to φ = 0.600 (top right), while the Blazhko phase ψ overlays theinterval 0.77–0.03. The exposure time is 900 s, but the time step is shorter. There is a line doubling which remains over 75% of the cycle. The vertical line indicates the sodium line laboratory wavelength.

In the text

	Fig. 7 Curves of velocity with respect to the stellar rest frame of the sodium and Hα over the pulsation phase φ = 0.8 to φ = 2.2. Shown are the blue components, V (NaB) (blue filled line) and V (HαB) (blue dot-dashed line), and their respective red components, V (NaR) (green dotted line) and V (HαR) (red dashed line). The mean uncertainty is ± 0.6 km s−1.
In the text

Fig. 8 Phased kinetic and radiative flux measurements with trend curve. Measurements points are the most probable values of |Vshock| and Fr (φ)∕Fr(1.04). Associated uncertainties are computed in Tables 8 and 6. Red dots are measured values from August 13, 2014 while the red triangle represents the night of October 12, 2013. Velocities from April 2017 are black squares with the corresponding trend curve (grey dotted line). Panel a: evolution of |Vshock| plotted with a cubic spline interpolation curve. In the phase interval φ = 0.90−0.92, Vshock decreases much more than in the phase interval φ = 0.92−1.04. Thus, for 40% of the period, the shock front velocity decreases by a factor of more than three. Panel b: evolution of maximum values of Fr(φ)∕Fr(1.04) ratio plotted with a cubic spline interpolation curve. In the phase interval φ = 0.90−0.92, the Fr ratio decreases much more than in the phase interval φ = 0.92−1.04. The flow of radiative losses increases rapidly with the speed of the shock front.

In the text