Georgina Dransfield

Amaury H M J Triaud

Tristan Guillot

Djamel Mekarnia

David Nesvorný

Nicolas Crouzet

Lyu Abe

Karim Agabi

Marco Buttu

Juan Cabrera

Davide Gandolfi

Maximilian N Günther

Florian Rodler

François-Xavier Schmider

Philippe Stee

Olga Suarez

Karen A Collins

Martín Dévora-Pajares

Steve B Howell

Elisabeth C Matthews

Matthew R Standing

Keivan G Stassun

Chris Stockdale

Samuel N Quinn

Carl Ziegler

Ian J M Crossfield

Jack J Lissauer

Andrew W Mann

Rachel Matson

Joshua Schlieder

George Zhou

HD 28109 hosts a trio of transiting Neptunian planets including a near

Keywords: planets and satellites, detection -planets and satellites, dynamical evolution and stability -planets and satellites, fundamental parameters

come

INTRODUCTION

The discovery of over 4 500 extrasolar planets has in no way diminished our curiosity regarding our place in the universe; if anything, we have more questions than ever.

Following early results by Doppler surveys (e.g. Mayor et al. 2011;Howard et al. 2010b), the Kepler survey (Borucki et al. 2010) revealed an abundant population of planets comparable in size to Neptune (Fressin et al. 2013), the majority occupying periods shorter than a few hundred days (Howard et al. 2010a). The formation of these Neptune-size worlds remains debated; one possible pathway is that 'pebbles' (large dust particles) separate from the gas in the planet-forming disc and drift inward (Johansen & Lambrechts 2017). Another possibility is that the pebbles are accreted into cores before migrating inward (Coleman & Nelson 2016); thus both models predict that most of the mass growth takes place in situ (e.g. Lambrechts et al. 2019;Bean et al. 2021).

Both the 'drift' and 'migration' formation pathways predict the formation of resonant chains, where adjacent planets find themselves in a mean-motion resonance with their neighbours (Bean et al. 2021). However most planets are not found with commensurate orbital periods (Fabrycky et al. 2014). Studies show that in ∼95% of cases the resonance chain becomes unstable following the dispersal of the gas disc (Matsumura et al. 2010;Pu & Wu 2015;Izidoro et al. 2021). Alternatively, the migratory pathway might be more chaotic (Rein 2012;Lee & Chiang 2015).

Statistical studies of the large samples of known planets have shown that multiplicity is common, i.e. many planets are in multiplanet systems (Fabrycky et al. 2014). Those planets that find themselves close to a mean-motion resonance with their neighbours allow for masses to be constrained without spectroscopic follow-up, by monitoring transit timing variations (TTVs) in the system (e.g. Agol et al. 2005;Holman & Murray 2005;Holman et al. 2010;Lithwick et al. 2012).

Neptune-sized planets orbiting bright stars provide excellent dynamical laboratories to further investigate these exciting systems. TESS, the Transiting Exoplanet Survey Satellite (Ricker et al. 2015) launched in 2018 to search for exoplanet candidates; one of its primary mission aims is to deliver to the community 50 planets smaller than 4 R ⊕ , with masses measured by high precision radial velocity (RV). Not all systems are suitable for RV mass measurements, but it is predicted that TESS will also observe ∼ 30 systems with TTVs, with approximately a third of those being suitable for dynamical mass measurements (Hadden et al. 2019). Another of TESS's goals is the identification of planets orbiting bright stars, to permit detailed investigations into the dynamics and chemical composition of exoplanetary atmospheres. The discovery and analysis of HD 28109 (TOI-282) fulfils two of those primary objectives.

Stars within TESS's northern and southern continuous viewing zones (CVZ) are continuously observed for a year, allowing for the discovery of long period planets (> 20 days). These planets are likely to be favourable candidates for JWST as its CVZ falls within TESS's; however they are challenging to confirm with ground-based photometry, especially if TTVs introduce large timing uncertainties.

Unless, of course, one is observing from Antarctica. ASTEP (Antarctic Search for Transiting ExoPlanets Guillot et al. 2015;[START_REF] Mékarnia | [END_REF]) is a 40 cm telescope located on the Antarctic Plateau at a latitude of -75.1 deg an elevation of 3, 233 m. This unique location has a thin atmosphere, completely dry air, and virtually uninterrupted observing between late May and late July (Crouzet et al. 2010a,b). The visibility covers TESS's CVZ and the area where Low-Earth Orbit telescopes have trouble reaching (e.g. CHEOPS Broeg et al. 2013). For this reason, ASTEP is the observatory (ground or space) best suited for the photometric confirmation and follow-up of long-period transiting systems.

In this work we present the discovery, validation and characterisation of three sub-Neptunes orbiting HD 28109. We begin by characterising the host in Section 2, followed by Sections 3 and 4 where we describe the identification and validation of planetary candidates in the system. In Section 5 we describe our global analysis of all available photometric and radial velocity data to characterise the planets. We then place the three planets in the context of the field, including potential for atmospheric characterisation of the system and further high precision radial velocity follow-up in Section 6. Finally, we summarise our findings and conclude in Section 7.

STELLAR CHARACTERISATION

HD 28109 (TOI-282) is a bright (V=9.42; J = 8.476) main-sequence star, of spectral type F8/G0V. Its right ascension and declination are 04:20:57.19 -68:06:09.68, and it has a parallax of 7.13 mas (Gaia Collaboration et al. 2021, 2016), placing it at a distance of 140 pc (Bailer-Jones et al. 2021).

All planetary measurements are derived from our knowledge of the star; we therefore begin by characterising the radius, mass, effective temperature and spectral type of the star.

Reconnaissance Spectroscopy

We observed TOI-282 on three nights spanning 2019 February 14 to 2021 January 07 with the CHIRON spectrograph on the 1.5 m SMARTS telescope. CHIRON is a high resolution echelle spectrograph fed by an image slicer and a fiber bundle, located at Cerro Tololo Inter-American Observatory (CTIO), Chile. The spectra have a resolution of R = 80000 with a wavelength coverage of 4500 -8900 Å (Tokovinin et al. 2013). The wavelength solution is provided by bracketing Thorium-Argon cathode-ray lamp observations, and the spectra are extracted and wavelength calibrated with the official CHIRON pipeline (Paredes et al. 2021). We extracted the RVs by fitting the spectral line profiles, which were measured via least-squares deconvolution of the observed spectra using synthetic templates (Donati et al. 1997). The three RVs exhibit a root mean square scatter of 58 𝑚/𝑠, which is not significant given the 44 𝑚/𝑠 mean per-point uncertainty. With these data, we cannot place strong limits on the presence of additional companions, but the two-year time span of the observations does allow us to measure long-term trends that might be indicative of very massive outer companions. The best-fit linear trend (0.076 ± 0.091 𝑚/𝑠/𝑑𝑎𝑦) is consistent with zero, which suggests it is unlikely that there is a stellar companion in the system, but the precise radial velocities described in Section 4.3 can provide tighter constraints on the masses of the transiting planets and limits on the presence of other bound companions.

We also use the CHIRON spectra to determine the effective temperature, surface gravity, and metallicity of the host star by matching against observed spectra that have previously been classified using SPC (Buchhave et al. 2012). Interpolation to the final parameters is performed with a gradient-boosting regressor implemented in the python module. We measure 𝑇 eff = 6120 ± 50 K, log 𝑔 ★ = 4.13 ± 0.10, and [Fe/H] = 0.0 ± 0.10. Following [START_REF] Gray | The Observation and Analysis of Stellar Photospheres[END_REF] and Zhou et al. (2018), we derive 𝑣 sin 𝑖 ★ by fitting broadening kernels to the instrumental, macroturbulent, and rotational line profiles, and estimate 𝑣 sin 𝑖 ★ = 7.7 ± 0.5 𝑘𝑚 𝑠 -1 .

Spectral Energy Distribution

As an independent check on the derived stellar parameters, and in order to determine an estimate for stellar age, we perform an analysis of the broadband spectral energy distribution (SED). Together with the Gaia EDR3 parallax, we determine an empirical measurement of the stellar radius following the procedures described in Stassun & Torres (2016); Stassun et al. (2017Stassun et al. ( , 2018)). We pulled the 𝐵 𝑇 𝑉 𝑇 magnitudes from Tycho-2, the 𝐽𝐻𝐾 𝑆 magnitudes from 2MASS, the W1-W4 magnitudes from WISE, and the 𝐺𝐺 RP 𝐺 BP magnitudes from Gaia. In addition, we pulled the FUV and NUV fluxes from GALEX in order to assess the level of chromospheric activity, if any. Together, the available photometry spans the full stellar SED over the wavelength range 0.2-22 𝜇m (see Figure 1).

We perform a fit using Kurucz stellar atmosphere models (Kurucz 1979), with the effective temperature (𝑇 eff ) and metallicity ([Fe/H]) adopted from the spectroscopic analysis. The remaining free parameter is the extinction, 𝐴 𝑉 , which was limited to the maximum line-of-sight value from the dust maps of Schlegel et al. (1998). The resulting fit is very good (Figure 1) with a best-fit 𝐴 𝑉 = 0.01 ± 0.01 and a reduced 𝜒 2 of 2.8 (excluding the GALEX FUV flux, which is consistent with a modest level of chromospheric activity; see below). Integrating the (unreddened) model SED gives the bolometric flux at Earth of 𝐹 bol = 4.37 ± 0.15 × 10 -9 erg s -1 cm -2 . Taking the 𝐹 bol and 𝑇 eff together with the Gaia EDR3 parallax, with no systematic offset applied (see, e.g., Stassun & Torres 2021), gives the stellar radius as 𝑅 ★ = 1.446 ± 0.035 R . Finally, estimating the stellar mass from the empirical relations of Torres et al. (2010) and a 6% error from the empirical relation itself gives 𝑀 ★ = 1.26 ± 0.08 M , consistent with the mass estimated empirically from the stellar radius together with the spectroscopic log 𝑔 which gives 𝑀 ★ = 1.03 ± 0.24 M .

We can also estimate the stellar age by taking advantage of the observed modest chromospheric activity in the UV together with empirical age-activity-rotation relations. For example, we can estimate the log 𝑅 HK = -4.56 ± 0.05 via the GALEX FUV excess and the empirical relations of Findeisen et al. (2011). That implies an age of 𝜏 = 0.9 ± 0.3 Gyr via the empirical activity-age relations of Mamajek & Hillenbrand (2008). Finally, we can further corroborate the activity-based age estimate by also using empirical relations to predict the stellar rotation period from the activity. For example, the empirical relation between 𝑅 HK and rotation period from Mamajek & Hillenbrand (2008) predicts a rotation period for this star of 7.2 ± 1.1 d, which is compatible with the (projected) rotation period 𝑃 rot /sin 𝑖 ★ = 9.51 ± 0.62 d inferred from the stellar radius above together with the spectroscopic 𝑣 sin 𝑖 ★ .

We present the adopted stellar parameters in Table 1 3 IDENTIFICATION OF PLANETARY CANDIDATES HD 28109 (TOI-282) is in the southern continuous viewing zone (CVZ) for TESS which means it is observed in all southern sectors, the only exception being sector 32 where the target was off the edge of the CCD by 3 pixels. At the time of writing, the star has been observed in 25 sectors at two-minute cadence.

In the following sections we describe the identification of candidates by TESS, as well as our own search for further candidates in the data.

Note on Nomenclature

Throughout Section 3 we refer to the host (HD 28109) using its TESS alias: TOI-282, where TOI is TESS Object of Interest. When referring to candidate planets, we also use the TESS nomenclature and add numeric suffixes in the order that candidates were identified. Once the candidates are confirmed we change their names. Thus, TOI-282.01 becomes HD 28109 c; TOI-282.03 becomes HD 28109 d; and TOI-282.04 is HD 28109 b, with letters indicating the order from closest to farthest from the host star.

TESS data

All 2 minute cadence TESS data is reduced by the SPOC (Science Processing Operations Center) pipeline [START_REF] Jenkins | Society of Photo-Optical Instrumentation Engineers[END_REF]) and the photometry data products are delivered in the form of simple aperture photometry (SAP) or Presearch Data Conditioning Simple Aperture Photometry (PDCSAP), where the latter has been corrected for instrument systematics. The lightcurves are then searched for transit-like signals; candidates that have a signal-to-noise ratio ≥ 7 are reported as threshold crossing events (TCEs).

The first TCE on TOI-282 was detected in Sector 1, and a second was detected in Sector 2. The multi-sector data validation report for Sectors 1-2 showed this as a planetary candidate with a period of ∼18 days. It wasn't until Sectors 1-6 were available that the individual candidates could be disentangled. The 1-6 multi-sector validation report reported three planet candidates on TOI-282: candidate .01 with a period of ≈56 days, .02 with a period of ≈31 days, and finally .03 with a period of ≈84 days. Candidate .02 was later retired as a false alarm1 .

The PDCSAP lightcurves for Sectors 1-13 and Sectors 27-39 can be found in Figure 2, where transits of TOI-282.01 are highlighted in yellow, while the transits of TOI-282.03 are in green.

Search For Additional Transiting Candidates

We made use of the custom S2 pipeline presented in Pozuelos et al. ( 2020) and Demory et al. (2020) to perform a search for additional transiting candidates in the TESS data. S uses the L package (Lightkurve Collaboration et al. 2018) to download TESS PDCSAP data from NASA Mikulski Archive for Space Telescope (MAST); any outliers, defined as any points over 3𝜎 above the running mean, are then removed. The pipeline then uses W (Hippke et al. 2019) to detrend the data using the biweight method testing several window lengths. In this case, we used six window sizes between 0.64 and 3.55 days in order to optimise signal detection efficiency (SDE). To aid in the detection of low SNR signals, S also optionally applies a a Savitzky-Golay filter (Savitzky & Golay 1964) to smooth the data and increase precision.

The search for periodic signals is carried out using the T L S (TLS) package (Hippke & Heller 2019) as it is optimised to search for periodic signals with transit-like shapes. We searched a wide parameter space, with periods ranging from 5 to 100 days setting the minimum signal-to-noise ratio (SNR) to 5. This is in contrast to the SPOC pipeline, which sets the minimum SNR to 7.

We recover both candidates .01 and .03 in the first instance and find that the SDE is consistent across all window sizes, although the SNR and transit depths are at their highest using a window size of 2.8259 days. We also find a third significant signal at a period of 22.89±0.01 days. This signal has an SNR of 23.8 and SDE of 13.1 using the same window-size as for the first two candidates; its depth of just 0.14ppt (parts per thousand) would make it extremely challenging to detect from the ground. Nevertheless, if this is indeed a third planet in the system, its size would be classed as a mini-Neptune just beyond the so-called radius valley (Fulton et al. 2017;Van Eylen et al. 2018), and its period would place it within ∼ 5% of a 5:2 third order resonance with candidate .01. The TESS lightcurve folded on this signal is presented in Figure 3, along with the Lomb-Scargle periodogram.

While this signal does not match any of the TOIs on ExoFOP, further scrutiny reveals that this might have been the period originally associated with TOI-282.02. Dietrich & Apai (2020) include the candidates of TOI-282 in their sample, and cite the orbital period for the innermost candidate as ∼ 22.89 days. Additionally, when we scrutinise some of the first multi-sector data validation reports produced by the SPOC pipeline, we find that Candidate 2 has a period of ∼ 22.89 days 3We now adopt this re-identified candidate as TOI-282.04 and indicate in dark pink the detected transits on Figure 2.

Two other periodic signals with low significance were also recovered which did not pass our vetting tests and are therefore attributed to instrument systematics.

VETTING AND VALIDATION OF CANDIDATES

In this section we describe first the ground-based follow-up observations we conducted, staring with high-resolution imaging of the host to search for blended companions. We then describe the photometric observations conducted from several sites between January 2019 and August 2021. In the first instance, the purpose of ground-based follow-up is to validate the planetary nature of the candidates; in Section 4.4 we describe how we used the available data on this system to rule out all feasible false positive scenarios.

We note that while the validation tests described here allowed us to confirm the planetary nature of TOI-282.01, .03 and .04 initially, the analysis we describe in Section 5.2 reveals anti-correlated transit timing variations between planets .01 and .03. Therefore, many of the tests here described are not necessary to validate these two objects as planets since TTVs are not known to have astrophysical false-positives. We nevertheless kept the validation steps here for completeness.

High Resolution Imaging

If an exoplanet star has a close companion (bound or line of sight), "third-light" flux from the companion can lead to an underestimated planetary radius if not accounted for in the transit modelling (Ciardi et al. 2015;Furlan & Howell 2017;Matson et al. 2018).

Gemini South Telescope

To search for close-in (bound) companions unresolved in other follow-up observations, we obtained high-resolution imaging observations from Gemini South's Zorro speckle interferometric instrument4 . TOI-282 was observed using Zorro on UT 2020 December 23. Zorro provides simultaneous high-resolution speckle imaging in two optical bands, 562/54 and 832/40 nm, with output data products including a reconstructed image, and robust limits on companion detections (Howell et al. 2011(Howell et al. , 2021)). Figure 4 (upper panel) shows our 5𝜎 detection limit contrast curves and the corresponding reconstructed speckle image in 832 nm. We find that TOI-282 is a single star with no companions detected in our observations down to a contrast level of 5 to 8 magnitudes from the Gemini 8-m diffraction limit out to 1.2". At the distance of TOI-282 (d=140 pc) these angular limits correspond to spatial limits of 2.8 to 168 AU.

Southern Astrophysical Research Telescope

We also searched for stellar companions to TOI-282 with speckle imaging on the 4.1-m Southern Astrophysical Research (SOAR) telescope (Tokovinin 2018) on UT 2019 February 2, observing in Cousins I-band, a similar visible bandpass as TESS. This observation was sensitive to a 7.3-magnitude fainter star at an angular distance of 1 arcsec from the target. More details of the observation are available in Ziegler et al. (2020). The 5𝜎 detection sensitivity and speckle auto-correlation functions from the observations are shown in Figure 4 (middle panel). No nearby stars were detected within 3 of TOI-282 in the SOAR observations.

VLT NaCo

We collected high-resolution AO images of TOI-282 with VLT/NaCo [START_REF] Lenzen | Instrument Design and Performance for Optical/Infrared Ground-based Telescopes[END_REF][START_REF] Rousset | Society of Photo-Optical Instrumentation Engineers[END_REF] on UT 2019 January 25, using the 𝐾 𝑠 filter. We collected 9 frames, each with exposure time 20s, and dithered the telescope by 2" in a grid-like pattern between each frame, and constructed a sky background frame by median combining the science images. We processed the data using a custom pipeline which corrects bad pixels, subtracts the sky frame and applies a flat-field correction, and then aligns the stellar position between each frame and co-adds the images. To test the sensitivity of these images, we injected faint fake PSFs into the data, and scaled them such that they could be redetected at 5𝜎. Significances are averaged radially and presented in Figure 4 (lower panel); the data are sensitive to companions 5mag fainter than the star beyond 230mas, and to companions 7mag fainter than the star in the back-ground limited regime. We searched for companions in the reduced images by eye. Detector persistence at the dither positions causes a faint point source 2" directly to the south of the companion, but by inspecting the individual images we confirm that this is not a true visual companion. Apart from this persistence, no point sources are seen anywhere in the field of view, and we confirm that the star is single to the limit of our resolution.

Photometric Follow-up

Both candidates identified by TESS (TOI-282.01 and .03) have long periods (∼ 56 d and ∼ 84 d), long transit durations (> 8 h), and their transits are shallow. These constraints combined made TOI-282 a challenging system to follow-up from the ground, accounting for why the .01 only has seven photometric follow-up observations to date with three detections, while the .03 candidate has two observations with no transit detections. Below we describe first the observations carried out by Hazelwood, followed by the observations taken by ASTEP and LCO.

Hazelwood

The Hazelwood Observatory is a backyard observatory with a 0.32 m Planewave CDK telescope working at f/8, a SBIG STT3200 2.2k × 1.5k CCD, giving a 20 × 13 field of view and 0.55 per pixel. The camera is equipped with 𝐵, 𝑉, 𝑅𝑐, 𝐼𝑐, 𝑔 , 𝑟 , 𝑖 and 𝑧 filters. Typical FWHM is 2.2" to 2.7". The Hazelwood Observatory, operated by Chris Stockdale in Victoria, Australia, observed an ingress of TOI-282.01 in Ic on UT 2019 January 21; the transit was below the detection threshold, resulting in a flat lightcurve, but a nearby eclipsing binary (NEB) check cleared the field within 2.5 .

Coincidentally, this observation spanned a rare double transit that included the innermost candidate (TOI-282.04) as well. The successful NEB check therefore also applies to this candidate, ruling out nearby stars as the source of the event.

ASTEP

ASTEP is a custom 0.4 m Newtonian telescope equipped with a 5-lens Wynne coma corrector and a 4k × 4k front-illuminated FLI Proline KAF-16801E CCD. The camera has an image scale of 0."93 pixel -1 resulting in a 1 𝑜 × 1 𝑜 corrected field of view. The focal instrument dichroic plate splits the beam into a blue wavelength channel for guiding, and a non-filtered red science channel roughly matching an 𝑅 𝑐 transmission curve (Abe et al. 2013;Guillot et al. 2015). The telescope is automated or remotely operated when needed. Due to the extremely low data transmission rate at the Concordia Station, the data are processed on-site using an automated IDL-based aperture photometry pipeline [START_REF] Mékarnia | [END_REF]). The calibrated light curve is reported via email and the raw light curves of about 1,000 stars of the field are transferred to Europe on a server in Roma, Italy and are then available for deeper analysis. These data files contain each star's flux computed through 10 fixed circular apertures radii, so that optimal light curves can be extracted.

The ASTEP location at Dome C, on the Antarctic plateau, means that it enjoys excellent photometric conditions and exceptionally long nights (Crouzet et al. 2018). ASTEP is therefore well placed to observe targets with long orbital periods and transit durations.

We observed a transit of the candidate TOI-282.01 with ASTEP on the night of UT 2020 June 08 following our interest to confirm planets with long transits that we were particularly well able to observe during the austral winter. TOI-282.01 was identified by searching the TESS T F [START_REF] Jensen | Tapir: A web interface for transit/eclipse observability[END_REF]. Our detection of the full transit, the first of its kind as part of the photometric follow-up effort of TFOP (TESS Follow-up Observing Program) Sub-Group 1, prompted us to follow the system and look more closely at the available TESS data.

LCO

We observed TOI-282 in Pan-STARRS 𝑧-short band on UTC 2021 August 30 from the Las Cumbres Observatory Global Telescope (LCOGT; Brown et al. 2013) 1.0 m network node at Siding Spring Observatory. The 1 m telescopes are equipped with 4096 × 4096 SINISTRO cameras having an image scale of 0. 389 per pixel, resulting in a 26 × 26 field of view. The images were calibrated by the standard LCOGT BANZAI pipeline [START_REF] Mccully | Proc. SPIE. p[END_REF], and photometric data were extracted using A I J (Collins et al. 2017). The target star was saturated to provide higher photometric precision for fainter stars nears within 2. 5 of the target star. The observations were scheduled to include ingress coverage for a transit of TOI-282.03 with more than ±3𝜎 ephemeris coverage according to a public ephemeris from the SPOC pipeline. However, according to the linear ephemeris extracted in Section 5.1 of this work, the observations covered from the time of nominal ingress to 90 minutes later. According to the TTV analysis of Section 5.2, the ingress would have happened 45 minutes earlier than the linear ephemeris from this work on this epoch, resulting in all in-transit coverage. Although these data do not rule out nearby eclipsing binaries due to the ingress timing, we do not find an obvious deep event in any star within 2. 5 of the target star that could be the source of the TESS signal.

Radial Velocity Follow-up

We acquired 8 high-resolution spectra (R≈140 000) with the Echelle SPectrograph for Rocky Exoplanets and Stable Spectroscopic Observations (ESPRESSO; Pepe et al. 2021) on the 8.2 m ESO Very Large Telescope (VLT; Paranal observatory, Chile). The observations were performed in service mode between UT 2019 October 15 and UT 2019 December 27 as part of our program 0104.C-0003 (PI: Rodler). The exposure time was set to 1115 s, leading to a mean signal-to-noise (S/N) ratio of ∼160 per pixel at 550 nm. We reduced the ESPRESSO data and extracted the radial velocity measurements by employing ESO's ESPRESSO data reduction pipeline, version 2.2.15 , reaching a mean RV precision of 1.1 m/s.

We also gathered 7 high-resolution spectra (R≈115 000) with the High Accuracy Radial Velocity Planet Searcher (HARPS; Mayor et al. 2003) spectrograph mounted at the 3.6 m ESO telescope (La Silla observatory, Chile), as part of ESO programs 60.A-9700 and 60.A-9707. The spectra were collected between UT 2019 May 24 and UT 2021 January 29, setting the exposure time to 1800 s, which led to a mean S/N ratio of ∼100 per pixel at 550 nm. We reduced the spectra using the dedicated HARPS data reduction software (DRS) and extracted the radial velocity (RV) measurements by cross-correlating the Echelle spectra with a G2 numerical mask (Baranne et al. 1996;Pepe et al. 2002;Lovis & Pepe 2007). The resulting radial velocity measurements have a mean precision of 3.58 m/s

Validation

In this section we make use of data validation reports combined with our ground-based follow-up observations to validate the planetary nature of the three candidates orbiting TOI-282. The tests we describe here are standard procedure for the discovery of TESS planets, but we note again that for the two outer candidates the presence of anticorrelated TTVs render these steps somewhat redundant.

Vetting Tests

The TESS data validation report contains several tests to check the validity of the candidates (Twicken et al. 2018). These tests include checking for a difference in depth between odd and even transits, the presence of a potential secondary eclipse implying possible binarity, and any centroid offsets to ensure the event is on target. Both TOI-282.01 and .03 passed these tests successfully and were given false alarm probabilities of 2.14 × 10 -162 and 2.21 × 10 -34 respectively, indicating that these events are highly unlikely to be instrumental in nature. Having passed these tests, these candidates were ripe for follow-up observations to validate their planetary nature.

The false alarm probability for TOI-282.04 was calculated by S as 8 × 10 -5 ; S also provides a vetting stage for promising signals, where Field Of View (FOV) plots are generated through the internal usage of TPFP (Figure 5) (Aller et al. 2020) to check for possible contamination sources; it also injects the search results into L (Eisner et al. 2020), which is used as the main vetting engine. S also provides folded lightcurves to check for even/odd transits and discard the detection of harmonics or sub-harmonics of the real signal. We scrutinise the data validation report produced by L and note no centroid offsets at the time of transits, nor were there any sharp changes in background flux.

Ruling Out False Positives

Nearby eclipsing binaries (NEB) scenarios were ruled out for all three candidates by the photometric observations, as described in Section 4.2.1. Subsequent observations of TOI-282.01 by ASTEP later further confirmed the events on target with clear detections of the transits.

Additionally, our high-resolution observations described in Section 4.1 rule out any blended bound companions down to a projected orbital separation of 2.8AU. A 0.1 M companion at this orbital separation would produce an RV variation of ∼ 1.6 km/s. The very low scatter in our radial velocity measurements described in Section 4.3 rules out the presence of any bound stellar-mass companions in the system.

Statistical Validation with S and T

S uses the statistical validation package T presented in Giacalone et al. (2021) to calculate the false positive probability (FPP) for each of the candidates. T queries the TESS Input Catalog (TIC) for stars within 10 pixels of the target and calculates the probability that the transits detected are caused by false positive scenarios. In all, 18 scenarios are tested, including that the signal is caused by a transiting planet on the target star with the input orbital period. S automates the preparation stage for T by preparing the lightcurves and apertures; the statistical validation is then carried out five times to produce mean probabilities for each scenario.

The criteria for statistical validation of a planetary candidate is stated as FPP < 0.015; we compute FPPs of 0.0075, 0.035 and 0.203 for TOI-282.01, .03 and .04 respectively. This immediately places TOI-282.01 well below the threshold for validation, although the probability remains too high for TOI-282.03 and .04.

However, we are able to add further to this body of evidence by returning to the fact that this is a multi-planet system. Lissauer et al. (2012) derived an expression for a 'multiplicity boost', where the presence of additional planets in the same aperture increases the probability that the signals are planetary in nature:

𝑃 3 ≈ 50 𝑃 1 50 𝑃 1 + (1 -𝑃 1 ) (1)
Where 𝑃 1 is the probability of planethood without taking the additional planets into account, and P 3 if the probability with the 3-planet multiplicity boost. This decreased the FPP for TOI-282.01, .03 and .04 to 0.00015, 0.0007 and 0.005 respectively. All three planets now meet the criteria for statistical validation showing that overwhelmingly, the planet hypothesis is the most favoured. We note that while the multiplicity boost described about was developed for Kepler candidates, a similar boost for TESS candidates is described in [START_REF] Guerrero | [END_REF] which for small planets (R < 6R ⊕ ) decreases the FPP by a factor of ≈ 54. We choose to apply the boost described in Lissauer et al. (2012) as it is more conservative.

Validation Conclusions

Having ruled out all astrophysical false positive scenarios, calculated false positive probabilities and taken into account multiplicity we can now consider the system validated.

From this point on we refer to the three validated planets by their formal names in order of increasing semi-major axis:

ANALYSIS

With the planetary nature of the signals validated, we now present our analysis of all available data for HD 28109 b, c and d. We begin by describing our global analysis of the available photometry for the system in Section 5.1, followed by our dynamical analysis in Section 5.2. In Section 5.3 we make use of priors from our photometric and dynamical fits to fit the available radial velocity data.

Global Photometric Analysis

We carried out a global photometric analysis of the datasets described in Section 4 using A (Günther & Daylan 2021, 2019). A is a publicly available versatile inference package capable of jointly fitting photometric and radial velocity datasets from different instruments using (Maxted 2016) for lightcurve and RV models, (Foreman-Mackey et al. 2013) for MCMC, D (Speagle 2020) for dynamic nested sampling, and (Foreman-Mackey et al. 2017) for Gaussian Process (GP) models. We chose to use this package as it can simultaneously fit for white (uncorrelated) noise and red (correlated) noise for different instruments, while offering a wide selection of models to choose from and compare evidence.

We adopted the signal parameters found in Section 3.3 as uniform priors, and the stellar parameters derived in Section 2 as Gaussian priors and we fit using the nested sampling algorithm. Nested sampling [START_REF] Skilling | Bayesian Inference and Maximum Entropy Methods in Science and Engineering: 24th International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering[END_REF]) works by drawing a number of live points from the prior distribution, and then removing the point with the lowest likelihood. Another point is then drawn from the prior while requiring that the overall likelihood now be higher than before. This process then repeats until the change in Bayesian evidence, parametrised as Δ ln 𝑍, falls below a certain threshold. We begin the algorithm with 1500 live points and set the threshold to Δ ln 𝑍 < 0.01 as recommended by the A documentation. We computed quadratic limb darkening coefficients 𝑢 1 and 𝑢 2 for HD 28109 in TESS and 𝑅 𝑐 filters using P (Parviainen & Aigrain 2015) and the PHOENIX stellar atmosphere library provided by (Husser et al. 2013). The resulting coefficients are presented in Table 3; these are also input to our global fit as Gaussian priors.

We initially fit for a two-planet model using only HD 28109 c and HD 28109 d, and then for a three-planet model incorporating HD 28109 b. One of the key advantages of nested sampling over MCMC (Markov Chain Monte Carlo) is the calculation of the Bayesian evidence at each stage of the algorithm; this allows us to determine which model is statistically favoured by calculating the Bayes' factor (Kass & Raftery 1995). In both cases, we fit simultaneously for all transit parameters, quadratic stellar limb-darkening coefficients in the Kipping (2013) parametrisation, as well as the baseline model and the white noise scaling for each instrument. These latter two ensure that uncertainties arising from instrument systematics are appropriately propagated to derived physical parameters.

In Section 3.3 we found that all three planets had the highest SNR in TESS when the flux was detrended; to detrend the data for semiperiodic stellar variability (red noise) present in all instruments, we simultaneously fit a GP using a Simple Harmonic Oscillating (SHO) kernel. The parametrisation used in the fit is summarised in Table 3.

We find that the 3-planet model is vastly statistically favoured over the 2-planet model, with a Bayes Factor in excess of 60 000. A 𝐵𝐹 > 150 is usually considered strong enough statistical evidence to confirm a more complex model (Kass & Raftery 1995). The fit's result further confirms the existence of HD 28109 b in our data.

In both the 2-and 3-planet model fits, eccentricity and argument of periastron were set as free parameters. We note that the resulting derived orbital parameters suggest that the two outer planets could have small to moderate eccentricities (𝑒 c = 0.116 +0.30 -0.079 and 𝑒 d = 0.19 +0.34 -0.15 ) and yield calculated host densities within 1𝜎 of the prior, while the inner planet yields a very large eccentricity of

𝜌 ★,obs ≡ 3𝜋(𝑎/𝑅 ★ ) 3 𝑜𝑏𝑠 𝐺𝑃 2 (2)
where 𝜌 ★,obs is the host density derived from transit parameters assuming a Keplerian circular orbit. The manner and extent of the discrepancy between 𝜌 ★,𝑜𝑏𝑠 and the external prior on stellar density can be indicative of erroneous assumptions made regarding the orbits (Kipping 2014). In particular, for such an eccentric orbit we would expect the ratio 𝜌 𝑜𝑏𝑠 /𝜌 𝑡𝑟 𝑢𝑒 to be much greater than unity.

We therefore repeat the 2-and 3-planet fits constraining eccentricities to zero and then compare the evidence for each of the models to see which is statistically most likely. We find that all 3-planet models are preferred over all 2-planet models, with Bayes Factors in excess of 50 000. We also find that all circular models are preferred over models with non-zero eccentricities. In the case of the 3-planet circular model vs. the 3-planet eccentric model we find a Bayes Factor of almost 3 000, indicating that a more complex lightcurve modelling with eccentric orbits is not preferred.

The results of the global fit, including derived parameters, are presented in Table 3. In Figure 6 we present the three individual ASTEP transits for planet c along with the best fitting model, and in Figure 7 we present the TESS phase folded lightcurves for each planet along with the best fitting model.

Dynamical Analysis

In this section we consider a non-Keplerian model for the planets of HD 28109; we first fit individual transit times and then use these timings to infer dynamical masses for the planets.

Transit Timing Variations

Given the proximity of the two outer planets to a first order 3:2 meanmotion resonance (within ∼ 0.3%), we would expect these planets to experience some mutual gravitational influence leading to transit timing variations (TTVs). This is supported by visual inspection of the individual transit fits resulting from the nested sampling, as several appear to have noticeable offsets in time between the model and the data.

We therefore re-fit the photometry once again using A , this time allowing the midtime of each individual transit to vary about a linear ephemeris for all three planets. Fixed Lin. Ephem. Fit ln 𝜎 jitter;espresso (ln km/s) -5.99 +0.38 -0.36 gp ln Q(flux) -4.9 +1.9 -2.2

Fixed Lin. Ephem. Fit gp ln 𝜔 0 (flux) 6.0 +2.2 -1.9

Fixed Lin. Ephem. Fit

Derived Parameters

Companion radius; R p (R ⊕ ) 2.199 When allowing for TTVs, we fit for all the same parameters as in the linear fit, except that we now fix the linear ephemeris using the values derived from the most statistically favoured fit in Section 5.1. We also constrain eccentricities to zero for all three planets and once again apply Gaussian priors on the stellar density and quadratic limb darkening coefficients.

In preparation for a TTV fit, A attempts to guess the location of the flux minimum during each transit present in the data, with these guesses then used as uniform priors for the fit. However, we find that almost all of the guesses for planet b are affected by local stellar variability or nearby transits of other planets biasing the fit; instead, for each transit's mid-time, we impose a uniform prior defined as the linear predictions ±120 minutes for all planets. Once again we simultaneously fit a GP with a SHO kernel to account for stellar variability.

The results of this fit are presented in Figure 8. HD 28109 c and d present very clear and significant anti-correlated TTVs with peak to peak amplitudes of ∼ 50 mins and ∼ 100 mins respectively. While the shifts in timing do appear to show a sinusoidal variation between early and late, it is most likely that what we see here is the shorter term 'chopping' part of the signal (Deck & Agol 2015) given that the 'Super-Period' for these planets is expected to be in excess of 9 000 days.

For HD 28109 b the errors on the timings are very large due to how shallow individual events are (≈ 0.19ppt) and therefore our results are consistent with no TTVs. We note that for some individual events the errors were significantly smaller than average. Visual inspection of each transit fit showed that this is probably caused by artifacts in the lightcurves. We therefore amplified the errors on all timings of planet b to be at least the mean timing uncertainty.

Mass Estimates From TTVs

We perform dynamical analysis of the observed transit times with a symplectic 𝑁-body integrator code described in Nesvorný et al. (2012). We instruct the code to simultaneously fit all measured transit times, using MULTINEST (Feroz & Hobson 2008) to perform the regression, setting the integration time-step to 0.5 d, roughly 1/50 of the planet b's orbital period. We consider dynamical models with 2, 3 and 4 planets. The model with 3 planets has 20 parameters: planet-to-star mass ratios, orbital periods, eccentricities, longitudes of periapsis, impact parameters, difference in nodal longitudes, and reference epochs between a reference time (2458337 BJD UTC) and the first observed transit of each planet. The nodal longitude of the innermost planet in each fit is fixed at 270 • , and therefore the 2 and 4 planet models have 13 and 27 parameters, respectively. We use uniform priors for all parameters except for the impact parameter, since our light-curve fits provide strong constraints which are used here. The host star mass was fixed at 𝑀 ★ = 1.26 M . The uncertainties discussed below therefore do not account for the uncertainty in the stellar mass.

We find that the existing measurements are not good enough to uniquely determine the planet properties. The first problem arises because there is a large radial separation between the inner planet b and the two outer planets: the measured TTVs of c and d thus most likely reflect their mutual interaction (rather than those of planet b). In addition, given the relatively large measurement errors of b's TTVs it is not clear whether b's TTVs contain any useful information about planets c and d. Thus, as b can be decoupled from the rest of the system it is difficult to meaningfully constrain its properties. Indeed, when we instruct MULTINEST to perform a 3-planet fit, the mass of planet b always reaches to the upper edge of the mass prior even if this is clearly nonphysical (i.e., mass exceeding that of Jupiter for the estimated radius 𝑅 𝑏 = 2.199 𝑅 ⊕ ). We therefore limit the uniform prior on the mass of planet b to M b = 11M ⊕ .

All fits performed here confirm that the masses of c and d are in the planetary range. Our best fit indicates m c /M * = (1.9 +1.0 -0.7 ) × 10 -5 and m d /M * = (1.4 +0.6 -0.5 ) × 10 -5 , suggesting that both planets have low densities. The orbital eccentricities of c and d are most likely low (e < 0.1), but some 2-planet fits identify modes with eccentric orbits (e > 0.1). In Table 3 we present the 2𝜎 upper limits for the eccentricities of all three planets as this parameter remains poorly constrained; the mean values of the posterior distributions for the eccentricities are e b = 0.1519 +0.0894 -0.0944 , e c = 0.0391 +0.0406 -0.0231 , and e d = 0.0238 +0.0313 -0.0157 . The posteriors for the longitude of periastron are not included in Table 3 as they are unconstrained by our fits.

Given the wide orbital separation between planets b and c, we perform an additional fit including a putative fourth planet. This planet may be non-transiting or too small to be detected. We used a uniform prior between 28 d to 44 d for the 4th planet orbital period. The 4-planet fits give masses of c and d that are consistent with the masses reported above (and low orbital eccentricities). The best parameters for the hypothetical 4th planet are m 4 /M * 7.6 × 10 -6 and P 4 39.3 d. With this orbit, the two inner planets would be wide of the 3:2 resonance. The Bayesian evidence for this fit is significantly lower than for the 3-planet fit, indicating that the 3planet fit is preferred.

The eccentricities of planets c and d are low and when that is the case, the resonant and near-resonant dynamics can be studied analytically following Nesvorný & Vokrouhlický (2016). There are three variables to consider, each of them being a combination of orbital elements. Constant 𝛿 is an orbital invariant that defines the position of the system relative to the 3:2 resonance, the resonant angle 𝜓 is a function of the usual 3:2 resonant angles 𝜎 1 and 𝜎 2 , and variable Ψ is a combination of planetary masses, semi-major axes and eccentricities (see Nesvorný & Vokrouhlický (2016) for details). The resonant librations of 𝜓 can only happen for 𝛿 > 0.945. Figure 9 shows position of planets c and d in the terms of the dynamical variables 𝛿, 𝜓, and Ψ. About two thirds of the posterior samples are non-resonant with the two orbits being slightly wide of the exact resonance. About one third of the sample, however, show resonant librations with a generally low libration amplitude. Additional observations will be needed to better constrain the location of the c and d planets relative to the 3:2 resonance.

We predict that the amplitudes of the TTVs will be very large (∼ 1000 min for planet c and ∼ 2000 min for planet d) once one whole TTV period has been sampled. In Figure 10 we present the predicted TTVs of the system for the next 11 000 days. Continued monitoring of the system will allow us to both improve timing predictions and refine the dynamical mass measurements of planets c and d.

Radial Velocity Analysis

The global photometric analysis yielded strong priors on all orbital parameters; additionally, we can place priors on the masses of the two outer planets from the dynamical fit. We therefore fit the radial velocity data set to gain an estimate on the mass of the innermost planet using A . We adopt Gaussian priors on the all orbital parameters for the three planets (P, T 0 , , R p /R * , (R * +R p )/a, cosi,), as well as the radial velocity semi-amplitudes of the outer planets (K).

We place a wide uniform prior between 0.1-6.0 m/s (corresponding to 0.5 -30 M ⊕ ) on the semi-amplitude of the inner planet and fit for all three periods, reference times and semi-amplitudes. In our models, we vary the RV offset and RV jitter for each individual dataset. The jitter parameter appears as an additional RV noise term added in quadrature to the nominal RV uncertainties and subsums any instrumental or stellar effects (such as activity) not considered by the formal RV errors.

In Figure 11 we present the phased radial velocity curves. The radial velocity semi-amplitude posterior for the inner planet yields a mass estimate of M b = 18.50 +9.12 -7.61 M ⊕ ; for the outer planets we find mass estimates consistent with the estimates from the dynamical analysis: M c = 7.99 +2.13 -2.03 M ⊕ and M d = 5.53 +0.98 -0.94 M ⊕ . We do however note that to precisely constrain the masses of the system far more measurements are needed, and therefore this result is just a first estimate.

SYSTEM DISCUSSION

HD 28109 is host to at least three planets, with two very close to a first order mean-motion resonance. So far, this is the brightest TESS star, and the third brightest star overall, known to host TTV planets. In Figure 12 we present the current sample of planets known to show transit timing variations in multi-planetary systems, plotted as planet radius and orbital period vs. host V magnitude, highlighting the position of the three planets of the HD 28109 system 6 . Planets plotted as dark grey circles have measured masses (by RV or TTV) while grey circles have masses estimated by mass-radius relations. The only host brighter than HD 28109 known to host TTV planets with measured masses is WASP-18, but as the outer planet is non-transiting full characterisation of the system is not yet possible (Pearson 2019).

In this section we discuss the system in the context of other known exoplanetary systems, and current and near-future instrumentation capabilities.

The planets of HD 28109 on the Mass-Radius Diagram

In Figure 13 we present a Mass-Radius diagram showing the positions of HD 28109 b, c, and d in relation to other known exoplanets with masses measured using TTVs or precise radial velocities. The mean densities of Earth, Neptune, Saturn and pure iron are also shown.

The masses of HD 28109 c and d estimated from TTVs indicate that they are under-dense compared with most other sub-Neptunian mass planets. Given the paucity of planets in this size and mass regime with precisely constrained parameters at orbital separations comparable to HD 28109 b and c, we cannot determine at this time whether or not their densities are truly anomalous.

The mass estimate of HD 28109 b is rather more puzzling. First estimates from radial velocities in this work suggest that this planet is significantly denser than typical planets of its size. Atmospheric escape in sub-Neptunes is thought to be driven by photoevaporation (Owen & Wu 2013) and/or core-powered mass loss (Ginzburg et al. 2018). At its orbital separation, the gravity of planet b would be sufficient to retain a H/He dominated envelope if stellar irradiation is the main driver [START_REF] Mordasini | [END_REF]. Should the large mass be 6 Data for this plot was retrieved from the NASA Exoplanet Archive (Akeson et al. 2013) on 2021 October 20: https: //exoplanetarchive.ipac.caltech.edu/cgi-bin/TblView/ nph-tblView?app=ExoTbls&config=PSCompPars. confirmed with radial velocities, it could favour core-powered mass loss as the main driver atmospheric escape in this scenario, as the heat source is the planet's own residual formation energy (Gupta & Schlichting 2020).

Given the two outer planets' relatively large orbital radii (> 50 days), it is unlikely that they suffered major mass loss (Gupta & Schlichting 2020). Preliminary models of the contraction of the planetary envelope based on the approach used by Guillot et al. (2006) and Valencia et al. (2010) and neglecting the core luminosity indicate that matching their mass, radius and age require a hydrogen-helium envelope which consists of 19% of the planetary radius for planet c and 10% for planet d. Such a relatively significant envelope mass fraction may result from a formation at large orbital distances (Lee & Chiang 2016). These values are to be considered with caution however because including core luminosity and tidal heating that may be caused by the resonances would yield smaller envelopes (Millholland 2019).

Planets with vastly different densities within the same system, which is likely to be the case for TOI-282, provide invaluable laboratories to explore the various outcomes of planet formation with the same conditions.

Prospects For Precise Mass Measurements With Radial Velocities

Stellar activity in HD 28109 poses a challenge on RV follow-up observations given the small RV semi-amplitudes of these planets.

In our RV analysis, we find that the RV jitter terms in the HARPS and ESPRESSO measurements are on the order of 2.4 m/s and 0.4 m/s , respectively (Table 3). Given the TTV masses of planets c and d, their expected RV semi-amplitudes are K c ≈ 1.14 m/s and K d ≈ 0.71 m/s. For planet b, the range is broad, spanning from the Chen & Kipping (2017) estimate to the mass of a pure iron planet, K ≈ 0.91 -2.4 m/s. Given the length of the period of the outermost planet, sampling the full orbit will be a considerable undertaking.

With the goal to further characterize the system, we are currently conducting a radial velocity survey with ESPRESSO; once concluded, the results will be presented in a follow-up paper. These high-precision radial velocity measurements will help to disentangle the eccentricity of planet b, as well as shed light on the possible presence of additional non-transiting planets in the system. In addition, these measurements will finally put precise constraints on the masses of these planets, which will be crucial for atmospheric follow-up observations.

Dynamical Exploration of Planet b's Parameters

Our dynamical and photometric fits provided us with TTV mass estimates for the two outermost planets, along with strong evidence for circular orbits. However, the mass and eccentricity of planet b remain largely unconstrained, and a dynamical analysis might be able to restrict plausible parameters of planet b. The first 3-planet photometric fit in Section 5.1 found an eccentricity of ∼0.8, while the dynamical fits presented in Section 5.2.2 favour an eccentricity of ∼0.17.

We use the Mean Exponential Growth factor of Nearby Orbits chaos index (𝑌 (𝑡), MEGNO; Cincotta & Simó 2000) to assess the stability of the HD 28109 system with the full range of eccentricity and mass estimates for planet b. MEGNO is used to evaluate whether a body's trajectory will be stable following perturbations of its initial conditions. The time-averaged MEGNO 𝑌 (𝑡) will tend to 2 for 𝑡 → ∞ if the motion of the body is quasi-periodic, while chaotic behaviour will cause 𝑌 (𝑡) to tend to infinity instead.

We made a 50 × 50 mass-eccentricity grid for planet b with eccentricity values evenly spaced between 0 -0.8, and masses evenly spaced between 5 -30 M ⊕ . We used the MEGNO implementation of the N-body integrator R (Rein & Liu 2012) which in turn uses the Wisdom-Holman WHF code (Rein & Tamayo 2015). We used steps of 1/50th of the period of planet b (0.5 days) and integrated for 10 000 orbits of the outermost planet (∼ 2 300 years) to construct a two-dimensional MEGNO map. We find that the full mass range is non-chaotic but only for eccentricities smaller than ∼0.4; eccentricities up to ∼0.5 can be stable but only for the largest masses in the range we tested. This indicates that a wide range of mass values are possible for planet b, in line with the wide mass prior used for our dynamical and RV fits.

Prospects For Atmospheric Characterisation with HST and JWST

To assess the suitability of HD 28109 b, c and d for atmospheric characterisation using transmission and emission spectroscopy, we calculated for each planet the Transmission and Emission Spectroscopy Metrics following Kempton et al. (2018). The TSMs for each planet We produced model transmission spectra for all three planets using E _T (Kempton et al. 2017) assuming solar metallicity and C/O ratio for the host, and cloud-free H 2 dominated atmospheres for the planets. Using these spectra we then used P E (Batalha et al. 2017) to simulate Hubble Space Telescope (HST hereafter) observations of their atmospheres using transmission spectroscopy. Our simulations indicate that for the two outermost planets a single HST transit would be sufficient to carry out reconnaissance spectroscopy (as in de Wit et al. 2018de Wit et al. , 2016)), to assess the presence of clouds in their atmospheres. In the event of a cloud-free atmosphere, we could detect the presence of water in transmission to 8𝜎 and 3.5𝜎 respectively. To reach a 3.5𝜎 detection for planet b could require a large number of transits, but with its mass unconstrained for the time being we cannot make a fair estimate. As an example, should the mass be consistent with the Chen & Kipping (2017) estimate, we would require 6 transits, each spanning 8 orbits of HST to reach a 3.5𝜎 detection. We present the model spectra along with simulated HST observations in Figure 14.

We do also note that the simulations described above assume the TTV masses for planets c and d, and the mass estimated from massradius relations for planet b. Should any mass be significantly larger the scale height of the atmosphere will decrease, leading to lower significance detections.

CONCLUSIONS

In this work we have presented the discovery and confirmation of three sub-Neptune sized planets orbiting HD 28109. Given the brightness of the host, the planets are well-suited to precise radial velocity measurements, which will allow comparison of TTV masses with RV masses. This system is also the brightest TESS system discovered to date known to host TTV, making it an exquisite laboratory for planet formation theory. Additionally, we have shown that the outer two planets are well suited to atmospheric reconnaissance with current instrumentation, and they will be high-priority planets for follow-up with JWST and Ariel.

This work has also demonstrated the capabilities of ASTEP as a observatory specialising in planets with long periods, very long transits, and large TTVs. The uninterrupted night enjoyed from Antarctica during two months of the year enables observations of systems that cannot be feasibly followed from anywhere else in the world. Continued selection of these candidates for follow-up will help the community fill out the most sparsely populated parameter spaces in the ever-growing field of exoplanets.

DATA AVAILABILITY

TESS data products are available via the MAST portal at https://mast.stsci.edu/portal/Mashup/Clients/ Mast/Portal.html. Follow-up photometry and high resolution imaging data for HD 28109 are available on Exo-FOP at https://exofop.ipac.caltech.edu/tess/target. php?id=29781292. These data are freely accessible to ExoFOP members immediately and are publicly available following a oneyear proprietary period. The radial velocity data used in this paper are available via the ESO Public Archive at http://archive.eso. org/eso/eso_archive_main.html.

Figure 1 .

 1 Figure 1. Spectral energy distribution (SED) of TOI-282. Red symbols represent the observed photometric measurements, where the horizontal bars represent the effective width of the passband. Blue symbols are the model fluxes from the best-fit Kurucz atmosphere model (black).

Figure 2 .

 2 Figure 2. PDCSAP lightcurves of TOI-282, as produced by the SPOC pipeline. The top panel shows the photometry for sectors 1-13, while the bottom panel contained sectors 27-37, with a gap where sector 32 was missing. Binned data with a cadence of 30 min are shown in purple.

Figure 3 .

 3 Figure 3. Results of our search for additional candidates in the data, using T L S as implemented by S . Upper panel: Lomb-Scargle periodogram showing the detected period and its harmonics. Lower panel: TESS lightcurve phase-folded on this period, with a transit model overplotted.

Figure 4 .

 4 Figure 4. Contrast curves for HD 28109/TOI-282. The top panel shows the results from Gemini South, the middle panel shows the results from SOAR, and the bottom panel shows results from VLT/NaCo. All three panels are inset with final reconstructed images; none of the observations revealed bound companions within 5 magnitudes of the target.

•

  TOI-282.04 → HD 28109 b • TOI-282.01 → HD 28109 c • TOI-282.03 → HD 28109 d

Figure 5 .

 5 Figure 5. Images of the field of TOI-282/HD 28109, showing that the star does not have any bright contaminating sources nearby. Left Panel: Target pixel files showing TOI-282 as a red star and other nearby sources as orange circles. The size of the circles is representative of the source's brightness. Middle Panel: Sloan Digital Sky Survey (SDSS York et al. 2000) image of the surrounding field of TOI-282 (indicated by a red cross). Right Panel: TESS target pixel file of one of the 25 sectors in which TOI-282 was observed, created with TPFP (Aller et al. 2020). The pipeline aperture is shown by the white shaded region, and all nearby sources from Gaia (Gaia Collaboration et al. 2018) down to a magnitude difference ΔM = 8 are overplotted as red circles.

Figure 6 .

 6 Figure 6. Detrended transits observed with ASTEP on the nights of 2020 June 08, 2021 March 15 and 2021 May 10. The green points are the detrended flux, while the white circles are binned. Curves superimposed are the best-fitting models for each transit.

Figure 7 .

 7 Figure 7. Global fit of the TESS photometric data for each of the planets of HD 28109. Points in the background are the phase-folded 2-minute cadence data, while white circles are the binned points. Curves are 20 samples drawn from the posterior.

Figure 8 .

 8 Figure 8. Transit timing variations (in minutes) for each of the planets of HD 28109 spanning two years of TESS data plus ground-based photometric monitoring of planet c. The points show the difference between the fitted midtimes and the linear predictions, while the lines show the model predictions from our dynamical analysis in Section 5.2.2.

Figure 9 .

 9 Figure 9. The 3:2 resonance structure diagram following Nesvorný & Vokrouhlický (2016). Systems with 𝛼 = 𝑎 c /𝑎 d < 𝛼 res (𝛼 > 𝛼 res ), where 𝛼 = 0.763 corresponds to the exact resonance, have orbits just wide (narrow) of the resonance. The resonance region where 𝜓 librates is shaded. The separatrices and stable point are solid. Three planetary systems are plotted: HD 28109 (this work), TOI-216 (Dawson et al. 2021; Nesvorny et al. 2021) and KOI-142(Nesvorný et al. 2013). TOI-216 is firmly in the resonance and KOI-142 is firmly outside the resonance. For HD 28109, the blue dots show a sample of 500 posterior solutions from our dynamical 3-planet fit, of which roughly two-thirds are non-resonant and roughly one-third are resonant.

Figure 10 .

 10 Figure 10. Model predictions of the transit timings variations for the planets of HD 28109 for the next 11 000 days, spanning more than one full TTV period for the two outer planets.

Figure 11 .Figure 12 .

 1112 Figure 11. Phased radial velocity measurements for each of the three planets. Triangles are points from ESPRESSO, while circles are points from HARPS; the curves are 20 model samples from the posterior. The data uncertainties include the RV jitter.

Figure 13 .

 13 Figure 13. Mass-radius diagram showing the positions of HD 28109 b, c and d. Planets c and d are plotted as squares with their masses are estimated from TTVs.As the mass of planet b is still poorly constrained, we plot its position using its RV estimated mass as a hexagon, while the average mass for planets of its radius (estimated mass following Chen & Kipping (2017)) is plotted as a cross.

Figure 14 .

 14 Figure 14. Model transmission spectra of HD 28109 b, c and d with simulated HST observations. Each point corresponds to a single transit depth with no offset; the curves represent the simulated transmission spectra.

Table 1 .

 1 Stellar parameters adopted for this work.

	Designations	HD 28109, TIC 29781292, TOI-282, 2MASS J04205712-
		6806095, UCAC4 110-003794, WISE J042057.17-
		680609.5, Gaia DR2 4668163021600295552, HIP 20295,
		TYC 9154-01248-1
	Parameter	Value	Source
	T mag	8.9387±0.006	Stassun et al. (2019)
	B mag	9.91±0.03	Høg et al. (2000)
	V mag	9.38±0.02	Høg et al. (2000)
	G mag	9.3063±0.0028	Gaia Collaboration et al.
			(2021)
	J mag	8.476±0.020	Cutri et al. (2003)
	H mag	8.256±0.024	Cutri et al. (2003)
	K mag	8.175±0.023	Cutri et al. (2003)
	W1 mag	8.136±0.023	Wright et al. (2010)
	W2 mag	8.175±0.020	Wright et al. (2010)
	W3 mag	8.155±0.020	Wright et al. (2010)
	W4 mag	7.997±0.161	Wright et al. (2010)
	Distance	140.087±0.194	Bailer-Jones et al. (2021)
	SpT	F8/G0V	Houk & Cowley (1975)
	𝑅 ★	1.446 ± 0.035𝑅	This work
	𝑀 ★	1.26 ± 0.08𝑀	This work
	Age	1.1±0.1 Gyr	This work
	T eff	6120±50 K	This work
	log 𝑔 ★	4.13±0.10	This work
	[Fe/H]	0.0±0.1	This work
	𝑣 sin 𝑖 ★	7.7±0.5 km/s	This work

Table 2 .

 2 Summary of ground-based follow-up observations carried out for the validation of .03 and .04.

			Follow-up Observations
			High Resolution Imaging
	Observatory	Filter	Date	Sensitivity Limit	Result
	VLT	𝐾 𝑠	2019 January 25	Δ𝑚 = 5 at 0.23 -4.5"	No sources detected
	SOAR	𝐼 𝑐	2019 February 02	Δ𝑚 = 7.3 at 1 -3"	No sources detected
	Gemini South	562 nm	2020 December 23	Δ𝑚 = 5 at 0.1 -1.2"	No sources detected
	Gemini South	832 nm	2020 December 23	Δ𝑚 = 8 at 0.1 -1.2"	No sources detected
			Photometric Follow-up
	Observatory	Filter	Date	Coverage	Result
			TOI-282.01 / HD 28109 c
	Hazelwood	𝐼 𝑐	2019 January 19	Ingress +50%	Field cleared of NEBs within 2.5"
	ASTEP	𝑅 𝑐	2020 June 08	Full	Detection
	ASTEP	𝑅 𝑐	2021 March 15	Ingress +50%	Detection
	ASTEP	𝑅 𝑐	2021 May 10	Ingress +90%	Detection
	ASTEP	𝑅 𝑐	2021 July 05	Full	Non-detection
	ASTEP	𝑅 𝑐	2021 August 30	Full	Non-detection
	LCO	𝑧 𝑠	2021 August 30	Full	Field cleared of NEBs within 2.5"
			TOI-282.03 / HD 28109 d
	ASTEP	𝑅 𝑐	2021 March 08	Ingress +30%	Field cleared of NEBs within 2.5"
	ASTEP	𝑅 𝑐	2021 August 23	Ingress +40%	Field cleared of NEBs within 2.5"
			TOI-282.04 / HD 28109 b
	Hazelwood	𝐼 𝑐	2019 January 19	Ingress +50%	Field cleared of NEBs within 2.5"
			Spectroscopic Observations
	Instrument	Wavelength Range	Date Range	Number of Spectra	Use
	CHIRON	450 -890 nm	2019 February 14 -2021 January 07	3	Stellar characterisation
	ESPRESSO	380 -788 nm	UT 2019 October 15 -UT 2019 December 27	8	Rule out stellar companion
	HARPS	383 -693 nm	UT 2019 May 24 -UT 2021 January 29	7	Rule out stellar companion

  TOI-282.01 was observed on four further occasions by ASTEP (UT 2021 March 15, UT 2021 May 10, UT 2021 July 05, UT 2021 August 30) resulting in two partial transits and two non-detections due to poor weather. TOI-282.03 was observed on the UT 2021 March 08 and again on the UT 2021 August 23. Both observations resulted in non-detections of the shallow event, but the field was cleared of NEBs within 2.5".

Table 3 .

 3 Priors used in our fit, along with fitted and derived parameters. Uniform priors are indicated as U (lower bound, upper bound) and normal priors are indicated as N (mean, standard deviation).

		HD 28109 b	HD 28109 c	HD 28109 d	Source
		Fit Parametrisation and Priors
	Transit Depth R p /R ★	U (0.008, 0.015)	U (0.023, 0.035)	U (0.017, 0.037)
	Inverse Scaled Semi-major Axis (R ★ + R p )/a	U (0.05, 0.06)	U (0.015, 0.035)	U (0.015, 0.03)
	Orbital Inclination cos i	U (0.000, 0.04)	U (0.000, 0.04)	U (0.000, 0.04)
	Transit Epoch T 0	U (2458344.7, 2458344.9)	U (2458337.12, 2458337.32)	U (2458355.58, 2458355.98)
	Planet Mass M p /M ★	U (0.0, 0.000026)	U (0.0, 0.0005)	U (0.0, 0.0005)
	Period P	U (22.8, 23.0)	U (55.6, 56.6)	U (84.2, 86.0)
	Eccentricity e	U (0.0, 0.5)	U (0.0, 0.5)	U (0.0, 0.5)
	Longitude of Periapsis 𝜛	U (0, 360)	U (0, 360)	U (0, 360)
		Limb Darkening Coefficients
	TESS u 1	0.4041 ± 0.0007	ASTEP u 1	0.4758 ± 0.0009
	TESS u 2	0.1324 ± 0.0018	ASTEP u 2	0.1392 ± 0.0022
	TESS q 1	N (0.2878, 0.0024)	ASTEP q 1	N (0.3782, 0.0034)
	TESS q 2	N (0.3766, 0.0013)	ASTEP q 2	N (0.3868, 0.0014)
	External Priors		GP Priors
	Stellar Mass M ★	N (1.26, 0.08)	Amplitude Scale gp ln S 0 (flux)	U (-15, 5)
	Stellar Radius R ★	N (1.446, 0.035)	Damping gp ln Q(flux)	U (-10, 10)
	Stellar Effective Temperature T eff	N (6120, 50)	Frequency gp ln 𝜔 0 (flux)	U (-10, 10)
		Fitted Parameters
	R p /R ★ (R ★ + R p )/a cos i T 0 (BJD) P (d)	0.01441 +0.00031 -0.00037 0.05017 +0.00022 -0.00012 0.03964 +0.00024 -0.00040 2458344.81772 ± 0.00757 22.89104 +0.00035 -0.00036	0.02632 ± 0.00027 0.02231 +0.00073 -0.00066 0.0081 +0.0018 -0.0021 2458377.80109 +0.00724 -0.00733 56.00819 +0.00194 -0.00202	0.01999 +0.00035 -0.00037 0.01688 +0.00053 -0.00047 0.0064 +0.0013 -0.0015 2458355.67324 ± 0.00432 84.25999 +0.00744 -0.00662	Fixed Lin. Ephem. Fit Fixed Lin. Ephem. Fit Fixed Lin. Ephem. Fit Dynamical Fit Dynamical Fit
	e	< 0.3307	< 0.1203	< 0.0864	Dynamical Fit
	Fitted RV Jitter		Fitted GP Hyperparameters
	ln 𝜎 jitter;HARPS (ln km/s)	-7.8 +1.6 -2.8	gp ln S 0 (flux)	-19.72 +0.20 -0.19

MNRAS 000,1-19 (2022)

The 31-day transit event was caused by a combination of noise and confusion with the 84-day transit events in early TESS data

. 2 Sis publicly available at https://github.com/franpoz/ SHERLOCK.MNRAS 000,1-19 (2022)

All data validation reports for TOI-282 are publicly available at https: //exofop.ipac.caltech.edu/tess/target.php?id=29781292.

https://www.gemini.edu/sciops/instruments/alopeke-zorro/ MNRAS 000,1-19 (2022)

https://www.eso.org/sci/software/pipelines/espresso/ espresso-pipe-recipes.html MNRAS 000,1-19 (2022)

http://www.astropy.org MNRAS 000,1-19 (2022)

This paper has been typeset from a T E X/L A T E X file prepared by the author.MNRAS 000,1-19 (2022)

ACKNOWLEDGEMENTS

The authors wish the thank the reviewer for their comments and feedback which helped to improve and clarify the manuscript. This work makes use of observations from the ASTEP telescope. ASTEP benefited from the support of the French and Italian polar agencies IPEV and PNRA in the framework of the Concordia station program, from OCA, INSU, and from Idex UCAJEDI (ANR-15-IDEX-01). The authors would also like the winterover staff of Concordia station; their work makes the continued operations of ASTEP possible. Some of the observations in the paper made use of the High-Resolution Imaging instrument Zorro obtained under Gemini LLP Proposal Number: GN/S-2021A-LP-105. Zorro was funded by


Jean-Luc Bertrand-Krajewski
email: jean-luc.bertrand-krajewski@insa-lyon.fr

Lôide Angelini-Sobrinha

Vincent Pons

Jean-Christophe Grimard

How Green Roof Performance Depends on Climatic Conditions and Climate Change? First Results of Simulation Scenarios

Keywords: Hydrological simulation, Sizing, Retention, Vegetated roof, Source control, Stormwater management

Hydrological simulations were used to analyse the performance of an advanced green roof with additional underlying storage. Simulations were done with 20-year data series of precipitation and evapotranspiration or temperature for the cities of Lyon, Melbourne, and Seoul, considering a current climate scenario (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019). Additionally, two future scenarios for the city of Lyon (2035-2044 and 2090-2099) were simulated. Simulations were done for substrate depths ranging from 60 to 240 mm. The advanced green roof has an additional underlying storage layer, which was simulated with storage capacities from 52 to 150 mm. This abstract shows some of the first results regarding retention performance. There is a significant variability between the three cities (climate context) and, in Lyon, between past (2000(climate context) and, in Lyon, between past ( -2019) ) and future (2090-2099) climates.

Introduction

Event and long-term hydrological performance indicators of green roofs (retention, actual evapotranspiration, peak flow attenuation, outflow, etc.) depend on the meteorological and climatic conditions prevailing at their locations. For a given green roof with a particular design (e.g. substrate depth, additional storage volume, controlled outflow, etc.), its performance will vary according to its geographical location and to the meteorological time series. In addition, for a given location, the performance will also vary in the future due to climate change. The Olympluie project aims to evaluate and compare the performance of various green roofs for different cities around the world, and, for some of these locations, estimate how the performance will change by using future time series derived from the IPCC RP 8.5 scenario. In this abstract, some of the first results obtained for Lyon, Melbourne and Seoul are presented.

Material and methods

Green roofs

The advanced green roofs simulated in this work include i) a substrate depth ranging 60 to 240 mm, and ii) an additional underlying storage layer below the substrate ranging from 52 to 150 mm (or L/m 2 ) [START_REF] Bournique | Modélisation hydrologique d'une solution de toiture bleue -verte de gestion des eaux pluviales[END_REF]. The green roofs are vegetated with sedums and equipped with outflow controllers. The advanced green roof (commercial name: OASIS) is manufactured by the French company Le Prieuré and will be installed on the roofs of the Athletes' Village building managed by the SOLIDEO company for the Paris 2024 Olympic Games (hence the project name "Olympluie", i.e. Olympic rain).

Historical time series

Rainfall intensity (I, mm/h) and potential evapotranspiration (PET, mm/h) time series have been created from raw data sets collected for various cities around the world with time steps ranging from 1 to 60 minutes. Intensity raw data on the one hand, and PET or air temperature raw data on the other hand, were available with time steps from 1 minute to 1 hour (resp. 1 day). In case PET data were not available, PET was estimated from the measured or calculated mean daily temperature by means of the [START_REF] Oudin | Which potential evapotranspiration input for a lumped rainfall-runoff model? Part 2 -Towards a simple and efficient potential evapotranspiration model for rainfall-runoff modelling[END_REF] formula [START_REF] Oudin | Which potential evapotranspiration input for a lumped rainfall-runoff model? Part 2 -Towards a simple and efficient potential evapotranspiration model for rainfall-runoff modelling[END_REF]. The resulting time series used for modelling scenarios are summarized in Table1. Other cities are presently considered to extend the list of cases. Future time series For some cities, IPCC RCP 8.5 scenario regional modelling results are available. In this case, these data, combined with historical time series, have been used to create future time series by downscaling IPCC data by means of a method based on the MRC (multiplicative random cascade) approach [START_REF] Pons | Forecasting green roof detention performance by temporal downscaling of precipitation time-series projections[END_REF]. In this abstract, only future time series for Lyon are shown.

Green roof model

The model used for the simulations is the Gepeto model developed by INSA Lyon in partnership with Le Prieuré. The model was previously tested and calibrated with experimental data [START_REF] Bournique | Modélisation hydrologique d'une solution de toiture bleue -verte de gestion des eaux pluviales[END_REF][START_REF] Principato | Test of a vegetated roof hydrological model with two experimental data sets[END_REF].

In the following simulations, the green roof settings are default values, the most important ones being the substrate depth (80 mm), the substrate storage capacity (52 mm), the cells capacity just below the substrate (4 mm), the additional storage layer (52 mm), the controlled outflow (10 L/s/ha), the controlled outflow starting depth in the storage layer (9.83 mm), and two outflow controllers.

Simulation scenarios

Simulations are carried out with historical and future time series given in Table1. In the first phase of work presented in this abstract, two design elements were varying: i) the thickness of the substrate, from 60 to 240 mm, and ii) the rainwater underlying storage capacity, from 52 to 150 mm. All simulations are run with a 6-min time step.

Results

Some examples of simulation results are given hereafter for Lyon, Melbourne, and Seoul. The selected performance indicator is the interception efficiency, i.e. the percentage of rainwater volume that is retained and/or evapotranspirated by the green roof. An interception efficiency of e.g. 60% means that 60% of the rain falling over the green roof is stored and/or evapotranspirated. Consequently, 40% of the rainfall volume is discharged by the green roof either as controlled outflow or as overflow in case the green roof maximum storage capacity (in both substrate and storage layers) is reached. Other results and performance indicators (peak flows, overflows, etc.) for various green roof configurations will be presented during the conference.

In the case of Lyon, the global interception efficiency for the period 2000-2019 for a substrate depth of 80 mm and a storage layer depth of 52 mm is 66.7%. The interception efficiency is only 60.0% for the future climate period 2090-2099, showing the effect of the climate change on the performance. The box plot of the annual interception efficiencies for both periods is shown in Figure 9. For the period 2000-2019, the annual interception efficiency varies from 58.3% to 87.6%. For the period 2090-2099, it varies from 51.8% to 70.8%. Figure 10 shows the comparison of annual interception efficiencies for the same roof as in the case of Lyon simulated in Melbourne and Seoul. In Melbourne, the annual interception efficiency ranges from 65.0% to 96.2%, which is better that in the case of Lyon due to differences in climatic conditions (not detailed here). In Seoul, the annual interception efficiency ranges from 34.5% to 71.3%: this lower performance, compared to Lyon and Melbourne, is due to the very strong rainfall events which are observed during the monsoon period from mid-June to mid-September. In this case, the green roof design must be adapted, with higher substrate depth, or with a higher storage layer depth, or a combination of both (results not shown here) to reach a given level of performance.

Conclusion

The Olympluie project simulates various green roof design scenarios for different climatic contexts, including future climate derived from IPCC data sets. Many performance indicators will be evaluated for a detailed comparison of scenarios, allowing to adapt the green roof design to each geographical and climatic context.

Figure 9 .

 9 Figure 9. Box plot of annual interception efficiency in Lyon for the periods 2000-2019 (left, n = 20) and 2090-2099 (right, n = 10). Hsubstrate = 80 mm, Hstorage = 52 mm.

Figure 10 .

 10 Figure 10. Box plot of annual interception efficiency in Melbourne (left) and Seoul (right) for the period 2000-2019 (n = 20). Hsubstrate = 80 mm, Hstorage = 52 mm. Simulations for Melbourne and Seoul over the period 2000-2019 gives global interception efficiencies equal to 79.9% and 47.1% respectively. Figure10shows the comparison of annual interception efficiencies for the same roof as in the case of Lyon simulated in Melbourne and Seoul. In Melbourne, the annual interception efficiency ranges from 65.0% to 96.2%, which is better that in the case of Lyon due to differences in climatic conditions (not detailed here). In Seoul, the annual interception efficiency ranges from 34.5% to 71.3%: this lower performance, compared to Lyon and Melbourne, is due to the very strong rainfall events which are observed during the monsoon period from mid-June to mid-September. In this case, the green roof design must be adapted, with higher substrate depth, or with a higher storage layer depth, or a combination of both (results not shown here) to reach a given level of performance.

Table 1 .

 1 Cities and time series used for simulation scenarios. Due to limited space, only cities marked with * are illustrated in this abstract.

	City, Country	Historical I and PET	Future I and PET	Simulation time step of
		time series	time series	I and PET time series
	Lyon*, France	2000-2019	2035-2044, 2090-2099	6 min
	Melbourne*, Australia	2000-2019		6 min
	Seoul*, South Korea	2000-2019		1, 6 and 15 min
	Montreal, Canada	2005-2021		1, 6 and 15 min
	Los Angeles, USA	2000-2019		1 hour
	Dourados, Brazil	2000-2019		15 min and 1 hour
	Cape Town, South Africa	2000-2019		1 hour

Acknowledgment

This work was realised within the Graduate School H2O'Lyon (ANR-17-EURE-0018) and Université de Lyon (UdL), as part of the pro-gramme 'Investissements d'Avenir' run by Agence Nationale de la Recherche (ANR). Thanks to all colleagues in Melbourne, Seoul, Montreal, Los Angeles, Dourados and Cape Town who shared and/or helped to get data in their respective cities.

Biography

The presenting author is Professor at the University of Lyon, France, in INSA Lyon (National Institute of Applied Sciences). He was head of the laboratory DEEP (Wastes, Water, Environment, Pollutions) in the period 2014-2022 and is co-president of the multidisciplinary OTHU -Field Observatory for Urban Hydrology (www.othu.org) since 2023. His main research topics are related to urban hydrology and urban water systems, especially hydraulic and pollutants transport modelling in sewer systems, field monitoring and metrology, stormwater management, and integrated approach for urban water systems. More information and publications are available on his website at jlbkpro.free.fr