Red Dots: A temperate 1.5 Earth-mass planet in a compact multi-terrestrial planet system around GJ1061

TL;DR

This study reports the discovery of a three-planet system around the nearby star GJ 1061, including a potentially habitable Earth-mass planet, using high-precision radial velocity measurements from the HARPS spectrograph.

Contribution

The paper presents the detection of three low-mass planets around GJ 1061, with detailed analysis confirming their stability and potential habitability, expanding knowledge of nearby multi-planet systems.

Findings

Three planet candidates with periods close to 1:2:4 ratio.

One planet lies within the star's habitable zone.

The system is dynamically stable over long timescales.

Abstract

Small low-mass stars are favourable targets for the detection of rocky habitable planets. In particular, planetary systems in the solar neighbourhood are interesting and suitable for precise characterisation. The Red Dots campaigns seek to discover rocky planets orbiting nearby low-mass stars. The 2018 campaign targeted GJ 1061, which is the 20$^{\rm th}$ nearest star to the Sun. For three consecutive months we obtained nightly, high-precision radial velocity measurements with the HARPS spectrograph. We analysed these data together with archival HARPS data. We report the detection of three planet candidates with periods of $3.204\pm 0.001$, $6.689\pm 0.005$ and $13.03\pm 0.03$ days, which is close to 1:2:4 period commensurability. After several considerations related to the properties of the noise and sampling, we conclude that a 4$^{\rm th}$ signal is most likely explained by stellar rotation, although it may be due to a planet. The proposed three-planet system (and the potential four-planet solution) is long-term dynamically stable. Planet-planet gravitational interactions are below our current detection threshold. The minimum masses of the three planets range from $1.4\pm 0.2$ to $1.8\pm 0.3$ M$_\oplus$. Planet d, with $m \sin i = 1.68\pm0.25$ M$_\oplus$, receives a similar amount of energy as Earth receives from the Sun. Consequently it lies within the liquid-water habitable zone of the star and has a similar equilibrium temperature to Earth. GJ 1061 has very similar properties to Proxima Centauri but activity indices point to lower levels of stellar activity.

Tables

Table 1. Table 1: Stellar parameters for GJ 1061.

Parameter	Value	Ref.
Alias name	L 372-58
$\alpha$	03 35 59.700	Gaia
$\delta$	$-$44 30 45.725	Gaia
${\mu }_{\alpha }\cos \delta$ [mas/yr]	$745.286\pm 0.118$	Gaia
${\mu }_{\delta }$ [mas/yr]	-$373.673\pm 0.138$	Gaia
$\pi$ [mas]	$272.2446\pm 0.0661$	Gaia
$V$ [mag]	$13.06\pm 0.07$	G14
$J$ [mag]	$7.523\pm 0.02$	2MASS
Sp. type	M5.5V	H02
$T_{\mathrm{eff}}$ [K]	$2953\pm 98$, $2999\pm 41$	S18, G14
$L$ [10${}_{}{}^{-3}L_{\odot }^{}$]	$1.7\pm 0.1$, $3$	S18, G14
$R$ [$R_{\odot }$]	$0.156\pm 0.005$, $0.19$	S18, G14
$M$ [$M_{\odot }$]	$0.12\pm 0.01$, $0.14$	S18, G14
[Fe/H] [dex]	$-0.08\pm 0.08$	N14
$\gamma$ [km/s]	1.06 $\pm$0.01	this work
$v\sin i$ [km/s]		this work
$\log L_{H\alpha }/L_{\mathrm{bol}}^{\ast }$	4.88 (inactive)	this work
$\log L_X/{10}^{-7}$ W	26.07	S04
Age [Gyr]	$>7.0\pm 0.5$	W08

Table 2. Table 2: Properties of the photometric data sets a

Data set	Season	$\Delta T$	$N_{\mathrm{obs}}$	$N_n$	rms
		[d]			[mmag]
ASH2 V	2018	97	426	32	6.8
ASH2 R	2018	97	426	32	8.6
MONET-S	2018	96	4320	43	21.2
AAVSO	2018	123	6987	103	11.8
TESS	2018	53	27602	44	1.5
MEarth T11	2017	92	1031	31	10.1
MEarth T13	2016/17	243	10607	171	8.8
ASAS-SN	2014/18	1558	1054	407	60.5

Table 3. Table 3: Model and reference statistical parameters from the simultaneous fit of HARPS pre and post upgrade data for the 3-Keplerian signals plus SHO model. Although we fit full Keplerian orbits, they are indistinguishable from circular ones (Fig.  7 - 8 ). The derived values for semi-major axis and planetary mass take the stellar mass uncertainty listed in Table  1 into account. Due to aliasing, we list two solutions for planet d (solutions with P = d 12.4 {}_{\mathrm{d}}=12.4 denoted with *. Δ ​ ln ⁡ L Δ 𝐿 \Delta\ln L indicates the fit improvement with respect to the model with one less planet. False Alarm Probabilities (FAPs) are calculated from the improvement of the log-likelihood statistic following Baluev ( 2009 ) .

Keplerian parameters	Planet b	Planet c	Planet d	Planet d*
Detection order	#1	#2	#3	#3*
$P$ [d]	$3.204\pm 0.001$	$6.689\pm 0.005$	${13.031}_{-0.032}^{+0.025}$	${12.434}_{-0.023}^{+0.031}$
$K$ [m/s]	$2.43\pm 0.24$	${2.48}_{-0.29}^{+0.28}$	$1.86\pm 0.25$
$K$* [m/s]	$2.54\pm 0.25$	$2.39\pm 0.24$		$1.76\pm 0.28$
$e$ [-]a
$\omega$ [deg]	${145}_{-65}^{+81}$	${88}_{-85}^{+95}$	${157}_{-71}^{+88}$	${278}_{-214}^{+49}$
$t_{\mathrm{peri}}$ [d]b	${0.61}_{-0.72}^{+0.56}$	${0.2}_{-1.5}^{+1.9}$	${6.4}_{-2.5}^{+3.1}$	${8.0}_{-2.4}^{+2.0}$
Derived parameters
$a$ [au]	$0.021\pm 0.001$	$0.035\pm 0.001$	$0.054\pm 0.001$	$0.052\pm 0.001$
$m\sin i$ [$M_{\oplus }$]	${1.38}_{-0.15}^{+0.16}$	$1.75\pm 0.23$	${1.68}_{-0.24}^{+0.25}$
$m\sin i$* [$M_{\oplus }$]	${1.44}_{-0.16}^{+0.17}$	$1.74\pm 0.20$		${1.57}_{-0.25}^{+0.27}$
$\lambda$ [deg]b,c	$66.0\pm 5.7$	$80.0\pm 7.6$	${335.7}_{-8.0}^{+8.4}$	$72.2\pm 8.9$
$t_c$ [d]b,d	$3.4\pm 0.1$	$6.9\pm 0.2$	$4.2\pm 0.7$	$13.1\pm 0.7$
$F$ [$S_{\oplus }$]c	$3.8\pm 0.7$	$1.4\pm 0.2$	$0.6\pm 0.1$	$0.6\pm 0.1$
HARPS pre [m/s]	${1.0}^d$
HARPS post & Red Dots [m/s]	${1.06}_{-0.17}^{+0.18}$
$\Delta \ln L$	16.8	19.1	23.3
FAP [10-5]	4.5	1.6	0.12

Table 4. Table 4: HARPS pre-upgrade, pre-Red Dots data obtained by SERVAL .

BJD	RV	${\sigma }_{\mathrm{RV}}$
	[m/s]	[m/s]
2452985.714000	-35.96	1.82
2452996.738000	-25.68	2.64
2453337.750000	-42.07	1.64
2454341.869000	-15.33	2.04
2455545.668000	-33.95	1.63
2455612.524000	-13.61	1.83
2455998.513000	-9.65	1.72

Table 5. Table 5: HARPS post-upgrade data obtained by TERRA and SERVAL (dLW). Outliers are flagged (RV or dLW) in the last column. Data are given in 1-day averages and 1 m/s has been quadratically added to the RV uncertainties. Full table in online data.

BJD	RV	${\sigma }_{\mathrm{RV}}$	dLW	${\sigma }_{\mathrm{dLW}}$	flag
	[m/s]	[m/s]	[$m^2/s^2$]	[$m^2/s^2$]
2458020.814834	1.48	1.69	-12.1	3.6
2458040.854023	-6.77	1.06	2.9	1.4
2458043.783336	2.56	1.66	0.9	1.4
2458052.765154	-1.54	1.59	2.6	1.4
2458053.866459	-6.16	1.05	7.7	1.4
2458054.699823	-7.11	1.26	-0.3	1.2

Table 6. Table 6: Photometric observing facilities.

Acronym	Location	Tel	FOV	CCD	Scale	Band(s)
		(m)	(arcmin2)		(arcsec pix${}_{}{}^{-}1$)
ASH2	SPACEOBS, Chile	0.40	54.0$\times$82.0	2.7k$\times$4k	1.20	V,R
MONET-S	SAAO, South Africa	1.20	12.6$\times$12.6	2k$\times$2k	0.37	R
AAVSO	SPACEOBS, Chile	0.40	47$\times$47	4k$\times$4k	0.69	V
TESS	TESS satellite	4$\times$0.105	4$\times$(24${}_{}{}^{\circ }\times$24∘)	4k$\times$4k	21.1	TESS
MEarth-S	CTIO, Chile	0.40	26.0$\times$26.0	2k$\times$2k	0.76	RG715
ASAS-SN	worldwide	4$\times$0.14	4.5 deg2	2k$\times$2k	8.0	V

Table 7. Table 7: Ground-based photometry, full table in online data.

BJD	rel. mag	$\sigma$ rel. mag
ASH2 R
2458309.830903	0.0062	0.0016
2458310.834277	0.0048	0.0024
ASH2 V
2458309.829251	0.0001	0.0020
2458310.835848	0.0028	0.0040
MONET-S R
2458308.622806	0.0022	0.0010
2458318.563727	0.0160	0.0011
AAVSO R
2458308.901969	0.0119	0.0017
2458309.899256	0.0137	0.0015
MEarth RG715
2457905.943122	-0.0043	0.0029
2457906.939185	-0.0067	0.0073

Table 8. Table 8: Equivalent widths obtained by TERRA in the HARPS post-upgrade, data. The flag 0 0 indicates rejected data due to possible flares. Full table in online data.

BJD	H$\alpha$	Na D1	Na D2	flag
	Å	Å	Å
2458020.814834	-0.12	1.17	0.95	1
2458040.854023	-0.16	1.17	0.90	1
2458043.783336	-0.02	1.23	0.97	1
2458052.765154	-0.03	1.18	0.93	1

Equations

$$k(\tau )=a{e}^{-c\tau }$$

$$S(\omega )=\sqrt{\frac{2}{\pi }}\frac{S_0{\omega }_0^4}{({\omega }^2-{{\omega }_0}^2{)}^2+{{\omega }_0}^2{\omega }^2/Q^2}$$

Full text

††pubyear: 2019††pagerange: Red Dots: A temperate 1.5 Earth-mass planet in a compact multi-terrestrial planet system around GJ 1061††thanks: This research has made use of the services of the ESO Science Archive Facility. Based on observations collected at the European Southern Observatory under ESO programmes 072.C-0488(E), 183.C-0437(A), 0101.C-0516(A), and 198.C-0838(A).–F

1 Introduction

The Red Dots collaboration111https://reddots.space/goals/ is an effort to detect exoplanets orbiting our nearest stellar neighbours, i.e., stars within 5 pc, by concentrating the observational efforts on one star at a time. The Red Dots observing cadence ensures the good sampling of signals only partially constrained from previous campaigns. This strategy was decisive in the discovery of Proxima b, orbiting our closest stellar neighbour Proxima Centauri (Anglada-Escudé et al., 2016); it contributed significantly to the discovery of Barnard’s Star b (Ribas et al., 2018), and revealed the paucity of terrestrial planet candidates in Barnard’s Star.

Within 5pc, there are 15 exoplanet systems known, mostly around M-stars. Three of them are high multiplicity systems: YZ Ceti (Astudillo-Defru et al., 2017c), Wolf 1061 (Wright et al., 2016; Astudillo-Defru et al., 2017b), and GJ 876 (Rivera et al., 2010). For planetary systems around low-mass stars similar to our target GJ 1061 (a late type M-star with low mass of $M=0.12$ M⊙), state-of-the-art planet formation simulations show that planets form outside the water ice-line. They migrate into warm orbits close to the central star when they attain masses of  0.5–1 ${\rm M_{\oplus}}$ (Coleman et al., 2017; Alibert & Benz, 2017). Typically, these planets migrate in resonant chains and continue to accrete solid material. This means that multi-planet systems of hot and warm Earth-mass and super-Earth planets is a common outcome, with the majority of these systems remaining in or close to resonant configurations. Observed planetary systems around stars similar to GJ 1061, such as YZ Ceti, Proxima, TRAPPIST-1 (Gillon et al., 2017; Grimm et al., 2018), as well as Teegarden’s star (Zechmeister et al., 2019) closely match the simulations from recent work which uses both pebble and planetesimal accretion scenarios to form the planetary systems (Coleman et al., 2019).

While other programmes such as HIRES/Keck (Butler et al., 2017) or the HARPS survey (e.g. Astudillo-Defru et al., 2017c) have detected numerous planets by working on large samples of stars, our strategy is different. We perform uniform, nightly cadence spectroscopic and photometric monitoring of a very small number of targets (up to three per season) to cover the periods of temperate orbits several times in each season (up to three months). This strategy has proven optimal (Anglada-Escudé et al., 2016) for detection of low-mass planets in hot to temperate orbits around nearby red-dwarfs, and for identifying and mitigating intrinsic stellar variability and spurious correlations as also demonstrated by Damasso & Del Sordo (2017).

In this paper we first introduce GJ 1061 with its basic stellar parameters (Sect. 2) followed by a presentation of the spectroscopic and photometric data (Sect. 3). In the analysis (Sect. 4), we first describe the signal detection in the radial velocity data and the detailed analysis of the signals (4.1). In order to assess their origin being due to planetary orbits or stellar activity we then analyse the spectroscopic activity indicies and the photometry (4.2). The results are finally discussed in Sect. 5.

2 GJ 1061

2.1 General properties of GJ 1061

The stellar parameters, including activity indices, are summarised in Table 1. At a distance of $\sim$3.67 pc (Gaia Collaboration et al., 2018), GJ 1061 is the 20th nearest star to the Sun222http://www.recons.org/ and has an effective temperature $T_{\rm eff}$ slightly below 3000 K (Gaidos et al., 2014; Stassun et al., 2018). Two different $T_{\rm eff}$ values are reported by Stassun et al. (2018) and Gaidos et al. (2014). While both are consistent, the stellar radius, luminosity, and mass differ between both studies. Since this affects the derived planetary parameters we explicitly indicate in Sect. 4 and 5 which parameters were used.

2.2 Activity, age, and rotation period estimates

Mid-M dwarf stars generally show significant magnetic activity (e.g., West et al., 2008; Reiners et al., 2012; Barnes et al., 2014; Jeffers et al., 2018b, and references therein). GJ 1061 is relatively magnetically inactive with little or no signs of H$\alpha$ emission (Barnes et al., 2014) and a so far undetermined rotation rate. Schmitt & Liefke (2004) reported a very low detection of X-ray emission at the level of $\log{L_{\rm X}}\sim 26.07$. Given its very low activity GJ 1061 has been used by Barnes et al. (2017) as a slowly rotating inactive reference for Doppler imaging studies.

The chromospheric Ca ii index of GJ 1061 is one of the lowest ($\log{R^{\prime}_{\rm HK}}=-$5.32) among several thousand stars investigated by Boro Saikia et al. (2018, see Fig. 3). This is in agreement with the value $\log{R^{\prime}_{\rm HK}}=-$5.754 obtained by Astudillo-Defru et al. (2017a). These $\log{R^{\prime}_{\rm HK}}$ values suggest a rotation period ranging from 50 to 130 days according to the relations by Suárez Mascareño et al. (2015) and Astudillo-Defru et al. (2017a) and are consistent with the inferred rotation period in the range of 50–200 d from the X-ray activity-rotation relation by Pizzolato et al. (2003), Wright et al. (2011), and Reiners et al. (2014). These ranges of potential rotation period are consistent with the absence of H$\alpha$ emission (Newton et al., 2017) and is consisten with $v\sin i<5$ km s*-1* reported by Barnes et al. (2014) who as well reported $-6.74\leq{\rm log}{L_{\rm H\alpha}/L_{\rm bol}}\leq 6.12$ and very low amplitude activity. From the lack of H$\alpha$ emission we can infer that the star is older than the so-called activity lifetime of an M5 or M6 dwarf, which is 7.0$\pm$0.5 Gyr (West et al., 2008). Neither circumstellar material nor stellar multiplicity has been detected for GJ 1061 (Avenhaus et al., 2012; Rodriguez et al., 2015) using the WISE and HERSCHEL satellite. Thus, though the star’s bulk properties are very similar to those of Proxima Centauri (Anglada-Escudé et al., 2016), GJ 1061 seems to be less active. These findings are further supported by the pilot precision radial velocity survey of a sample of 15 M5V - M9V stars conducted over a six night span by Barnes et al. (2014); GJ 1061 showed the lowest variability with r.m.s. variability $<5$ ms*-1*.

3 Observations

3.1 HARPS Red Dots

Observations during 54 nights were made using the HARPS spectrograph at the ESO 3.6 m at La Silla (Mayor et al., 2003) from July to September 2018 with exposure times of 1800 s and an average signal-to-noise ratio of $\sim 20$ at 650 nm. We also analyse public archival data (HARPS public333ESO programme 198.C-0838(A), PI Bonfils) treating the seven pre-fiber upgrade and 107 post-fiber upgrade data as separate sets (Lo Curto et al., 2015). All reduced spectra were processed with both TERRA (Anglada-Escudé & Butler, 2012) and SERVAL (Zechmeister et al., 2018) codes, yielding radial velocities (RVs), pseudo equivalent widths (pEW) of H$\alpha$ and Na i D, and differential line widths (dLW, see Sect. 3.2 of Zechmeister et al.). We use nightly averages in cases of multiple exposures per night resulting in 98 post-upgrade measurements.

3.2 Photometry

Contemporaneous photometry from ASH2444Astrograph for South Hemisphere II, MONET-South555MOnitoring NEtwork of Telescopes, and AAVSO666American Association of Variable Star Observers was obtained to determine the stellar rotation period of GJ 1061. Moreover, the TESS mission (Ricker et al., 2015b) observed this target in continuous mode for $\sim$50 consecutive days during our Red Dots campaign. We also used archival photometry from MEarth (Berta et al., 2012) and ASAS-SN (Kochanek et al., 2017). These data sets are summarised in Table 2. The corresponding observing facilities are described in Appendix E and Table 6, the photometric data are listed in Table 7.

All new ground-based CCD measurements were obtained by the method of synthetic aperture photometry. Each CCD frame was corrected in a standard way for bias and/or dark, and flat field by instrument specific pipelines. From a number of nearby and relatively bright stars within the frames, the best sets were selected as reference stars.

4 Data analysis

4.1 Signal detection in the radial velocities, and model parameters

The RV analysis incorporates all measurements listed in Tables 4 and 5, except for one outlier777The excluded RV point is a 5.5 $\sigma$ outlier from JD=2458376.8825. The spectrum shows no anomaly except the +20.45 m/s RV deviation. It is listed in Table 5 for completeness.. The first step of the analysis is to detect the most likely periodic signals in the data and assess their significance using likelihood periodograms (Baluev, 2009). We use an iterative sequential procedure consisting of: i) computing a likelihood periodogram to identify the period of the (next) most probable signal, ii) including the signal in the model and repeating the process until no significant signals remain. Using our favoured model described below, this procedure is illustrated in Fig. 1, showing four detected signals. Sometimes a new period is significant, but ambiguous due to aliasing. To ensure a secure and unique period determination, we additionally require that the maximum likelihood solution (highest peak) must have a $\Delta$ ln L $>$ 5 compared to the second highest peak. In a Bayesian sense, this is equivalent to requiring that the model probability (assuming uniform priors) of the best solution must be $e^{5}\sim 150$ times more probable than any other solution with an alternative period and the same number of parameters. In addition to Keplerian signals, the model also includes a zero-point offset and a white noise jitter parameter for each set, a global linear trend, and a model for the correlated noise (see below). We also used the Systemic (Meschiari et al., 2009) interface to visually verify that the solution was indeed unique.

In the second step, i.e. the detailed analysis, we deployed a python script using celerite (Foreman-Mackey et al., 2017) and emcee (Foreman-Mackey et al., 2013) to characterize the credibility intervals for all the parameters in the model. For an increasing number of planets ($N=0,1,2,3,4$) the parameters for each Keplerian model are the orbital period $P$, the radial velocity semi-amplitude $K$, the eccentricity $e$, the periastron longitude $\omega$ as well as the periastron passage time $t_{\rm peri}$888The eccentricities of the first three Keplerian signals are small. Allowing eccentric Keplerian solutions for them does not significantly improve upon circular orbit fits.. The 4th signal was alternatively modelled using the stochastically driven damped harmonic oscillator kernel (SHO) provided by celerite to account for correlations between data points using a Gaussian Process (GP) framework (see Appendix B for more details). Model parameters for this kernel are an amplitude $S_{0}$, a resonance frequency $\omega_{0}$ (we convert it in a period $P_{0}$), and the quality factor $Q$ of the oscillator converted to a damping time scale, $\tau_{\rm d}$. Like in the first step, the model also includes a zero-point offset and a white noise jitter parameter999The jitter parameter for the HARPS pre-upgrade data result in an unrealistic high value of about 10 m/s which is likely due to the fact that this data set only contains seven data points scattered over about 10 years. We therefore fixed this jitter parameter to 1 m/s. This has no impact on the overall results. for each instrument, and a global linear trend.

The first and second signals (planets b and c) are unique with no severe aliasing ambiguities. However, for the third signal (planet d), the period could be either 13.0 d or 12.4 d given that the $\Delta\ln L$ between the two solutions is only about $2.5$. In frequency space, the difference between $1/13.0\,{\rm d}^{-1}$ and $1/12.4\,{\rm d}^{-1}$ is $1/270$, which corresponds to the highest secondary lobe of the window function for this time series, thus the confusion is probably caused by aliasing. The bulk properties of the planet candidates (period and mean longitude) do not change dramatically between solutions but the amplitude and therefore mass is affected, although within error bars. Further observations can to pin-down the precise period and orbital elements. Both solutions are listed in Table 3.

A similar, but more problematic, situation occurs with the 4th signal. From the periodograms we cannot robustly distinguish between P=53 d and P=130 d because these two solutions only differ by a $\Delta\ln L$ of $\sim$3; including either in the model removes the other, indicating an aliasing between the two. Treating the 4th signal as a Keplerian orbit results in a rather high eccentricity of e=0.34. When modelled with the SHO-kernel, the fit favours a period of P${}_{0}=53$ d rather than 130 d and an unconstrained long damping time $\tau_{\rm d}>130\,d$. Since the period of the 4th signal is overlapping with the probable range for the rotation period, we discuss the interpretation of that signal as due to activity or planetary orbit in Section 4.2.

We also tested whether a $N=4$ model with $P_{4}=53$ d including an SHO kernel would result in a statistically significantly better fit. This is not the case since the log-likelihood increases only by 3. By comparison of the log-likelihood of the different models ($N=0\ldots 4$, SHO)and taking potential activity induced RV signals into consideration (see below) we identified our favoured model as the one with three planets and a SHO kernel to describe the correlated noise.

Table 3 summarizes the most probable parameter values, their credibility intervals, and the statistical significance of each signal for our favoured model with three Keplerian signals and the SHO kernel. The improvement of the fit compared to $N-1$ signals is indicated as $\Delta\ln L$. Statistical significances are given through the False Alarm Probability estimates derived from the improvement in the log-likelihood statistic (Baluev, 2009). Fig.  2 shows the best fit for all post-upgrade HARPS data. Fig.  3 shows the final phase-folding with the model specified in Table 3. The periodograms of the data with these four signals subsequently removed are shown in Fig. 1, the periodograms of the corresponding models are shown in Fig. 10. The posterior distribution of the parameters are shown in the Appendix in Figs. 6 to 8.

As a final check, we investigated the phase coherence using data from 2017 (HARPS public) and 2018 (HARPS Red Dots) separately. The first three signals are detected in both sub-sets with parameters matching the full analysis, however, larger uncertainties (about factor 10 in period and factor 2 in mass). Neither season is sufficiently long for a clear detection of the fourth signal. This is further illustrated with Fig.  11 in Appendix D. The signals from the three planets are coherent over two seasons, i.e. these detections are robust.

4.2 Activity or planets?

After the signal detection, we assess the validity of the interpretation of the signal as being due to planetary orbits or due to stellar activity. We employ spectroscopic activity indices as well as photometry in the following.

4.2.1 Spectroscopic indices

We measured pseudo-equivalent widths (EW) of the two Na i D lines and H$\alpha$ chromospheric lines as obtained with TERRA (Table 8) to search for evidence of periodicities in the stellar activity (Fig.   15). For this investigation, we excluded epochs with EW measurements showing significant deviations from the quiescent state of the chromosphere of the star (occasional small flares, cosmic rays, and background contamination can all cause spurious measurements of EW) through a two-step $3\sigma$ clipping. These excluded measurements are flagged with [math] in the online version of the data tables. The H$\alpha$ index shows a signal at about 130 d which we can well model (Fig.  14) with a Gaussian Process using the SHO kernel described above.

Changes in the width of the average spectral line are known to correlate with spurious Doppler shifts. We also investigated differential change in the line width (or dLW, Table 5) from SERVAL (Zechmeister et al., 2018). The periodogram of this index (Fig.  17) show two possible signals at 27 d and 13.5 d, but at rather low significance (FAP$\sim$1-10%). The latter is close to, but not exactly at the period of one of the RV signals, so we examined the issue further. We fitted the dLW-data with three models: a Keplerian model accounting for strictly periodic but not necessarily sinusoidal variations; and with two different kernels provided by celerite to model correlations between data points. The REAL kernel describes correlation through an exponential excitation or decay at a characteristic time scale, the SHO kernel is as described above (Sect. 4.1). The fit with a periodic Keplerian-like signal yields 26.07$\pm 0.04$ d, where the eccentricity also accounts for the signal at 13.5 d (first harmonic). However, this signal is not coherent over all data, but only for about one observing season of about 100 d length. The fit of the dLW (Fig.  16) using both kernels results very short characteristic and damping timescales (of few days), which further supports that the coherence of the signal is short lived. On the other hand, the third RV signal is very coherent throughout the season (see colour coded phase folded plots Fig. 3), so we consider it very unlikely that the dLW and RV signals are related at a significant level. The dLW signal does not match the photometric period either (Sect. 4.2.2) but it is about half of the 53 d period present in the RV data. This raises doubt on the planetary origin of the 53 d signal.

It should be noted that the dLW may also be due to instrumental or observational effects. Even when the star is very inactive, spurious signals can arise from varying Moon contamination (periodicity close to the detected $\sim 27$ d), varying background and therefore spectral line contrast, or small instrumental variations such as small focus changes.

4.2.2 Photometry analysis

TESS (Ricker et al., 2015a) observed GJ 1061 in Sector 3 and 4 (TIC79611981). With the planetary radii derived in Sect. 5 we expect transit depths between about 3 and 40 mmag, depending on the planetary composition. Using the Mikulski Archive for Space Telescopes (MAST)101010https://mast.stsci.edu/, we inspected the TESS data for transits, finding no transit signal at any of the RV periods. TESS data has a too short time-span ($\sim$ 60 days) for a detection of the potential long-term rotation modulation in GJ 1061.

We combined most ground-based time series (omitting ASAS-SN due to large rms and ASH2 V) that we could find on the star, including the new Red Dots photometric data obtained (quasi-)simultaneously to the HARPS observations (Fig.  12). The combined analysis using an SHO kernel as well as individual offsets and jitter terms shows several strong and clear signals in the range of 60–150 d (Fig.  13, top). We are aware of the fact that a combination of ground-based photometry using different filters as well as comparison stars is prone to artefacts. We therefore also investigate the data sets individually. ASH2 R as well as AAVSO show a long-term periodicity which is however unconstrained due to insufficient duration. In the MEarth T13 data set periodicity at 130$\pm 5$ d is seen most clearly. The 130 d signal is not coherent over the full photometric data set; this is consistent with an origin in relatively short-lived active regions on the slowly rotating stellar surface. The MEarth T13 data set also shows power in the 5 d and 2.5 d range at low amplitude. These, however, disappear in the residuals (Fig.  13, bottom).

Using the SHO kernel, we estimate a damping time $\tau_{\rm d}\sim 20$ d, i.e. significantly shorter than the underlying oscillator period of $\sim$130 d further supporting this picture. Alternatively, the REAL kernel yields $\tau_{\rm d}$ $\sim 44$ d. Correlation kernels following an over-damped oscillator produce signals in a broad period range as discussed in Ribas et al. (2018)111111As a word of caution we like to note that the parameters of the GP models may not have a physical meaning but due to instrumental systematics or due to the combination of inhomogeneous data like in this case..

We therefore assessed the chances of the 4th signal being a false positive induced by correlated noise with the characteristics of those of the star activity traced by the photometry. Since the period of 130 d being detectable in the RV, photometric, and Hα activity index data, we concentrate on the 53 d signal in the following. We generated 50 000 realisations of correlated noise using the REAL kernel with a $\tau_{\rm d}\sim 44$ d at the sampling of the RV data. For each of those we obtained the highest peak from a GLS periodogram (Zechmeister & Kürster, 2009). As expected, we find a broad distribution of periods among those peaks with a broad maximum in the range 100 to 400 d. A signal with a period of between 51 and 57 d (the full range of periods from the $N=4$ planet fit) is found with 16 % chance. For the shorter orbital periods a false-alarm probability is negligible. While a planetary origin for the 53 d signal cannot be ruled out, we consider stellar activity is the most likely explanation for this 4th signal.

Summarizing, both, the 53 d and the 130 d signal are within the expected range for the stellar rotation period as mentioned in Sect. 2.2. Both have possible counter parts in the activity indices (twice the period of the dLW) or photometry. The common detection of the 130 d signal in RV, photometric, and activity index data justify to identify the 130 d period as the stellar rotation period.

4.3 A stable compact planetary system

We tested the dynamical properties of the model using mercury6 (Chambers, 1999) and Systemic, and assuming four planets as the dynamically most challenging case. Consistent with expectations based on separations compared to mutual Hill radii, and the period commensurabilities (close to 1:2:4:16), the system survives at least $10^{8}$ years. This is, however, requiring small eccentricities for planets b to d, especially planet c needs an eccentricity near zero in order to keep the system long-term stable. We found that, when adopting the minimum planet masses, the system is not in a mean-motion resonance since the libration angles rotate. Significant exchanges of angular momentum are nevertheless present and manifest as regularly oscillating eccentricities. A dynamical fit results in a very subtle shift in orbital periods within the error bars which, over the duration of observation, causes a slight phase shift without corresponding to a statistically better model. The planetary system remains stable even if viewed almost face-on, at an inclination of $i=10^{\circ}$. This dynamical robustness precludes strong constraints on the system’s inclination. For randomly oriented systems, geometrical arguments disfavour small values of $i$ and thus masses larger than 3-4 times the reported minimum masses. Probabilistic arguments based on the Kepler mission’s detection statistics and comparison to RV detected samples (e.g. Bixel & Apai, 2017) makes it a priori very unlikely that these planets have masses larger than 3 ${\rm M_{\oplus}}$.

5 Discussion and conclusions

We find three RV signals that can be attributed to planets with minimum masses close to the mass of the Earth matching the expectation from state-of-the-art planet formation simulations (Coleman et al., 2017; Alibert & Benz, 2017). All three signals are coherent over the observational time span $\sim$1 year. A fourth signal is also significant, but considering all available data we find that it is mostly likely caused by correlated noise induced by stellar activity and rotation; we find $P_{\rm rot}\sim 130\,$d, since this signal is detected in the RV, photometric as well as Hα index data. Whether the signal at 53 d is only due to correlated noise or planetary origin would require longer-term observations with dense coverage over about five months in order to check its coherence as well as amplitude stability.

Using the stellar parameters (Table 1) and planetary parameters (Table 3) we show the three planets in relation to the star’s habitable zone (Fig. 4), assuming the minimum mass as the actual mass. Fig. 4 shows their insolation, that of other Earth-like planets, and the limits of the optimistic and conservative habitable zone (Kopparapu et al., 2014). The two sets of stellar parameters result in different insolation values for the planets. We propagate the uncertainties in $T_{\rm eff}$ into the uncertainties in insolation values. With an orbital period of around 13 d, planet d is the most interesting regarding potential habitability as it lies well within the classic liquid water habitable zone for both sets of stellar parameters. The range of possible sizes for planet d are shown in Fig. 5 for various planetary composition ranging from a cold H+He to a pure Fe composition (Zeng et al., 2016).

This study illustrates the efficiency of the Red Dots observing strategy. Our homogeneously sampled data over one season has been essential to the robust detection of the three short period planets. Additional planets in longer-period orbits may also be present in the system. The non-zero linear trend of 1.8$\pm 0.4$ m s*-1* yr*-1*(Fig. 9) is indicative of a possible long-period companion, but a much longer observational baseline is needed to explore this. The sequence of period ratios close to 2 for the inner three planets may hint at an additional undiscovered planet exterior to planet d. Sub-Earth mass planets in periods longer than 20 d are below the current detection threshold, but could be revealed with further high precision RV measurements.

Acknowledgements

S.D., S.J., A.R., L.S., and M.Z. acknowledge support from the Deutsche Forschungsgemeinschaft under Research Unit FOR2544 “Blue Planets around Red Stars”, project no. DR 281/32-1. S.J. acknowledges project JE 701/3-1 and DFG priority program SPP 1992 “Exploring the Diversity of Extrasolar Planets” (RE 1664/18). E.R., C.R.-L., M.J.L.-G. and N.M. acknowledge financial support from the Spanish Agencia Estatal de Investigación through projects AYA2016-79425-C3-3-P, ESP2017-87676-C5-2-R, ESP2017-87143-R and the Centre of Excellence “Severo Ochoa” Instituto de Astrofísica de Andalucía (SEV-2017-0709). J.R.B. and C.A.H. are supported by STFC under grant ST/P000584/1. Z.M.B acknowledges CONICYT-FONDECYT/Chile Postdoctorado 3180405. E.P. and I.R. acknowledge support from the Spanish Ministry for Science, Innovation and Universities (MCIU) and the Fondo Europeo de Desarrollo Regional (FEDER) through grants ESP2016-80435-C2-1-R and ESP2016-80435-C2-2-R, as well as the support of the Generalitat de Catalunya/CERCA programme. G.A.-E. research is supported by STFC Consolidated Grant ST/P000592/1.

Data were partly obtained with the MONET/South telescope of the MOnitoring NEtwork of Telescopes, funded by the Alfried Krupp von Bohlen und Halbach Foundation, Essen, and operated by the Georg-August-Universität Göttingen, the McDonald Observatory of the University of Texas at Austin, and the South African Astronomical Observatory. Data were partly collected with the robotic 40-cm telescope ASH2 at the SPACEOBS observatory (San Pedro de Atacama, Chile) operated by the Instituto de Astrofífica de Andalucía (IAA). We acknowledge the effort of various pro-am to obtain photometric monitoring during the Red Dots 2018 campaign through the AAVSO.

Appendix A Radial Velocity data

In Table 4 and 5 we list the all radial velocity data for GJ 1061. Radial velocities are obtained using SERVAL and TERRA. Note that the full Table 5 is available online.

Appendix B Modeling details

We start our python-based detailed analysis using scipy.optimize.minimize. For each Keplerian signal we use the parameters provided from the GLS analysis as starting values (amplitude, period, periastron time). The best fit and its uncertainties are used to then initialize the MCMC run with 500 walkers and 3000 steps. As priors for all parameters, we use a normal distribution with the best fit as mean and 100 times the uncertainty as the standard deviation. After the run we check that the parameter distribution was not affected by the initialisation. Boundaries are set only for numerical reasons, e.g. to avoid numerical problems with too large or too small values for the kernel parameters. To avoid non-physical parameter combinations, the absolute value of the eccentricity is used when calculating the Keplerian models. The absolute value of the velocity amplitude is used in order to avoid degenerate solutions. In case the final best solution from the MCMC chains is significantly better than the starting parameters, this procedure is repeated with the new best parameter set.

The correlation of the data is modeled using the Gaussian Process Regression framework. We use the celerite package (Foreman-Mackey et al., 2017) since it only scales linearly with the number of data points and not cubical. The REAL kernel of celerite represents an exponential decaying coherence of the data and has two parameters, the amplitude $a$ and the inverse damping time scale $c$ (to avoid negative values for $a$ and $c$, the logarithm of the values are used):

[TABLE]

The SHO kernel represents a stochastically driven, damped harmonic oscillator. It is defined through its power spectral density (PSD) using three parameters, an amplitude $S_{0}$, the quality factor $Q$ describing the damping and the oscillator frequency $\omega_{0}$. Like in the REAL kernel, the logarithms of the parameters are used. The quality factor is related to the inverse of the damping time. The PSD of this term is:

[TABLE]

Appendix C MCMC corner plots

We show the posterior distribution from 500 walkers and 30 000 steps (the first half rejected as burn-in) are as corner plots using corner (Foreman-Mackey, 2016) in Figs. 6 to 9, first for the full fit parameter set and then parameter subsets for the individual planets (Figs. 7 and 8) as well as for instrumental parameters (Fig.  9). For the parameter subsets we show absolute values for the eccentricities, the longitude of periastron and the periastron passage time are converted modulo $2\pi$ and orbital period, respectively. The subset corner plots also contain derived parameters, i.e. the semi-major axis, the planetary minimum mass as well as the longitude of periastron. For the former two we propagate the uncertainties in the stellar mass. The 4th signal is modeled using the SHO-kernel.

Appendix D Periodograms

The subsequent subtraction of Keplerian signals from the radial velocity data is shown in Fig.1 as GLS-periodograms (Zechmeister & Kürster, 2009). In addition, the periodograms of the signals are shown in Fig.10.

We also use the Stacked Bayesian Generalized Lomb-Scargle (SBGLS) periodogram (Mortier et al., 2015; Mortier & Collier Cameron, 2017) in order to check the development of the signals with increasing number of data points (Fig.  11). The three planetary signals are indeed stable in period and their probability is increasing with the number of data points as expected for coherent signals.

Appendix E Photometry

The used photometric facilities are summarized in Table 6. They are:

ASH2. The ASH2 (Astrograph for South Hemisphere II) telescope) is a robotic 40 cm telescope with a CCD camera STL11000 2.7k$\times$4k, FOV 54$\times$82 arcmin2. The telescope is at SPACEOBS121212http://www.spaceobs.com/en (San Pedro de Atacama Celestial Explorations Observatory), at 2450 m above the sea, located in the northern Atacama Desert in Chile. During the present work, only subframes with 40% of the total FOV were used, resulting in a useful FOV of 21.6$\times$32.8 arcmin2. Data are obtained as differential magnitude relative to a main reference star. A set of five check stars were selected and used for checking purposes. In total, 32 epochs in V and R bands were obtained in the period between July and October, 2018, with a time span of 97 days.

MONET-S. The 1.2 m MONET-South131313https://monet.uni-goettingen.de/ telescope (MOnitoring NEtwork of Telescopes) is located at the SAAO (South African Astronomical Observatory), South Africa. It is equipped with a Finger Lakes ProLine 2k$\times$2k e2v CCD, FOV 12.6$\times$12.6 arcmin2. We performed aperture photometry with the AstroImageJ package using a set of about five comparison stars. 43 epochs in the R band were collected during the Red Dots campaign.

AAVSO. Within the AAVSO collaboration, the Remote Observatory Atacama Desert (ROAD) (Hambsch, 2012) was used. This is also located at SPACEOBS. It is equipped with a 40 cm f/6.8 Optimized Dall-Kirkham (ODK) reflector from Orion Optics, UK, and a CCD camera ML16803 4k$\times$4k from FLI, USA, equipped with Astrodon UBVRI photometric filters. Data were reduced using the LesvePhotometry package with two comparison stars, one for reference and another one for checking purposes. These observations were also collected simultaneously with the Red Dots campaign.

TESS. At the end of Red Dots campaign, GJ 1061 was also observed by the TESS141414http://heasarc.gsfc.nasa.gov/docs/tess/ satellite (Ricker et al., 2015b). This was carried out during two consecutive sectors of 27 days each, with a effective total time span of about 53 days between September and November, 2018, with the main aim of investigating possible transit signals.

MEarth-S. The MEarth project151515https://www.cfa.harvard.edu/MEarth/ consists of two robotically controlled observatories dedicated to monitoring thousands of M-dwarf stars (Berta et al., 2012): MEarth-N (North), at Fred Lawrence Whipple Observatory (FLWO), USA, since 2008 and MEarth-S (South), at Cerro Tololo Inter-American Observatory (CTIO), Chile, since 2014. Each array consists of eight identical 40 cm robotic telescopes (f/9 Ritchey-Chrétien Cassegrain), each equipped with a 2k$\times$2k CCD camera, FOV 26$\times$26 arcmin2, sensitive to red optical and near-infrared light. MEarth project generally uses an RG715 long-pass filter, except for the 2010-2011 season when an I715-895 interference filter was choosen. In the case of GJ 1061, two time series are available from the MEarth-S project, which have been collected with telescopes number 11 (T11) (year 2017) and 13 (T13) (2016-2017).

ASAS-SN. The ASAS-SN project161616http://www.astronomy.ohio-state.edu/asassn/ (All-Sky Automated Survey for Supernovae) (Kochanek et al., 2017) currently consists of 24 telescopes distributed on 6 units around the Globe, located in Hawaii (Haleakala Observatory), Chile (CTIO, 2 units), Texas (McDonald Observatory), South Africa (SAAO) and China. Each station consists of four 14 cm aperture Nikon telephoto lenses, each with a thermo-electrically cooled, back-illuminated, 2k$\times$2k, Finger Lakes Instruments, ProLine CCD camera. The field of view of each camera is roughly 4.5 deg2 with a pixel scale of 8.0 arcsec.

The combined photometric data are listed in Table 7 and shown in Fig.  12 together with a model using the SHO-kernel from celerite. The signal at about 130 d is well described with this quasi-periodic oscillations as visible from the periodograms in Fig.  13. The analysis is presented in the main text in Sect. 4.2.2.

Appendix F Activity indicators

TERRA and SERVAL provide various spectroscopic activity indicators, among those the H$\alpha$ as well as NaD equivalent width and the differential line width. These measurements are listed in Tables 8 and 5. The spectroscopic activity indicators have been investigated for signals with periods at the potential planetary periods to avoid false positive planet detections as well as for the detection of the stellar rotation. In Fig.  16 we show the fit of the differential line width, in Fig.  14 for the Hα equivalent width. The periodograms of the data, signals and residuals are shown in Figs. 17 and 15. The analysis is presented in the main text in Sect. 4.2.1.