Metallicity Dependence of Giant Planets around M Dwarfs

The field M dwarf sample is constructed based on the SpeX spectroscopic survey carried out by Terrien et al. (2012) between 2011 and 2013. The whole stellar sample consisted of 886 nearby M dwarfs mostly selected according to their proper motions (Lépine & Shara, 2005; Lépine & Gaidos, 2011) with a small set of planet hosts as well as wide binaries that serve as abundance calibration stars (Rojas-Ayala et al., 2010). All targets were observed with SpeX in the short cross-dispersed (SXD) mode at a spectral resolving power of $\approx 2000$ covering the $JHK$ bands simultaneously. Later work from Terrien et al. (2015) provided the [Fe/H] and [M/H] measurements for most stars based on different calibrations (see Section 4). We refer the readers to Terrien et al. (2015) for more details about the survey description, observations, and data reduction.

To obtain precise stellar properties, we crossmatch the full catalog with Gaia DR3 (Gaia Collaboration et al., 2023) through their 2MASS identifiers (Cutri et al., 2003; Skrutskie et al., 2006) to retrieve the stellar kinematic properties. We restrict the sample to stars with parallax and proper motion measurements from Gaia, which excluded 109 stars. We do not make a specific cut on the significance of these measurements as most stars are nearby and have well-determined parallax and proper motions. We find that 24 stars in Terrien et al. (2015) do not have metallicity measurements, thus we remove them from our sample. To minimize the effect of giant planets in the field M dwarf sample, we filter out and exclude 7 known giant planet M dwarf hosts, all of which are included in our planet sample below, thus the final field M star sample ends with 746 stars. The full catalog is available in Table LABEL:field_M_star_sample. We note that we do not have prior information on whether the rest of the stars in the field M dwarf sample host giant planets since most of them do not have radial velocity follow-up observations. Therefore, the field M dwarf group could be regarded as a mixed sample, containing stars without gas giants as well as a small fraction of stars hosting giant planets. Here, we demonstrate that such a contamination will not affect the final conclusion. According to the occurrence rate measurements, we conservatively estimate that there are about two M dwarfs hosting hot Jupiters and twenty M stars harboring cold Jupiters in the field M dwarf sample (Sabotta et al., 2021; Gan et al., 2023a; Bryant et al., 2023). We randomly choose 22 stars in the field M dwarf sample, assuming they host giant planets and have the highest metallicity (0.7 dex) as in our planet sample. We repeat the simulation above for 5000 times and we find the metallicity distribution almost does not change with the median [Fe/H] varying between -0.07 and -0.06 dex for the $H$-band measurements, 0.07 and 0.08 dex for the $K$-band measurements as well as [M/H] ranging between -0.04 and -0.02 dex for the $H$-band results, 0.0 and 0.02 dex for the $K$-band results, all of which are within $1\sigma$.

We constructed the M dwarf giant planet sample with the help of information from the NASA Exoplanet Archive (Akeson et al., 2013) accessed on 2023 April 1st. We accepted all M dwarfs with effective temperature $T_{\rm eff}<4000$ K and stellar mass $M_{\ast}<0.65\ M_{\odot}$ that have confirmed giant planets with mass $0.2\ M_{\rm Jup}\leq M_{p}<13.6\ M_{\mathrm{Jup}}$ detected by either transit or radial velocity surveys. The raw selection during this step is based on stellar parameters from the literature, thus from different sources (i.e., spectroscopic facilities and methodologies). Both stellar and planet parameters are refined in a homogeneous way (see Section 4). After excluding objects that are beyond the sky coverage of SpeX (Dec$<-30^{\circ}$), a total of 22 M dwarfs including 5 multi-giant-planet systems (GJ 1148, GJ 317, GJ 3512, GJ 849 and GJ 876) are left. Figure 1 shows the color-magnitude diagram and the refined stellar mass versus effective temperature distribution of the two samples. Our planet sample eventually contains these 22 M stars with 27 giant planets.

We emphasize that our planet sample includes planetary systems detected by different surveys, and the potential bias between different selection functions is challenging to characterize. Therefore, we are not able to probe the giant planet occurrence rates in different stellar metallicity bins as other works (e.g., Fischer & Valenti, 2005). Instead, we treat the M dwarf planet sample and the field M dwarf group as independent populations and test whether their metallicities could be drawn from a single distribution. The reason is that M dwarfs with confirmed giant planets are rare. For a specific transit or RV survey, the total number of detected systems is limited, prohibiting attempts to compare the occurrence rates in several metallicity bins due to large Poisson noise.

We first determined the spectral types for the planet sample by comparing the 1d merged spectra to the IRTF spectral library (Cushing et al., 2005; Rayner et al., 2009). We utilized the SPLAT code (Burgasser & Splat Development Team, 2017) to normalize the data and find the best matches through $\chi^{2}$ minimization. In the match, we exclude the regions close to 1.1–1.2, 1.3–1.5 and 1.75–2.0 $\mu m$ with strong telluric absorption that may dominate the residuals. We show an example of our spectrum match in Figure 3.

The spectral types of the field M dwarf sample come from Terrien et al. (2015), which instead applied the NIR calibration developed in Newton et al. (2014) based on the $\rm H_{2}O$-$\rm K2$ index defined by Rojas-Ayala et al. (2012). We examined the consistency of spectral type (SpT) between two methods with the 7 overlapped giant planet host stars and we find that $\Delta_{\rm SpT}\leq 1.0$. We thus neglected the difference and used them as inputs from the same source in the following steps. The spectral type of our planet sample ranges from M0V to M5V while the field M dwarf sample spans M0V to M5.5V.

Building on the derived spectral types, we measured the iron abundance [Fe/H] and overall metallicity [M/H] following the same method as in Terrien et al. (2015). We made use of the metal algorithm (Mann et al., 2013) to determine the metallicity based on the $H$-band (Terrien et al., 2012) and $K$-band (Rojas-Ayala et al., 2012) relations. The methodology is calibrated with a sample of wide binaries containing a Sun-like primary and a low-mass M star secondary, by determining the correlation between the metallicity of the central star and the features in the companion dwarf spectra (Mann et al., 2013). Since the $J$-band calibration has a larger scatter, we chose not to use it in this work. We note that the $K$-band spectra are expected to provide more accurate metallicity results because they are less affected by the telluric contamination than $H$-band spectra (Dressing et al., 2019). The systematic errors of $H$- and $K$-band calibrations are 0.09 and 0.08 dex, which are considered in the following analysis.

Figure 4 presents the comparisons of [Fe/H] and [M/H] outputs from $H$ and $K$ band spectra with subscripts of M13,H and M13,K. Overall, we find that the results agree well with each other although the $K$-band outputs appear to be slightly higher than those from $H$-band toward the high metallicity end. The discrepancy might be due to the wavelength calibration or the data reduction procedure but the exact reason for the offset is still unclear. Therefore, we chose to independently investigate the metallicity distribution based on the results from two bands. We caution the readers these calibrations were untested for M dwarfs beyond M5 and metallicity greater than 0.56 dex and the corresponding outputs could be interpreted as extrapolations since the wide binary sample used in Mann et al. (2013) did not have stars outside this parameter space. The majority of stars in both our planet sample and field M dwarf sample are within these ranges (see Figure 4) so we consider the outliers are not expected to significantly impact the final results.

In this section, we obtain the stellar properties of our planet sample in a homogeneous way to refine the planet parameters, which enables the investigation of the relation between stellar metallicity and other planet properties. In short, we apply the empirical relations, calibrated with nearby binaries, between stellar physical parameters and absolute magnitude as well as color.

Combining the parallax from Gaia DR3 (Gaia Collaboration et al., 2023) and $m_{K}$ magnitude from 2MASS (Cutri et al., 2003; Skrutskie et al., 2006), we first computed the absolute magnitude $M_{K}$ of each M dwarf in our sample. The stellar masses $M_{\ast}$ were then estimated through the $M_{\ast}$-$M_{K}$ relation²²2https://github.com/awmann/M_-M_K- derived by Mann et al. (2019), which is suitable for stars spanning a mass range of $0.075\leq M_{\ast}\leq 0.7\ M_{\odot}$ with a systematic uncertainty of about 3%. Since the effect of [Fe/H] on the $M_{\ast}$-$M_{K}$ relation is weak, about $0.0\pm 2.2$% change in mass per dex change in [Fe/H] for stars in the solar neighborhood (Mann et al., 2019), we do not take metallicity into account here. In terms of stellar radius $R_{\ast}$, we employed the third-order polynomial relation between $R_{\ast}$ and $M_{K}$ (see Eq. 4 in Mann et al., 2015). This relation is valid for stars that have spectral types ranging from K7 to M7 with a precision of 3%.

Furthermore, we estimated the stellar effective temperature using the $T_{\rm eff}$-color relation reported in Mann et al. (2015). We retrieved the $V$-band magnitudes from APASS (Henden et al., 2016) and $JH$-magnitudes from 2MASS (Cutri et al., 2003; Skrutskie et al., 2006). We opted to use the relation that includes two stellar color terms $V-J$ and $J-H$, which gives a scatter of about 50 K. Combined with the typical spectroscopic error, we conservatively adopt a systematic uncertainty of 70 K. The final uncertainties of $M_{\ast}$, $R_{\ast}$, and $T_{\rm eff}$ are determined through error propagation after accounting for the aforementioned systematic errors of the methodology. We provide the information of the field M star catalog in Table LABEL:field_M_star_sample with each column described in Table LABEL:table_column. The stellar properties of our planet sample are listed in Table LABEL:giant_planet_M_star_sample.

Finally, we recalculated the planetary mass $M_{p}$, semi-major axis $a$, and equilibrium temperature $T_{\rm eq}$ using the updated stellar properties above together with the orbital results from the literature. For each planet, we randomly draw a set of parameters including effective temperature $T_{\rm eff}$, stellar mass $M_{\ast}$ and radius $R_{\ast}$, period $P$, eccentricity $e$, and radial velocity semi-amplitude $K$ from normal distributions, centered at the best-fits with $\sigma$ adopted as the higher value of their lower and upper uncertainty to be conservative. We fixed the eccentricity at zero if the published result is 1) smaller than 0.1 and consistent with 0 within $3\sigma$; or 2) only an upper limit. We repeated the sampling process 5000 times, and recorded the median and standard deviation of each distribution as the final value and uncertainty of planetary parameters. Here, we do not induce the constraint on the scaled semi-major axis $a/R_{\ast}$ from the light curve so the uncertainty is higher for transiting systems compared with the value in literature. Table LABEL:giant_planet_properties summarizes the refined planet properties of our planet sample.

With the [Fe/H] and [M/H] homogeneously measured based on spectral data from the same instrument through the same methodology, we investigate the metallicity distribution of the two samples we constructed. Since the behaviors of [Fe/H] and [M/H] are almost identical, we mainly focus on [Fe/H] in the following sections and put all results of [M/H] in the Appendix.

Figure 5 illustrates the iron abundance distribution and the corresponding cumulative function of field M dwarfs (black line) and M stars with confirmed giant planets (orange line) from two bands. Each two-giant-planet system is treated as two single-planet systems independently. We also repeat the same analysis below but treat those M dwarfs hosting two giant planets as single stars, which gives similar results. By visual inspection, it appears that giant planets were mostly found around metal-rich M dwarfs, suggesting a correlation between their formation and stellar metallicity, similar to FGK counterparts (Santos et al., 2004; Fischer & Valenti, 2005; Johnson et al., 2010; Sousa et al., 2011; Maldonado et al., 2020; Osborn & Bayliss, 2020). In addition, the iron abundance of M dwarfs harboring giant planets tends to show a bimodal distribution with one peak roughly located at the median [Fe/H] of field M dwarfs and the other at about 0.4 dex. However, the two-peak phenomenon is not seen in [M/H] (see Figure 9). Since the planet number in each [Fe/H] bin is small, we attribute the feature to Poisson noise and do not over-interpret it.

To statistically examine whether the field M dwarfs and M stars hosting giant planets have an identical underlying [Fe/H] distribution, we carry out Kolmogorov-Smirnov (K-S; Hodges, 1958) and Anderson-Darling (A-D; Scholz & Stephens, 1987) tests. As the sample sizes are small, we apply the permutation method for the A-D test to obtain accurate $p$-values with the number of resamples set to $10^{6}$. We calculate the K-S and A-D statistic between two samples using the functions embedded in SciPy (Virtanen et al., 2020), which yield $p$-values of $4.9\times 10^{-7}$ and $9.9\times 10^{-7}$ for the $H$-band as well as $6.3\times 10^{-5}$ and $4.2\times 10^{-5}$ for the $K$-band metallicities, corresponding to $5\sigma$ and $4\sigma$ significance, respectively. To account for the uncertainty of each metallicity measurement, we resample the $\rm[Fe/H]_{M13,H}$ and $\rm[Fe/H]_{M13,K}$ by randomly drawing samples, assuming Gaussian distributions $\mathcal{N}{\rm([Fe/H],\sigma_{[Fe/H]}^{2})}$. Since the metallicity uncertainty of the field star data set is unavailable in Terrien et al. (2015), we assign 0.09 dex and 0.08 dex systematical errors³³3These systematical errors of the methodology dominate the final uncertainty according to our findings based on the planet sample (see Table LABEL:giant_planet_M_star_sample). to the $H$ and $K$ band measurements before randomizing. We loop the procedure 10000 times, record the $p$-values, and compute the fraction of trials with $p\leq 0.003$, corresponding to the $3\sigma$ significance criterion. We list all details of the K-S along with the A-D tests in Table 5. To sum up, we find that nearly all trials lead to $p$-values smaller than the 0.003 threshold, hence indicating a strong significance. Therefore, we conclude that the null hypothesis, the iron abundance distributions of field M dwarfs and M stars with giant planets are the same, could be rejected. Our findings suggest that giant planets favor metal-rich M stars, which agree with the conclusions from previous studies (Johnson & Apps, 2009; Rojas-Ayala et al., 2010).

Some initial efforts found that hot and warm Jupiters around M dwarfs may have different preferences of stellar metallicity (Gan et al., 2022, 2023b), hinting at different formation histories (Pollack et al., 1996; Boss, 2002). However, the large uncertainty of [Fe/H] as well as the systematic biases from heterogeneous instruments and methodologies that are challenging to characterize impede these works to draw a robust conclusion. Here, we revisit the puzzle using our samples characterized through a uniform pathway.

We divide the full giant planet sample into two subclasses: hot and warm Jupiters based on the refined planet properties. Following Gan et al. (2023b), we designate hot Jupiters (HJ) as planets having $M_{p}\geq 0.2\ M_{\mathrm{Jup}}$ and $a/R_{\ast}\leq 20$ and warm Jupiters (WJ) as planets within the same mass range but $a/R_{\ast}>20$. In this way, we find 10 HJs and 17 WJs, including 5 multi-warm Jupiter systems. Among two groups, there are one HJ and one WJ whose host stars plausibly have sub-solar metallicities: NGTS-1 ($\rm[Fe/H]_{M13,H}=-0.55\pm 0.17$ dex, $\rm[Fe/H]_{M13,K}=-0.32\pm 0.11$ dex)⁴⁴4For NGTS-1, Bayliss et al. (2018) fixed an overall metallicity of $\rm[M/H]=0$ dex in the analysis but did not have an estimate on $\rm[Fe/H]$. and TYC 2187-512-1 ($\rm[Fe/H]_{M13,H}=-0.21\pm 0.09$ dex, $\rm[Fe/H]_{M13,K}=-0.21\pm 0.08$ dex)⁵⁵5For TYC 2187-512-1, Quirrenbach et al. (2022) reported two $\rm[Fe/H]$ values: $0.08\pm 0.19$ dex from Passegger et al. (2019) and $-0.18\pm 0.08$ dex from Marfil et al. (2021).. We present the metallicity distributions of two groups in Figure 5 (HJ: red line, WJ: blue line). We examine whether the HJ and WJ populations are similar to each other through the same method as in Section 5.1. We find that the $p$-values of both the K-S and A-D tests are approximately $10^{-1}$. In fact, almost none of randomly drawn samples have $p$-values falling below 0.003. All results can be found in Table 5. Similar features are also seen for [M/H] (see Table 6). Consequently, we are not able to reject the null hypothesis that two distributions are identical. At this point, we conclude that HJs and WJs present analogous preferences on the metallicity of their host M dwarfs. As an independent inspection, we vary the $a/R_{\ast}$ boundary between 15 and 50 to define HJs and WJs and rerun the analysis. We find that the conclusion does not depend on the exact choice of $a/R_{\ast}$.