THE NEAR-ULTRAVIOLET LUMINOSITY FUNCTION OF YOUNG, EARLY M-TYPE DWARF STARS

We cross-correlated the LG11 catalog of 8889 nearby (d ≲ 60 pc), bright (J < 10), K7–M5 stars with the final GALEX data release. We searched the GALEX All-sky Imaging Survey (AIS), which contains both NUV and FUV sources, but kept only NUV matches, as GALEX was much more sensitive in this bandpass (e.g., a preliminary cross-correlation by LG11 found five times more matches to NUV than FUV sources).

Our cross-correlation included a correction for stellar proper motion between the J2000 epoch of LG11 and the epochs of the individual GALEX AIS observations. We used the GalexView online tool¹³ to perform a preliminary cross-correlation between the LG11 and GALEX AIS catalogs with a 1 arcmin match radius. This search radius accounted for the resolution of GALEX (∼5 arcsec) plus the maximum proper motion of any LG11 star over 10 years (roughly the time between the J2000 epoch and the last GALEX AIS observation). This preliminary search returned multiple GALEX AIS source matches for most LG11 stars, as expected from the large search radius. For each LG11 star we calculated its expected position on the sky at the observation date of each of its GALEX AIS matches. We then reperformed the cross-correlation, using the adjusted LG11 positions and a reduced matching criterion of 5 arcsec (the resolution of GALEX) to identify the final source matches.

There were 1251 LG11 stars with multiple matches to the GALEX AIS catalog, even after correcting for proper motion and applying the stricter 5 arcsec search radius. For these we simply took the closest match. This was justified because almost all of these multiple matches were due to repeated GALEX AIS observations of a given area of sky: the matches had different GALEX AIS tile numbers and/or different exposure times and are thus likely the same star. However, there were 49 multiple matches that had the same GALEX tile number and exposure time, therefore representing the case when two objects are close enough on the sky to be confused by GALEX. The NUV magnitude differences for these 49 multiple matches were only ∼0.30 on average, and only eight of these sources were ultimately used in the NUVLF. Thus the choice of which source to match should not significantly affect our derived NUVLF. Figure 2 compares the cross-correlation results before and after proper motion correction: the proper motion correction significantly improved the cross-correlation results, adding ∼300 new matches compared to the number of matches prior to proper motion correction and significantly decreasing the residual angular separations of the matches, as illustrated by the histogram shift toward smaller offsets.

Zoom In Zoom Out Reset image size

This cross-correlation process identified 5267 LG11 stars with NUV counterparts in the GALEX AIS catalog. Figure 3 plots m_NUV − K_S, a distant independent measure of NUV luminosity, versus V − J, a proxy for stellar effective temperature, for this sample. m_NUV is from the GALEX AIS catalog; J and K_S are from the Two Micron All Sky Survey (Skrutskie et al. 2006); and V is from the AAVSO Photometric All-Sky Survey (Henden et al. 2012), the Tycho-2 and Hipparcos catalogs (Perryman & ESA 1997), or generated from USNO-B magnitudes (Lépine & Shara 2005).

Zoom In Zoom Out Reset image size

Most sources in Figure 3 fall along a locus (designated by gray points according to the criterion described below), which we interpret as the basal level of NUV emission for early M dwarfs. To describe this basal NUV emission as a function of stellar effective temperature, we fit a line to median values of m_NUV − K_S versus V − J, using only stars with errors <10% in all bands and V − J colors <3.5. We iteratively removed outlier stars and reperformed the fit until the remaining stars were Gaussian-distributed in m_NUV − K_S about the median-fit line. This gave a best-fit (designated by the black dashed line in Figure 3) with the following parameters:

$$\begin{equation} m_{\rm NUV}-K_{\rm S} = 7.72 + 1.66 (V-J) . \end{equation} \tag{ 1 }$$

This basal NUV locus is presumably an extension of the locus of inactive solar-type stars, which reaches m_NUV − K_S ∼ 10 at V − J ∼ 2 (see Figure 2 in Findeisen et al. 2011). However, our basal locus is much bluer than predicted by PHOENIX model spectra (Allard et al. 2013; Rajpurohit et al. 2013) using solar abundances (Caffau et al. 2011) and log g = 5.0. Thus this locus likely does not exclusively represent photospheric NUV emission. It also cannot be an artifact of a constant GALEX flux limit, as the lower right-hand region of Figure 3 is populated, whereas the lower left-hand region is not. Instead, this locus likely indicates persistent NUV line emission from a higher-temperature upper chromosphere. Stelzer et al. (2013) previously noted that all M dwarfs appear to exhibit NUV emission in excess of their expected photospheric value from stellar atmospheric models.

The distribution in Figure 3 also features a smaller population of stars with NUV emission significantly in excess of the empirically determined basal value. We interpret these NUV-luminous sources as being mostly young stars exhibiting heightened activity but also including FPs that appear active for reasons other than stellar youth. We identified the 1210 NUV-luminous stars from the 5267 NUV-detected stars as those with m_NUV − K_S colors at least 2.5σ (∼1.12 mag) bluer than their expected basal value given by Equation (1). We identify these stars as blue points in Figure 3.

We only consider early M dwarfs (M0–M3) in the remainder of this work. This is because we identify our sample of young M dwarfs by their heightened UV luminosity (see Section 2.3), a property that has been observed to evolve significantly with stellar age for early M dwarfs and at a much slower rate for late M dwarfs (M5–M9). Namely, observations indicate that early M dwarfs are rapidly rotating and active for only ∼1 Gyr, whereas late M dwarfs remain in this state for up to ∼8 Gyr (West et al. 2008; Rebassa-Mansergas et al. 2013). This difference may be due to the transition to fully convective interiors around the M4 spectral subtype (Reiners 2012; Gastine et al. 2013). We therefore removed stars with spectral subtypes later than M3 as a first step to ensuring our measured NUV excess is due to stellar youth rather than a longer duration of stellar activity.

To identify late M dwarfs we used the empirical relation between V − J color and spectral subtype derived in LG11 and then revised in Lépine et al. (2013; see their Equation (12)). These color-based assignments provide only rough estimates (±1 spectral subtype); however, they also provide a uniform method of spectral classification, which is important for statistical studies such as this. Moreover, a rough estimate is sufficient because we are using these spectral subtypes to identify late M dwarfs, and the spectral subtype boundary between early and late M dwarfs is not well defined (i.e., somewhere between M3 and M5).

Removing the late M dwarfs from our original sample of NUV-detected M dwarfs (see Section 2.1) resulted in a final sample of 4805 NUV-detected early M dwarfs, of which 797 were identified as NUV-luminous (see Section 2.3). All references hereafter to our "sample" are to this final selection that includes only early M dwarfs, unless explicitly stated otherwise. Parameters of the NUV-luminous sources are presented in Table 4.

Figure 3 suggests that the locus of stars exhibiting basal NUV emission (gray points; see Section 2.1) has an intrinsic width. We estimated this intrinsic width by first taking the standard deviation of distances in m_NUV − K_S from the locus median (black dashed line) to all stars in the locus with errors <10% in all bands. We accounted for measurement errors by subtracting in quadrature the median measurement errors for each band: the m_NUV errors were taken from the GALEX AIS catalog, the V-band errors were taken from the respective catalogs (see Section 2.1), and the K_S- and J-band errors were considered negligible at <1%. Errors in V − J were translated to m_NUV − K_S, using the slope of our median fit to the locus (see Equation (1)) before subtracting in quadrature. We found that the locus of stars exhibiting basal NUV emission in Figure 3 has an intrinsic width of ∼0.50 mag in m_NUV − K_S after accounting for measurement errors.

The locus width could be the product of several factors including stellar variability, interstellar extinction, unresolved binaries, and metallicity variations. We first investigated stellar variability using the 1202 LG11 stars in our sample with multiple matches to the GALEX AIS catalog that were likely repeated observations of the same star at different epochs rather than source confusion (see Section 2.1). We found the maximum difference in NUV magnitude for each multiple match and then took the standard deviation as an estimate of the error introduced by stellar variability. We then subtracted this value (∼0.24) in quadrature along with the measurement errors (see above), which reduced the estimated intrinsic locus width by only ∼0.07 mag. This shows that stellar variability is not a significant contributor to the locus width.

We also confirmed that the locus width was not due to interstellar extinction, A(V). UV sources are subject to significant interstellar extinction, as evidenced by the significant drop in GALEX detections near the Galactic plane (see Figure 1 in Bianchi et al. 2011). However, interstellar extinction should be negligible for our sample, as most LG11 stars reside within 60 pc and are therefore contained within a low-density (∼0.005 atoms cm⁻³) region known as the "Local Bubble" (Cox & Reynolds 1987). Nevertheless, we tested whether interstellar extinction could account for the locus scatter by searching for a minimum locus width as a function of assumed extinction per parsec, as described below. We found distances to each star using the J-band photometric distance, where M_J was obtained from V − J color (although LG11 computes photometric distances, we recompute them here using an updated color magnitude relation given by Equation (22) in Lépine et al. 2013). We used extinction coefficients from Yuan et al. (2013), which produced reddening corrections of E_{(NUV − K)} = 2.24 × A(V) and E_{(V − J)} = 0.77 × A(V). We tested A(V) values ranging from 0 to 0.001 mag pc⁻¹ at a cadence of 1.5 × 10⁻⁶ mag pc⁻¹. This encompassed values well beyond the expected interstellar extinction within 60 pc; assuming A(V) ∼ 1 mag kpc⁻¹ along the Galactic plane and conservatively assuming a Local Bubble that is 10% the typical density of the interstellar medium, the expected interstellar extinction within 60 pc is A(V) ∼ 0.0001 mag pc⁻¹. We applied reddening corrections to each star on the basis of their individual distances and then remeasured the locus width for each A(V) test value, finding no local minimum in the locus scatter. Rather, attempting to correct for extinction only increased the scatter of the locus.

We therefore concluded that stellar variability and interstellar extinction do not contribute significantly to the locus width. However, the locus width may be due to unresolved binaries, metallicity-dependent stellar colors, and continued variation of the basal NUV emission level with age. Unfortunately, our data set did not allow us to investigate these possibilities in detail.

We checked for FUV and X-ray counterparts to our sample. To obtain FUV counterparts we simply took the FUV sources associated with our NUV matches to the GALEX AIS catalog. To identify stars with X-ray counterparts, we cross-correlated our sample with the ROSAT All-Sky Survey Bright Source Catalog (Voges et al. 1999) and Faint Source Catalog (Voges et al. 2000). We used a 25 arcsec search radius around the LG11 coordinates, corresponding to the 2σ ROSAT positional uncertainty determined by Voges et al. (1999). We converted the PSPC detector count rate into an X-ray flux, F_X, using the conversion factor from Schmitt et al. (1995): CF = (5.30HR + 8.31)10⁻¹² erg cm⁻² count⁻¹, where HR is the first hardness ratio from the ROSAT catalog. We did not correct for proper motion for this cross-correlation, due to the large positional uncertainty of ROSAT compared to GALEX. However, we checked for mismatches with background galaxies and quasars by plotting F_X/$F_{K_{\rm S}}$ as a function of J − K_S but found no significant outliers, i.e., sources with F_X/$F_{K_{\rm S}} \ge 0.1$ (Kouzuma & Yamaoka 2010).

Only ∼8% of our sample (387 of 4805 sources) had detectable flux in all 3 wavelength bands (NUV, FUV, and X-ray). However, ∼85% of this multiwavelength subsample (328 of 387 sources) was also selected as NUV-luminous in Section 2.3, which means ∼40% of our NUV-luminous subsample (328 of 797 sources) was detected in all 3 bands. Figure 4 shows F_FUV/F_NUV versus F_X/F_NUV for the 387 multiwavelength sources. Sources with "hard" spectra are located in the upper right of the figure, whereas those with "soft" spectra are located in the lower left. Sources with high FUV but low X-ray emission could be M dwarfs with white dwarf companions (MD+WD pairs), as white dwarfs emit strongly in the FUV due to their hot photospheres but lack coronae from which X-ray emission typically originates. Figure 4 highlights a known WD+MD pair and several MD+WD candidates identified by their high F_NUV/F_J ratios (e.g., see Figure 3 in Shkolnik et al. 2011).

Zoom In Zoom Out Reset image size

We obtained medium-resolution (λ/Δλ ∼ 1000) optical spectra for 2128 out of the 4805 M dwarfs in our sample. The majority of these (1307 spectra) were acquired using the Super-Nova Integral Field Spectrograph (SNIFS; Aldering et al. 2002; Lantz et al. 2004) mounted on the University of Hawaii 2.2 m telescope atop Mauna Kea. SNIFS uses a dichroic mirror to separate incoming light into blue (3200–5200 Å) and red (5100–9700 Å) spectrographic channels. We only used spectra from the red channel, as M dwarfs have very low signal in the blue channel. All spectra had signal-to-noise ratios (${\rm {\rm S/N}}$s) ≳ 80 per resolution element in the red channel while avoiding the nonlinear regime of the detector. Details of our SNIFS data reduction method can be found in Mann et al. (2012) and Lépine et al. (2013). SNIFS is an integral field spectrograph and therefore also provides limited spatial information in the form of image cubes. A SNIFS image cube covers 6 × 6 arcsec at 0.4 arcsec per pixel.

The remaining 821 spectra were obtained using four instruments on three different telescopes: the Mark III spectrograph and the Boller and Chivens CCD spectrograph on the 1.3 m McGraw–Hill telescope at the MDM Observatory on Kitt Peak (564 spectra); the RC spectrograph on the 1.9 m Radcliffe telescope at the South African Astronomical Observatory (SAAO) in South Africa (67 spectra); and the REOSC spectrograph on the 2.15 m Jorge Sahade telescope at the Complejo Astronómico El Leoncito Observatory (CASLEO) in Argentina (190 spectra). Details of the data reduction methods for these spectra are in Gaidos et al. (2014b).

We observed 193 M dwarfs in our sample with the Robo-AO laser adaptive optics and imaging system (Baranec et al. 2013, 2014) mounted on the Palomar Observatory 1.5 m telescope. These observations were taken from 2013 August 13 to 2014 May 25 (UT). Robo-AO has a field of view of 44 × 44 arcsec at 43.10 mas per pixel. Typical point-spread function (PSF) widths achieved at red-visible wavelengths are in the range of 0.12–0.15 arcsec, and companions down to six magnitudes fainter than the primary can be detected (Law et al. 2014). We used a Sloan i'-band filter (York et al. 2000) for stars with V < 13 and a long-pass filter cutting on at 600 nm (hereafter LP600) for fainter stars. Observations consisted of a sequence of full-frame-transfer detector readouts of the electron-multiplying CCD camera at the maximum rate of 8.6 Hz for a total of 90 s of integration time. The individual images were corrected for detector bias and flat-fielding effects before being combined through post facto shift-and-add processing that used the target as the tip-tilt reference star with 100% frame selection.

We queried SIMBAD for our entire sample of NUV-luminous stars in order to identify any known FPs. We first searched for tight binaries, including spectroscopic binaries (SBs), eclipsing binaries (EBs), and RS Canum Venaticorum (RS CVn) binaries. We followed up these candidates in the literature to confirm that they had orbital periods <10 days, at which point tidal interactions between companions likely result in synchronized orbits and therefore enhanced, persistent NUV emission beyond stellar youth (e.g., see Meibom et al. 2006 and references therein). We then searched for close binaries with separations <5 arcsec (i.e., unresolved by GALEX) and white dwarf or late M secondary components, which also emit persistent NUV emission at older ages. We did not remove systems containing secondary M dwarf components with unknown spectral subtypes. We also searched for sources in our sample with background NUV objects within 5 arcsec and then inspected each of them individually to confirm that they were not double entries in the SIMBAD database.

We found seven EBs, seven RS CVn systems, and eight SBs in our sample. We also found 18 close binaries with late M or white dwarf components within 5 arcsec. There were nine systems with background sources that were the probable source of NUV emission, rather than the M dwarf. These background FPs included a variety of sources: one contact EB, two SBs, one RS CVn star, two late M dwarfs, one white dwarf, and two flare stars.

4.2.1. Robo-AO: Late M Companions

We searched our 193 Robo-AO images (see Section 3.2) for close binaries whose secondary components may be causing FP NUV emission. We first searched for binaries using a principal component analysis (to flag elongated sources) and a Gaussian source finder (to flag multiple sources). We followed up these candidate binaries with manual (by-eye) checks to confirm the existence of clearly resolved binaries. We then developed a method that utilizes binary contrast ratios to identify systems where the secondary component is likely to be a late M dwarf with elevated NUV emission (see Section 2.4).

For the principal component analysis, we calculated for each star an elongation factor, ε (the ratio between the longer principal axis and the shorter principal axis), and an orientation angle, θ (the angle between the longer principal axis and the vertical image axis; the vertical axis of Robo-AO images is 239 right of north). Only sources that were >10σ above the noise were considered, and the noise was calculated using an outlier-resistant estimate of the dispersion in an area of empty sky around the source. In theory, single stars should have uniformly distributed θ with ε = 1, whereas binary stars should also have uniformly distributed θ, but ε > 1 due to a companion skewing the otherwise symmetric distribution. However, Figure 5 illustrates that point sources in our data set tend to have slight elongation (ε ∼ 1.05) and positive tilt (θ ∼ 10°). To account for this, in the principal component analysis we ignored sources with these systematic effects (i.e., all stars inside the black box in Figure 5) and flagged all other sources as potential binaries due to their significant elongation. The Gaussian source finder was used to search the remaining parameter space (i.e., inside the black box in Figure 5) by looking for positive brightness perturbations that were >3σ above the noise. We flagged all images with multiple positive brightness perturbations as potential binaries. We then used by-eye checks on all candidate binaries to identify only the clearly resolved systems for further analysis. These binaries and their calculated separations are listed in Table 2.

Zoom In Zoom Out Reset image size

Table 2. Robo-AO Binaries^a

Star	Date	ρ	P.A.	Δm
	(UT)	(arcsec)	(deg)	(mag)
PM_I00234+2418	2012 Oct 9	2.15	255	3.1
PM_I00235+0947S	2013 Oct 25	3.75	140	0.3
PM_I00235+2014	2013 Oct 23	1.72	139	1.1
PM_I00505+2449S	2013 Oct 23	1.02	320	0.8
PM_I00531+4829	2013 Oct 24	1.21	52	1.6
PM_I00574+3736	2013 Jul 24	4.65	289	5.4
PM_I01133+5855	2013 Oct 25	2.00	214	3.4
PM_I01146+2057	2013 Jul 24	1.40	329	2.1
PM_I01376+1835	2013 Oct 23	1.59	23	0.1
PM_I01410+5308E	2013 Oct 24	3.95	78	1.2
PM_I01480+4652	2013 Oct 22	1.09	248	1.5
PM_I01491+0624	2013 Oct 25	0.69	358	1.0
PM_I02024+1034	2013 Oct 24	0.85	34	0.4
PM_I02208+3320	2013 Oct 23	1.45	337	1.9
PM_I02347+1251	2013 Oct 25	1.17	39	1.4
PM_I02408+4452	2013 Oct 22	5.92	342	5.0
PM_I02560+1220	2013 Oct 25	0.91	287	0.8
PM_I03053+2131	2013 Oct 23	0.63	297	0.6
PM_I04284+1741	2013 Oct 23	1.63	79	1.8
PM_I04310+3647	2013 Oct 22	0.77	215	0.2
PM_I04333+2359	2013 Aug 14	0.75	314	0.6
PM_I04453+1334	2013 Oct 25	5.58	169	4.1
PM_I04499+2341E	2013 Aug 14	2.39	92	0.8
PM_I04540+2200	2013 Oct 24	3.90	320	3.8
PM_I05228+2016	2013 Oct 24	3.57	90	1.8
PM_I05341+4732	2013 Aug 15	2.33	14	1.0
PM_I06088+4257	2013 Oct 24	1.24	352	0.8
PM_I06212+4414	2013 Oct 24	1.30	203	3.2
PM_I06268+4202	2013 Oct 22	0.86	175	1.2
PM_I11505+2903S	2013 Apr 22	0.47	347	0.6
PM_I17038+3211	2012 Aug 07	1.34	147	1.8
PM_I20514+3104	2013 Oct 25	1.44	135	1.2
PM_I21010+2615	2013 Oct 25	0.44	108	0.2
PM_I21221+2255	2013 Aug 14	5.31	306	4.8
PM_I21410+3504	2013 Jul 24	4.42	72	0.3
PM_I22006+2715	2013 Oct 24	5.30	356	5.1
PM_I22234+3227	2013 Oct 24	1.30	248	0.5
PM_I23045+4014	2013 Oct 25	0.84	334	0.6
PM_I23063+1236	2013 Oct 25	0.42	327	0.4
PM_I23300+1643	2013 Oct 25	0.95	147	0.2
PM_I23318+1956Wn	2013 Oct 25	5.29	81	1.8
PM_I23450+1458	2013 Oct 25	1.15	178	0.3
PM_I23535+1206S	2013 Oct 25	5.71	165	0.8
PM_I23578+3837	2012 Aug 28	0.50	243	1.6

Note. ^aRaw pixel positions were converted to on-sky separations (ρ) and position angles (P.A.) by applying a distortion solution derived from Robo-AO data. Typical uncertainties are conservatively estimated to be σ_ρ ∼ 0.05 arcsec, σ_PA ∼ 20, and σ_Δmag ∼ 0.1 mag.

Download table as:  ASCIITypeset image

Despite the high resolution of Robo-AO, even the closest binaries resolved by this instrument are too far apart (i.e., several AU) to be tidally locked. However, Robo-AO can easily resolve binaries with separations <5 arcsec and thus unresolved with GALEX. Therefore contrast ratios derived from Robo-AO images can be used to identify the binaries in our sample with late M secondary components likely causing FP activity. To identify such systems, we utilized the empirical relation between absolute Sloan i-band magnitude (M_i) and spectral subtype derived in Hawley et al. (2002). By calculating the difference between M_i at M3.5 and M_i at all other spectral subtypes, we derived the maximum allowable contrast ratio (Δm_lim) as a function of primary spectral subtype such that both components are early M dwarfs. We confirmed that Δm_lim was applicable to images taken through either Robo-AO filter (Sloan i-band or LP600) by measuring the contrast ratios of a binary in our sample that was imaged in both bands; the difference in the measured contrast ratios was only ∼0.1 mag.

For each Robo-AO binary, we measured the contrast ratio by performing aperture photometry on each companion using the IRAF phot routine and a 20 pixel circular aperture. This aperture size was based on when the curve of growth for a typical single star reached an asymptotic value. Closer binaries required smaller aperture radii (5–15 pixels) to avoid contamination from companions. For sky subtraction we used median-combined, manually sampled patches of nearby empty sky around each system to ensure proper sky measurements. Figure 6 shows measured contrast ratios for all of our Robo-AO binaries with separations ≲ 5 arcsec. The black dashed line indicates our calculated Δm_lim as a function of primary spectral subtype. We found 26 binaries with primary spectral subtypes ⩽M3 but secondary components likely to be >M3, making them potential FPs.

Zoom In Zoom Out Reset image size

The observational completeness of our Robo-AO FP search was limited by two factors: (1) the maximum contrast ratio that Robo-AO can detect and (2) the probability that a late M companion was unresolved because its projected distance from the primary at the time of observation was too small. To estimate the completeness due to (1), we used the Cruz et al. (2007) J-band luminosity function, n(M_J), to calculate the number of M dwarf companions expected to have spectral subtypes ⩾M3.5 but earlier than the limit imposed by the maximum contrast ratio that Robo-AO can detect, Δm_max. We then compared that number to the total number of M dwarf companions expected to have spectral subtypes ⩾M3.5:

$$\begin{equation} C =\frac{\sum _{i=1}^{{\rm binaries}} \int _{M_{ J, 3.5}}^{M_{ J, SpT_{i}}+\Delta m_{{\rm max}, { i}}} n(M_{ J})dM_{ J} }{\sum _{i=1}^{{\rm binaries}} \int _{M_{ J, 3.5}}^{M_{J, 9.0}} n(M_{ J})dM_{ J}}. \end{equation} \tag{ 2 }$$

We adopted Δm_max = 6.0, which is the maximum contrast ratio that Robo-AO can achieve at high "contrast performance" (i.e., PSF sizes >0.15 arcsec) for separations >1.0 arcsec in the Sloan i-band (see Figure 5 in Law et al. 2014). The median PSF size and separation for the binaries in our Robo-AO sample were ∼0.17 arcsec and ∼1.41 arcsec, respectively, justifying this choice of Δm_max. We used relations from Hawley et al. (2002) to obtain i − J color as a function of M_J, which allowed us to translate our Robo-AO contrast ratios from i to J band for use in Equation (2). We calculated each primary's M_J, using its V − J color according to Equation (22) from Lépine et al. (2013). Evaluating Equation (2) gave a completeness factor of ∼0.97.

To estimate completeness due to (2), we compared the median PSF size in our Robo-AO images (∼0.17 arcsec) to the GALEX PSF size (∼5.0 arcsec) and assumed a uniform source distribution on the sky (equivalent to an isotropic, uniform distribution of orbits with a log orbital period; Duquennoy & Mayor 1991). This gave a ∼3% probability of an unresolved companion and thus a ∼0.97 completeness factor. The product of the completenesses from (1) and (2) gave a final observational completeness of C ∼ 0.94 for our Robo-AO FP search.

4.2.2. Optical Spectra: Missing Hα Emission

Emission in the Balmer α line of hydrogen (Hα; λ ≈ 6563 Å) is another indicator of stellar activity and is strongly correlated with NUV emission (e.g., Lépine et al. 2013). FPs resulting from GALEX source confusion can therefore potentially be identified by the absence of Hα emission despite significant amounts of NUV emission (i.e., NUV emission at least 2.5σ above the basal level; see Section 2.3). This is because sources that are unresolved by GALEX's ∼5 arcsec beam can be resolved by spectroscopy with ∼1 arcsec resolution. We therefore computed Hα equivalent widths (EW_Hα) for the 2128 stars in our sample with medium-resolution optical spectra (see Section 3.1).

To calculate EW_Hα we first shifted each spectrum to its rest frame by applying wavelength offsets found by matching our observed spectra to PHOENIX model atmospheres. We used the BT-SETTL version of the PHOENIX atmospheric model code (Allard et al. 2013; Rajpurohit et al. 2013) with the CIFIST grid and Caffau et al. (2011) abundances for the Sun. Following Lépine et al. (2013), we measured the flux within a 14 Å-wide spectral region (6557.61–6571.61 Å, in air) relative to pseudo-continuum regions (6500–6550 Å and 6575–6625 Å, in air). Errors were calculated using a Monte Carlo method that assumed Gaussian-distributed noise and random wavelength calibration errors of 0.5 Å. As noted by Lépine et al. (2013), this choice of continuum region systematically underestimates EW_Hα values due to differences between our pseudo-continuum and the true spectral continuum. We therefore applied a small offset (0.3 Å) so that stars with basal NUV emission had a mean EW_Hα ≈ 0. The correction did not affect our FP analysis because the offset was applied to the entire population.

Figure 7 shows our measured EW_Hα values as a function of our selection cutoff for NUV-luminous stars (see Section 2.3). As expected, the stars with basal NUV emission (gray points, corresponding to those in Figure 3) are located in a band with negligible EW_Hα. The NUV-luminous stars (blue points, corresponding to those in Figure 3) form a distinct locus where increasing NUV emission corresponds to increasing EW_Hα values. However, a subset of NUV-luminous stars have lower-than-expected EW_Hα values (i.e., they lie significantly below the NUV-luminous locus) and are presumably FPs due to GALEX source confusion. We fit a line to median values of EW_Hα versus relative NUV emission for the NUV-luminous population (black dashed line in Figure 7) to determine the expected EW_Hα value for a young M dwarf with a given NUV emission level. We identified 37 stars with EW_Hα values more than 3σ below the expected value, making them likely FPs.

Zoom In Zoom Out Reset image size

Interestingly, all known WD+MD binaries in our sample, as well as the new candidate WD+MD systems identified in this work (see Section 2.6), were flagged as FPs using this detection method. As shown in Figure 7, these WD+MD pairs have anomalously high NUV emission relative to their EW_Hα values, as expected (see Section 2.6). It is also important to note the distinct lack of stars with high EW_Hα values but basal NUV emission levels (upper left quadrant in Figure 7); because the NUV and Hα measurements were taken at different epochs, this may illustrate consistent levels of activity for the vast majority of stars in our sample. This would suggest that our sample contains very few flare stars and thus a low probability of misidentified FPs due to stellar variability. Of course, the same reasoning could be applied to the significant population of sources with high NUV emission but low EW_Hα values to argue for evidence of flaring; however, as discussed above, this population also contains FPs, which complicates the interpretation.

We estimated the observational completeness for this FP detection method using an injection and recovery method. We replaced the region of Hα emission in each star's spectrum with the median of the surrounding continuum flux and then injected a synthetic Hα signal with an EW equal to the expected value for that star's NUV emission level (using the black dashed line in Figure 7 as a guide). We also added random Gaussian noise scaled to the noise in the surrounding continuum regions. We remeasured the EW_Hα values, using the same method as above, repeating 100 times and taking the average of the resulting EWs for each star. We then checked whether this EW_Hα value was within 3σ of the expected value. In all cases the signal was recovered, implying C ∼ 1.0 for our EW_Hα FP detection method.

4.2.3. SNIFS Integral Field Spectra: Hα-emitting Companions

SNIFS image cubes provide both spatial and wavelength dimensions that can be used in conjunction to search for FPs. In particular, they can be used to find binary systems appearing young due to unresolved late M companions with persistent activity despite being old. In such cases, the early M primary exhibits basal NUV emission but dominates the continuum signal (making the system appear as a single early M dwarf), whereas the unresolved late M companion is the source of the stronger Hα emission and presumably also the NUV flux. This configuration is detectable as a shift between the centroid of a white-light image versus the centroid of an Hα image. We used this method of detecting FPs in addition to our Robo-AO image analysis (Section 4.2.1) because Hα traces activity and we have more SNIFS image cubes than Robo-AO images.

We performed this analysis on the 242 stars in our NUV-luminous sample with SNIFS image cubes. To create the white-light image we summed the SNIFS image cube over all wavelengths covered by the spectrum. We identified the source centroid location by employing a principal component analysis, using only points that were >10σ above the noise. The noise was calculated using an outlier-resistant estimate of the dispersion around the outer edge of the image; we used the outer edge due to the small size of SNIFS images (14 × 14 pixels). To create the Hα image, we summed the SNIFS image cube across only wavelengths covering the Hα spectral line (see Section 4.2.2 for our Hα line parameters), then subtracted the continuum and divided by the noise. The continuum was estimated using the median value of surrounding wavelength regions multiplied by the number of wavelength elements covered by the spectral line, and the noise was found by taking the standard deviation of the continuum regions. To identify the Hα centroid location, we again used a principal component analysis but first applied a mask to consider only pixels that were also used to calculate the white-light centroid. We discarded systems lacking any significant Hα emission.

Figure 8 shows the distribution of offsets between the white-light and Hα centroids. We fit a Rayleigh function to the distribution using the IDL routine mpfit (Markwardt 2009), which returned a mean offset of μ ≈ 0.37 pixels and a dispersion of σ ≈ 0.20 pixels. We flagged systems with centroid offsets >3σ above the mean as potential FPs, then used by-eye checks to confirm 25 systems with clear shifts in their image centroids. We estimated the observational completeness as the ratio of the area over which FPs could be detected (i.e., the area of the annulus from r = 3σ to r = r_max, where r_max is the maximum centroid offset in our sample) to the total survey area (i.e., the area of a circle with r = r_max). This resulted in an observational completeness of C ∼ 0.96 for our centroid offset detection method.

Zoom In Zoom Out Reset image size

4.2.4. SuperWASP Light Curves: Tidally Locked Binaries

Time-series photometry from the SuperWASP exoplanet survey (Pollacco et al. 2006) can be repurposed to identify very short period, interacting stellar binaries (e.g., Norton et al. 2011). These systems of tidally locked, synchronously rotating stars remain magnetically active and rapidly rotating due to the transfer of angular momentum from their orbits to their spins. They exhibit photometric variability because of transits and/or fixed patterns of spots established by the interacting magnetic fields of the companions. Even though a young, single M dwarf can also have a surplus of star spots due to elevated activity, these spots tend to be uniformly distributed across the surface, which likely dampens any induced light curve variability (see Barnes et al. 2011, and references therein). Moreover, spots on single stars tend to migrate, causing the phases of their light curve signals to change over time.

We therefore cross-referenced our NUV-luminous sample with the SuperWASP database, which is available online.¹⁴ We used a 3 arcsec search radius and only considered light curves with more than 1000 data points due to the limited photometric precision of SuperWASP. This resulted in a sample of 312 NUV-luminous sources with SuperWASP light curves. We inspected each of these light curves for stellar variability by computing their Lomb–Scargle periodogram (Scargle 1987), which estimates a frequency spectrum based on a least-squares fit to sinusoids. We only considered signals with a false-alarm probability of <0.1% and ignored periods within 2% of 1 day and fractions thereof (1/4, 1/3, 1/5, etc.), as those are likely artifacts due to observing schedules (Gaidos et al. 2014a). We used a stricter 10% filter around periods with stronger systematics (0.5, 1.0, 1.5, and 2.0 days).

To identify candidate tidally interacting binaries, we selected light curves that featured large amplitude variations (>2.6%) at short periods (<10 days) with phases that appeared perfectly Keplerian (i.e., stable over many cycles). The period criterion was based on the orbital period at which tidal interactions between companions begin to synchronize orbits and therefore enhance stellar activity (see Meibom et al. 2006, and references therein). The amplitude criterion was obtained from running the Lomb–Scargle periodogram on all known EBs in LG11 with SuperWASP light curves, then taking the minimum amplitude recovered as our lower limit. Figure 9 shows this limit in the context of all SuperWASP light curve amplitudes derived from our sample. From the preliminary list of candidate FPs generated by applying the above criteria, we used by-eye checks to identify 15 light curves with clear sinusoidal signals consistent with those of known tidally locked binaries. We then divided each of these light curves into halves and reran the analysis on both sections to check that the period remained unchanged, indicating variability due to regular eclipses rather than varying star spot patterns.

Zoom In Zoom Out Reset image size

We estimated observational completeness using a method of injection and recovery of artificial sinusoidal signals. We randomly selected 100 SuperWASP light curves from our sample and randomized the fluxes of the data points. We then injected signals with randomly selected periods of 0.1–10 days (i.e., our FP period search criteria) and amplitudes of 2.6–6.7% (i.e., our FP amplitude search criteria limited by the maximum amplitude found in our sample). We then reperformed our Lomb–Scargle periodogram search using the same period filters and false-alarm probability as before to test whether we would have recovered the injected signals as candidate FPs. This produced an observational completeness of C ∼ 0.83 for our SuperWASP FP detection method.

Construction of an accurate NUVLF requires the identification and consideration of FPs, i.e., systems that appear NUV-luminous for reasons other than stellar youth. Our approach was to (1) estimate the overall FP rate among our NUV-luminous sample based on FP detection methods for which we could determine observational completenesses (Section 4.2); (2) remove FPs identified using these FP detection methods as well as known FPs from the literature (Section 4.1); and then (3) use our derived overall FP rate to statistically correct the remaining NUV-luminous stars in our sample not yet flagged as FPs (Section 5.2). The last step was necessary to account for the fact that not all NUV-luminous sources were observed with all FP detection methods and also because our FP detection methods cannot each detect all possible types of FPs (e.g., the Robo-AO method cannot resolve tidally locked binaries).

We used a maximum likelihood approach to estimate the overall FP rate in our NUV-luminous sample. When screening a source with our FP detection methods in Section 4.2, the source is either identified as an FP, or it is not. We can describe the likelihood of these two possible outcomes, f(p) and g(p), respectively, in terms of the overall FP rate (p) and the completenesses of the different FP detection methods that were applied to the source (C_j):

$$\begin{equation} g(p) = (1-p) + p \prod _{j}^{{\rm missed}} (1-C_{j}) \end{equation} \tag{ 3 }$$

and

$$\begin{equation} f(p) = p \left[1-\prod _{j}^{{\rm missed}}(1-C_{j})\right ]. \end{equation} \tag{ 4 }$$

Equation (3) is the probability that a source was not identified as an FP: the first term is the probability that the source was not an FP, whereas the second term is the probability that the source was an FP but that it was missed by all the FP detection methods that were applied to it. Equation (4) is the probability that a source was identified as an FP: the term in the brackets is the probability that the source was identified as an FP by the detection methods that were applied to it, and this term is then multiplied by the actual FP rate. This approach assumes that there are no "false negatives" (i.e., it assumes that our FP detection methods cannot wrongly classify a source as an FP). It also assumes that a source is an FP if just one of our methods detected it as an FP.

We found the value of p that maximizes the likelihood:

$$\begin{equation} \ln {L} = N_{\rm FP} \ln p + \sum _{k}^{\rm NFP} \ln \left[ 1- p \left(1 - \prod _{m} (1-C_{m}) \right) \right], \end{equation} \tag{ 5 }$$

where N_FP is the total number of FPs found using our FP detection methods and C_m is the observational completeness of each FP detection method applied to a given star that was not found to be an FP (NFP). Solutions to Equation (5) for all possible values of p are shown in Figure 10, indicating a most likely FP rate of p ∼ 0.16 ± 0.02. We estimated the errors by fitting an inverted parabola around the peak in the log likelihood, then using the uncertainty on the curvature of the fit (i.e., $\sigma _{x}=1/\sqrt{-2c}$ for a parabola described by y = a + bx + cx²).

Zoom In Zoom Out Reset image size

We defined a fractional NUV luminosity to describe the distribution of NUV luminosities with respect to a basal value:

$$\begin{equation} R^{\prime }_{\rm NUV} = \frac{L_{\rm NUV}-L_{\rm basal}}{L_{\rm bol}}, \end{equation} \tag{ 6 }$$

where L_NUV is the total NUV luminosity, L_basal is the basal NUV luminosity, and L_bol is the bolometric luminosity. To obtain these parameters from the observed quantities shown in Figure 3, we rearranged the color relation $m_{\rm NUV}-K_{\rm S}=M_{\rm NUV}-{M}_{K_{\rm S}}$ to obtain an expression for absolute NUV magnitude: ${M}_{\rm NUV}={M}_{K_{\rm S}}+(m_{\rm NUV}-K_{\rm S})$. We calculated $M_{K_{\rm S}}$ by translating V − J color into M_J, using Equation (22) from Lépine et al. (2013), then converting M_J to $M_{K_{\rm S}}$, using J − K_S = 0.854 (J − K_S color varies little for early M dwarfs; here we used the median of our sample). To translate M_NUV into L_NUV, we used the GALEX zero points¹⁵ to convert magnitude into flux density (erg⁻¹ s⁻¹ cm⁻² Å⁻¹), then multiplied by the effective GALEX NUV bandwidth (Δλ = 732 Å) to get the NUV flux, F_NUV. We then input F_NUV into L_NUV = 4πd²F_NUV, where d = 10 pc (in accordance with the definition of absolute magnitude). We also applied this method to the median fit shown in Figure 3 (black dashed line) to obtain L_base as a function of V − J color. To calculate L_bol, we used the standard equation $L_{\rm bol}={L}_{\odot }\times 10^{-0.4({M}_{{\rm bol}}-{M}_{{\rm bol}, \odot })}$, where M_{bol, ☉} = 4.7554 mag and L_☉ = 3.8270 × 10³³ erg⁻¹ s⁻¹ (Mamajek 2012). To obtain M_bol, we used the K_S-band bolometric correction from Leggett et al. (2001).

The "1/V_max" method (Schmidt 1968) is used to construct luminosity functions by accounting for the bias of flux-limited surveys toward intrinsically bright sources. This is done by inversely weighting sources by the volume of space over which they could have been detected by the survey. Thus bright sources are assigned smaller weights to correct for being detectable to larger distances. The luminosity function is then calculated by summing the weights (rather than the number of stars) in each luminosity bin, giving units of stars pc⁻³.

To apply the 1/V_max method to our sample, for each star we calculated the maximum distance at which it would have been included in LG11 and also detected by GALEX. Thus the limiting detection distance of each star was determined by one of three factors: (1) the J-band magnitude limit (J < 10) of the LG11 catalog; (2) the proper motion limit (≳ 40 mas in the north and ≳ 100 mas in the south) of the LG11 catalog; and (3) the GALEX sensitivity limit. Because GALEX sensitivity varies across the sky due to varying tile exposure times, we estimated the limiting NUV magnitude for each star in our sample (see Section 5.3) to create a map of limiting NUV magnitude across the sky. Then for each star we compared its two LG11 limiting detection distances to each GALEX limiting detection distance across the sky, recording the smallest value in each case. The average of the cube of these smallest detection distances was then used to calculate V_max for that star. The star's contribution to the NUVLF was then determined by its weight, 1/V_max.

To account for FPs, we first removed all FPs found in the literature (Section 4.1) or with our FP detection methods (Section 4.2). We then statistically accounted for FPs in the remaining stars by multiplying their 1/V_max weights by 1 − p (i.e., the probability of not being an FP). For NUV-luminous stars we used p = 0.16 (see Section 4.3), and for all other stars we used p = 0. There were also 92 sources in our sample with limiting detection distances that were found to be smaller than their actual distance. This was mostly due to anomalies in survey sensitivities (e.g., 67 of these sources had declinations < − 20°, where LG11 proper motion limits can vary). We removed these sources from our sample before constructing the NUVLF; however, only nine of these sources were members of our NUV-luminous sample (these are flagged in Table 4).

The resulting NUVLF is shown in the main panel of Figure 11. By construction, the peak at $R^{\prime }_{\rm NUV}\approx 0$ consists of stars with basal levels of NUV emission (designated by gray points in Figure 3), whereas the extended tail toward higher $R^{\prime }_{\rm NUV}$ values consists of the NUV-luminous stars (designated by blue points in Figure 3). The negative values in the distribution result from the subtraction of a median-fit basal level when calculating $R^{\prime }_{\rm NUV}$ (see Equation (6)). We applied a small offset ($R^{\prime }_{\rm NUV}\sim 1.5\times 10^{-8}$) to shift the distribution peak to zero. This was essentially a correction to our median basal fit (black dashed line in Figure 3), which was likely skewed to higher $R^{\prime }_{\rm NUV}$ values because (1) we were not considering upper limits and (2) we were only considering stars with the lowest NUV flux errors, which typically have the highest NUV fluxes.

Zoom In Zoom Out Reset image size

The shape of the peak in the main panel of Figure 11 is dictated by photometry errors and the intrinsic width of the basal NUV locus (see Section 2.5). Thus we needed to extract the extended tail of NUV-luminous stars for our analysis of the NUVLF of young early M dwarfs. We did this by reflecting the negative side of the distribution about the ordinate and subtracting it from the positive side of the distribution. This essentially removed the population of stars with basal NUV emission from the distribution. This assumes that the distribution of $R^{\prime }_{\rm NUV}$ values in the basal population is symmetric about zero, but does not assume any distribution in particular. The median GALEX counts for basal NUV sources (∼50 counts) were sufficiently high that the distribution due to GALEX Poisson errors should be fairly symmetric. This extraction process gave a residual distribution that represents the NUVLF of young early M dwarfs, shown as the blue dashed line in the inset of Figure 11 and tabulated in Table 3. Errors were found using a standard bootstrap method. We sampled with replacement from the original set of $R^{\prime }_{\rm NUV}$ values until we obtained a bootstrap sample that contained the same number of data points as the original sample. We then reconstructed the NUVLF by following the same steps as above but instead using the bootstrap sample. We repeated this 100 times and used the standard deviation in each bin as an estimate of our errors.

Table 3. Near-ultraviolet Luminosity Functions

log (R')	$\rho _{V_{{\rm max}}}$a	ρAvnib
	(10−5 stars pc−3 dex−1)	(10−5 stars pc−3 dex−1)
−7.70	−29.6 ± 37.4	−1.5 ± 17.3
−7.50	36.1 ± 41.6	14.8 ± 28.4
−7.30	55.1 ± 52.0	94.9 ± 35.9
−7.10	77.3 ± 78.7	86.6 ± 40.3
−6.90	91.6 ± 74.9	134.6 ± 38.6
−6.70	200.1 ± 65.6	144.2 ± 42.6
−6.50	155.2 ± 58.8	104.9 ± 39.0
−6.30	163.9 ± 52.1	76.3 ± 33.7
−6.10	154.4 ± 37.5	126.6 ± 24.3
−5.90	208.3 ± 33.3	205.1 ± 19.1
−5.70	176.1 ± 27.5	176.1 ± 16.5
−5.50	57.9 ± 21.4	57.9 ± 10.3
−5.30	4.2 ± 5.4	4.2 ± 4.3

Notes. ^aNUVLF from Section 5.2 (no upper limits). ^bNUVLF from Section 5.3 (with upper limits).

Download table as:  ASCIITypeset image

Most of the LG11 stars without matches in the GALEX AIS catalog (2226 out of 3622) are either late M dwarfs (and therefore not considered in this study) or stars within 20° of the Galactic plane (where GALEX AIS coverage is sparse; Bianchi et al. 2014). To determine which of the remaining 1396 LG11 stars were true nondetections (i.e., stars located in an area of sky observed by GALEX AIS but too faint to be detected), we searched for GALEX AIS tiles covering their coordinates, using the aforementioned GalexView online tool (see Section 2.1). We only considered GALEX AIS tiles with centers less than 05 from the LG11 coordinates of the candidate nondetections. This separation limit ensured that the stars would have actually been located on the tile but not on the tile edge, where photometry can be significantly degraded (Bianchi et al. 2014). This search resulted in the identification of 638 nondetections, for which we estimated upper limits to incorporate into our NUVLF.

To calculate upper limits, we derived a relation between GALEX exposure time (t_NUV) and limiting NUV magnitude (m_lim), using the online GALEX Exposure Time Calculator.¹⁶ We queried the tool for exposure times, given various m_NUV values, and required S/N ∼ 2. We used a hypothetical star with T_eff = 3500 K and coordinates α = 240° and δ = −11° (i.e., an early M dwarf with median LG11 declination and located away from the Galactic plane). To check the relation, we plotted m_NUV versus t_NUV for the 5267 GALEX-detected LG11 stars; as expected, the derived relation between t_NUV and m_lim corresponded to the lower (NUV-dim) bound of the NUV-detected population. We then obtained m_lim values for the 638 nondetections by using the t_NUV values of their associated GALEX AIS tiles with our derived relation.

We considered these upper limits in our NUVLF, using the method of Avni et al. (1980). This method employs a nonparametric, recursive approach to statistically account for upper limits when constructing a luminosity function (see Equation (6) in Avni et al. 1980). As before, we replaced the number of stars in each $R^{\prime }_{\rm NUV}$ bin with the sum of their 1/V_max weights in order to account for survey biases toward brighter stars. However, because we were considering GALEX upper limits with this method, we relaxed any GALEX constraints on V_max by only considering the detection distance limits imposed by LG11 when calculating V_max. We again had to apply a small offset ($R^{\prime }_{\rm NUV}\sim 6.5\times 10^{-9}$) to shift the distribution peak to zero; this correction was much smaller than before, most likely because we are now taking into account upper limits. The resulting NUVLF is shown as the red line in the inset of Figure 11 and tabulated in Table 3 with errors calculated using the same bootstrap method described in Section 5.2.

The two principal uncertainties in our derivation of the NUVLF are (1) the overall FP rate, p, which we used to statistically correct for FPs when constructing the NUVLF (Section 4.3); and (2) the distances used to compute the limiting detection volumes, V_max, which we used to weight each star's contribution to the NUVLF (Section 5.2). We address these two issues below.

Our derived overall FP rate of p ∼ 16% is consistent with the ∼16% SB rate among nearby X-ray luminous M dwarfs found by Shkolnik et al. (2009). One might expect our FP rate to be higher, as our definition of an FP encompasses additional, wider binaries. However Shkolnik et al. (2009) selected their sample based on X-ray fluxes from the ROSAT All Sky Survey, which was less sensitive to active stars than the GALEX AIS. Thus their sample was more biased toward the most active objects, which likely have higher FP rates (see below). Agreement between our FP rate and that of Shkolnik et al. (2009) may therefore simply be a coincidence resulting from several factors.

The NUV luminosities in our sample span several dex (see Figure 11). Thus a single, overall FP rate may be insufficient to describe the entire population, as the most NUV-luminous sources may have significantly higher FP rates. We tested this by progressively removing the dimmest NUV sources from our sample, then repeating the maximum likelihood estimation in Section 4.3 to rederive p using the cropped sample. Results are given in Figure 12, which shows a constant FP rate of p ∼ 16% (i.e., our overall FP rate) until m_NUV − K_S ∼ 11. At higher NUV luminosities, the FP rate steadily increases, reaching p ∼ 80% by m_NUV − K_S ∼ 10. To test the implications of this result, we rederived our NUVLF, using this varying FP rate (instead of the constant p ∼ 16% value) when multiplying the 1/V_max weights by 1 − p to statistically account for FPs (see Section 5.2). We did not expect significant changes to the NUVLF, as there are few stars in our sample with high enough NUV luminosities to require FP rates that are significantly larger than the overall FP rate. The results are compared to the original NUVLF in the inset of Figure 11. The consistency between the two NUVLFs suggests that our derivation is not particularly sensitive to this varying FP rate. Still, our use of a single FP rate is a gross simplification of reality, where there are multiple, sometimes unrelated sources of FPs, each of which cannot be detected by all our methods applied in Section 4.2. A more rigorous approach would use a modified version of Equation (5) to estimate multiple p values, one for each FP source, and also account for the completenesses of each method for each FP source. However, the outcome of the sensitivity analysis described above suggests that our results would not change significantly.

Zoom In Zoom Out Reset image size

Another source of uncertainty in our derivation of the NUVLF is our estimation of stellar distances, which we used to calculate V_max and thus the weighted contribution of each star to the NUVLF. For most sources we used a J-band photometric distance, where M_J was estimated from V − J color using Equation (22) in Lépine et al. (2013). However, we substituted more accurate parallax distances for the 561 stars in our sample that also had trigonometric parallax measurement with errors <10%. We used the average difference between these photometric and parallax distances to estimate fractional errors as a function of stellar distance. These ranged from ∼50% for distances of ∼3 pc (the minimum distance in our sample) to ∼10% for distances of ∼40 pc (the distance containing 95% of our sample). We then tested how these distance uncertainties impacted our derived NUVLF, using a Monte Carlo approach. We perturbed each distance by a random Gaussian deviate scaled to our estimated fractional errors, then reran our derivation of the NUVLF. We repeated this 100 times and took the standard deviation in each luminosity bin as an estimate of the impact of our distance uncertainties on our NUVLF. The impact appeared to be negligible: the variation was only ∼14% of the original NUVLF in each bin, well within the errors estimated by bootstrap sampling (see Section 5.2).

Relations between stellar age and emission at high-energy wavelengths are typically expressed as two-parameter power laws of the form F_λ = αt^β (c.f., Ribas et al. 2005), where t is the age of the star and the two parameters, α and β, are the zero point and slope of the power law, respectively. Stelzer et al. (2013) derived an age-activity relation of this form for early M dwarfs by combining their sample of 159 nearby field M dwarfs (assuming an age of ∼3 Gyr) with an additional sample of young (∼1 Myr) M dwarfs from the TW Hya association. They found that β = −0.84 ± 0.08 at GALEX NUV wavelengths.

Because stellar ages for our sample are mostly unknown, we first attempted to infer an age-activity relation by fitting our observed NUVLF to model NUVLFs constructed from assumed power-law age-activity relations. To create model populations, we assumed a constant star-formation rate (and thus a uniform age distribution) with a maximum age of ∼10 Gyr (i.e., the approximate age of the Galactic disk at the present solar radius; Bergemann et al. 2014). We then used a two-parameter power law, described above, to assign $R^{\prime }_{\rm NUV}$ values to each synthetic star based on its model age. We created the model NUVLF by binning the synthetic population according to the same $R^{\prime }_{\rm NUV}$ bins as in the inset of Figure 11, then normalizing the model distribution such that the integral under the model function equaled that of the real function. We searched for best-fit parameters by minimizing reduced χ². We calculated reduced χ² by taking the difference between the model and observed NUVLFs at each $R^{\prime }_{\rm NUV}$ bin, applying the errors shown in the inset of Figure 11 to the observations, and then dividing by the number of bins minus the number of power-law parameters. We found best-fit parameters of β ≈ −1.29 and α ≈ 2.4 × 10⁻⁶ with reduced χ² ≈ 6.3. The best-fit model is compared to the observed NUVLF in the inset of Figure 11. Clearly, this simple model is unable to account for the observed NUVLF. The discrepancy may be due to two key assumptions in our simplified model: (1) our neglect of a constant or "saturated" level of NUV emission at very young ages and/or (2) our assumption of a constant star-formation rate. We discuss (1) below and (2) in Section 6.3.

Saturated (i.e., constant) NUV emission in young early M dwarfs was recently reported by Shkolnik & Barman (2014). They used early M members of nearby YMGs to derive an NUV age-activity relation that showed NUV emission remaining constant for ≈300 Myr, then declining as a power law with β = −0.84 ± 0.09. This power-law index agrees with the previous findings of Stelzer et al. (2013), who did not consider saturated NUV emission due to the lack of stars younger than a few hundred Myr in their sample. We therefore used the YMG members in our sample to empirically derive a power-law age-activity relation that included a saturation component. We first searched the literature to identify 32 candidate YMG members in our sample and then required high (⩾95%) membership probabilities, which we obtained using the Bayesian Analysis for Nearby Young AssociatioNs; (BANYAN; Malo et al. 2013) online tool.¹⁷ For the Hyades, which is not included in BANYAN, we used the approach of Shkolnik & Barman (2014) by requiring a kinematic link to the YMG. We also removed any FPs identified in Section 4. This resulted in a final set of 20 YMG members, which are flagged in Table 4.

Table 4. Parameters of UV-luminous early M Dwarfs

Star	V	J	mNUV − KS	FNUV	FFUV	FX	EWHα	$\log (R^{\prime }_{\rm NUV})$	Robo-AOa	Low Hαb	Shifted Hαc	SuperWASPd	FP?e
				(μJy)	(μJy)	(erg cm−2 s−1)	(Å)
PM_I00001+6943	13.52	9.70	10.95	44.04	99.99	2.95E-13	99.99	−5.67	–	–	–	–	N
PM_I00024-4601	12.43	9.18	11.98	27.05	5.78	99.99	99.99	−6.22	–	–	–	No	N
PM_I00059+4129	12.96	9.40	11.56	33.08	10.22	99.99	1.14	−7.48	No	No	No	Yes	D
PM_I00072-1036	11.85	8.83	11.29	70.18	11.78	4.48E-13	0.80	−5.89	–	No	No	–	N
PM_I00107-2039	13.39	9.48	11.18	43.77	10.82	9.27E-13	99.99	−5.77	–	–	–	–	N
PM_I00117-1139	12.67	9.85	11.05	34.89	99.99	99.99	−0.04	−5.83	–	No	–	–	N
PM_I00118-5521	13.34	9.32	11.50	36.21	18.78	99.99	3.18	−5.86	–	No	–	–	N
PM_I00162+1951W	11.79	7.88	11.64	117.08	22.29	2.39E-12	6.16	−5.90	–	No	–	–	N
PM_I00166+3000	12.79	9.11	12.42	20.21	3.56	1.02E-13	99.99	−6.38	–	–	–	–	N
PM_I00197-2233	13.77	9.88	12.44	9.59	99.99	99.99	99.99	−6.33	–	–	–	–	N
PM_I00206-5340	13.80	9.68	13.30	5.18	99.99	99.99	99.99	−6.71	–	–	–	–	N
PM_I00211+4456	13.62	9.60	10.89	52.93	99.99	1.07E-12	5.20	−5.67	No	No	No	–	N
PM_I00216-4605	12.21	8.32	11.14	133.11	28.10	7.46E-12	1.94	−5.76	–	No	–	No	N
PM_I00234+2418	13.03	9.75	11.01	40.74	99.99	99.99	−0.31	−7.48	Yes	No	–	No	D
PM_I00235+0947S	12.75	9.79	10.59	56.81	10.92	99.99	−0.77	−7.48	No	No	Yes	–	D
PM_I00235+2014	11.00	8.14	10.84	193.75	99.99	1.35E-12	0.30	−7.48	No	No	Yes	No	D
PM_I00241-6211	11.33	8.38	11.00	144.75	32.66	1.04E-12	0.94	−7.48	–	No	–	–	L
PM_I00245-2522	13.62	9.84	11.29	29.08	5.05	3.07E-13	2.68	−5.84	–	No	No	No	N
PM_I00250-3646	12.46	8.64	12.00	42.38	12.20	9.64E-13	1.12	−7.48	–	No	–	No	L

Notes. ^aYes = found as FP by Robo-AO FP detection technique (Section 4.2.1); No = not found as FP by Robo-AO FP detection technique; – = not observed by Robo-AO. ^bSame as (a) but for missing Hα FP detection technique (Section 4.2.2). ^cSame as (a) but for shifted Hα FP detection technique (Section 4.2.3). ^dSame as (a) but for SuperWASP FP detection technique (Section 4.2.4). ^eL = found as FP in literature (Section 4.1); D = found as FP, using detection techniques (Section 4.2); N = not found as FP in literature or by FP detection techniques; R = removed from NUVLF (Section 5.2);\break; Y = YMG member used in age-activity relation derivation (Section 6.2).

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

In order to extend our age-activity relation to older field stars, as well as make our derived age-activity relation comparable to previous works, we redefined our fractional NUV luminosity ($R^{\prime }_{\rm NUV}$; see Equation (6)). Instead of subtracting the observed basal NUV luminosity (L_basal), we subtracted the photospheric NUV luminosity predicted by PHOENIX models (L_phot) to obtain $R^{*}_{\rm NUV}$ (the predicted values of L_phot as a function of V − J are shown by the dash-dotted line in Figure 3). The $R^{*}_{\rm NUV}$ values for our 20 YMG members are plotted as a function of their age in Figure 13. Also shown are 28 older field stars from our sample, which were identified by having space motions >1σ from the mean of active/young stars in at least two spatial dimensions (space motion information was obtained from Lépine et al. 2013). Similar to the findings of Shkolnik & Barman (2014), we observed a saturated level of NUV emission lasting a couple Myr, followed by a power-law decline.

Zoom In Zoom Out Reset image size

To derive our age-activity relation, we fit a broken power law to the median $R^{*}_{\rm NUV}$ values of each age group. To obtain the best-fit parameters and associated errors, we constructed 100 bootstrap samples from our data, found the model parameters associated with the minimum reduced χ² for each bootstrap sample, then took the mean and standard deviation. The results are shown in Figure 13, indicating a saturated NUV emission level at $R^{*}_{\rm NUV}=2.3\pm 0.19\times 10^{-6}$ until 160 ± 70 Myr of age, after which NUV emission declines as a power law with slope β = −0.83 ± 0.11. Our β value agrees well with those found by Stelzer et al. (2013) and Shkolnik & Barman (2014). Our saturation timescale is also consistent with that of Shkolnik & Barman (2014), although their value of 300 Myr is slightly higher, possibly due to our lack of data points between 100 and 600 Myr.

We investigated possible star-formation rate histories by deriving an age distribution for our sample, using a Monte Carlo approach. For each star we perturbed its $R^{*}_{\rm NUV}$ value by a random Gaussian deviate scaled to its error on $R^{*}_{\rm NUV}$. We also perturbed each parameter of the age-activity relation derived in Section 6.2 in an analogous manner. We used these perturbed values to estimate an age for each star and then constructed an age distribution by summing the 1/V_max weights of the stars in each age bin. We repeated this 100 times, then summed the normalized distributions to create the final age distribution shown in Figure 14. We found that at young ages our derived age distribution varied greatly depending on the input parameters, resulting in large errors and thus an uncertain distribution at young ages.

Zoom In Zoom Out Reset image size

There has been much discussion on the star-formation rate history of the Solar neighborhood. Gizis et al. (2002) used a spectroscopic survey of 676 nearby M dwarfs to infer a constant star-formation rate over the last 4 Gyr. However, there have been several studies (which do not utilize M dwarfs) that indicate elevated star-formation rates in recent history. Hernandez et al. (2000) used Hipparcos data to claim rapidly fluctuating star-formation rates with frequencies of ∼0.5 Gyr; however, our data do not have sufficient time resolution to be compared to their work. Bonatto & Bica (2011) used a sample of 442 star clusters within 1 kpc to show a recent (220–600 Myr) local burst in star formation that is twice the average star-formation rate. Tremblay et al. (2014) used the luminosity function of white dwarfs within 20 pc to show enhanced star-formation rates within the last 5 Gyr, compared to that of 5–10 Gyr. The latter two studies appear to be most consistent with our Figure 14, if the slope in our age distribution can be considered significant.

Our analysis indicates an era of saturated NUV emission for young M dwarfs lasting ∼100–200 Myr (although possibly up to ∼300 Myr; see Shkolnik & Barman 2014). Roughly 120 stars in our sample have sufficiently high NUV luminosities ($R^{*}_{\rm NUV}\gtrsim 2.3\times 10^{-6}$ or m_NUV − K_S ≲ 10.9) to place them in this saturation age interval. Correcting for an FP rate of ∼25% (see Figure 12) reduces this saturated sample to ∼90 stars or ∼2% of our sample. Any planets in the habitable zones of these young M dwarfs will be exposed to persistent, elevated NUV irradiation. Because the dissociation energies of several key atmospheric molecules are in the NUV (e.g., H₂O at 2398 Å, CO₂ at 2247 Å, CH₄ at 2722 Å), the atmospheres of these planets can be significantly altered by photodissociation. Although detailed studies of these processes are just beginning, recent results suggest the implications may be significant (e.g., see Miguel et al. 2015 for the effects of Lyman α radiation on the atmosphere of mini-Neptune GJ 436b, which orbits an M3 star).

The saturation timescale of a few hundred Myr also roughly corresponds to the final "giant impact" phase predicted by terrestrial planet formation models. This interval scales with the orbital period, which is proportionally shorter for planets in the compact habitable zones of M dwarfs (Morbidelli et al. 2012). This implies that the early atmospheres of planets around young M dwarfs are subject to erosion via heat injection from impactors in addition to NUV irradiation by their host stars.

We also find that the vast majority of stars in our saturated sample (∼70% after FP correction) have F_FUV/F_NUV ⩾ 0.1, which is at least two orders of magnitude above the solar value of ∼0.001 (see also France et al. 2013). This raises the potential for high rates of abiotic atmospheric O₂ and O₃ (produced from CO₂)—two molecules that have been proposed as biosignatures on Earth-like planets (see Tian et al. 2014, and references therein). Moreover, even some of the older M dwarfs in our sample, which exhibit basal levels of NUV emission above model-predicted photospheric values (see Figure 3), have F_FUV/F_NUV ⩾ 0.1. Thus very blue light may remain an important consideration for habitable-zone planets around these very red stars.