A Likely Magnetic Activity Cycle for the Exoplanet Host M Dwarf GJ 3512

For the analysis of the radial velocity curve, we used an adaptive importance sampling (AIS) method. In particular, we implemented a version of the algorithm described in Table 2 of Bugallo et al. (2017), with adaptation subject to the quadratic distance between the observations and the model. In particular, the proposal pdf is adapted at some step only if a sample with a lower quadratic distance is found. This distance is used to determine the noise variance that is included in the likelihood function in the next step. In this sense, the scheme acts as a maximum-likelihood algorithm. AIS methods have been used in other research fields but they have been rarely applied by the astrophysical community (Wraith et al. 2009; Lewis & Bridle 2011). Compared to MCMC methods, AIS techniques offer a more natural solution to the problem of model selection because they yield direct estimates of the model evidence without the post-processing needed in MCMC (Stokes et al. 2017; Speagle 2020).

The implementation of the AIS we used for this work is parallelized. This parallelization enables the reduction of the computational run-time and the use of a large number of samples at each step. Note that in most of the MCMC methods commonly applied in the astrophysical literature, one sample is generated at each step (e.g., Gregory 2011; Foreman-Mackey et al. 2013). For our study, we generated 10⁵ samples and performed 1000 iterations. At the end of the process, a total of 10⁸ samples were generated for each configuration, i.e., for the model with a single planet and that of two planets. The Monte Carlo sample that maximizes the likelihood was determined together with the noise variance (γがんま parameter) and the Bayesian evidence. A resampling step was performed at the end of the process to determine confidence levels for each orbital parameter.

We searched for the signature of up to two planets in the radial velocity data of GJ 3512. This selection is supported by the frequency analysis of the radial velocity curve (Section 3). Only Keplerian orbits were considered. Commonly, each single orbit is represented by five parameters: period (P), amplitude (K), eccentricity (e), argument of periastron (ωおめが), and last periastron passage (τたう). The mean velocity of the star (v₀) is a common factor to all the orbits. Therefore, we inverted the following equation:

$$\begin{eqnarray}&&{\boldsymbol{v}}={v}_{0}+\displaystyle \sum _{i=1}^{n}{K}_{i}\left[\cos \left({{\boldsymbol{\nu }}}_{i}+{\omega }_{i}\right)+{e}_{i}\cos ({\omega }_{i})\right]+{\boldsymbol{s}},\end{eqnarray} \tag{ 1 }$$

where

$$\begin{eqnarray}&&{{\boldsymbol{M}}}_{i}=\displaystyle \frac{2\pi }{{P}_{i}}\left({\boldsymbol{t}}-{\tau }_{i}\right),\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{{\boldsymbol{M}}}_{i}={{\boldsymbol{E}}}_{i}-{e}_{i}\sin \left({{\boldsymbol{E}}}_{i}\right),\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\tan \displaystyle \frac{{{\boldsymbol{\nu }}}_{i}}{2}=\sqrt{\displaystyle \frac{1+{e}_{i}}{1-{e}_{i}}}\tan \displaystyle \frac{{{\boldsymbol{E}}}_{i}}{2},\end{eqnarray} \tag{ 4 }$$

and where n is the number of planets we consider and ${\boldsymbol{t}}$ is the time vector of the observations. ${{\boldsymbol{M}}}_{i}$ is the mean anomaly of the ith component, ${{\boldsymbol{E}}}_{i}$ is its eccentric anomaly, and ${{\boldsymbol{\nu }}}_{i}$ is the true anomaly. The vector ${\boldsymbol{s}}$ is a random variable with zero mean and covariance matrix ${\boldsymbol{\Sigma }}$ and it represents observational noise. The matrix ${\boldsymbol{\Sigma }}$ can be expressed as

$$\begin{eqnarray}&&{\boldsymbol{\Sigma }}=\gamma {{\boldsymbol{\Sigma }}}_{\mathrm{obs}},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{{\boldsymbol{\Sigma }}}_{\mathrm{obs}}={{\boldsymbol{\sigma }}}_{\mathrm{obs}}^{\top }{\boldsymbol{I}}{{\boldsymbol{\sigma }}}_{\mathrm{obs}}.\end{eqnarray} \tag{ 6 }$$

Here, ${\boldsymbol{I}}$ is the identity matrix, ${{\boldsymbol{\sigma }}}_{\mathrm{obs}}$ is the standard deviation of the radial velocity determined at ${\boldsymbol{t}}$, ${{\boldsymbol{\sigma }}}_{\mathrm{obs}}^{\top }$ is the transpose of the standard deviation, and γがんま is the variance of the noise to be inferred by optimization (see Section 4.1). The scalar ${\gamma }^{1/2}$ is commonly interpreted as a jitter induced by the stellar magnetic activity (e.g., Oshagh et al. 2017). However, it contains also information of the model truncation.

Prior and proposal pdfs were defined for each parameter. Uniform density functions were selected for the prior pdfs. For the proposal pdfs, we considered uniform density functions only for e and ωおめが, while we used normal pdfs for the remaining parameters. The choice of normal proposal pdfs permits easy adaptation in the algorithm. The mean and variance of those normal pdfs were chosen according to the data. For instance, the semi-amplitude of the radial velocity curve is 70–80 m s⁻¹ and, thus, the mean of the normal pdf for K was chosen inside that range. The standard deviation for the same parameter was chosen large enough to permit the algorithm to jump to other regions of the parameter space. Similar arguments are applied to the remaining proposal pdfs. The algorithm is adaptive and will explore the parameter space conveniently. Note that we did not fix the periods selected from the frequency analysis (Section 3). Windowing also causes leakage (Harris 1978; Phillips et al. 2008), which shifts the peaks in the periodogram if the time window does not cover exactly a multiple of the period of the signal. Therefore, the peaks in a periodogram are not always located at the actual frequency of the signal. AIS algorithms adapt the proposal pdf at each step. In this work, the proposal pdfs were adapted according to the Monte Carlo sample that maximizes the likelihood. Table 1 lists the prior and initial proposal pdfs for the two components considered in this work. Note that we adapted only the Gaussian proposals. The uniform proposals are preserved at each step.

Table 1.  Prior and Proposal pdfs Used in the Analysis of the Radial Velocity Curve of GJ 3512

	Planet 1			Planet 2
Parameter	Prior	Proposal		Prior	Proposal
P (day)	${ \mathcal U }(150,250)$	${ \mathcal N }(202,20)$		${ \mathcal U }(700,3000)$	${ \mathcal N }(1300,100)$
K (m s−1)	${ \mathcal U }(60,75)$	${ \mathcal N }(70,10)$		${ \mathcal U }(10,20)$	${ \mathcal N }(15,5)$
e	${ \mathcal U }(0,0.6)$	${ \mathcal U }(0,0.6)$		${ \mathcal U }(0,0.1)$	${ \mathcal U }(0,0.1)$
ωおめが (rad)	${ \mathcal U }(0,2\pi )$	${ \mathcal U }(0,2\pi )$		${ \mathcal U }(0,2\pi )$	${ \mathcal U }(0,2\pi )$
τたう (day)	${ \mathcal U }(0,250)$	${ \mathcal N }(0,20)$		${ \mathcal U }(0,2000)$	${ \mathcal N }(0,20)$
Common to the Configurations with 1 and 2 Planets
v0 (m s−1)	${ \mathcal U }(-10,10)$ ${ \mathcal N }(0,5)$

Download table as:  ASCIITypeset image

The Bayesian (or model) evidence Z was determined for each model. We calculated Z ratios to settle the most probable number of planets. The results are shown in Table 2. The model with one planet (hereafter M1) is clearly less probable than the model with two planets (M2). Note that the ratios in Table 2 are in logarithmic scale. The presence of the less massive planet is also supported by the value of ${\gamma }_{{\rm{M}}2}$, which is lower than ${\gamma }_{{\rm{M}}1}$. Another indicator for model selection is the Bayesian information criterion (BIC) of Schwarz (1978). The BIC is a logarithmic magnitude that accounts for the increase in the number of parameters. It is defined as

$$\begin{eqnarray}&&\mathrm{BIC}=k\mathrm{ln}n-2\mathrm{ln}\widehat{L},\end{eqnarray} \tag{ 7 }$$

where k is the number of dimensions (the number of parameters of the model to be estimated), n is the number of observations (data points), and $\widehat{L}$ is the maximum value of the likelihood function for the candidate model. If the model errors are assumed independent and normally distributed, the BIC can be determined as

$$\begin{eqnarray}&&\mathrm{BIC}=n\mathrm{ln}\widehat{{\sigma }^{2}}+k\mathrm{ln}n,\end{eqnarray} \tag{ 8 }$$

where $\widehat{{\sigma }^{2}}$ is the error variance for $\widehat{L}$. The method chooses among all the (model) variances explored, the one which maximizes the likelihood. The comparison of the value of the BIC for M1 and M2 shows that

$$\begin{eqnarray}&&{\rm{\Delta }}\mathrm{BIC}={\mathrm{BIC}}_{{\rm{M}}1}-{\mathrm{BIC}}_{{\rm{M}}2}=342,\end{eqnarray} \tag{ 9 }$$

which indicates that M2 is much more probable than M1 (the model with the lowest BIC value is preferred).

Table 2.  Model Probability Tests

Model M1		Model M2
$\mathrm{log}\tfrac{{Z}_{{\rm{M}}1}}{{Z}_{{\rm{M}}1}}$	${\gamma }_{{\rm{M}}1}$	$\mathrm{log}\tfrac{{Z}_{{\rm{M}}2}}{{Z}_{{\rm{M}}1}}$	${\gamma }_{{\rm{M}}2}$
0	11.09	$1.04\times {10}^{3}$	1.32

Note. M1 is for the configuration with one planet, M2 is for a configuration with two planets

Download table as:  ASCIITypeset image

The orbital parameters we inferred for both models are listed in Table 3, together with the masses of the planets (we assumed ${M}_{\star }=0.123$ ${\text{}}{M}_{\odot };$ Morales et al. 2019). The best fit (M2) is plotted in Figure 4. For any of the models (M1 and M2), the orbit and the mass of GJ 3512b determined in this work are essentially the same as those obtained by Morales et al. (2019). In model M2, the orbital parameters inferred for GJ 3512c are in agreement with the ranges proposed by Morales et al. (2019) for their suggested second planet. We inferred a period of $P\approx 1600$ days, a semi-amplitude of K = 14.1 m s⁻¹, and an eccentricity of $e\approx 0$ for this planet. The minimum mass derived for GJ 3512c, $m\sin i=0.20\,{M}_{J}$ ($\sim 50\,{M}_{\mathrm{Earth}}$), is in the range of the super-Neptune planets. This configuration perfectly reproduces the observed periodogram. Figure 5 shows the periodogram of the data (continuous, black line) and that for the simulated curve of the two planets. Contrarily to the simulated curve for a single planet, this configuration is able to reproduce the periodogram at low frequencies (small window in the two figures). As a result, the presence of other planets in the system cannot be inferred from the frequency analysis of the current data. We notice that the small uncertainty in the period of GJ 3512c is a product of the joint posterior distribution of these data, which is very narrow. The AIS algorithm draws many samples and it assigns weights to all them. Then, it adapts the mean values and the variances of the normal proposal pdfs according to the weighted means and variances of the drawn samples and repeat the process until the requested number of iterations is completed. At the end of the process, hundreds of millions of samples along with their corresponding weights are generated. A resampling process is performed with those weights to determine the posterior distribution. For this data set, the marginal posterior distribution of the period of GJ 3512c is very narrow. Nevertheless, other observations could give slightly different results. Note that this procedure is different from the MCMC methods, where all the accepted samples are used to construct the marginal posterior distributions without weighting. This result is intuitive in some sense and also extremely important, in our opinion. We believe that this point deserves further studies. Indeed, without considering this result, one may believe that the chain produced by an MCMC algorithm or the weighted samples generated by an IS scheme fully represent the posterior pdf but, actually, they just represent the tails of the posterior distribution (i.e., meaningless regions). Namely, it is extremely important to discover the narrow region (a peak) of very high values of the posterior.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3.  Orbital Parameters and Mass of the Planets Revealed with Each Model in Table 2

Parameter	GJ 3512b	GJ 3512c
Model M1
v0 (m s−1)	1.89 ± 0.23
P (day)	${204.46}_{-0.05}^{+0.06}$	⋯
K (m s−1)	${67.9}_{-0.1}^{+0.1}$	⋯
e	${0.37}_{-0.02}^{+0.03}$	⋯
ωおめが (rad)	${2.32}_{-0.01}^{+0.01}$	⋯
${\tau }^{1}$ (day)	${199.6}_{-0.4}^{+0.3}$	⋯
a (au)	${0.338}_{-0.002}^{+0.001}$	⋯
$m\sin i$ (MJ)	${0.45}_{-0.01}^{+0.02}$	⋯
Model M2
v0 (m s−1)	5.52 ± 0.23
P (day)	${203.69}_{-0.02}^{+0.09}$	${1599.6}_{-0.8}^{+1.1}$
K (m s−1)	${71.8}_{-0.1}^{+0.2}$	${15.1}_{-0.2}^{+0.2}$
e	${0.44}_{-0.01}^{+0.01}$	${0.0183}_{-0.0001}^{+0.0001}$
ωおめが (rad)	${2.195}_{-0.012}^{+0.006}$	${0.46}_{-0.01}^{+0.03}$
${\tau }^{1}$ (day)	${196.5}_{-0.4}^{+0.5}$	${1119.2}_{-0.3}^{+0.4}$
a (au)	${0.337}_{-0.001}^{+0.001}$	${1.292}_{-0.003}^{+0.003}$
$m\sin i$ (MJ)	${0.46}_{-0.01}^{+0.02}$	${0.20}_{-0.01}^{+0.01}$

Note. Last periastron passage prior to the first observation.

Download table as:  ASCIITypeset image

The process of Bayesian inference was repeated with the measures from the VIS and NIR channels separately, for the model M2. The value determined for each orbital parameter and each planet is shown in Table 4. The results obtained previously for the entire VIS+NIR data set are included in this table for comparison. Table 4 shows that the inference method obtains very similar results for the orbital parameters of GJ 3512b with any of the three combinations: VIS, NIR, and VIS+NIR. Likewise, the orbital parameters of GJ 3512c inferred with the VIS channel data are very similar to those determined for the entire data set. Contrarily, the value of the period of GJ 3512c inferred with the NIR data is very different. Nevertheless, this value is very uncertain. This result shows that the marginal posterior pdf has a strong dependence on the data, especially if the period of the planet is not fully covered.

Table 4.  Inferred Orbital Parameters for Model M2 Using Data from the VIS and NIR Channels, Compared to the Results for the Combination of the Two Channels (Table 3)

	VIS Channel		NIR Channel		VIS + NIR Channels
Parameter	GJ 3512b	GJ 3512c	GJ 3512b	GJ 3512c	GJ 3512b	GJ 3512c
v0 (m s−1)	7.01 ± 0.37		5.56 ± 0.65		5.52 ± 0.23
P (day)	${203.90}_{-0.23}^{+0.02}$	${1594.3}_{-4.8}^{+5.1}$	${203.44}_{-0.34}^{+0.26}$	${1395.5}_{-93.2}^{+110.5}$	${203.69}_{-0.02}^{+0.09}$	${1599.6}_{-0.8}^{+1.1}$
K (m s−1)	${72.73}_{-0.76}^{+0.18}$	${15.6}_{-0.4}^{+0.3}$	${73.69}_{-1.1}^{+1.3}$	${14.5}_{-0.7}^{+1.1}$	${71.8}_{-0.1}^{+0.2}$	${15.1}_{-0.2}^{+0.2}$
e	${0.43}_{-0.01}^{+0.02}$	${0.0043}_{-0.0015}^{+0.0060}$	${0.45}_{-0.01}^{+0.01}$	${0.0078}_{-0.0078}^{+0.0604}$	${0.44}_{-0.01}^{+0.01}$	${0.0183}_{-0.0001}^{+0.0001}$
ωおめが (rad)	${2.182}_{-0.082}^{+0.025}$	${0.0007}_{-0.0007}^{+0.0023}$	${2.183}_{-0.022}^{+0.017}$	${0.0008}_{-0.0008}^{+0.0213}$	${2.195}_{-0.012}^{+0.006}$	${0.46}_{-0.01}^{+0.03}$
${\tau }^{1}$ (day)	${196.3}_{-0.1}^{+0.3}$	${1016.5}_{-0.9}^{+0.7}$	${197.1}_{-1.6}^{+1.4}$	${942.6}_{-150.3}^{+191.8}$	${196.5}_{-0.4}^{+0.5}$	${1119.2}_{-0.3}^{+0.4}$

Note. Last periastron passage prior to the first observation.

Download table as:  ASCIITypeset image

GJ 3512 is a moderately active M5.5 dwarf. Morales et al. (2019) determined a rotation period of 87 ± 5 days. The authors found peaks in the periodograms of photometric data of the star at different bands consistent with previous determinations by Newton et al. (2016). However, signatures of this rotation period are not clearly detected in the magnetic activity indices. Figure 6 shows the GLS periodograms for those indices (H${}_{\alpha }$ line width, differential line width, and chromatic index) compared to the periodogram of the R-band photometric data. The dashed line delimits the false-alarm probability of 0.01%. Only the photometric data show clearly the signature caused by the rotation period of the star. In the R band, the peak in the periodogram is observed at ∼86 days. In the periodogram of the ${{\rm{H}}}_{\alpha }$ line and differential line profile measures, a feature can be observed around that frequency but with very low power. This result is not surprising because the ${{\rm{H}}}_{\alpha }$ line width of M dwarfs is a tracer of the stochastic activity of the star (e.g., Crespo-Chacón et al. 2006).

Zoom In Zoom Out Reset image size

To produce the periodogram of the R-band photometric observations in Figure 6, we used all the available data. An offset was applied to the values of the first data set like in Morales et al. (2019). This offset compensates zero-point differences between the two TJO observation campaigns. The value of the offset was determined by comparing the means of the observations in the first observation period to those of the second observation period. The corrected light curve is shown in Figure 7. The first set of observations shows a decrease in the magnitude of the star with a duration of approximately 47 days. After that, it remains nearly flat. Contrarily, the R-band photometric data acquired later shows more variability. The period of photometric observations coincides with the last two sets of radial velocity observations (see Figure 4). Therefore, we constructed a new GLS periodogram with only R-band photometric data from the last data set, which expands for approximately 300 days. The result is shown in Figure 8. In this periodogram, two intense signals are detected. The signal at a lower frequency corresponds to a period of 86 days. The second, and more intense, peak is the first harmonic of the previous signal and it is located at a period of ∼42 days. These two signals were noticed by Morales et al. (2019, supplementary material) but the authors obtained lower significance possibly due to the mixture of data from different instruments. The frequency ratio of both peaks is 2.02. The presence of harmonics in the periodogram of photometric data from stars was investigated by Reinhold & Arlt (2015). The authors simulated the light curves of stars with a different number of spots at distinct latitudes and obtained some relations between the heights of the peaks of the harmonics and the latitude of the spots. In particular, they concluded that low-latitude spots cause intense peaks of the harmonics in the periodogram and less sine-shaped light curves. In Figure 8, the intensities of the fundamental period and the first harmonic are similar and the light curve is far from a sine shape. If the argument of Reinhold & Arlt (2015) is valid, GJ 3512 would contain a large spot at a very low latitude.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Although the hypothesis of the low-latitude single spot is attractive, a very intense first harmonic peak as that detected by us was not derived from the simulations of Reinhold & Arlt (2015). Neither did Santos et al. (2017) obtain very intense first harmonics in their simulations of single spots at distinct latitudes. High first harmonic intensities like that observed in GJ 3512 were detected by Balona (2013) in several A-type stars. The author attributed the intensity of the harmonics to the location of the spots and the rotation axis (see also Balona et al. 2015). Another plausible explanation for the intense first harmonic in the periodogram of GJ 3512 is the presence of two spots at longitudes separated by ∼180°. A similar configuration was observed in the M4 dwarf GJ 1243 (Davenport et al. 2015). In that star, the primary spot was stable over years, while the secondary spot decayed in 100–500 days. These timescales are compatible with the observations of GJ 3512, where the fundamental period of 84–90 days seems stable but the first harmonic is detected only in the last data set. The photometric curve of GJ 3512 is also compatible with the observation of a flip-flop effect (Berdyugina 2005; Hackman et al. 2013). This effect consists in the inversion in longitude of the location of the spots on the star. During the inversion process, spots of similar size coexist at longitudes separated by approximately 180° (e.g., Berdyugina et al. 1999). This pattern may produce the first harmonic in the periodogram like in GJ 3512 (Figure 8).

GJ 3512 shows an emission ${{\rm{H}}}_{\alpha }$ line (Gizis et al. 2002), with $\mathrm{log}{L}_{{{\rm{H}}}_{\alpha }}/{L}_{\mathrm{bol}}=-4.4$ dex (Reiners et al. 2018). The latter is common to many old M dwarfs (López-Santiago et al. 2010). No magnetic activity saturation is expected for this star (see Figure 7 in Martínez-Arnáiz et al. 2011). Therefore, it is likely that GJ 3512 has a long activity cycle similar to that of the Sun. Magnetic activity cycles have been observed in several low-mass stars. Lorente & Montesinos (2005) studied the relation between the cycle length and the rotation period of a number of active stars. According to their results, GJ 3512 should have a magnetic activity cycle of ${P}_{\mathrm{cycle}}\gt 20$ yr. However, Böhm-Vitense (2007) showed that inactive stars have shorter activity cycles than those observed in active stars with the same rotation rate. From that study, a cycle length of approximately 15 yr can be inferred for GJ 3512. This value is similar to the activity period found for GJ 551, another M5.5 inactive dwarf (Suárez Mascareño et al. 2016).

In view of the previous discussion, we investigated the presence of the magnetic activity cycle in the radial velocity data. The signature of the 87 day period is not observed in the periodogram of the radial velocity curve (Figure 1), but it is common that a signal passes unnoticed in a periodogram if it is not persistent in time. Note that the power of the signal generated by GJ 3512b is very high and it can mask the remaining signals (Section 3). Therefore, an amplitude modulated sinusoid was fitted to the residuals using AIS. Our aim was to find the magnetic activity cycle of GJ 3512. The model used was

$$\begin{eqnarray}&&{\boldsymbol{v}}=a\sin \left(\displaystyle \frac{2\pi }{{P}_{\mathrm{cycle}}}{\boldsymbol{t}}+{\phi }_{1}\right)\sin \left(\displaystyle \frac{2\pi }{{P}_{\mathrm{rot}}}{\boldsymbol{t}}+{\phi }_{2}\right)+b,\end{eqnarray} \tag{ 10 }$$

where a is the amplitude of the sinusoid, P_cycle is the period of the activity cycle, P_rot is the rotation period, ${\phi }_{1}$ and ${\phi }_{2}$ are phases, and b is an offset (for a similar analysis using solar data, see Edmonds 2016). The results are shown in Table 5 and the best fit is plotted in Figure 9. The residuals of the fit have a standard deviation of 5.66 m s⁻¹. This value is close to the mean uncertainty of the VIS channel and it is below the mean uncertainty of the NIR channel. It is also lower than the standard deviation of the residuals from the two-planet fit, which is close to 8 m s⁻¹. Note, however, that the residual dispersion is dominated by the points in the second observation period, where the curve becomes flat. The rotation period of the star is detected, with a value of ${P}_{\mathrm{rot}}\approx 88$ days, similar to the value determined previously (see Section 5.1). The result of the fit also suggests an activity cycle of approximately 14 yr. This value would be in agreement with the relation between cycle length and rotation period determined by Böhm-Vitense (2007). Nevertheless, this parameter is not well constrained due to the length of the observations and a longer cycle would be possible.

Zoom In Zoom Out Reset image size

Table 5.  Parameters of the Amplitude Modulated Sinusoid Fitted to the Residuals of the Radial Velocity Curve

Parameter	Value
a (m s−1)	${4.7}_{-2.1}^{+3.2}$
Pcycle (day)	${4950}_{-2500}^{+2200}$
${\phi }_{1}$ (rad)	${5.8}_{-0.6}^{+0.4}$
Prot (day)	${87.9}_{-1.0}^{+4.2}$
${\phi }_{2}$ (rad)	${2.2}_{-0.5}^{+1.5}$
b (m s−1)	$-{0.04}_{-0.43}^{+0.21}$

Download table as:  ASCIITypeset image