NGTS-7Ab: An ultra-short period brown dwarf transiting a tidally-locked and active M dwarf

TL;DR

This paper reports the discovery of NGTS-7Ab, the shortest period transiting brown dwarf around an active M dwarf, highlighting tidal interactions, magnetic activity, and the system's hierarchical triple nature.

Contribution

It presents the first detection of an ultra-short period brown dwarf transiting an active M dwarf, with detailed analysis of tidal synchronization and system dynamics.

Findings

Shortest period transiting brown dwarf around a star to date

Host star is magnetically active with multiple flares

System likely to experience orbital decay within 5-10 Myr

Abstract

We present the discovery of NGTS-7Ab, a high mass brown dwarf transiting an M dwarf with a period of 16.2 hours, discovered as part of the Next Generation Transit Survey (NGTS). This is the shortest period transiting brown dwarf around a main or pre-main sequence star to date. The M star host (NGTS-7A) has an age of roughly 55 Myr and is in a state of spin-orbit synchronisation, which we attribute to tidal interaction with the brown dwarf acting to spin up the star. The host star is magnetically active and shows multiple flares across the NGTS and follow up lightcurves, which we use to probe the flare-starspot phase relation. The host star also has an M star companion at a separation of 1.13 arcseconds with very similar proper motion and systemic velocity, suggesting the NGTS-7 system is a hierarchical triple. The combination of tidal synchronisation and magnetic braking is expected to drive ongoing decay of the brown dwarf orbit, with a remaining lifetime of only 5-10 Myr.

Tables

Table 1. Table 1: Stellar properties for each star. We have listed the photometry used in our SED fitting. We show the parallax and proper motions for reference, however do not use them all in our analysis for the reasons outlined in Sect.  3.2.2 . The references are: 1. Gaia Collaboration et al. ( 2018a ) , 2. Henden & Munari ( 2014 ) , 3. Skrutskie et al. ( 2006 ) , 4. Cutri & et al. ( 2014 ) .

Property	NGTS-7A	NGTS-7B	Source
R.A [$°$]	352.5216665551376	352.52202473338	1
Dec [$°$]	-38.96992064512876	-38.97006605140	1
Gaia Source ID	6538398353024629888	6538398353024172032	1
${\mu }_{R.A.}$ (mas yr-1)	$-27.003\pm 0.112$	$-28.601\pm 0.112$	1
${\mu }_{Dec}$ (mas yr-1)	$-16.225\pm 0.178$	$-14.776\pm 0.364$	1
Parallax (mas)	$7.2497\pm 0.1203$	$6.5232\pm 0.0787$	1
$B$	$17.091\pm 0.072$		2
$V$	$15.502\pm 0.028$		2
$g^{\prime }$	$16.187\pm 0.044$		2
$r^{\prime }$	$14.940\pm 0.010$		2
$i^{\prime }$	$13.822\pm 0.127$		2
Gaia $G$	$14.9154\pm 0.0020$	$15.5134\pm 0.0012$	1
$J$	$11.832\pm 0.030$		3
$H$	$11.145\pm 0.026$		3
$K_s$	$10.870\pm 0.019$		3
$W1$	$10.740\pm 0.022$		4
$W2$	$10.660\pm 0.020$		4

Table 2. Table 2: Parameters from our fitting of NGTS-7AB system for the scenarios defined in Sect.  3.2.3 . In scenario (i) we have placed both stars at the distance of NGTS-7B and fit for the radius, whereas in scenario (ii) we assumed both stars were on the main sequence and fit for both radius and distance. The bold values for scenario (i) are to indicate it is our favoured scenario, as discussed in Sect.  4.1 . Here we report the median of the 1D distribution for each parameter, along with the errors determined from the 16th and 84th percentiles. Limb darkening parameters with asterisks had priors applied when fitting (Sect.  3.3 ).

Scenario	(i)	(ii)
SED Fitting
TeffA (K)	${\mathrm{𝟑𝟑𝟓𝟗}}_{-\mathrm{𝟖𝟗}}^{+\mathrm{𝟏𝟎𝟔}}$	${3393}_{-31}^{+30}$
TeffB (K)	${\mathrm{𝟑𝟑𝟓𝟒}}_{-\mathrm{𝟏𝟒𝟕}}^{+\mathrm{𝟏𝟕𝟐}}$	${3300}_{-42}^{+44}$
$\log g$A	${4.89}_{-0.28}^{+0.40}$	${4.82}_{-0.23}^{+0.39}$
$\log g$B	${4.98}_{-0.34}^{+0.37}$	${4.99}_{-0.34}^{+0.36}$
$R_A$ ($R_{\odot }$)	${0.61}_{-0.06}^{+0.06}$	${0.34}_{-0.04}^{+0.04}$
$R_B$ ($R_{\odot }$)	${0.46}_{-0.07}^{+0.08}$	${0.28}_{-0.03}^{+0.03}$
$D_A$ (pc)	${152.67}_{-2.01}^{+2.01}$	${88.04}_{-8.79}^{+8.91}$
$D_B$ (pc)	${152.70}_{-1.99}^{+1.98}$	${88.04}_{-8.79}^{+8.91}$
Transit Parameters
Period (hours)	${16.2237952}_{-0.0000018}^{+0.0000026}$	16.2237957${}_{-0.0000021}{}^{}{}_{}{}^{+0.0000024}$
Time of transit centre (days) $T_{centre}$ (HJD - 2456658.5)	${1050.053304}_{-0.0000055}^{+0.0000017}$	1050.053311${}_{-0.0000125}{}^{}{}_{}{}^{+0.0000099}$
$R_A/a$	${0.20213}_{-0.00257}^{+0.00310}$	${0.20215}_{-0.00258}^{+0.00366}$
$R_{BD}/a$	${0.04710}_{-0.00061}^{+0.00093}$	${0.04725}_{-0.00062}^{+0.00121}$
$a$ (AU)	${0.0139}_{-0.0014}^{+0.0013}$	${0.0078}_{-0.0008}^{+0.0009}$
$i$ (°)	${88.43520}_{-1.10843}^{+0.98314}$	${88.43124}_{-1.29644}^{+1.01065}$
Surface brightness ratio	${0.03620}_{-0.01198}^{+0.01148}$	${0.03763}_{-0.01225}^{+0.01296}$
SAAO LD1*	${0.24872}_{-0.02002}^{+0.02043}$	${0.25023}_{-0.02080}^{+0.02006}$
SAAO LD2	${0.06045}_{-0.12362}^{+0.12719}$	${0.06297}_{-0.14502}^{+0.12759}$
EulerCam LD1*	${0.53550}_{-0.01732}^{+0.01645}$	${0.53480}_{-0.01690}^{+0.01760}$
EulerCam LD2	${0.15415}_{-0.24625}^{+0.19985}$	${0.17269}_{-0.22253}^{+0.17810}$
NGTS LD1*	${0.36273}_{-0.05013}^{+0.02752}$	${0.36208}_{-0.05178}^{+0.02798}$
NGTS LD2	${0.38254}_{-0.11373}^{+0.12664}$	${0.36759}_{-0.12993}^{+0.14530}$
Spot Parameters
Spot 1 $l$ ($°$)	${74.68895}_{-3.21393}^{+3.82344}$	${75.22438}_{-3.48933}^{+4.62037}$
Spot 1 $b$ ($°$)	${50.01602}_{-11.89659}^{+8.70891}$	${49.48639}_{-12.49710}^{+9.24656}$
Spot 1 size ($°$)	${13.87737}_{-2.67395}^{+3.60749}$	${13.76827}_{-2.56000}^{+3.05474}$
Spot 1 brightness factor	${0.48236}_{-0.25854}^{+0.18490}$	${0.46430}_{-0.24309}^{+0.17512}$
Spot 2 $l$ ($°$)	${176.06974}_{-3.76279}^{+4.63372}$	${176.58879}_{-4.08154}^{+5.62414}$
Spot 2 $b$ ($°$)	${77.97929}_{-2.10508}^{+1.81747}$	${77.47726}_{-2.56784}^{+1.98424}$
Spot 2 size ($°$)	${30.25273}_{-3.82615}^{+3.62982}$	${30.22503}_{-4.27825}^{+3.72268}$
Spot 2 brightness factor	${0.27168}_{-0.17487}^{+0.16630}$	${0.30954}_{-0.19153}^{+0.15558}$

Table 3. Table 3: \textcolor blackHARPS radial velocities for NGTS-7A and NGTS-7B from our analysis in Sect.  3.4 . The radial velocity of NGTS-7B is fixed to be constant during our analysis. The signal-to-noise ratios correspond to the spectral order 66 centered at 653 nm.

${\mathrm{BJD}}_{\mathrm{TDB}}$ (-2,450,000)	${\mathrm{RV}}_A$ (km s-1)	${\mathrm{RV}}_A$ error (km s-1)	${\mathrm{Contrast}}_A$ (per cent)	${\mathrm{RV}}_B$ (km s-1)	${\mathrm{RV}}_B$ error (km s-1)	${\mathrm{Contrast}}_B$ (per cent)	S/N
8364.50765417	-25.369	1.513	2.584	-7.751	0.054	9.261	6.1
8364.52902662	-20.596	1.640	2.673	-7.751	0.054	11.787	6.7
8373.47768215	14.435	2.326	2.703	-7.751	0.054	12.936	4.5
8373.49931132	17.142	3.578	1.744	-7.751	0.054	16.270	5.8
8373.52064613	19.385	1.159	2.608	-7.751	0.054	6.286	7.7

Equations

$$\delta_{filter}=\bigg{(}\frac{R_{BD}}{R_{A}}\bigg{)}^{2}\bigg{(}1+\bigg{(}\frac{F_{B}}{F_{A}}\bigg{)}\bigg{)}^{-1}$$

$$\delta_{eclipse}=\bigg{(}\frac{R_{BD}}{R_{A}}\bigg{)}^{2}\frac{\int F_{BD}(T_{BD})\,S(\lambda)\,d\lambda}{\int F_{A}\,S(\lambda)\,d\lambda}+A_{g}\bigg{(}\frac{R_{BD}}{a}\bigg{)}^{2}$$

$${\tau }_a\approx \frac{1}{13}\frac{h}{{\alpha }_{mb}{C}_{\ast }{\Omega }^3}$$

Full text

††pubyear: 2019††pagerange: NGTS-7Ab: An ultra-short period brown dwarf transiting a tidally-locked and active M dwarf–A

1 Introduction

The discovery of brown dwarfs in transiting exoplanet surveys provides a unique opportunity to probe the parameters of these substellar objects. With radii similar to Jupiter and masses between 13 and $\sim$78 $M_{\mathrm{J}}$ (e.g. Chabrier et al., 2000; Halbwachs et al., 2000), brown dwarfs are believed to form through molecular cloud fragmentation or gravitational instability, as opposed to the core accretion process that is commonly thought to form giant planets (e.g. Chabrier et al., 2014). Although the youngest substellar objects can have radii similar to early M stars (e.g. Stassun et al., 2006), as they age they undergo gravitational contraction (Lissauer, 2004). As brown dwarfs age their luminosity and temperature also decreases, resulting in their spectral energy distribution shifting towards longer wavelengths. As such, lone brown dwarfs can be identified in photometric surveys from their colours (e.g. Pinfield et al., 2008; Folkes et al., 2012; Reylé, 2018). However, as close companions to pre-main or main sequence stars, such identification is not possible and we must rely on their effects on the host star.

The large masses of brown dwarfs should provide easily detectable signatures in radial velocity measurements relative to those of exoplanets (e.g. km s*-1* instead of m s*-1* Brahm et al., 2016; Carmichael et al., 2019). Despite this, the number of transiting brown dwarfs relative to exoplanets remains low, with currently only 19 known to date.

The paucity of brown dwarfs on short periods around main sequence stars has previously been termed the “brown dwarf desert”, from radial velocity and transit observations (e.g. Campbell et al., 1988; Marcy & Butler, 2000; Grether & Lineweaver, 2006). This driving factor for this desert is typically attributed to the different formation mechanisms of low and high mass brown dwarfs in binary systems (e.g. Ma & Ge, 2014). High mass brown dwarfs ($\gtrapprox$ 43 $M_{\mathrm{J}}$) are believed to form through molecular cloud fragmentation, whereas their lower mass counterparts form within the protoplanetary disc. However, along with their formation pathways, a contributing element for the brown dwarf desert may be inward orbital migration of the brown dwarf (e.g. Armitage & Bonnell, 2002). One way of driving this is thought to be through tidal interactions between brown dwarfs and their host stars (e.g. Pätzold & Rauer, 2002; Damiani & Díaz, 2016), along with the effect of the magnetic braking of the host star (e.g. Barker & Ogilvie, 2009; Brown et al., 2011). If the companion is close enough, tidal interactions can decay its orbit, moving the companion inwards. The angular momentum lost from this orbit is expected to be transferred to the spin of the host star (e.g. Bolmont et al., 2012), eventually resulting in a state of spin-orbit synchronisation. In this state, the orbital and spin periods are equal. Such synchronisation has been detected in transiting brown dwarf systems before, for example in CoRoT-15b (Bouchy et al., 2011), a 63$M_{\mathrm{J}}$ brown dwarf orbiting an F7V star with a period of 3.06 days. Along with this, brown dwarf systems have shown behaviour close to synchronisation (e.g. WASP-128b Hodžić et al., 2018). However, during this process, magnetic braking will remove angular momentum from the system (Barker & Ogilvie, 2009). This acts to spin down the star, which in turn exacerbates the orbital decay of the companion. As such, even though (pseudo) spin-orbit synchronisation may be achieved, for active stars the magnetic braking can still drive the decay of the companion orbit. The combination of these effects eventually results in the engulfment of the brown dwarf by the host star. The timescale of this orbital decay is dependent on a number of factors, notably the stellar radius (Damiani & Díaz, 2016). Consequently, the decay timescale is expected to be shortest for brown dwarfs around G and K type stars (e.g. Guillot et al., 2014), making brown dwarf companions rarer around these stars (as noted by Hodžić et al., 2018) and contributing to the desert.

For M stars the orbital decay timescale is expected to be longer than G and K stars, due to the strong dependence of tides on stellar radius (e.g. Damiani & Díaz, 2016). This is in spite of the strong magnetic activity of M stars, which can manifest itself as both saturated quiescent X-ray emission and transient activity such as stellar flares (e.g. Hilton, 2011; Jackman et al., 2019). Of the 19 transiting brown dwarfs known to date, 4 brown dwarfs have been identified transiting M stars.Two of these systems are hierarchical triples consisting of two M dwarfs and a brown dwarf (NLTT41135 B, LHS 6343C; Irwin et al., 2010; Johnson et al., 2011), with the two M dwarfs in close proximity on the sky (2.4″, 55 AU and 0.55″, 20 AU respectively). Both these systems are believed to have ages greater than 1 Gyr and be in stable configurations. The third system, AD 3116 (Gillen et al., 2017), is a M+BD system discovered in the Praesepe open cluster and has an age of $\sim$ 700 Myr. This age makes it one of the younger transiting brown dwarf systems and useful for testing brown dwarf models with age. The fourth M+BD system is LP 261-75 (Irwin et al., 2018), a M+BD transiting pair with a distant visual brown dwarf companion (Reid & Walkowicz, 2006). LP 261-75 is expected by Irwin et al. (2018) to have an age of several Gyrs despite the high activity of the M dwarf primary, which instead suggests an age in the 130-200 Myr range (e.g. Reid & Walkowicz, 2006). This strong activity instead is associated with tides from interactions between the brown dwarf and the host star. These four systems show the range of ages and configurations these systems can have, highlighting how further observations of transiting brown dwarfs are required to understand their formation and evolution. In particular, the discovery of unstable systems are needed in order to test evolutionary scenarios.

In this paper we report the discovery of NGTS-7Ab, a brown dwarf transiting an active M star on a 16.2 hour orbital period. The host star’s rotation period is locked to the orbit of the brown dwarf, posing questions about the formation and evolution of such systems. We present our detection with NGTS, along with follow up photometric and spectroscopic measurements to constrain the radius and mass of the brown dwarf and M star host. We also present a detection of the secondary eclipse with NGTS, which we use to measure the temperature of NGTS-7Ab. This system is heavily diluted by a possibly associated nearby source. We describe the steps taken to account for this, along with presenting different scenarios based on the assumptions taken. We also discuss the possible formation scenarios of this system and outline how it may evolve in the future.

2 Observations

2.1 Photometry

2.1.1 NGTS

NGTS-7 was observed with NGTS for 130 nights between 2016 May 04 and 2017 Jan 11, using a single camera. The phase folded lightcurve is shown in Fig. 1. Observations were obtained in the custom NGTS filter (520-890nm) with a cadence of 13 seconds. For a full description of the NGTS instrument and pipeline processing see Wheatley et al. (2018). The NGTS lightcurves were detrended using a version of the sysrem algorithm, as done for previous NGTS discoveries (e.g. Bayliss et al., 2018; Raynard et al., 2018; West et al., 2018).

This star was originally identified as an object of interest due to the detection of flares as part of the NGTS flare survey (e.g. Jackman et al., 2018, 2019). We subsequently identified a 16.2 hour periodicity. We then noted transit events of 4.3 per cent depth occurring on the same period.

Gaia DR2 resolves two stars with a separation of 1.13″, while all other catalogues list it as a single source. \textcolorblackThe catalogue photometry and astrometry is given in Tab. 1. To confirm the source of the transits we perform centroiding using the vetting procedure described by Günther et al. (2017). We describe this analysis in Sect. 3.1 and refer to the two sources as NGTS-7A and NGTS-7B  where NGTS-7A is the transit source. The two stars have Gaia G magnitudes of 14.9 (NGTS-7A) and 15.5 (NGTS-7B), meaning that there is non-negligible dilution present in our photometry, something we discuss and account for in Sect. 3.3.

2.1.2 SAAO

Follow up photometry of NGTS-7 was obtained at the South African Astronomical Observatory (SAAO) on 2018 Aug 08 (I band, secondary eclipse), 2018 Aug 11 (I band, primary transit), and again on 2018 Oct 04 (I band, secondary eclipse) using the 1.0m Elizabeth telescope and “shocnawe”, one of the SHOC high speed CCD cameras (Coppejans et al., 2013). On each occasion, sky conditions were clear throughout the observations, with the seeing around 2 arcseconds. The data were reduced with the local SAAO SHOC pipeline developed by Marissa Kotze, which is driven by python scripts running iraf tasks (pyfits and pyraf), and incorporating the usual bias and flat-field calibrations. Aperture photometry was performed using the Starlink package autophotom. We used a 5 pixel radius aperture that maximised the signal to noise ratio, and the background was measured in an annulus surrounding this aperture. One bright comparison star in the $2.85\times 2.85$ arcminute field of view was then used to perform differential photometry on the target. The two stars identified by Gaia DR2 coincident with the position of NGTS-7 were not resolved in these data. Figure 2 shows the primary transit observed on 2018 Aug 11. A stellar flare can be clearly seen shortly before transit ingress.

2.1.3 EulerCam

One transit of NGTS-7 was observed with EulerCam on the 1.2m Euler Telescope at La Silla Observatory (Lendl et al., 2012). These observations were obtained on the night of 2018 Sept 01, in the V band filter and are shown in Fig. 2. The data were bias and flat field corrected then reduced using the PyRAF implementation of the “phot" routine. An aperture radius and ensemble of comparison stars were used such that the scatter in the out of transit portion of the lightcurve was minimised.

2.1.4 TESS

NGTS-7 was observed at a 30 minute cadence with the NASA Transiting Exoplanet Survey Satellite (TESS) (Ricker et al., 2015) between 2018 Aug 27 and 2018 Sep 19, in Sector 2. A 15 $\times$ 15 pixel (5.25’ $\times$ 5.25’) cutout was obtained from the TESS full-frame image stacks using the tesscut routine 111https://github.com/spacetelescope/tesscut. This cutout is shown in Fig. 3. Aperture masks were chosen by-eye to exclude nearby bright sources up to 2.5’ away. The 21" pixel-scale of TESS creates a PSF of NGTS-7 which is blended with at least 3 significantly bright stars ($\Delta G<3.5$ mag.) As it is not possible to completely exclude the flux from these blended stars in TESS we chose our aperture to enclose them, with the knowledge the TESS lightcurve will be diluted. We estimated the per-pixel background contribution by selecting 8 pixels West of the aperture that do not include any stars brighter than $G=18.4$ (3.5 magnitudes fainter than NGTS-7A). This region is shown as the magenta box in Fig. 3. This was subtracted from the aperture-summed flux to create a background-corrected light curve.

The transit seen in the TESS light curve is both shallower and more V-shaped than that from NGTS, despite the similar bandpasses of NGTS and TESS. This is due to a combination of additional dilution in the TESS data (from the neighbouring sources) and the 30 min cadence which smears out the transit (which has a duration of only 1.3 hours, e.g. Smith et al., 2018). Due to these effects we do not use the TESS lightcurve in our transit fitting (Sect. 3.3). \textcolorblackHowever, we do use it in Sect. 3.7 where we discuss the phase of the out-of-transit variations of NGTS-7.

2.2 Spectroscopy

2.2.1 HARPS

We obtained high-resolution spectroscopy for NGTS-7A with the HARPS spectrograph on the ESO 3.6m telescope (Mayor et al., 2003). Five measurements \textcolorblackwith an exposure time of 1800s were taken on the nights beginning 2018 Sept 02 and 2018 Sept 11 as part of programme ID 0101.C-0889(A). Due to the relative faintness of the source we used the high efficiency fibre link (EGGS), with a fibre size of 1.4″ instead of the usual 1.0″ mode. Consequently, these spectra contain light from both NGTS-7A and NGTS-7B \textcolorblackand we see a narrow and a broad peak in the Cross Correlation Functions (CCFs) shown in Fig. 5. \textcolorblackThe RVs of NGTS-7A and NGTS-7B along with the respective contrasts from our analysis in Sect. 3.4 are given in Tab. 3.

2.2.2 SAAO

Follow up spectroscopy of NGTS-7 was also obtained from SAAO on the 1.9m telescope using the SpUpNIC instrument (Crause et al., 2016) between the dates 2018 Sept 09 and 2018 Sept 11. 14 spectra with a resolution of R=2500 were obtained in total, with a wavelength range of 3860-5060Å. We have combined these spectra to give the average spectrum shown in Fig. 6. Observations were performed with a slit width of 1.8″ and average seeing of 2″, once again meaning both NGTS-7A and NGTS-7B are present in our data. Figure 6 shows clear TiO and CaI absorption features expected for M dwarf spectra. Along with this, we see several emission lines from the Balmer series, as well as He I and Ca II, showing at least one of the stars is chromospherically active.

3 Analysis

The observations of Sect. 2 were combined with available catalogue photometric and astrometric information. We use this information to \textcolorblackconfirm the source of the transits and characterise both NGTS-7A and NGTS-7B.

3.1 Identifying the Source of the Transit

\textcolor

blackIn Sect. 2 we noted that Gaia DR2 resolves two sources with a separation of 1.13″ at the position of NGTS-7. To confirm which source our transit signal is coming from, we performed centroiding using the vetting procedure described by Günther et al. (2017). We identify that the transit and the out-of-transit modulation comes from Gaia DR2 6538398353024629888, the brighter of the two sources. Figure 7 shows the phase folded transit and X centroid position, showing how the shape of the phase-folded centroid data follows the shape of the phase-folded lightcurve. While individual NGTS pixels are 5” across, the NGTS centroiding procedure is able to identify centroid shifts below 1” in size, meaning we are confident that we have identified the correct host star and now refer to this star as the primary star, or NGTS-7A. We refer to the neighbouring source as NGTS-7B and discuss it further in Sect. 3.2.3.

\textcolor

blackOut of transit modulation on the orbital period can be due to either ellipsoidal variation (e.g. Drake, 2003; Welsh et al., 2010) or reflection effects (e.g. Armstrong et al., 2016; Eigmüller et al., 2018). However, neither of these could adequately explain the number or position in phase of the maxima seen in Fig. 1 (just before the primary transit). The most natural explanation is that this out-of-transit modulation is due to starspots on the host star and that the spin period of NGTS-7A is locked to the orbital period of the transiting body. This places NGTS-7A in a state of spin-orbit synchronisation (e.g. Ogilvie, 2014). The change in the out of transit modulation in the TESS data can be explained by the evolution of starspots in the interval between the NGTS and TESS observations.

\textcolor

blackThe 16.2 hour period rotation of NGTS-7A will result in its observed CCF in our HARPS spectra being rotationally broadened. This broadened peak will also move around with a 16.2 hour period. In Sect. 2.2.1 we noted that our HARPS spectra contain light from NGTS-7A and NGTS-7B and the presence of a narrow and broad peak in our HARPS CCFs, seen in Fig. 5. As we will discuss in Sect. 3.4 we find the broad peak moves on a 16.2 hour period, as we might expect if NGTS-7A has a transiting body and itself is in a state of spin-orbit synchronisation. The rapid rotation of NGTS-7A and it being chromospherically active (as evidenced by the observed starspots) presumably means NGTS-7A is the source of the multiple stellar flares in the NGTS and SAAO lightcurves (e.g. Hawley et al., 2014). Along with this, NGTS-7A is likely the dominant source of the observed emission lines in our SAAO spectra.

\textcolor

blackBased on our observations and the evidence presented here we are confident that NGTS-7A is the source of the observed transits. Along with this we believe NGTS-7A is in a state of spin-orbit synchronisation with its companion, which will have spun up NGTS-7A to keep it at the observed period.

3.2 Stellar Parameters

Throughout this paper, all of our photometry is measured in apertures which contain the light from both NGTS-7A and NGTS-7B. Consequently, in order to obtain accurate parameters for NGTS-7Ab we need to estimate the dilution from NGTS-7B. We have done this through fitting the Spectral Energy Distribution (SED) of both NGTS-7A and NGTS-7B, using a combination of the information available from blended catalogue photometry and Gaia photometric and astrometric data. This information is listed in Tab. 1.

3.2.1 Gaia Photometry

While both sources have a Gaia G magnitude, only the primary star has BP and RP photometry. The Gaia G magnitudes for both stars are published in Gaia DR2, and are derived from fitting the line spread function (LSF) of each star from windows which are approximately 0.7$\times$2.1 arcsec2 in the along scan and across scan directions respectively (Gaia Collaboration et al., 2016). We used the Gaia Observation Scheduling Tool (GOST222https://gaia.esac.esa.int/gost/) to check the scans of NGTS-7A and NGTS-7B used for Gaia DR2. We obtained 35 scans, which are plotted over a SkyMapper i band image (Wolf et al., 2018) in Fig. 8. From Fig. 8 we can see that over 75 per cent of the scans which went into Gaia DR2 fall along (or close to) the position angle separating the two stars. \textcolorblackGiven the 1.13″ separation of the stars and the ability of the LSF to resolve sources of this separation (e.g. Fig. 7 in Fabricius et al., 2016), we expect only minimal contamination between the stars in the Gaia G band photometry. Consequently, we use the Gaia G band photometry in our analysis.

On the other hand, the BP and RP photometry is measured from the total flux in a 3.5x2.1 arcsec2 region (Evans et al., 2018). An example of this region is shown in Fig. 8, showing that the BP and RP photometry will be of both NGTS-7A and NGTS-7B combined. This is reflected in the BP-RP excess factor of 2.054 for NGTS-7A. The BP-RP excess factor is the sum of light from the BP and RP bands compared to the G band, and should ideally be around 1 for a single, non-contaminated, star. A value around 2 suggests the BP and RP photometry is comprised of flux from two similar stars. Indeed, we find that the BP and RP photometry of NGTS-7A fails the filter from Arenou et al. (2018), which is used to remove contaminated stars from their analysis. Consequently we do not use the Gaia BP and RP photometry of NGTS-7A in our analysis.

3.2.2 Gaia Astrometry

For both sources we initially test the quality of the Gaia astrometry by calculating both the Unit Weight Error (UWE) and the Renormalised UWE (RUWE). We compare the UWE against the filter specified by Lindegren et al. (2018) and check whether the RUWE is below the recommended value of 1.4 for a clean astrometric sample. We found that NGTS-7A suffers from significant astrometric excess noise (astrometic excess noise sig=71.6, RUWE=3.4), resulting in it failing both filters. NGTS-7B, while having non-zero astrometric excess noise (astrometic excess noise sig=4.4, RUWE=1.3), passes both filters. When calculating the astrometric solution of each star, Gaia DR2 assumes a single object. The astrometric excess noise is the extra noise that is required by the single source solution to fit the observed behaviour. High levels of astrometric excess noise \textcolorblackare a sign that the single source solution has failed, possibly due to unresolved binarity (e.g. Gaia Collaboration et al., 2018b). We also check each star further by comparing them against sources of similar magnitude, colour and parallax in the full Gaia DR2 sample. Both stars are outliers from the main sample in terms of their astrometric quality. We note in particular that each has a correlation between their parallax and proper motion components. One possibility for the low quality of the astrometric parameters for NGTS-7A and NGTS-7B may be levels of blending due to their proximity. Lindegren et al. (2018) has noted that during scanning of close sources the components can become confused, through a changing photocentre.

\textcolor

blackDue to it failing the recommended astrometry filters we have decided not to use the astrometric solution of NGTS-7A in our analysis. As we explain in Sect. 3.2.3 we consider two scenarios. The first of these uses only the astrometric solution of NGTS-7B and fixes both NGTS-7A and NGTS-7B at the distance of NGTS-7B, while the second doesn’t use Gaia parallaxes and assumes both sources are on the main sequence.

3.2.3 A possible wide binary

A scenario mentioned in Sect. 3.2.2 which may be responsible for the low quality astrometry of NGTS-7A and NGTS-7B is that the two sources are a wide binary. If they are a wide binary, then we would expect them to be at the same distance. NGTS-7A and NGTS-7B have very similar proper motions, which supports this assumption, which are shown for reference in Tab. 1. However, as the proper motions are measured as part of the Gaia astrometry and may have levels of contamination, we have decided to seek out additional evidence. Wide binaries have previously been identified in both TGAS (e.g. Andrews et al., 2017) and Gaia DR2, with Andrews et al. (2018) finding that, as expected, real binaries will have similar systemic velocities, whereas chance alignments will not. From our radial velocity analysis in \textcolorblackSect. 3.4 we found that NGTS-7A and NGTS-7B have systemic velocities of $-4.2\pm 0.8$ km s*-1* and $-7.7\pm 0.1$ km s*-1* respectively. Using the distance of NGTS-7B results in a projected separation of 173 AU. This projected separation and the difference in systemic velocities places NGTS-7A and NGTS-7B well within the Andrews et al. (2018) sample of genuine wide binaries, instead of being a chance alignment on the sky. Consequently, it is very likely that NGTS-7A and NGTS-7B  are in fact a wide binary and are at the same distance. If so, this would provide a way of constraining the distance to NGTS-7Ab, along with placing it in a hierarchical triple system. \textcolorblackChecking for possible memberships of known associations using the BANYAN $\Sigma$ online tool 333http://www.exoplanetes.umontreal.ca/banyan/banyansigma.php reveals no likely associations (Gagné et al., 2018).

Following this, we have devised two separate scenarios on the assumption that NGTS-7A and NGTS-7B are in a wide binary. These are as follows:

We fix both sources at the distance of NGTS-7B, assuming they are a wide binary. 2. 2.

We believe neither Gaia DR2 parallax, instead fixing them at the same distance and assuming they are on the main sequence.

These scenarios both avoid using the poor astrometric solution of NGTS-7A.

3.2.4 SED Fitting

To determine the SED of both stars we have fitted two separate components simultaneously using a custom SED fitting process which utilises the PHOENIX v2 grid of models (Husser et al., 2013), following a similar method to Gillen et al. (2017). Initially we generated a grid of bandpass fluxes and spectra in Teff-$\log g$ space, which allowed us to interpolate across these parameters. We fit for Teff, $\log g$, along with the radius, $R$, and distance, $D$, of each star. We have chosen to fix the metallicity at the Solar value. Prior to fitting we inflated the errors of catalogue photometry by 2.5 per cent to account for the observed variability in the NGTS lightcurve. During fitting we compare the combination of fluxes from each star to the observed values, for all filters in Tab. 1 except Gaia G (which is used as a prior to normalise the respective SEDs). To explore the full posterior parameter space we use emcee (Foreman-Mackey et al., 2013) to generate an MCMC process, using 200 walkers for 50,000 steps, disregarding the first 25,000 as a burn in.

We have used a range of physically motivated priors in our modelling which we outline here. Firstly, the radii and distances are used in our model to scale the flux from each star by $(R/D)^{2}$. For scenario (i) (Sect. 3.2.3) we have placed a Gaussian prior on the distance of each star, using the value from Bailer-Jones et al. (2018) for NGTS-7B, $152.7\pm 1.9$ pc. In this scenario the fitted radius of each star is allowed to vary freely. For scenario (ii) we fit for the distance, which we also force to be the same for the two stars. We have placed a Gaussian prior on the fitted radius for each star, using the Mann et al. (2015) Teff-radius relation. For this prior, we have used the 13.4 per cent error given by Mann et al. (2015) as the standard deviation of the Gaussian prior to allow some variation. In both scenarios we have placed a prior on the synthetic Gaia G band flux for each star, using the observed flux values. This was done to anchor each star to observations.

Table 2 gives the results of each fit. For both scenarios we retrieve two stars with temperatures corresponding to M3-M4 spectral type (e.g. Pecaut & Mamajek, 2013). This similarity in spectral type matches what we would expect from the Gaia BP-RP excess factor, as discussed in Sect. 3.2.1. However, for scenarios (i) and (ii) we measure very different stellar radii for NGTS-7A and NGTS-7B. Investigating the posterior distribution of our SED fit reveals a strong correlation between the Teff and radius values of NGTS-7A and NGTS-7B. An example of this for Teff is shown in Fig. 9, with the full corner plot shown in Fig. 15. The full corner plot for scenario (ii) is shown in Fig. 16. The correlation between Teff and radius arises from the similarity of the two sources in spectral type, along with the availability of only the Gaia G magnitude to separate them. This correlation needs to be taken into account when determining the uncertainties in the age and mass of NGTS-7A. To incorporate these correlations we fit the 2D posterior distributions from our SED fitting with ellipses covering 68 per cent of our distribution. We have used these ellipses to probe the extremes of parameter space and incorporate the observed correlations into our analysis (Sect. 3.2.5). For each parameter we also report the 16th, 50th and 84th percentile of the marginalised 1D distributions in Tab. 2.

3.2.5 Primary Mass

When we fix both stars to the Gaia distance for NGTS-7B (scenario (i)), the median radius of the primary star is approximately 75 per cent oversized in radius compared to that of a main sequence star of the median Teff. One possible reason for this is that NGTS-7A and NGTS-7B are pre-main sequence stars and as such both have a larger than expected radius (e.g Jackman et al., 2019). In order to estimate the mass of the primary star we compared each source to the PARSEC isochrones (Bressan et al., 2012), assuming that NGTS-7A and NGTS-7B are the same age (reasonable if we assume they are bound). From comparing the median radius and Teff of NGTS-7B to the PARSEC isochrones we obtained an age estimate of 55 Myr and a mass of 0.35 $M_{\odot}$ for NGTS-7B. However, using this age estimate with the fitted parameters of NGTS-7A results in different mass estimates based on whether we use the median Teff (0.35 $M_{\odot}$) or the radius (0.55 $M_{\odot}$).

A potential reason for this discrepancy is the effect of starspots on NGTS-7A. For both main and pre-main sequence stars, modelling of starspots has shown they can act to both increase the stellar radius and decrease Teff (Jackson & Jeffries, 2014; Somers & Pinsonneault, 2016). The combined effect of these changes can be a diminished stellar luminosity (Jackson & Jeffries, 2014), which results in discrepancies when comparing to unspotted stellar models.

To correct for the effects of spots on our mass estimate for a given age, we used the PARSEC models to identify which unspotted models give a luminosity equal to or up to 10 per cent greater than the current value (this is approximately the change in luminosity caused by the sudden appearance of spots simulated by Jackson & Jeffries, 2014). This was done for the median Teff and radius values of each parameter, as given in Tab. 1. In this analysis we used the 1$\sigma$ extremes determined from the posterior distribution error ellipses from Sect. 3.2, in order to take correlations into account. This resulted in an age of $55^{+80}_{-30}$ Myr for NGTS-7B. We calculate the luminosity of NGTS-7A based on the median and the error ellipse also, to give a range of possible masses depending on the obtained primary parameters and age of NGTS-7B. Using the possible ages of the neighbour to determine the unspotted model we estimate the mass of the primary as 0.48^{+0.03}_{-0.12}$$M_{\odot}. For the errors we have combined the extremes from the age of the neighbour and whether the luminosity is altered by the appearance of spots.

Based on the age estimate of $55^{+80}_{-30}$ Myr for this system we have also searched for signs of Li 6708Å  absorption in our HARPS spectra. Primordial lithium is quickly depleted within the interiors of M stars (e.g. Chabrier et al., 1996) and is typically removed from their photosphere within 45-50 Myr (see Fig. 4 of Murphy et al., 2018). We do not find any sign of Li 6708Å  absorption in our HARPS spectra, consistent with our estimate of $55^{+80}_{-30}$ Myr and it suggests that the system cannot be much younger than this if scenario (i) is true.

We note we have assumed in this section that the companion star does not also suffer from spots, which may be unlikely for a young system. The presence of spots would alter the inferred age and hence mass of the primary star. However, we do not identify any significant modulation in either the NGTS or TESS lightcurves which could be attributed to spots on the companion.

We also note that while there exist empirical relations to attempt to correct for the effects of magnetic activity on measured Teff and radii (Stassun et al., 2012), using the ratio of the H$\alpha$ and bolometric luminosity, $\log L_{\mathrm{H{\alpha}}}/L_{\mathrm{Bol}}$. These relations are used to bring the Teff and radius values closer to expected model values, which can then be used to calculate the age and mass of NGTS-7A. Unfortunately it is likely that our measurements of H$\alpha$ luminosity for NGTS-7A are contaminated by NGTS-7B to an uncertain degree (from H$\alpha$ emission of its own). Consequently, we have chosen not to use these relations to adjust our fitted values here, but do discuss this further in Sect. 4.4.

For the second scenario where we have assumed both stars are drawn from the Mann et al. (2015) Teff-radius relation, we calculate a distance of $88.04^{+8.91}_{-8.79}$pc, given in Tab.2. To calculate the mass in this scenario we use the empirical mass relation of Benedict et al. (2016) for main sequence M stars. We have calculated the value of $M_{\mathrm{K_{s}}}$ for NGTS-7A using the best fitting SED model and the fitted distance. Using this relation with calculated distance of $88.04^{+8.91}_{-8.79}$pc for the NGTS-7A we calculate the primary mass $M_{\mathrm{A}}$ to be 0.24\pm 0.03$$M_{\odot}.

3.3 Transit and Spot Fitting

In order to model the transits of NGTS-7Ab we used the ELLC package (Maxted, 2016). ELLC is a binary star model that allows for multiple spots to be included on each star and as such can be used to model both transits and spot modulation at the same time.

We simultaneously fit the NGTS, SAAO and EulerCam lightcurves to ensure consistent transit parameters across our entire dataset. For the NGTS data we fit a transit model combined with a two spot model, to account for the out of transit modulation. We tested our fitting using both a single and double spot model, however we found a single spot was unable to match the average out of transit behaviour seen in Fig. 1. The transit in the TESS data is blurred by the 30 min cadence of the observations, and also suffers additional dilution from a number of blended sources (see Fig.3), and so we decided not to include the TESS light curve in our fit. We can use the TESS data to see spot modulation has changed between the NGTS observations and the SAAO and EulerCam follow up lightcurves (which were obtained at similar times to the TESS data). Consequently we did not use the NGTS spot model to fit the SAAO or EulerCam follow up lightcurves. However, the SAAO lightcurve of primary transit on 2018 Aug 11 does show evidence of the spot minimum during the single night, consistent with the TESS data. We incorporated this into our fitting as a quadratic term which we fit simultaneously with the transits.

The SAAO lightcurve also includes a flare just before ingress, which we masked out for our fitting but analyse in Sect. 3.8. For each bandpass we directly fitted independent limb darkening profiles. We used a quadratic limb darkening profile and generated our initial limb darkening parameters using the Limb Darkening Toolkit, (LDtk; Parviainen & Aigrain, 2015), using the best fitting SED from Sect. 3.2. During fitting we allowed each second limb darkening coefficient (LD2) to vary, while keeping the first (LD1) constant to reduce degeneracy in the fit. For each photometric band we also incorporate a dilution term, to account for the flux from the neighbouring star. For each band we use a Gaussian prior based on the expected dilution (and standard deviation) from our SED fitting. To estimate the expected dilution in a given bandpass we convolve the SED for each star with the specified filter curve and take the ratio of measured values. In order to take the observed correlations into account we sample the expected values for the Gaussian prior directly from the posterior distribution of the SED fits. For each filter we use the dilution term to correct the transit model as

[TABLE]

where $\delta_{filter}$, $R_{BD}$, $R_{A}$ are the transit depth in the chosen filter, radii of the companion and NGTS-7A respectively, while $F_{B}$ and $F_{A}$ are the fluxes of NGTS-7B and NGTS-7A in the specified bandpass. In the ideal scenario where $F_{B}$=0 we can see this becomes the usual transit depth equation. During our preliminary fitting we found the eccentricity to be consistent with zero when applying the Lucy & Sweeney (1971) criterion. Consequently for our final fitting we fixed the eccentricity at zero, i.e that the orbit has circularised. Due to the high time cadence of NGTS it is not feasible to fit the entire NGTS lightcurve for each step of the MCMC process. Instead we bin the lightcurve to 1000 bins in phase, using the period and epoch specified for that step. We chose 1000 bins in order to preserve the information in the ingress and egress. In order to sample the posterior parameter space we used emcee with 200 walkers for 50,000 steps and disregarding the first 25,000 as a burn in. We did this for both scenarios (i) and (ii), using the dilution values from the relevant SED model. The values of the best fitting parameters are shown in Tab. 2.

Using the results of our transit fitting for scenarios (i) and (ii), we measure the radius of NGTS-7Ab to be $1.38^{+0.13}_{-0.14}$ RJ and $0.77\pm 0.08$ RJ respectively. As brown dwarfs are expected to shrink with age (e.g. Baraffe et al., 2003), scenario (i) would imply a younger brown dwarf than scenario (ii), consistent with our age estimation from Sect. 3.2.5. The single period in our fitting is able to model both the orbital and spin periods, supporting our conclusion in Sect. 3.1 that the system is in a state of spin-orbit synchronisation.

Our best fitting spot model suggests the presence of two spot regions with a large size and a low brightness factor. Each region can be interpreted either as a single large spot of constant brightness, or as a series of smaller, darker, spots spread over a similar area. As we only fit for the dominant spots, our model is unable to rule out the presence of spots elsewhere on the star. It is most likely that these are smaller than our fitted regions however, as large spots elsewhere could act to decrease the observed variability (e.g. Rackham et al., 2018).

3.4 Radial Velocity

When analysing the HARPS data to measure the radial velocity shifts due to NGTS-7Ab we used the standard HARPS data reduction software (DRS) to obtain \textcolorblackour measured CCFs. Initial analysis with the M2 mask showed no variation between phases in the CCFs, with a seemingly unchanging peak at $-7.7\pm 0.1$ km s*-1*. It was realised that due to the fast rotation of NGTS-7A, the spectral lines were too broadened for the M2 mask (which uses a fine grid of molecular lines), resulting in a low signal-to-noise ratio CCF. It was found that analysing with the K5 mask (which uses fewer lines and is less susceptible to the fast rotation) showed both a CCF peak due to the background source and a shallow wide peak due to the motion of NGTS-7A, shown in Fig. 5. The increased width of this peak is due to the fast rotation of NGTS-7A. We confirmed both peaks were also present when using earlier spectral type masks, albeit at a lower signal to noise. With the CCFs from the K5 mask we simultaneously fit all our HARPS CCFs with two Gaussians plus an additional linear background term. Each Gaussian is allowed to vary in amplitude and midpoint, but is required to have a constant width. To perform our simultaneous fitting we once again use an MCMC process with emcee, with 200 walkers for 20,000 steps. We use the final 5000 steps to calculate our parameters and the results of our fitting are shown in Fig. 5. We folded the measured CCF peak midpoints in phase using the orbital period from Sect. 3.3 and fit a sinusoidal signal, shown in Fig. 12. \textcolorblackWe also list the measured midpoints and amplitudes in Tab. 3. As the orbit of NGTS-7Ab to have circularised (Sect. 3.3), we fitted the RV data using a single sinusoid. We fixed the period and epoch of this sinusoid to the values measured from our transit fitting. \textcolorblackWith this fit we measure a systemic velocity of $-4.2\pm 0.8$ km s*-1* and a semi-amplitude of $25.9\pm 0.9$ km s*-1* for NGTS-7A. We measure a systemic velocity of $-7.7\pm 0.1$ km s*-1* for NGTS-7B. Combining our measurement of the semi-amplitude for the radial velocity curve with the mass of NGTS-7A we calculate a mass of $75.5^{+3.0}_{-13.7}$ MJ for the transiting source for scenario (i). For scenario (ii) we obtain a value of $48.5\pm 4.3$ MJ.

The measured $48.5\pm 4.3$ MJ mass for scenario (ii) places NGTS-7Ab within the brown dwarf regime, making the system an brown dwarf transiting a main sequence M star. Our result for scenario (i) places NGTS-7Ab at the upper end of the brown dwarf regime, near the hydrogen-burning mass limit of $\sim$ 78$M_{\mathrm{J}}$ (Chabrier et al., 2000).

3.5 Rotational Broadening

We can also use our HARPS data to investigate the level of rotational broadening for NGTS-7A \textcolorblackand in turn constrain our radius measurement. We can construct a lower limit by assuming the profile of NGTS-7B is non-rotating and assuming a Gray (2005) profile to artificially broaden it to match the profile of NGTS-7A. We have assumed a limb darkening coefficient for the rotational profile of 0.55. Artificially broadening the CCF of NGTS-7B gives a lower limit of 31 km s*-1* for $v\sin i$, equivalent to a radius of 0.41$R_{\odot}$. This value is greater than the measured radius for NGTS-7A for scenario (ii) (a main sequence system with $R_{A}$=0.34 $R_{\odot}$) and is only consistent with scenario (i) (a pre-main sequence system with $R_{A}$=0.61 $R_{\odot}$).

3.6 Secondary Eclipse and Brown Dwarf temperature

As part of our fitting of the NGTS lightcurve we have identified evidence of a secondary eclipse for NGTS-7Ab, shown in Fig. 1. The presence of a secondary eclipse by its very nature implies non-negligible levels of flux from the brown dwarf itself. To estimate the temperature of NGTS-7Ab we equate the depth of the secondary eclipse to the ratio of fluxes in the NGTS bandpass,

[TABLE]

where $F_{\mathrm{BD}}(T_{\mathrm{BD}})$ and $F_{\mathrm{A}}$ are the SEDs of the brown dwarf (with temperature $T_{\mathrm{BD}}$) and NGTS-7A respectively, $S(\lambda)$ is the transmission curve of the NGTS filter (Wheatley et al., 2018) and $A_{\mathrm{g}}$ is the geometric albedo. For the SED of the primary star we use the results from our SED fitting. To generate the spectrum of the brown dwarf we have used the BT-Settl models (Allard et al., 2012), since the PHOENIX v2 models do not cover the full range of temperatures we wish to probe. For each model spectrum we have renormalised it to the distance of the primary star and to the expected brown dwarf radius. We opted to use these models instead of a blackbody due to the strong absorption features expected in the brown dwarf spectrum (e.g. Martín et al., 1999). We measured $\delta_{\mathrm{eclipse}}$ from the best fitting transit and spot model, making sure to correct for the effect of dilution in the NGTS bandpass. By including $A_{\mathrm{g}}$ we can also account for the effects of reflection. We have solved Eq. 2 in two limiting cases. These are $A_{\mathrm{g}}$=0 (no light is reflected) and $A_{\mathrm{g}}$=0.5. Iterating $T_{\mathrm{BD}}$ between 1200 and 3500 K returns estimates of 2880 K ($A_{\mathrm{g}}$=0.5) and 3200 K ($A_{\mathrm{g}}$=0) for scenarios (i) and (ii).

Comparing these temperatures to the Baraffe et al. (2015) models for an isolated 75.5 $M_{\mathrm{J}}$ brown dwarf results in ages up to 80 Myr, depending on the chosen value of $A_{g}$. This is in agreement with our estimate of 55 Myr for the age of this system assuming our scenario (i) in which the system is located at the distance implied by the Gaia DR2 parallax of the companion NGTS-7B (Sect. 3.2.3). In contrast, the Baraffe et al. (2015) models for a 48.5 $M_{\mathrm{J}}$ brown dwarf is not able to match the measured temperature range at any age. This high temperature of the brown dwarf, heavily disfavours and effectively rules out scenario (ii), in which both M stars were assumed be on the main sequence and hence at a smaller distance. Note that in scenario (ii) the brown dwarf would have to have a mass that was well below the hydrogen burning limit (Sect. 3.4).

3.7 Starspots

As part of our analysis in \textcolorblackSect. 3.1 and Sect. 3.3 we identified starspots were present in the NGTS and TESS data. When fitting the NGTS data, we allowed for two starspots in our model and assumed they were representative of the average starspot behaviour of NGTS-7A. Another check for whether this modulation is due to starspots is to look for evolution throughout the NGTS lightcurve. As starspots form and dissipate they will alter both the level of lightcurve modulation and the phase at which it occurs (e.g. Davenport et al., 2015; Jackman et al., 2018). To search for such changes within the NGTS data we split our data into 20 day sections. Visual inspection of the phase folded lightcurve in these sections showed that the modulation was slowly changing with time, indicative of starspot evolution. To show this we have phase folded each section in bins of 0.04 in phase and plotted the flux of each phase folded lightcurve against time in Fig. 13, following the method of Davenport et al. (2015).

From Fig. 13 we can see the movement of a dominant starspot group from around phase 0.5 to 0.25 over approximately 70 days. The level of modulation from this group is not constant, as the amplitude decreases at around 60 days in Fig. 13 before increasing again. One reason for this may be the dissipation and then formation of starspots from a large active region, which would act to change the overall level of modulation. From Fig. 13 we cannot identify any regions in the lightcurve where the starspot modulation disappears completely, meaning we are unable to measure the full starspot lifetime. Observations of M stars with Kepler have shown that they can have starspots with lifetimes on the order of years (Davenport et al., 2015; Giles et al., 2017). Consequently, it is not unexpected that we do not observe drastically changing spot regions within the NGTS data alone.

A possibility for the apparent shift in starspot phase is that the stellar spin period is not exactly the same as the orbital period. A small enough offset may appear as a drift in phase without appearing as an anomaly in the phase folded data. We use Fig. 13 to estimate what this drift may be, by assuming that the starspot distribution remains constant and using the change in phase as an estimate of the period difference. From this we find a shift of -0.28 in phase over 100 days, approximately equal to a difference of 162 seconds per day. The starspots moving backwards in orbital phase would imply the star is spinning slightly faster than the orbital period ($P_{\mathrm{orb}}/P_{\mathrm{spin}}=1.002$). One check for this is to mask the transits out of the original NGTS lightcurve and search for periodicity in the remaining data using a generalised Lomb-Scargle periodogram, using the astropy LombScargle package (Astropy Collaboration et al., 2013). Doing this and using 20,000 steps between 0 and 2 days returns a best fitting period of 16.204 hours, implying that NGTS-7A is slightly supersynchronous, spinning roughly 1 minute faster than the orbital period. From a sample of Kepler eclipsing binaries Lurie et al. (2017) noted a subset of short period, slightly supersynchronous systems. It was suggested that the slight supersynchronous nature of these systems may be due to them having a non-zero eccenticity (yet too small to be measured), which may be the case for NGTS-7A.

If we assume the starspot drift is constant with time, we can calculate the expected shift during the approximately 620 day gap between the end of NGTS and the start of TESS observations. We expect the starspot minimum to have shifted to phase 0.5 during the TESS observations. However, as seen in Fig. 4 this is where the starspot maximum occurs. This discrepancy however does not rule out the slight period difference, as the original starspot group may have decayed and been replaced by a new one at a different phase (e.g. Jackman et al., 2018).

3.8 Magnetic Activity

Along with the presence of starspots, NGTS-7A shows other clear signatures of magnetic activity. For instance, this source was originally highlighted as part of the NGTS flare survey. To find flares in the NGTS data, lightcurves are searched night by night for consecutive outliers about a set threshold. Full information about our detection method can be found in Jackman et al. (2018, 2019). From this process we identified four flares in the NGTS lightcurve and we have also identified one from our SAAO follow up lightcurve, which can be seen in Fig. 2. To calculate the flare energies we follow the method of Shibayama et al. (2013) and have assumed the flare can be modelled as a 9000 K blackbody. When calculating the flare energy, we have corrected each lightcurve for the expected dilution in the respective bandpass using our best SED fits from Sect. 3.2. From this we calculated energies ranging between $7.7^{+2.4}_{-1.8}\times 10^{32}$erg and $3.3^{+1.0}_{-0.8}\times 10^{33}$erg  or $2.5^{+0.7}_{-0.6}\times 10^{32}$erg and $1.1^{+0.3}_{-0.2}\times 10^{32}$erg for scenarios (i) and (ii) respectively. Based on the total observing time in the NGTS and I band filters, we measure the rate of flares above the minimum measured energy for NGTS-7A as $72\pm 32$ per year. The high rate of flares is similar to that of other known active M stars, such as GJ 1243 (Ramsay et al., 2013; Hawley et al., 2014) and YZ CMi (Lacy et al., 1976).

In Sect. 2.2.2 we noted the presence of emission lines from the Balmer series, helium and calcium, as shown in Fig. 6. By co-adding our HARPS spectra we were also able to identify the presence of $H\alpha$ emission. All of these emission lines are persistent, i.e. they appear in each individual spectrum, making us confident they are not just the product of a flare. \textcolorblackIn Sect. 3.1 we attributed these strong emission lines to NGTS-7A and their presence during quiescence is a clear sign that NGTS-7A is chromospherically active (e.g. Reid et al., 1995; Walkowicz & Hawley, 2009). Active M stars are known to show high energy flares more frequently than their inactive counterparts (Hawley et al., 2014), fitting in with our observation of multiple flares across datasets.

For our NGTS and SAAO data we have also checked where the flares occur in starspot phase. We find that the flares occur in the NGTS data at phases 0.42, 0.43, 0.30 and 0.56. All of these phases are when the two dominant active regions are in view. From comparing to spot modulation in the TESS lightcurve we also know that the flare observed in the SAAO follow up lightcurve occurred when the spots were in view, close to the spot modulation minimum. Previous studies of the flare-starspot phase relation for M stars have found flares appear to occur with a uniform distribution in starspot phase (e.g. Hawley et al., 2014; Doyle et al., 2018). This uniform distribution has been explained as either flares occuring in small active region, which do not not cause detectable spot modulation, or flares occuring in permanently visible active regions.

Systems with known inclinations can constrain which latitudes are permanently visible, something not known for the majority of stars. As we believe NGTS-7A has been spun up by NGTS-7Ab and the system is not inclined relative to us, the only permanently visible active regions would be at the pole. The fact that none appear when the dominant starspots are not in view suggests the flares are associated with the starspots dominating the modulation, as opposed to a permanently visible polar region or smaller spots elsewhere.

3.8.1 X-ray Activity

To determine the X-ray luminosity of NGTS-7A we have searched through available archival X-ray catalogues. NGTS-7 was detected during the Einstein 2 sigma survey conducted with the IPC instrument (Moran et al., 1996). It has an upper limit entry in the XMM upper limit server444http://xmm2.esac.esa.int/UpperLimitsServer/ (from an 8 second exposure slew observation) and was not detected in the ROSAT All-Sky Survey. For our analysis we have chosen to use the Einstein 2 sigma entry, due to it being a detection as opposed to an upper limit. Given a count excess of 8.1 counts over an exposure time of 1223 seconds, we obtain an Einstein IPC count rate of $6.6\times 10^{-3}$ counts s*-1*, with a signal-to-noise ratio of 2.35. We use the WebPIMMS interface555https://heasarc.gsfc.nasa.gov/cgi-bin/Tools/w3pimms/w3pimms.pl to calculate the flux in the 0.2-12.0 keV energy range. When doing this we use a Galactic nH column density of $1.7\times 10^{20}$ and an APEC optically-thin plasma model with $\log T=6.5$. From this we estimate an unabsorbed flux of $1.66\times 10^{-13}$ erg cm*-2s-1* between 0.2 and 12.0 keV. For our two scenarios of Sect. 3.2 we estimate $L_{\mathrm{X}}$ and $L_{\mathrm{Bol}}$ using the parameters from our best fitting SED. From this we obtain $\log L_{\mathrm{X}}=29.2$ and $\log L_{X}/L_{\mathrm{Bol}}$= -2.54 and -2.53 respectively. While these values imply NGTS-7A is more X-ray active than stars which show saturated X-ray emission (\log L_{X}/L_{\mathrm{Bol}}$$\approx -3; Pizzolato et al., 2003; Wright et al., 2018), these values are within the scatter of the Wright et al. (2011) sample. However, one has to take into account that NGTS-7B is within the Einstein IPC aperture, which has a spatial resolution of only ca. $1^{\prime}$. The detected flux may therefore stem from both stars together. If both are equally X-ray bright, this would reduce the $\log L_{X}/L_{\mathrm{Bol}}$ level for NGTS-7A to -2.84. Another possibility is that the Einstein exposure covered a flare of one of the stars, therefore registering a higher X-ray flux level compared to the quiescent level. To check for very large flares and confirm our choice of parameters in WebPIMMS we calculated the expected count rates in XMM and the ROSAT All-Sky Survey for comparison. In both cases we find find that the expected counts for the existing exposure times of XMM and ROSAT are below or at the respective upper limits. While this does not completely rule out a flare during the Einstein observation, it makes less likely; we are therefore confident that NGTS-7A is indeed an X-ray saturated star, fitting with our observations of rapid spin and magnetic activity.

4 Discussion

4.1 The nature of NGTS-7Ab

With an orbital period of 16.2 hours, NGTS-7Ab is the shortest period transiting brown dwarf around a main or pre-main sequence star to date. It is also only the fifth known brown dwarf transiting an M star (Irwin et al., 2010; Johnson et al., 2011; Gillen et al., 2017; Irwin et al., 2018). The host star is magnetically active, showing starspot modulation and flaring activity in both the NGTS and follow up lightcurves.

In Sect. 3.2.3 we formulated two possible scenarios for the nature of the NGTS-7 system. Scenario (i) places both stars at the distance implied by the Gaia DR2 parallax of NGTS-7B and results in a pre-main sequence system of roughly 55 Myr (Sect. 3.2.5), while scenario (ii) assumes both stars are on the main sequence. These two scenarios resulted in brown dwarf masses of $75.5^{+3.0}_{-13.7}$ MJ and $48.5\pm 4.3$ MJ respectively. In Sect. 3.5 we measured the rotational broadening of NGTS-7A and obtained a value of 31 km s*-1*, a value too high for a main sequence M star rotating with a period of 16.2 hours. In Sect. 3.6 we used the detection of the secondary eclipse of NGTS-7Ab to measure its temperature. We measured temperatures between 2880 K and 3200 K, depending on the geometric albedo of NGTS-7Ab. We found these measured temperatures could not be explained by a 48.5$M_{\mathrm{J}}$ brown dwarf at any age, heavily disfavouring scenario (ii) once again. Based on these pieces of evidence we conclude that scenario (i) is the most likely scenario and that NGTS-7Ab is a $55^{+80}_{-30}$ Myr brown dwarf, transiting a tidally-locked chromospherically active pre-main sequence M dwarf in a state of spin-orbit synchronisation.

4.2 Formation of NGTS-7Ab

It has previously been suggested that brown dwarfs around stars fall into two separate mass regimes (Ma & Ge, 2014; Grieves et al., 2017), with a boundary at around 42.5$M_{\mathrm{J}}$. It was suggested by Ma & Ge (2014) that the two populations of companion brown dwarfs are related to their formation mechanism. Lower mass brown dwarfs (< 42.5$M_{\mathrm{J}}$) are thought to form in the protoplanetary disc, in a mechanism similar to giant planets. Whereas, higher mass brown dwarfs \textcolorblackmay follow a formation path similar to stellar binaries and form through molecular cloud fragmentation. One reason for this separate mechanism is the limited mass available in protoplanetary discs to form companions, especially for discs around M dwarfs (e.g. Andrews et al., 2013; Ansdell et al., 2017). \textcolorblackBased on this analysis and the mass of NGTS-7Ab, we might expect molecular cloud fragmentation to be a more likely pathway for the formation of NGTS-7Ab.

If we believe that the two Gaia sources are physically bound and that this is a hierarchical triple system, then NGTS-7Ab is similar to both NLTT41135 B (Irwin et al., 2010) and LHS 6343C (Johnson et al., 2011). These systems are both M+M visual binaries where one star hosts a transiting brown dwarf. Both systems are stable with ages greater than 1 Gyr, however the presence of an outer body could help explain how NGTS-7Ab found its way onto a close orbit. \textcolorblackOne explanation for the tight orbit of NGTS-7Ab is that is has been moved inwards via the Kozai-Lidov mechanism (Kozai, 1962; Lidov, 1962), where an outer body drives periodic oscillations between the inclination and eccentricity of the inner orbit. If the brown dwarf is driven into a highly eccentric orbit it may then circularise via tidal forces, resulting in both the observed tight orbit and the spin up of NGTS-7A (e.g. Bolmont et al., 2012). From a sample of 38 high mass (> 7 $M_{\mathrm{J}}$) exoplanets and brown dwarfs, Fontanive et al. (2019) found that companions with orbital periods less than 10 days have circularisation timescales consistent with the Kozai-Lidov mechanism. We have estimated the timescale of the Kozai-Lidov mechanism ($\tau_{KL}$) for the NGTS-7 system using the formalism from Kiseleva et al. (1998) and find $\tau_{KL}$<55 Myr for initial orbits beyond 0.1 AU. For the outer orbit (of NGTS-7A and NGTS-7B) we have assumed a period of 2500 years (see Sect. 4.5) and an eccentricity of 0.5 (e.g. Raghavan et al., 2010). This timescale decreases for wider orbits. For the closer separations where the Kozai-Lidov mechanism may not have had enough time to operate, an alternative explanation may be that dynamical interactions during formation may have scattered NGTS-7Ab inwards and onto an eccentric orbit which was was then circularised through tidal forces.

On the other hand, if the two Gaia sources are not physically bound then NGTS-7Ab would be more similar to AD 3116, a relatively young (sub-Gyr) brown dwarf orbiting an M star in the Praesepe open cluster (Gillen et al., 2017). AD 3116 does not show any sign of a nearby binary component and the brown dwarf is on a 1.98 day period. In this situation, the brown dwarf companion most likely formed close enough to its host star to migrate inwards to its current position through interactions with the primary itself (e.g. Armitage & Bonnell, 2002), rather than being driven to a close orbit by a third body.

As mentioned previously, one mechanism for migrating orbiting bodies inwards is through the combination of tidal forces and the magnetic wind of the host star (e.g. Damiani & Díaz, 2016). These forces act in conjunction to migrate brown dwarfs inwards by transferring angular momentum from the orbit to the spin of the host star, which is then lost via magnetic breaking. The process acts with varying efficiency for different spectral types. These interactions have been argued to be particularly efficient for G and K stars (Guillot et al., 2014), due to their radiative interiors and moderate magnetic winds. F stars however have a much weaker wind, and the low masses and radii of M stars result in reduced tidal forces (Damiani & Díaz, 2016). Both of these factors result in increased migration timescales for F and M stars. While this depends on the initial position and age of system, these interactions could provide a feasible mechanism for moving NGTS-7Ab into its current position.

4.3 Future Evolution of NGTS-7Ab

The remaining lifetime of NGTS-7Ab will be set by the combination of tidal dissipation and magnetic braking from the spin-down of the star which together act to remove angular momentum from the orbit of the brown dwarf. In the synchronised state, the torque on the star due to the stellar wind is equal to the tidal torque (e.g. Damiani & Lanza, 2015; Damiani & Díaz, 2016) and the orbit of NGTS-7Ab is expected to decay on a timescale set primarily by the magnetic braking of the host star (e.g. Barker & Ogilvie, 2009).

To estimate the in-spiral time $\tau_{\mathrm{a}}$ of the orbit we follow Damiani & Díaz (2016) and use

[TABLE]

where $h$ is the orbital angular momentum of the system, $\alpha_{\mathrm{mb}}=1.5\times 10^{-14}$ is the magnetic braking parameter (Dobbs-Dixon et al., 2004; Damiani & Díaz, 2016), $C_{*}$ is the primary star moment of inertia and $\Omega$ is the angular velocity of the star in the synchronised state. For our two scenarios we estimate $\tau_{\mathrm{a}}$ as 5 and 10 Myr respectively, implying that NGTS-7Ab will not remain in the current state for long and is very close to the end of its lifetime.

This short remaining lifetime strengthens our conclusion in Sect. 4.1 that NGTS-7 is a young system consisting of pre-main sequence stars and a hot brown dwarf with an age of only 55 Myr.

4.4 The mass of NGTS-7A

To account for the effects of starspots on our stellar mass estimate for NGTS-7A in Sect. 3.2.5, we corrected for the expected decrease in luminosity, using the age from NGTS-7B. This was then compared directly to the unspotted PARSEC models to estimate the mass. This method assumes a limiting drop in luminosity up to 10 per cent, however it may be possible that the change is greater than this. An alternative way of accounting for starspots is to use the empirical relations of Stassun et al. (2012). These relations, from observations of low mass stars and eclipsing binaries, can be used to estimate the difference between observations and models due to magnetic activity. These corrections can be utilised with either $\log L_{\mathrm{H{\alpha}}}/L_{\mathrm{Bol}}$ or $\log L_{X}/L_{\mathrm{Bol}}$. In Sect. 3.8.1 we estimated $\log L_{X}/L_{\mathrm{Bol}}$ for the primary star, assuming both that it was the sole X-ray emitter ($\log L_{X}/L_{\mathrm{Bol}}$=-2.54) and that both stars were equally X-ray bright ($\log L_{X}/L_{\mathrm{Bol}}$=-2.84). $\log L_{X}/L_{\mathrm{Bol}}$ was calculated using the best fitting SED of the primary star alone and should thus provide a more constrained estimate of the magnetic activity. Using an average of the two values with the relations for Teff and radius of Stassun et al. (2012), we obtain correction factors of -6.5 per cent and 17 per cent respectively. Applying these correction factors and comparing the new model Teff and radius estimates to the PARSEC models, we obtain an age of 65 Myr and a mass of 0.47 $M_{\odot}$. These values are consistent with the age and mass obtained in Sect. 3.2.5, supporting our conclusion that magnetic activity (starspots) may have altered the SED of the primary star.

4.5 The orbit of the wide binary NGTS-7AB

In Sect. 3.2.2 and Sect. 3.2.1 we discussed the issues present in both the Gaia astrometry and photometry. Investigation of the scan angles used in Gaia DR2 showed over 75 per cent of scans passed through both NGTS-7A and NGTS-7B. We determined that the 1.13″ separation of the two sources was not enough to result in significant contamination of the Gaia G band photometry, however would result in blended BP and RP photometry. We determined that the close proximity may be responsible for the perturbed astrometric solution of each source, due to a shifting photocentre between scans.

Something else which has been noted as perturbing the astrometry of sources in Gaia DR2 is orbital motion. Gaia DR2 uses measurements obtained over an approximately two year timespan and orbital motion of a similar period could significantly affect the measured proper motions and parallaxes (Gaia Collaboration et al., 2018a). To see whether orbital motion could affect the astrometric solutions for NGTS-7A and NGTS-7B in a similar manner we have estimated the orbital period of the system, assuming a circular orbit. Using the masses of 0.48 $M_{\odot}$ and 0.35 $M_{\odot}$ for NGTS-7A and NGTS-7B and a separation of 173 AU we estimate the period as 2500 years. Consequently, if these sources are on a circular orbit it is unlikely orbital motion dominates the astrometric issues. We have also estimated the astrometric motion of NGTS-7A due to NGTS-7Ab to see whether this could be contributing to the astrometric noise. We estimate astrometric shifts of 0.012 and 0.013 mas for scenarios (i) and (ii) respectively, meaning it is unlikely NGTS-7Ab is causing significant astrometric noise (see also the analysis of the GJ2069 system from Mann et al., 2018).

The third data release of Gaia is planned to include information about binarity (e.g. Lindegren et al., 2018; Gaia Collaboration et al., 2018a, b), meaning we will be able to constrain these scenarios further. Along with this, it is expected that blending between close sources will be improved upon. \textcolorblackAO-assisted photometric observations could also help improve the SED fitting, \textcolorblackbetter defining the parameters of the system.

5 Conclusions

We have reported the discovery of NGTS-7Ab, a high mass transiting brown dwarf orbiting an M star with an orbital period of 16.2 hours. This is the shortest period transiting brown dwarf around a pre-main or main sequence star known to date and only the fifth brown dwarf transiting an M star host. Through the detection of starspot modulation in the NGTS data we have identified that the M star is in a state of spin-orbit synchronisation. We estimated an in-spiral time of 5 to 10 Myr. The short in-spiral time fits with the system being young and NGTS-7A being pre-main sequence M dwarf with an age of $55^{+80}_{-30}$ Myr. If so, then NGTS-7Ab has a mass of $75.5^{+3.0}_{-13.7}$ MJ, placing it at the upper end of the brown dwarf regime. Through our analysis we identified that NGTS-7A is chromospherically active, showing emission lines in spectra, strong X-ray emission and exhibiting multiple flares in our photometry. These flares appear to occur more often when the starspots are in view, suggesting the two are related.

The host star NGTS-7A has a neighbouring source, NGTS-7B, of similar brightness and proper motion and systemic velocity 1.13 arcsec away. By accounting for both stars in our SED fitting, we determined the two stars to have similar temperatures. Their very similar kinematics and close proximity on the sky strongly suggest they constitute a bound binary system. If so, we believe NGTS-7Ab is part of a hierarchical triple system and the presence of NGTS-7B may have had a role in moving the brown dwarf into its close orbit. \textcolorblackGaia DR3 and AO-assisted observations will be valuable in determining the system parameters more precisely in the future.

Acknowledgements

This research is based on data collected under the NGTS project at the ESO La Silla Paranal Observatory. The NGTS facility is funded by a consortium of institutes consisting of the University of Warwick, the University of Leicester, Queen’s University Belfast, the University of Geneva, the Deutsches Zentrum für Luft- und Raumfahrt e.V. (DLR; under the ‘Großinvestition GI-NGTS’), the University of Cambridge, together with the UK Science and Technology Facilities Council (STFC; project reference ST/M001962/1 and ST/S002642/1). JAGJ is supported by STFC PhD studentship 1763096 and would like to thank coffee. PJW, SG, TL, BTG, DP and RGW are supported by STFC consolidated grant ST/P000495/1. SLC is supported by an STFC Ernest Rutherford fellowship. MNG acknowledges support from MIT’s Kavli Institute as a Torres postdoctoral fellow. JSJ acknowledges support by Fondecyt grant 1161218 and partial support by CATA-Basal (PB06, CONICYT). CAW acknowledges support from Science and Technology Facilities Council grant ST/P000312/1. EG gratefully acknowledges support from the David and Claudia Harding Foundation in the form of a Winton Exoplanet Fellowship.

This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This paper uses observations made at the South African Astronomical Observatory (SAAO).

.1 Affiliations

1Dept. of Physics, University of Warwick, Gibbet Hill Road, Coventry CV4 7AL, UK

2Centre for Exoplanets and Habitability, University of Warwick, Gibbet Hill Road, Coventry CV4 7AL, UK

3Institute of Astronomy, University of Cambridge, Madingley Rise, Cambridge CB3 0HA, UK

4Department of Physics and Astronomy, University of Leicester, University Road, Leicester, LE1 7RH

5Astrophysics Group, Cavendish Laboratory, J.J. Thomson Avenue, Cambridge CB3 0HE, UK

6Department of Physics, and Kavli Institute for Astrophysics and Space Research,

Massachusetts Institute of Technology, Cambridge, MA 02139, USA

7Juan Carlos Torres Fellow

8Geneva Observatory, University of Geneva, Chemin des Mailettes 51, 1290 Versoix, Switzerland

9Institute of Planetary Research, German Aerospace Center, Rutherfordstrasse 2, 12489 Berlin, Germany

10Center for Astronomy and Astrophysics, TU Berlin, Hardenbergstr. 36, D-10623 Berlin, Germany

11Astrophysics Research Centre, Queen’s University Belfast, 1 University Road, Belfast BT7 1NN, UK

12Departamento de Astronomía, Universidad de Chile, Casilla 36-D, Santiago, Chile

13Centro de Astrofísica y Tecnologías Afines (CATA), Casilla 36-D, Santiago, Chile

14Instituto de Astronomía, Universidad Católica del Norte, Angamos 0610, 1270709, Antofagasta, Chile

15Leibniz Institute for Astrophysics Potsdam (AIP), An der Sternwarte 16, 14482 Potsdam, Germany

16Institute of Geological Sciences, FU Berlin, Malteserstr. 74-100, D-12249 Berlin, Germany

$\ddagger$ Winton Fellow

Appendix A SED fitting corner plot