GJ 3929: High-precision Photometric and Doppler Characterization of an Exo-Venus and Its Hot, Mini-Neptune-mass Companion

We observed three transits of GJ 3929b on the nights of 2021 February 26, 2021 April 30, and 2021 September 21, using the Astrophysical Research Consortium (ARC) Telescope Imaging Camera (ARCTIC; Huehnerhoff et al. 2016) at the ARC 3.5 m Telescope at Apache Point Observatory (APO). To achieve precise photometry on nearby bright stars, we used the engineered diffuser described in Stefansson et al. (2017).

The air mass of GJ 3929 varied from 1.00 to 1.66 over the course of its observation on 2021 February 26. The observations were performed using a 30 nm wide narrowband Semrock filter centered at 857 nm (described in Stefansson et al. 2018) owing to moderate cloud coverage, with an exposure time of 33.1 s in the quad-readout mode with 2 × 2 on-chip binning. In the 2 × 2 binning mode, ARCTIC has a gain of 2 e/ADU, a plate scale of 0228 pixel⁻¹, and a readout time of 2.7 s. We reduced the raw data using AstroImageJ (Collins et al. 2017). We selected a photometric aperture of 31 pixels (707) and used an annulus with an inner radius of 70 pixels (1596) and an outer radius of 100 pixels (228).

We also observed a transit of GJ 3929b on 2021 April 30. The air mass during observations varied between 1.00 and 1.51. The observations were performed using the same Semrock filter as described previously, with an exposure time of 45 s in the quad-readout mode with 2 × 2 on-chip binning. For the final reduction, we selected a photometric aperture of 33 pixels (752) and used an annulus with an inner radius of 58 pixels (1322) and an outer radius of 87 pixels (1984).

We observed a final transit of GJ 3929b on 2021 September 21. The air mass during observations varied between 1.21 and 3.22, and the resulting scatter in data points was >3 times the values of either previous ARCTIC night (rms_20210226 = 1000 ppm; rms_20210430 = 910 ppm; rms_20210921 = 3400 ppm). Consequently, we chose not to use this final ARCTIC transit during analysis of planet b.

We checked for air-mass correlation on each night but found little evidence for any significant correlation. A plot of the ARCTIC transits used in our final analysis is visible in Figure 4.

Zoom In Zoom Out Reset image size

We additionally use publicly available data taken by the Las Cumbres Observatory Global Telescope Network (LCOGT; Brown et al. 2013). These data were obtained from the Exoplanet Follow-up Observing Program (ExoFOP) website. ³² Two transits of GJ 3929b were obtained using the LCOGT. The first transit was obtained on 2021 April 10. Data were taken by both the SINISTRO CCDs at the 1 m telescopes of the McDonald Observatory (McD) and the Cerro Tololo Inter-American Observatory (CTIO). Both instruments have a pixel scale of 0.00389 pixel⁻¹ and a field of view (FOV) of 260 × 260.

A second transit was obtained on 2021 April 15. These data were taken simultaneously in four different filters (g', ${i}^{{\prime} }$, r ', and ${z}_{s}^{{\prime} }$) with the Multi-color Simultaneous Camera for studying Atmospheres of Transiting exoplanets 3 camera (MuSCAT3; Narita et al. 2020) mounted on the 2 m Faulkes Telescope North at Haleakala Observatory (HAL). It has a pixel scale of 0.0027 pixel⁻¹, corresponding to an FOV of 9.01 × 9.01.

As outlined in Kemmer et al. (2022), high air mass caused the CTIO observations to exhibit higher scatter. In fact, both transits on 2021 April 10 exhibit much higher scatter (rms_CTIO = 3300 ppm; rms_MCD = 2200 ppm) than on 2021 April 15 (rms_gp = 1010 ppm; rms_ip = 850 ppm; rms_rp = 910 ppm; rms_zs = 920 ppm). Consequently, for the same reasons outlined in Section 2.2.1, we chose not to utilize either transit from 2021 April 10 in our final analysis.

The publicly available data were calibrated by the LCOGT BANZAI pipeline (McCully et al. 2018), and photometric data were extracted using AstroImageJ (Collins et al. 2017). The resulting photometric data are the same as those that were utilized in Kemmer et al. (2022).

We observed GJ 3929 using the ShARCS camera on the Shane 3 m telescope at Lick Observatory (Srinath et al. 2014). GJ 3929 was observed using the K_S and J filters on the night of 2021 February 26. Instrument repairs prevented our observations from benefiting from Laser Guide Star (LGS) mode. Fortunately, GJ 3929 is sufficiently bright that LGS mode is helpful but not necessary. Further instrument repairs prevented our observations from using a dither routine to create master-sky images of GJ 3929. Instead, after a series of observations, we shifted several arcseconds to an empty region of sky and took images with the same exposure time for purposes of sky subtraction.

The raw data are reduced using a custom pipeline developed by our team (described in Beard et al. 2022). Using algorithms from Espinoza et al. (2016), we then generate a 5σ contrast curve as the final part of our analysis (Figure 2). We detect no companions at a ΔK_s = 4.85 at 076 and ΔK_s = 9.75 at 835. Additionally, we detect no companions at ΔJ = 4.54 at 109 and ΔJ = 7.62 at 899.

We note that observing conditions on 26 February 2021 were marginal. As a result of overcast conditions and poor seeing, the FWHM of the centroid in each reduced AO image was fairly large (077 and 111). Consequently, our final constraints on nearby luminous companions are not as tight as they might have been. However, our high-contrast images were taken in redder wave bands than than the ${z}^{{\prime} }$ filter used in Kemmer et al. (2022), and so we provide additional sensitivity toward detecting redder, cooler companions. In tandem, our results and those outlined in Kemmer et al. (2022) are consistent: we detect no nearby luminous companions as an explanation for the observed transit event.

We obtained RVs of GJ 3929 using the new, ultraprecise NEID spectrometer (Schwab et al. 2016) on the WIYN 3.5 m telescope at KPNO. NEID is an environmentally stabilized (Stefansson et al. 2016; Robertson et al. 2019) fiber-fed spectrograph (Kanodia et al. 2018) with broad wavelength coverage (3800–9300 Å). We observed GJ 3929 in High Resolution (HR) mode with an average resolving power ${ \mathcal R }$ = 110,000. The default NEID pipeline utilizes the cross-correlation function (CCF; Baranne et al. 1996) method to produce RVs. However, this method tends to be less effective on M dwarfs (e.g., Anglada-Escudé et al. 2012), and so we use a modified version of the SpEctrum Radial Velocity AnaLyser pipeline (SERVAL; Zechmeister et al. 2018) as described in Stefansson et al. (2021). SERVAL shifts and combines all observed spectra into a master template and compares this template with known reference spectra. We then minimize the χ² statistic to determine the shifts of each observed spectrum. We mask telluric and sky emission lines during this process. A telluric mask is calculated based on their predicted locations using telfit (Gullikson et al. 2014), a Python wrapper to the Line-by-Line Radiative Transfer Model package (Clough et al. 2005).

We obtained 27 observations of GJ 3929 between 2021 January 6 and 2022 January 27. Our first two nights of observation for this system used three consecutive 900 s exposures, but we later changed our observation strategy to one 1800 s exposure per night. We obtained a median signal-to-noise ratio (S/N) of 44.8 in order 102 (λ = 4942 Å) of NEID for each unbinned observation. The median unbinned RV error bar is 1.18 m s⁻¹. The error bars are estimated from expected photon noise. A total of 23 nightly binned RVs were obtained, though four were discarded because the laser frequency comb calibrator was not available on those nights, resulting in a less precise instrument drift solution that is insufficient for precision RV analysis. This left us with 19 nightly binned NEID RVs that were used in the analysis.

We observed GJ 3929 with the Habitable Zone Planet Finder (HPF; Mahadevan et al. 2012, 2014), a near-infrared (NIR; 8080–12780 Å), high-precision RV spectrograph. HPF is located at the 10 m Hobby–Eberly Telescope (HET) in Texas. HET is a fixed-altitude telescope with a roving pupil design. Observations on the HET are queue scheduled, with all observations executed by the HET resident astronomers (Shetrone et al. 2007). HPF is fiber-fed, with separate science, sky, and simultaneous calibration fibers (Kanodia et al. 2018), and has precise, millikelvin-level thermal stability (Stefansson et al. 2016).

We extracted precise RVs with HPF using the modified version of SERVAL Zechmeister et al. (2018) optimized for use for HPF data as described in detail in Stefánsson et al. (2020). The RV reduction followed similar steps to those outlined in Section 2.4.1.

We obtained 18 observations of GJ 3929 with HPF over the course of six observing nights. These data were taken between 2021 August 27 and 2022 March 11. We obtained three consecutive exposures on each observing night, resulting in a median unbinned RV error of 7.15 m s⁻¹. Data taken on BJD = 2,459,649 were excluded from our analysis owing to poor weather conditions. Our data set then consists of five nightly binned HPF RVs. Due to the small quantity of the HPF data, we considered fits that did not utilize HPF data. We found that model results did not differ meaningfully whether HPF data were utilized or not, and we include them in our final model for completeness. HPF spectra were still used to derive stellar parameters, as outlined in Section 3.

Our RV modeling also utilizes CARMENES RVs published in Kemmer et al. (2022). Kemmer et al. (2022) published 78 high-precision RVs as a part of their study of GJ 3929 using the CARMENES spectrograph (Quirrenbach et al. 2014). CARMENES is a dual-channel spectrograph with visible and NIR arms (${{ \mathcal R }}_{\mathrm{VIS}}$= 94,600; ${{ \mathcal R }}_{\mathrm{NIR}}$ = 80,400). CARMENES is located at the Calar Alto Observatory in Almería, Spain. RVs of GJ 3929 were taken between 2020 July 30 and 2021 July 19. Each observation lasted 30 minutes, with a median observation S/N of 74. Five RVs were discarded owing to a missing drift correction in Kemmer et al. (2022), and we do so as well. This results in a final data set containing 73 RVs. These RVs were taken using the visible arm of CARMENES and have a median uncertainty of 1.9 m s⁻¹.

We modeled GJ 3929's photometry using the exoplanet software package (Foreman-Mackey et al. 2021a). First, we downloaded the TESS PDCSAP flux using lightkurve (Lightkurve Collaboration et al. 2018). We then performed a standard quality-flag filter, removing data points designated as of poor quality by the SPOC pipeline, and we median-normalized the TESS data. We then combined the TESS data with our normalized ARCTIC and LCOGT data for joint analysis.

Initial fits to ARCTIC and LCOGT data appeared to have a slight residual trend, and so in our adopted fit we detrended ARCTIC and LCOGT photometry before combining the data sets. We utilized the NumPy polyfit function to fit a line for purposes of detrending (Harris et al. 2020). This function performs a simple least-squares minimization to estimate the linear trend. This detrending was performed before modeling the data, as we found that including a detrending term in the model did not meaningfully improve our results, while increasing the complexity of our model.

We found it best to partition the photometric data into four regions of interest: the TESS data (which consist of two consecutive sectors), two different nights of ARCTIC data, and a night of LCOGT data. Due to the possibility of systematic offsets between nights and the distinct conditions during each night of ARCTIC observations, we choose to treat each ARCTIC night separately in our model. Furthermore, the LCOGT data were taken with four different filters. Consequently, we model each filter separately. For each instrument−filter combination, then, we adopt a unique mean and jitter term. The mean terms are additive offsets to account for potential systematic shifts between nights and are simply subtracted from all data points when fitting. The jitter terms are meant to model additional white noise not properly accounted for in the error bars of the data set and are added in quadrature with the error bars. Our model thus consists of seven total mean terms and seven jitter terms.

The physical transit model was generated using exoplanet functions and the starry light-curve package (Luger et al. 2019), which models the period, transit time, stellar radius, stellar mass, eccentricity, radius, and impact parameter to produce a simulated light curve. We adopt quadratic limb-darkening terms to account for the change in flux that occurs when a planet approaches the limb of a star (Kipping 2013). The two ARCTIC transits were taken using the same Semrock filter, and so we expect their limb-darkening behavior to be the same. Thus, we adopt the same limb-darkening parameters for each ARCTIC transit. We adopt distinct limb-darkening terms for the LCOGT data taken with the SDSS ${g}^{{\prime} }$, ${i}^{{\prime} }$, ${r}^{{\prime} }$, and ${z}_{s}^{{\prime} }$ filters. We note that this results in six pairs of limb-darkening terms, in contrast to seven separate jitter and mean terms, but is physically motivated.

Similar to Kemmer et al. (2022), we choose not to include a dilution term in our final model. GJ 3929 does not have many neighbors and is much brighter than all of them (Figure 3). GJ 3929 has an estimated contamination ratio of 0.000765, meaning that 0.08% of its flux is possibly from nearby sources (Stassun et al. 2019). This suggests that a dilution term is not necessary.

Zoom In Zoom Out Reset image size

After constructing a physical transit model using starry, we compare it to the data after it has been adjusted to account for offsets, and we add our jitter parameters in quadrature with the error bars during likelihood estimation. Each free parameter is given a broad prior to prevent any biasing of the model, and we summarize the priors used in Table 4. The model is then optimized using scipy.optimize.minimize (Virtanen et al. 2020), which utilizes the Powell optimization algorithm (Powell 1998). This optimization provides a starting guess for posterior inference. We then used a Markov Chain Monte Carlo (MCMC) sampler to explore the posterior space of each model parameter. exoplanet uses the Hamiltonian Monte Carlo (HMC) algorithm with a No U-Turn Sampler (NUTS) for increased sampling efficiency (Hoffman & Gelman 2011). We ran 10,000 tuning steps and 10,000 subsequent steps and assessed convergence criteria using the Gelman–Rubin (G-R) statistic (Ford 2006). We considered a chain well mixed if the G-R statistic was within 1% of unity. All the parameters in our model indicated convergence using this metric.

Our photometry-only fits are consistent with the joint fits adopted in Section 4.5. A final plot of the photometry, folded to the period of planet b, is visible in Figure 4.

Zoom In Zoom Out Reset image size

We first used a Generalized Lomb−Scargle (GLS) periodogram (Zechmeister & Kürster 2009) to analyze the RVs of GJ 3929 and to identify any periodic signals. We estimate the analytical false-alarm levels and normalize the periodogram following the steps outlined in Zechmeister & Kürster (2009), which assume Gaussian noise. With this assumption, we scale the sample variance (and false-alarm levels) by $\tfrac{N-1}{2}$ in order to reproduce the population variance, which is the quantity of interest in our analysis. Consistent with Kemmer et al. (2022), we detected significant periodicities between 14 and 16 days. In contrast to Kemmer et al. (2022), however, we find that when including the new, more precise NEID RVs (median CARMENES RV error ∼1.6× median NEID RV error), as well as our HPF RVs, the 15-day signal has grown in power relative to the 14-day signal, suggesting that it might be the true signal. Relative peak strengths of alias frequencies in a periodogram do not always indicate the true period, however, and we detail a more formal model comparison later in the section. A plot of the combined data set periodogram and periodograms on NEID and CARMENES only are visible in Figure 5. After the subtraction of the longer-period planet c, the signal of the 2.6-day planet b is clearly identifiable in the periodogram.

Zoom In Zoom Out Reset image size

We used the RadVel software package to analyze the RVs of GJ 3929 (Fulton et al. 2018). RadVel models an exoplanet's orbit by solving Kepler's equation using an iterative method outlined in Danby (1988). Each planetary orbit is then modeled by five fundamental parameters: the planet's orbital period (P), the planet's time of inferior conjunction (T_c ), the eccentricity of the orbit (e), the argument of periastron (ω), and the velocity semiamplitude (K). We additionally include instrumental terms, γ and σ, which account for systematic offsets between instruments and excess white noise.

We construct the RV model in a Bayesian context, encoding prior information about each parameter as a part of the model. Similar to the fits described in Section 4.1, we adopt broad priors on the free parameters of our model to prevent any bias in our results, the primary exception being that during RV-only fits we put tight priors on P_b and T_con,b , as these are much more tightly constrained by transits than by RV fits. We emphasize, however, that our final adopted fit is a joint fit between RVs and transits, detailed in Section 4.5. Detailed prior information is available in Table 4.

In order to estimate the posterior probability of our model, we used an MCMC sampler to explore the posterior parameter space. RadVel utilizes the MCMC sampler outlined in Foreman-Mackey et al. (2013). We first used the Powell optimization method to provide an initial starting guess for each parameter (Powell 1998). We then ran 150 independent chains and assessed convergence using the Gellman–Rubin (G-R) statistic (Ford 2006). The sampling was terminated when the chains were sufficiently mixed. Chains are considered well mixed when the G-R statistic for each parameter is <1.03, the minimum autocorrelation time factor is ≥75, the max relative change in autocorrelation time is ≤0.01, and there are ≥1000 independent draws. All of our considered models eventually satisfied these conditions.

We additionally considered the inclusion of a Gaussian process (GP; Ambikasaran et al. 2015) model to account for coherent stellar activity. Kemmer et al. (2022) identify a rotation period of ∼120 days for GJ 3929. This value is derived from a combination of long-term photometry taken using the Hungarian Automated Telescope Network (HATNet; Bakos et al. 2004), the All-Sky Automated Search for SuperNovae (ASAS-SN; Shappee et al. 2014), and Joan Oró Telescope (TJO; Colomé et al. 2010) and periodogram analysis of the CARMENES Hα values. We use the combined Hα values from CARMENES and NEID to expand on this, plotted in Figure 6. While the maximum power occurs at a slightly shorter period than observed in Kemmer et al. (2022), we note that rotational variability is often quasi-periodic in nature and periodograms can have trouble distinguishing longer periods (Lubin et al. 2021). Our value observed here is still consistent with the previously reported value, and we make no amendment to the system's rotation period.

Zoom In Zoom Out Reset image size

The >100-day rotation period of this system is consistent with a quiet, slowly rotating star, and we normally would not expect a large RV signal due to activity. However, Kemmer et al. (2022) found an RV fit that included a GP to be preferred to an RV-only fit, and so we proceed with a series of fits, some of which include a GP. Our GP fits utilize the quasi-periodic GP kernel owing to its flexibility and wide application in exoplanet astrophysics (e.g., Haywood et al. 2014; López-Morales et al. 2016).

We also compared fits with the GP kernel that was adopted in Kemmer et al. (2022). Kemmer et al. (2022) utilized a combination of two simple harmonic oscillator (SHO; Foreman-Mackey 2018) kernels, outlined in more detail in Kossakowski et al. (2021). We refer to this as the double simple harmonic oscillator kernel (dSHO).

In order to explore the possibility of an additional planet in the GJ 3929 system and the plausibility of stellar activity interfering with RV signals, we perform a model comparison. Model comparisons vary in the number of planets included, whether or not we include a GP to account for stellar noise and eccentric fits, and whether or not the second planet is modeled as the 14-day signal or the 15-day signal. A full table of our model results is provided in Table 3. Our analysis found that both the quasi-periodic and dSHO GPs perform similarly in model comparison, and so we only include the quasi-periodic results for brevity. When comparing models, we use the Bayesian information criterion (BIC; Kass & Raftery 1995) and the corrected Akaike information criterion (AICc; Hurvich & Tsai 1993). The BIC of each model can be used to estimate the Bayes factor (BF), a measure of preference for one model over another. Half the difference in BIC between two models is used to estimate the Schwarz criterion, which itself is an approximation of the log BF. The AICc is an approximation of the Kullback–Leibler information, another metric for ranking the quality of models (Hurvich & Tsai 1993).

Kass & Raftery (1995) suggest that a log₁₀ BF > 2 (ln BF > 4.6) is decisive evidence for one model over another. For GJ 3929, our two-planet (∼15 days) model is preferred over the next best model, a two-planet GP (∼15 days), with a BF of 5.86 (RadVel estimates likelihoods using ln), suggesting a strong preference for the no-GP case. The AICc simply prefers the model that minimizes the AICc, which is also the two-planet model (∼15 days). Both methods of estimation are only asymptotically correct, but they are preferred by a wide enough margin and agree with one another. Consequently, we use these comparisons to justify selecting the two-planet model (∼15 days) without a GP as our best model.

A phase folded diagram of the RV model taken from our final joint fit is included in Figure 7.