Characterization of the radial velocity signal induced by rotation in late-type dwarfs

Abstract

1 INTRODUCTION

High precision radial velocity (RV) measurements give astronomers the possibility of detecting small exoplanets, down to the mass of the Earth. Unfortunately, intrinsic variations of the magnetic regions on the stellar surface induce RV variations that, if stable over a few rotation periods, can mimic a planetary signal. Recognizing and characterizing them is key to disentangle the true planet induced signals (Queloz et al. 2001; Dumusque et al. 2012; Robertson et al. 2014; Santos et al. 2014).

Activity-induced signals depend mainly on the activity level and spectral type of a star. Highly active stars usually rotate faster, inducing shorter period signals with larger amplitudes. These amplitudes also depend on the strength of the magnetic field, which is a function of the spectral type as it depends strongly on the depth of the convective zone. These signals are produced by a variety of physical phenomena in the stellar surface. Among them, spots are particularly relevant for rapidly rotating late-type dwarfs (Noyes et al. 1984; Saar & Donahue 1997; Santos et al. 2000; Hatzes 2002; Desort et al. 2007). For slowly rotating stars, plages become another important, but not so well understood, source of variability in the time-scales of rotation. RV signals induced by plages and spots are not expected to be well correlated (Meunier, Desort & Lagrange 2010).

The amplitude of the so-called RV jitter has been a subject of study for many years (Saar, Butler & Marcy 1998; Santos et al. 2000; Paulson et al. 2002; Wright 2005; Isaacson & Fischer 2010; Martínez-Arnáiz et al. 2010). Taking advantage of the precision achieved by HARPS (better than 1 m s⁻¹) and the sampling rate of planet-hunting surveys, we present here a detailed analysis of this jitter measuring its period and amplitude in RV time series data of a sample of late-type dwarfs. We study these rotationally induced periodic signals, with particular attention to the harmonics of the stellar rotation period (see also Boisse et al. 2011), and to the relationships between the measured amplitudes and the level of chromospheric emission [|$\log _{10}(R^{\prime }_\textrm{HK})$|]. We also investigate the relation between the amplitude of the induced RV signals and the amplitude of the modulation of the chromospheric emission measured as the Mount Wilson S Index (Noyes et al. 1984).

2 STELLAR SAMPLE AND DATA

We selected a sample of bright southern stars with HARPS spectra available in the European Southern Observatory (ESO) public data base. Our initial sample consists of 55 low activity [|$\log _{10}(R^{\prime }_\textrm{HK})\ <-4.4$|] main-sequence stars covering from late-F to mid-M type. The selected stars are part of the planet-hunting programme using HARPS (Mayor et al. 2003) and their rotation periods and magnetic cycles have been investigated using spectroscopic and photometric time series (Suárez Mascareño et al. 2015; Suárez Mascareño, Rebolo & González Hernández 2016). All of them had been observed more than 20 individual nights as of 2016 May, giving a total number of more than 9000 spectra. Fig. 1 shows the distribution of spectral types for the sample of stars. Table 1 shows relevant data for stars where a reliable detection of an activity-induced RV signal was obtained in this work.

Figure 1.

Distribution of spectral types in our sample. M-dwarfs have been separated into two groups, with the fully convective ones (later than M3.5) in a separate category. Dark fill shows the number of detections of periodic RV signals induced by rotation.

Open in new tabDownload slide

2.1 Spectroscopic data

HARPS is a fibre-fed high-resolution echelle spectrograph installed at the 3.6 m ESO telescope in La Silla Observatory (Chile). The instrument has a resolving power R ∼ 115 000 over a spectral range from 378 to 691 nm. It has been designed to attain extreme long-term RV accuracy. It is contained inside a vacuum vessel to avoid spectral drifts due to changes in temperature and air pressure. HARPS comes with its own pipeline providing extracted and calibrated spectra, as well as RV measurements and other data products such as cross-correlation functions (CCFs) and their bisector profiles.

For the analysis of the spectral indicators, we use the extracted order-by-order wavelength-calibrated spectra produced by the HARPS pipeline. In order to minimize atmospheric effects, we create a spectral template for each star to correct the order-by-order flux ratios for the individual spectra. We correct each spectrum for the Earth's barycentric RV and the RV of the star using the measurements given by the pipeline. We finally re-bin them into a constant wavelength step.

2.2 Photometry

We also use the light curves provided by the All Sky Automated Survey (ASAS) public data base. ASAS (Pojmanski 1997) is an all sky survey in the V and I bands running since 1998 at Las Campanas Observatory, Chile. Best photometric results are achieved for stars with V ∼ 8–14, but this range can be extended implementing some quality control on the data. ASAS has produced light curves for around 10⁷ stars at δ < 28°.

The ASAS catalogue supplies ready-to-use light curves with flags indicating the quality of the data. For this analysis, we relied only on good-quality data (grade ‘A’ and ‘B’ in the internal flags). Even after this quality control, there are still some high dispersion measurements that cannot be explained by a regular stellar behaviour. We reject those measurements by de-trending the series and eliminating points deviating more than three times the standard deviation from the median value.

3 STELLAR ACTIVITY INDICATORS AND RVs

We compute the Ca II H&K index, the |$\log _{10}R^{\prime }_{\text{HK}}$| and the H α index following Noyes et al. (1984), Lovis et al. (2011) and Suárez Mascareño et al. (2015, 2016). We use the indices to re-measure the rotation period of the selected stars and to be able to compare the measured RV semi-amplitudes with the activity indicators.

3.1 Radial velocities

RVs are taken directly from the measurements of the standard HARPS pipeline, except for M-type stars. The RV measurements in the HARPS standard pipeline is determined by a Gaussian fit of the CCF of the spectrum with a digital mask (Pepe et al. 2000). In the case of M-dwarfs, due to the huge number of line blends, the CCF is not Gaussian resulting in a less precise RV measurement and in a loss of sensitivity to the changes in the full width at half-maximum (FWHM) of the CCF. For this type of stars, we opted to do a different modelling of the CCF, using the combination of a second-order polynomial with a Gaussian function over a 15 km s⁻¹ window centred at the minimum of the CCF following Suárez Mascareño et al. (2017). The centre of the Gaussian function is taken as our RV measurement. This allows us to improve the stability of the measured RV as well as to improve the sensitivity of the FWHM measurements. For high signal-to-noise measurements, the extracted value is virtually the same.

3.2 Quality control of the data

As the sampling rate of our data is not well suited for modelling fast events, such as flares, and their effect in the RV is not well understood, we identify and reject points likely affected by flares by searching for abnormally high measurements of the activity indicators and clear distortions in the Balmer series (Reiners 2009).

4 ANALYSIS

4.1 RV signals induced by stellar rotation

In order to identify the rotation-induced RV signals, we first analyse the time series of the Ca II H&K and H α activity indicators, the time series of the variations of the FWHM of the CCF, and the light curve (when available) to determine the rotation period of the star. Some of the rotation periods were previously measured in Suárez Mascareño et al. (2015, 2016). In these cases, we measure it again including the last available data. To search for periodic signals in the series, we compute the power spectrum using a generalized Lomb–Scargle periodogram (Zechmeister & Kürster 2009) and if there is any significant periodicity we fit the detected signal using the rvlin package (Wright & Howard 2012).

To evaluate the false alarm probability of any peak in the periodogram, we follow Cumming (2004) modification over the work by Horne & Baliunas (1986) to obtain the spectral density thresholds for a desired false alarm level. This means our false alarm probability is defined as FAP = 1 − [1 − P(z > z₀)]^M, where P(z > z₀) = exp(−z₀) is the probability of z being greater than z₀, with z the target spectral density, z₀ the measured spectral density and M the number of independent frequencies. We search for the power values corresponding to 10, 1 and 0.1 per cent false alarm probability.

The mean rotation period of the star is estimated by calculating the weighted average of the detections obtained with different indicators. Errors are the peak-to-peak variation divided by the square root of the number of detections, and the global false alarm probability is the combined false alarm probability given by a Fisher's combined probability test (Fisher 1925).

Different activity proxies provide information about different parts of the stellar atmosphere. Photometry gives information mainly about photospheric variability, while the S_MW and H α indices reflect the variability at different heights of the stellar chromosphere and sometimes even at different latitudes. Thus, differential rotation may lead to differences in the rotation period estimates from these three variability indicators that could be larger than the error bars of individual measurements. The period of the induced RV signal should be similar but not necessarily coincident with the average period determined from activity proxies.

Once we have a measurement of the rotation period of the star, we search for periodic signals in the RV time series following the same procedure. To do that, we first subtract long-term trends, if present, using linear fits or incomplete orbits fitting sinusoidal models when the change of slope is apparent. Then we compute the power spectrum using a generalized Lomb–Scargle periodogram (Zechmeister & Kürster 2009) and iteratively fit all the detected signals with Keplerian models using the rvlin package (Wright & Howard 2012) until there are no significant remaining signals. During the process, we were able to recover all the published known planetary signals for the analysed stars (using exoplanet.eu as source). The analysis of their signals is beyond the scope of this paper. After cleaning the series from planetary signals, we search for periodic signals compatible with the rotation period previously measured.

Fig. 2 shows the periodograms for the _MW index, H α index and RV time series for the G-type star HD 41248. The three periodograms have a common periodicity that we interpret as the rotation signature of the star in the spectroscopic data. We detect a mean rotation of 23.8 d and a rotation-induced RV signal of 25.6 d with a semi-amplitude of 2.43 m s⁻¹. For the K-type star HD 125595, Fig. 3 shows a very similar picture. There is a common periodicity across the different indicators that marks the rotational modulation. We measure a mean rotation period of 38.7 d and a rotation-induced RV signal of 36.7 d with a semi-amplitude of 2.11 m s⁻¹. Fig. 4 shows the same scenario for the case of the M-dwarf star GJ 514. There is a common periodicity in the periodograms for the S_MW index, H α index and RV time series. Once again, we interpret this periodicity as the rotation signal of the star. We measure a mean rotation period of 30.0 d and a rotation-induced RV signal of 30.8 d with a semi-amplitude of 2.06 m s⁻¹. Fig. 5 shows the phase folded fits of the three time series of the three stars using the detected periodicities.

Figure 2.

Periodograms for the S_MW index, H α index and RV time series for the G-type star HD 41248. Light red fill identifies the rotational modulation. Horizontal lines show the FAP levels. Red dotted line shows the 10 per cent FAP level, green dashed line the 1 per cent level and blue solid line the 0.1 per cent level.

Open in new tabDownload slide

Figure 3.

Periodograms for the S_MW index, H α index and RV time series for the K-type star HD 125595. Light red fill identifies the rotational modulation. Horizontal lines show the FAP levels. Red dotted line shows the 10 per cent FAP level, green dashed line the 1 per cent level and blue solid line the 0.1 per cent level.

Open in new tabDownload slide

Figure 4.

Periodograms for the S_MW index, H α index and RV time series for the M-dwarf GJ 514. Light red fill identifies the rotational modulation. Horizontal lines show the FAP levels. Red dotted line shows the 10 per cent FAP level, green dashed line the 1 per cent level and blue solid line the 0.1 per cent level.

Open in new tabDownload slide

Figure 5.

Phase folded curve using the rotational modulation for S_MW index (left), H α index (centre) and RV (right) for the G-type star HD 41248 (top panels), the K-type star HD 125595 (middle panels) and the M-dwarf GJ 514 (bottom panels). Grey dots are the raw measurements after subtracting the mean value. Red dots are the same points binned in phase with a bin size of 0.1. The error bar of a given bin is estimated using the weighted standard deviation of binned measurements divided by the square root of the number of measurements included in this bin. This estimation of the bin error bars assumes white noise, which is justified by the binning in phase, which regroups points that are uncorrelated in time.

Open in new tabDownload slide

By following this procedure for each star, we are able to find rotation-induced RV signals in 37 stars of our sample. Table 2 and Fig. 6 show the results.

Figure 6.

Period signal activity-induced RV signal versus the rotation period determined using the activity proxies. The grey dashed line shows the 1–1 relationship.

Open in new tabDownload slide

There could be multiple reasons why we cannot find the rotation-induced RV signals in the remaining 17 stars of the sample. In some situations – low activity stars, specially in the case of faint stars – the RV-induced signal might be below the noise level. In some others, the geometric pattern of activity features might cause signals that are not at the rotation period of the star, but at one of its harmonics. The inclination of the stellar axis related to our line of sight might also play a role in the relative proportion between harmonics and the shape of the RV-induced signal by altering the apparent geometric pattern. Even minimizing the amplitude to an undetectable degree in the case of stars aligned pole-on. Short-lived activity regions could also undermine the detectability of the signals. In order to measure a coherent signal, it is important for the active regions to last at least a few rotational periods. It has been seen that, in the case of some M-type stars, active regions are stable over long time spans (Robertson et al. 2015), but it does not have to be true for all stars.

4.2 RV signals induced at the harmonics of the rotation period

When measuring rotation-induced RV signals, it is expected to find also modulation with periods around the first harmonics of the rotation period. The relative strength of each of the present signals depends on the surface configuration and the inclination of the rotation axis related to our line of sight. For some particular configurations, the signals at the harmonics of the rotation are dominant (Boisse et al. 2011).

We performed an analysis equivalent to the previous search for RV signals at the harmonics of the rotation period. A signal at one of these harmonics would be considered an activity-induced signal if there is an equivalent signal in any of the available activity indicators or a strong correlation between the RV signal and at least one activity indicator.

We found activity-induced signals associated with harmonics of the rotation period in nine stars, eight of them with FAPs smaller than 1 per cent. Table 3 shows the periods and semi-amplitudes of the detected signals. Only one of the detections corresponds to a star not included in Table 2. The star GJ 581 shows an induced signal at half the rotation period but not at the rotation period. Of the eight stars where we detect a signal both at the rotation period and at one of its harmonics, only one star shows the signal larger for the harmonic. The star HD1461 shows a 0.49 m s⁻¹ signal at the rotation period and 0.72 m s⁻¹ at the harmonic. In the seven remaining stars, the RV signal at the harmonic is always smaller.

Fig. 7 shows the combined periodogram for the RV of all the stars from Table 1, after removing the planetary signals, with the time axis normalized at the rotation period of each individual star. The rotation signal is the dominant feature in this combined periodogram, being the only clear signal. There is a marginal structure at half the rotation period, telling us once more that some of the stars have signals at the first harmonic of the rotation. No other prominent features are present in our combined periodogram.

Figure 7.

Combined periodogram of the RV time series of all of our stars with detected rotational modulations, after subtracting all planetary signals, with the period axis normalized at the rotation period. Light red fill shows the position of the rotation period and the first harmonics.

Open in new tabDownload slide

4.3 Implications for planet-hosting candidates

In the process of our search for rotation-induced RV signals, we recover all the previously published RV signals attributed to extrasolar planets. For most of them, we agree on the planetary origin of the signal, but there are four cases where we think a rotation-induced signal is more likely to be causing the RV modulation. The cases of HD 40307 e, GJ 676 A e, GJ 667C d and GJ 832 c.

HD 40307 was initially claimed to be the host of three super-Earths at 4.3, 9.6 and 20.4 d (Mayor et al. 2009). Then the system was expanded to six planets (Tuomi et al. 2013), planet e being a small super-Earth at ∼34.6 d orbital period. This period coincides within the uncertainties with our stellar rotation detection (36.5 d) and therefore we ascribe to rotation the RV signal, as we previously pointed in Suárez Mascareño et al. (2015).

GJ 676 A has four detected planet candidates. Planet b detected at ∼1050 d (Forveille et al. 2011) and planets c at ∼4400 d, d at ∼3.6 d and e at ∼35.4 d (Anglada-Escudé & Tuomi 2012). We recover the signals attributed to planets b, d and e and the trend attributed to planet c. As suggested in Suárez Mascareño et al. (2015), the RV signal of planet e is most likely caused by stellar rotation.

GJ 667 C is a star that hosts six planet candidates. Planets b at ∼7.2 d and c at ∼28.1 d (Bonfils et al. 2013), and planets d at ∼91.6 d, e at ∼62.2 d, f at ∼39.0 d and g at ∼256 d (Gregory 2012; Anglada-Escudé et al. 2013). Robertson & Mahadevan (2014) and Feroz & Hobson (2014) claimed that planets d and e are artefacts induced by stellar rotation, while Delfosse et al. (2013) had already identified the 91 d signal as an alias of an ∼106 d activity signal. We support the idea of the signal at ∼91 d being an artefact caused by the stellar rotation. We agree on the ∼106 d rotation period and we also found a power excess at ∼91 d in the S_MW periodogram (Suárez Mascareño et al. 2015), supporting the idea that the signal might be the mark of rotation measurements at different latitudes.

GJ 832 hosts two potential planets. Planet b at ∼3416 d (Bailey et al. 2009) and planet c at ∼35.7 d (Wittenmyer et al. 2014), the latter with an orbital period very close to our rotation period measurement. As in Suárez Mascareño et al. (2015), we interpret the signal of planet c as the rotation-induced RV signal (see also Bonfils et al. 2013).

5 DISCUSSION

The previous analysis of RV signals induced by stellar rotation provided a collection of 35 stars where we were able to measure the rotation RV signal and/or one of its harmonics, for stars going from late-F to mid-M-type. For our sample, the distribution of RV-induced signals peaks at 1–2 m s⁻¹. Fig. 8 shows the distribution of the semi-amplitudes of the periodic RV signal induced by the rotation of the stars for the different spectral types.

Figure 8.

Distribution of semi-amplitudes for the rotational modulations of the stars in our sample.

Open in new tabDownload slide

Figure 9.

Semi-amplitude of the rotation-induced signal versus chromospheric activity level |$\log _{10}(R^{\prime }_\textrm{HK})$|⁠. Empty symbols correspond to our data, filled symbols the literature data. The grey dashed lines show the best fit to the data for the two groups of stars: FGK-type stars in one group and M-type stars in the other.

Open in new tabDownload slide

Table 5 shows the parameters for the different fits for the two groups of stars.

We see some discrepancy between our results in Corot-7 and GJ 674 and the results of Boisse et al. (2011) and Bonfils et al. (2007). Our measurement of the amplitudes for the activity-induced RV signals is much smaller. In both cases, we are using much longer observation baselines, averaging the amplitude over longer periods of time. In the case of Corot-7, Boisse et al. (2011) performed the analysis during a time when the activity level of the star was higher than the average level of our full data set, and therefore a larger RV amplitude for the activity-induced signal is expected. In the case of GJ 674, the sampling rate of the modern data is much better at the scale pertinent for rotation period measurements than in the original data set. The amplitude of the rotation-induced RV signal might have been overestimated. For HD 41248, α Cen B and GJ 176, our measurements are compatible with those previously published. We included the cases of HD 166435, the Sun, GJ 3998 and GJ 15 A as a comparison. The RV semi-amplitude of their rotation signals match the expected value quite well according to their mean activity levels. See Table 4.

Even in the case of slow rotators, the RV signals induced by stellar rotation have semi-amplitudes larger than those of terrestrial planets in the habitable zones. In the case of FGK stars, these signals can drop below the 1 m s⁻¹ threshold for stars that are more quiet than the Sun, reaching sometimes the HARPS precision limit. For M-dwarfs even in the case of very slow rotators, these signals are larger than 1 m s⁻¹. We have not detected any rotation-induced RV signal smaller than 1 m s⁻¹ in an M-dwarf star. Even in the case of the most quiet stars, the RV signal induced by the rotation is still larger than those of terrestrial planets. Modelling and removing these signals become necessary.

Figure 10.

Semi-amplitude of the rotation-induced signal versus the semi-amplitude of the rotation signal in the S_MW index. There are two distinct populations, with a break at K2.5-type stars.

Open in new tabDownload slide

Table 6 shows the measured values for the parameters in eq-uation (2).

Finally, we study the phase difference between the induced RV periodic signals and the periodic modulation in the activity indicators. Bonfils et al. (2007) and Santos et al. (2014) found that these signals are not necessarily in phase. In Fig. 11, we show the phase difference between the RV signal and the S_MW signal against the colour of the stars and the mean activity level. We restrict this analysis to those stars where the period estimates of different indicators were consistent within error bars. We see an apparent evolution of the phase shift with the RV signal lagging behind the S_MW signal towards redder stars. Starting with G-type stars, which show small dispersion close to shift zero, and gradually increasing towards K-type stars that reach a phase shift close to 360° (zero again) for mid/late K-type stars, a trend that seems to continue for M-type stars. We do not find any correlation between the phase shift and the chromospheric activity level of the stars. Further investigation on a larger sample will be needed to confirm this behaviour.

Figure 11.

Phase shift between the RV and the S_MW signals against the colour B − V (top panel) and the mean |$\log _{10}(R^{\prime }_\textrm{HK})$| (bottom panel).

Open in new tabDownload slide

Then we study the phase shift between the RV signal and the H α index signal finding no clear relationship.

We also study the phase shift between the H α index and the S_MW index and find a group of stars (HD 30495, HD 59468, Corot-7, HD 209100, GJ 676A, GJ 536, GJ 667C, GJ 674 and GJ 526) show a difference in phase close to 180°, while the rest of the stars show a difference in phase close to zero. We do not see a clear correlation between the phase shifts and the average level of activity of each of these two groups.

The distinct behaviour of the different spectral types seen in Figs 9 and 10, along with phase shift between the RV signal and the S_MW signal (Fig. 11), gives a clue about the nature of the dominant surface features causing the induced RV variations. The relationship between the Ca II H&K emission and the effective change in RV varies for the different spectral types. Ca II H&K emission is generated mainly in the stellar plages (Shine & Linsky 1974), which are the main source of stellar-induced RV variations in slowly rotating solar-type stars (Lockwood et al. 2007; Dumusque, Boisse & Santos 2014; Shapiro et al. 2014). In this scenario, and assuming a high convection level – as in the case of the Sun – it is expected that the RV signal and the activity signal to be more or less in phase, as we see in our G-type stars. In cooler stars, it seems that spots become more and more important. An RV variation cased solely by spots would be expected to show a 90° shift between the RV signal and the activity index signal. The gradual change in the phase might tell us that the equilibrium between the contributions of both types of activity regions gradually changes towards smaller stars, with more complex contributions causing a wide variety of phase shifts between the two signals.

6 CONCLUSIONS

Using data for 55 late-type stars (from late-F to mid-M), we have analysed the RV time series searching for periodic signals that match the stellar rotation periods. For 37 stars, we have clearly found the induced RV signal by the stellar rotation or one of its harmonics, and measured its semi-amplitude.

Our study supports previously reported doubts on the Keplerian origin of the periodic RV signals attributed to planets HD 40307 e, GJ 676 A e, GJ 667C d and GJ 832 c, which can also be explained as rotation-induced activity signals.

We have investigated the correlation between the level of chromospheric emission, represented by the |$\log _{10}(R^{\prime }_\textrm{HK})$| index, and the measured semi-amplitude, and obtained a specific linear relationship between this index and the logarithm of the amplitude for each spectral type.

We have also investigated the correlation between the amplitude of the RV signals and the amplitude in the activity signals, and we find two different correlations for G to mid-K stars and late-K to mid-M dwarfs.

We have studied the phase shift between the period RV-induced signal caused by stellar rotation and the period signal from the S_MW index. We find an apparent evolution of the phase shift and the colour B − V. We do not find a correlation between the phase shift and the activity level.

The systematic measurement and characterization of stellar activity-induced RV signals is a necessary step for a reliable identification of the RV signals produced by terrestrial planets with short orbits. For the same activity level, M-type stars show larger activity-induced RV signals than G- and K-type stars. Because of their lower stellar mass, both very low-activity M dwarfs and late-K dwarfs offer a very good opportunity for the detection of terrestrial planets.

Acknowledgments

This work has been financed by the Spanish Ministry project MINECO AYA2014-56359-P. JIGH acknowledges financial support from the Spanish MINECO under the 2013 Ramón y Cajal program MINECO RYC-2013-14875. This work is based on data obtained from the HARPS public data base at the ESO. This research has made extensive use of the SIMBAD data base, operated at CDS, Strasbourg, France and NASAs Astrophysics Data System. We are grateful to all the observers of the following ESO projects, whose data we are using: 60.A-9036, 072.C-0096, 073.C-0784, 073.D-0038, 073.D-0578, 074.C-0012, 074.C-0364, 074.D-0131, 075.D-0194, 076.C-0878, 076.D-0130, 076.C-0155, 077.C-0364, 077.C-0530, 078.C-0044, 078.C-0833, 078.D-0071, 079.C-0681, 079.C-0927, 079.D-0075, 080.D-0086, 081.C-0148, 081.D-0065, 082.C-0212, 082.C-0308, 082.C-0315, 082.C-0718, 083.C-1001, 083.D-0040, 084.C-0229, 085.C-0063, 085.C-0019, 085.C-0318, 086.C-0230, 086.C-0284, 087.C-0368, 087.C-0831, 087.C-0990, 088.C-0011, 088.C-0323, 088.C-0353, 088.C-0662, 089.C-0050, 089.C-0006, 090.C-0421, 089.C-0497, 089.C-0732, 090.C-0849, 091.C-0034, 091.C-0866, 091.C-0936, 091.D-0469, 180.C-0886, 183.C-0437, 183.C-0972, 188.C-0265, 190.C-0027, 191.C-0505, 191.C-0873 and 282.C-5036.

REFERENCES