Constraints on HD113337 fundamental parameters and planetary system. Combining long-base visible interferometry, disk imaging and high-contrast imaging

HD113337 is a Main-Sequence F6V field star more massive than the Sun, hosting one (possibly two) radial velocity (RV) giant planet(s) and a cold debris disk (marked by an infrared excess). We used the VEGA interferometer on the CHARA array to measure HD113337 angular diameter, and derived its linear radius using the Gaia parallax. We computed the bolometric flux to derive its effective temperature and luminosity, and we estimated its mass and age using evolutionary tracks. We used Herschel images to partially resolve the outer disk, and high-contrast images of HD113337 with the LBTI to probe the 10-80 au separation range. Finally, we combined the deduced contrast maps with previous RV of the star using the MESS2 software to bring upper mass limits on possible companions at all separations up to 80 au, taking advantage of the constraints on the age and inclination (brought by the fundamental parameter analysis and the disk imaging, respectively). We derive a limb-darkened angular diameter of 0.386 $\pm$ 0.009 mas that converts into a linear radius of 1.50 $\pm$ 0.04 solar radius. The fundamental parameter analysis leads to an effective temperature of 6774 $\pm$ 125 K, and to two possible age solutions: one young within 14-21 Myr and one old within 0.8-1.7 Gyr. We partially resolve the known outer debris disk and model its emission. Our best solution corresponds to a radius of 85 $\pm$ 20 au, an extension of 30 $\pm$ 20 au and an inclination within 10-30 degrees for the outer disk. The combination of imaging contrast limits, published RV, and our new age and inclination solutions leads to a first possible estimation of the true masses of the planetary companions: $\sim 7_{-2}^{+4}$ Jupiter masses for HD113337 b (confirmed companion), and $\sim 16_{-3}^{+10}$ Jupiter masses for HD113337 c (candidate). We also constrain possible additional companions at larger separations.


Introduction
Thousands of exoplanets have been discovered for more than twenty years, exhibiting a wide diversity of properties (mass, separation, eccentricity, etc). Each planet detection method allows to estimate only some of these parameters. Most of the planetary companions known so far have been detected indirectly, either with the transit or the radial velocity (RV) method. For both techniques, the main derived parameter (i.e., the planet radius from transits, and the planet minimal mass from RV) de-Send offprint requests to: simon.borgniet@obspm.fr Partly based on observations made with the VEGA/CHARA spectro-interferometer. pends directly on the values of the star parameters (i.e., the stellar radius and mass, respectively). Hence a better precision on these stellar parameters leads to better estimations of the planetary parameters , Stassun et al., 2017, White et al., 2018. In the case of transiting companions, the combined use of the stellar mass and radius even allows to determine the true companion mass and thus, its density and expected composition (Ligi et al., 2016, Crida et al., 2018. In the case of RV planetary systems, another major uncertainty lies in the inclination (i) between the orbital plane and the observer line-of-sight. Indeed, RV alone only provide the companion minimal mass (m p sin i) and not its true mass, which is decisive to infer the companion true nature. While astrometric S. Borgniet et al.: Constraints on the HD113337 system. measurements of the primary star used in combination with RV can help to derive the true companion mass for the most massive ones (brown dwarfs; see e.g. Sahlmann et al., 2011, Bouchy et al., 2016, one can also try to estimate the system inclination by looking at other proxies. Deriving the inclination of the star's rotation axis or the inclination of a resolved debris disk are two such possibilities, if assuming that they all rotate in the same plane (i.e. the orbital plane).
Estimating precisely the main stellar fundamental parameters, i.e. the stellar mass (M ), the stellar radius (R ), and the stellar age is far from being straightforward. Most estimates of these parameters are based on the use of evolutionary models with constraints brought by various observations. A direct and accurate way to obtain stellar radii is to use long-baseline interferometry to directly measure the stellar angular diameter, which allows to reach an unbiased precision of ∼3% on R (see e.g. Ligi et al., 2012, Boyajian et al., 2015, Ligi et al., 2016. When combined to the stellar bolometric flux and the parallax, such a measurement of the stellar radius allows to derive new (and potentially unbiased) estimations of the stellar luminosity (L ) and effective temperature (T eff ). Then, once placed on an Hertzprung-Russell diagram, the mass and age can be determined through the use of stellar evolutionary models and the interpolation of isochrones. A good knowledge of the stellar age is essential in the case of directly imaged substellar companions, as the companion mass is determined through the use of mass-luminosity model. Hence, understanding the true nature and the formation processes of such imaged companions strongly depend on a good estimation of the age of the primary star. Famous cases are e.g. the giant planet orbiting around β Pictoris (Bonnefoy et al., 2014b) or the companion to GJ504 (D'Orazi et al., 2017. A key challenge to develop the theory of planetary formation and evolution processes is to understand the respective influence of the different stellar characteristics (e.g. the stellar mass, metallicity, effective temperature, etc) on these processes. While stellar metallicity is well known to be positively correlated to the giant planet (GP) frequency (Fischer & Valenti, 2005), a correlation between the stellar mass and the GP frequency and/or mass is yet to be fully investigated for stars more massive than the Sun (see e.g. Borgniet et al., 2017, and references therein). Interactions between giant planets and debris disks are another key topic to investigate in the context of planetary evolution.
We make here a case study of a system of high interest. HD113337 is a Main-Sequence star more massive than the Sun, that hosts one (possibly two) RV-detected giant planet(s) (Borgniet et al., 2019, hereafter BO19+), as well as an unresolved debris disk. We use multi-technique observations to better understand and/or constrain the properties of the primary star (through optical interferometry), to resolve the debris disk, and to explore the system's outer environment (through deep imaging). This paper is structured as follows: in Sect. 2, we present the HD113337 system, looking at the star, the debris disk and the planetary system. Second, we review in Sect. 3 the different observations that we made and the data reduction processes that we used. We present and discuss our results in Sect. 4. We specifically show how the combination of these different techniques allow us to better understand and constrain the HD113337 system.

The star
HD113337 is a Main-Sequence F6V-type star located at a distance d = 36.2±0.2 pc from the Sun (based on the parallax given by the Gaia second Data Release or DR2, Gaia Collaboration et al., 2016, and see details in Sect. 4.1). Its stellar mass is consistently estimated to be ∼1.4 M in the literature, based on spectroscopic (Allende Prieto & Lambert, 1999, hereafter AP99+) or photometric analyses (Geneva-Copenhagen Survey -hereafter GCS III+, Casagrande et al., 2011). Its estimated effective temperature ranges from 6670 ± 80 K based on the photometry (GCS III+) to 6760 ± 160 K based on the spectroscopy (AP99+). The fit of the spectral energy distribution (SED) by Rhee et al. (2007) gives T eff = 7200 K and provides a stellar radius estimation of 1.5 ± 0.15 R . As it is a field dwarf star, the age of HD113337 is the most difficult stellar parameter to estimate. The typical isochronal age derived from the photometric T eff is 1.5 +0.43 −0.55 Gyr (GCS III+). It is in agreement with the age of 1.6 +2.2 −0.8 Gyr derived by David & Hillenbrand (2015) from Strömgren photometry. We conducted two different analyses to derive independently an age estimation of HD113337 in Borgniet et al. (2014). Briefly, we first estimated the age of the bound distant (projected separation of ∼120 as or ∼4400 au) Mtype companion 2M1301+6337 to be 100 +100 −50 Myr. Second, we measured HD113337 Lithium abundance and estimated the corresponding age to be > 160 Myr, leaving our analysis unconclusive (for more details, see Borgniet et al., 2014, and references therein). Finally, activity-and rotation-related age diagnostics such as the relations derived by Mamajek & Hillenbrand (2008) do not apply to such an early spectral type.

The debris disk
HD113337 exhibits a clear infrared (IR) excess from ∼20 µm up to 1200 µm with a L IR /L BOL = 10 −4 fractional luminosity. Based on data from the InfraRed Astronomical Satellite (IRAS), Rhee et al. (2007) estimated the dust temperature to be ∼100 K and radius to be 18 au. Using Spitzer data that provide a better coverage on the spectral energy distribution (SED) of the disk, Moór et al. (2011) concluded that the disk has a dust temperature of ∼53 K, suggesting a disk radius of 55 ± 3 au. A more recent study by Chen et al. (2014) found that HD113337 SED was best fitted by a two-belt model, with a first, warm (316 ± 10 K) dust ring located at 1.7 au from the star and a second, cold (54 ± 5 K) dust ring at 179 au. The system was observed with the James Clerk Maxwell Telescope (JCMT) for the SCUBA-2 Observations of Nearby Stars (SONS) survey at both 450 µm and 850 µm by Holland et al. (2017), who reported upper limits of <75 mJy (5σ) and <3.6 mJy (3σ), respectively.

The planetary system
From 2006 to 2016, a RV survey of 125 northern AF-type dwarf stars (including HD113337) was carried out with the SOPHIE spectrograph at Observatoire de Haute-Provence (France). The aim was to search for giant planets and brown dwarfs (BD) around Main-Sequence stars more massive than the Sun. Clear periodic variations of HD113337 RV were detected and attributed to the presence of a ∼3 M Jup GP orbiting around the star with a ∼320-day period (∼1 au, Borgniet et al., 2014). After monitoring the system with SOPHIE for three additional years, HD113337 RV were found to exhibit a second Columns 1 and 5 give the observation date and UT time (on the science target). Columns 2, 3 and 4 give the calibrator identifier in the HD catalog, its uniform-disk angular diameter (θ CAL UD ) in the R-band, and the 1σ error bar on the calibrator diameter (σ CAL θ ). We took the calibrator diameters from the JMMC Stellar Diameters Catalog Version 2 (JSDC, Bourges et al., 2017), while we kept the corresponding 1σ uncertainties from the JSDC Version 1 (Lafrasse et al., 2010, see text). Column 6 gives the calibration sequence followed for the observation and column 8 is the baseline used. Column 7 gives the central wavelength of the 20 nm-wide spectral bands in which we computed the V 2 . Columns 9 and 10 give the projected base length B P and its orientation PA. Column 11 gives the corresponding calibrated V 2 value. Column 12 gives the Fried parameter (estimation of the quality of the atmosphere) for each observation night. Note that three observation points are not displayed here as the results were fully rejected (see text). periodicity on a longer timescale. The possible sources for this RV long-term variability were investigated and it was concluded that a possible origin was the presence of a second GP with a m p sin i ∼7 M Jup minimal mass on a wider orbit (∼5 au, see BO19+).
The combined presence of a debris disk, giant planet(s), as well as an ill-constrained age, makes HD113337 an object of high interest.

Optical interferometry of HD113337
The Center for High Angular Resolution Astronomy array (hereafter CHARA, ten Brummelaar et al., 2005) is the main optical and near-infrared interferometric array in the northern hemisphere. It hosts six 1-m telescopes arranged by pairs in a Y shape and oriented to the west (W1 and W2), east (E1 and E2) and south (S1 and S2), allowing a wide range of baseline orientations. The corresponding baselines range from ∼30 m to ∼330 m (i.e. a maximal angular resolution of 0.2-0.3 mas in the visible). The Visible SpEctroGraph and polArimeter (hereafter VEGA, Mourard et al., 2009) is one of the instruments operating in the visible at the focus of the CHARA array. VEGA is a spectro-interferometer which allows to combine the light coming from 2 to 4 telescopes simul-taneously, at different spectral resolutions (R = 6000 and 30000).
We observed HD113337 with two different telescope triplets (E1E2W2 and E1E2W1) chosen to (partially) resolve its small expected angular diameter (∼0.4 mas given its ∼27 mas parallax and the 1.5 R radius from Rhee et al., 2007). For each observation point, we tried to follow a calibrator-target-calibrator sequence (C-S-C, see Table 1) with either 30 or 40 blocks of 2500 short (10 ms) exposures per star to ensure an instrumental transfer function stable enough to correctly calibrate the target squared visibilities (V 2 ). While it is not necessarily mandatory, observing the calibrator (C) star twice (e.g. before and after the science (S) target) allows to monitor and take into account possible variations of the transfer function during the observation time (see Mourard et al., 2012, and below). Furthermore and if possible, using two different calibrators with well-defined angular diameters (C 1 -S-C 2 sequence) instead of one (C-S-C sequence) reinforces the robustness of the target V 2 computation by reducing its dependency to the calibrator diameter value and uncertainty. We used the SearchCal software (Bonneau et al., 2006) to select adequate calibrators at the different observation epochs. We acquired ten observation points on five different nights from May 2013 to June 2015, using VEGA red spectral channel in the ∼700 to ∼750 nm range. We computed the V 2 values in either two or three 20 nm-wide spectral bands (depending on the observing conditions) with the standard VEGA reduction pipeline vegadrs . Due to the small angular diameters of both our target and calibrators, we had to discard a signif- icant part of the data that revealed themselves not robust enough. We first discarded three entire observation points (UT 06:37 on 2013-05-25, UT 04:59 on 2014-07-07 and UT 07:32 on 2015-06-01), for which either the calibrated target V 2 computation process did not converge or the signal-to-noise ratio (S/N) on the V 2 values was very low ( 0). It also happens that we were not able to obtain measurements on the second calibrator for these three observation points (C-S sequence), which may have hindered the estimation of the transfer function at the time of the observations. Second, we also discarded V 2 measurements with S/N <4 , while ensuring that this did not bias our results (as done by Perraut et al., 2013). The detail of our observations (after selection) is given in Table 1.
The JMMC Stellar Diameter Catalog (JSDC) that provides the calibrator angular diameters has recently been updated (Bourges et al., 2017), with a significant increase (by 6-12% here) on the diameters of early-type calibrators (θ CAL UD ) and diameter uncertainties (σ CAL θ ) smaller by ∼50%. Here we chose to use the θ CAL UD values from the more recent JSDC2 while conservatively keeping the σ CAL θ values from the JSDC1 ( Table 1). The θ CAL UD change from JSDC1 to JSDC2 translates into an average increase of ∼7% on our V 2 values lower than 0.6 (2013-05-25 data) and no significant change on our V 2 values above 0.6. We fitted the visibility measurements with the JMMC fitting engine LITpro (Tallon-Bosc et al., 2008) based on a modified Levenberg-Marquardt algorithm 1 , and derived a uniformdisk (UD) angular diameter θ UD = 0.371±0.009 mas. We display the best model of the visibility function derived with LITpro in Fig. 1. If using the θ CAL UD values from the JSDC1 instead, we obtain θ UD = 0.364 ± 0.009 mas (i.e. less than a 2% difference). Furthermore, using the smaller σ CAL θ values from JSDC2 would lead to a θ UD uncertainty of 0.007 mas instead of 0.009 mas. Our θ UD robustness is due to the strong constraints brought on our model by the five measurements made with the E2W1 intermediate baseline on 2013-05-25 (our best data, see Fig. 1  We used the linear limb-darkening coefficients in the R-band provided by Claret & Bloemen (2011) to derive the corresponding limb-darkened angular diameter θ LD . For a solar metallicity and a null microturbulent velocity, we computed the limb-darkened diameters for T eff in the 6500 to 7000 K range, and log g in the 4 to 4.5 range. These coefficients vary at the third decimal level within the considered T eff and log g ranges. We thus considered the average limb-darkening coefficient on our parameter space and obtained a limb-darkened angular diameter of θ LD = 0.386 ± 0.009 mas (hence a ∼2.4% precision).

High-contrast imaging
HD113337 was observed on January 7, 2015 with the LMIRCam near-infrared camera (Hinz et al., 2008, Skrutskie et al., 2010 at the Large Binocular Telescope Interferometer (LBTI). The LBTI was operated in double-aperture mode: the secondary deformable mirrors were used to record two side-by-side adaptiveoptics (AO) images of HD113337 recorded by LMIRCam at Lband (3.68 -3.88 µm). The telescope+instrument do not have a derotator. Therefore, it automatically enables for passive angular differential imaging (ADI, Marois et al., 2006). We obtained 597×4.95 second and 605×4.95 second AO exposures of the target for the right and left eyes of the telescope, respectively. Less AO exposures were obtained on the right eye because of open loops. The field orientation changed by 35.1 • and 37.2 • during that sequence of exposures on the right and left eyes, respectively. The core of the star's point-spread function was saturated over a diameter of ∼128 mas during the observations to increase the dynamics of the recorded images. We therefore had to acquire non-satured exposures of the star before and after the sequence of saturated exposures for calibrating the astrometry and photometry using a neutral density (attenuation factor of 9 × 10 −3 ).
The data were reduced using the MPIA ADI pipeline (Bonnefoy et al., 2014a). The pipeline carried out the basic cosmetic steps on the raw frames (de-trending of the raw frames, bad pixel interpolation, sky-background subtraction, flat field calibration). The star position was registered in the resulting frames using the mpfit2Dpeak.pro IDL function 2 which allowed for a bi-dimensional Moffat function to be fitted onto the PSF wings while masking the saturated core. The parallactic angles were computed at the time of the observations. We applied the LOCI algorithm (Lafrenière et al., 2007) to evaluate and subtract the stellar halo in the left and right eye data, independently. No point source could be detected following that step (Fig. 2).

Herschel observations of the cold outer disk
The Herschel data of HD113337 were obtained using the Photodetector Array Camera and Spectrograph (PACS) instrument (Poglitsch et al., 2010), in the mini scanmap mode with simultaneous observations at 70 and 160 µm under the open time program OT 2 ksu 3 (Fig. 3). The data were obtained on March 3, 2012 with OBSID 1342243344 and 1342243345. Herschel PACS data reduction were performed following the procedure published by Balog et al. (2014) for calibration stars.
We fitted a 2-D Gaussian function to estimate the source's Full Width at Half Maximum (FWHM). At 70 µm, the measured FWHM is 7. 44×7. 15, ∼1.3 times larger than the FWHM of typical point sources (5. 76×5. 58). The measured FWHM at 160 µm is slightly larger than the typical value for point sources, but less significant. Because the source is marginally resolved at PACS wavelengths, we used aperture sizes of 12" and 22" to measure photometry with a sky annulus of 35"-45" at 70 µm and 160 µm, respectively. Including 7% of abso-

Determination of the stellar fundamental parameters
Linear radius -We used the Gaia (Gaia Collaboration et al., 2016) DR2 parallax π P = 27.61 ± 0.04 mas (Gaia Collaboration et al., 2018). According to the Gaia documentation, the published DR2 parallax uncertainties may be underestimated by a factor of ∼10% for bright stars such as HD113337 (G 5.9 mag). Furthermore, there are potential systematic errors on the DR2 parallaxes such as a global zero-point offset (Lindegren et al., 2018). We used the formula given by Lindegren et al. (2018) for bright stars to recompute the error on HD113337 DR2 parallax: i.e. we quadratically summed the published parallax uncertainty scaled by a factor 1.08 with an additional uncertainty of 0.021 mas. We thus obtained π P = 27.61 ± 0.05 mas. The corresponding distance is 36.2 ± 0.2 pc, i.e. in good agreement with the HD113337 distance of 36.18 ± 0.06 pc derived by Bailer-Jones et al. (2018) based on a Bayesian analysis of Gaia DR2 parallaxes.
We used our limb-darkened angular diameter θ LD in combination with the above Gaia parallax to derive the stellar radius and its error through a Monte-Carlo simulation R ± δR = θ LD + δθ LD 9.305 × (π P + δπ P ) . (1) We obtained R = 1.50 ± 0.04 R (precision better than 2.7%).
Bolometric flux -We collected spectroscopic and photometric data to compute the bolometric flux f bol . We determined HD113337 SED, and then computed the area under the curve of that distribution (Fig. 4). The flux distribution of HD113337 was determined by concatenating several flux distributions: for the ultraviolet region, the rebinned spectrum from the Sky Survey Telescope obtained by the IUE 'Newly Extracted Spectra'(INES) data archive 3 . Based on the quality flag listed in the IUE spectra (Garhart et al., 1997), we removed all bad pixels from the data as well as measurements with negative flux; for the visible and red regions, the STELIB spectrum (Le Borgne et al., 2003); in the near-infrared range, the J, H, and K 2MASS magnitudes. We used the formula in Cohen et al. (2003) to convert 2MASS magnitudes in fluxes.
At the shortest wavelengths, we performed a linear interpolation on logarithmic scale between 912 Å and 1842 Å, considering zero flux at 912 Å. At the longest wavelengths, we performed a linear interpolation on logarithmic scale using the 2MASS magnitudes and assuming zero flux at 1.6 × 10 6 Å. We estimated the uncertainty associated to the bolometric flux by considering the following conservative uncertainties, i.e., 3% uncertainty on the flux computed from the STELIB spectrum, 10% on the flux computed from the combined IUE spectra, and 15% on the flux derived from interpolations. Finally, we obtained a bolometric flux f bol = 1.05 ± 0.06 × 10 −7 erg/cm 2 /s.
Luminosity and effective temperature -From the parallax and the bolometric flux we derived the luminosity through a Monte Carlo method where C is the conversion from parsecs to cm (3.086×10 18 ), and π p the parallax in arcseconds. We found L = 4.29 ± 0.25 L , where the error bar is dominated by that of the bolometric flux.
Position in the Hertzsprung-Russell diagram -From the determined fundamental parameters, we set our target in the Hertzsprung-Russell diagram. We used the isochrone tool CMD 2.7 4 to derive the mass and the age of HD113337. We considered a metallicity [M/H] = 0.07 ± 0.02 based on different spectroscopic analyses (Boesgaard & Tripicco, 1986, Soubiran et al., 2016, thus Z spanning from 0.0169 to 0.0185, and Y spanning from 0.279 to 0.282. We obtained two solutions in agreement with our 1-σ error box: a young solution corresponding to an age of 15 +6 −1 Myr and a mass of 1.48 ± 0.08 M , and an old solution corresponding to an age of 1.25 ± 0.45 Gyr and a mass of 1.40 +0.03 −0.05 M (Fig. 5).
The fundamental parameters (R , T eff ) that we derived are generally in close agreement with previous determinations (Table 2). Regarding the stellar age and mass, finding such a degeneracy between a young and an old solution appears to be a typical result when carrying out this approach (see e.g. Ligi et al., 2016. The mass and age from our old solution agree with the mass and age ranges determined by Casagrande et al. (2011) and AP99+. However, our two age solutions for the adopted metallicity range do not fit the age of 150 +100 −50 Myr that we adopted in Borgniet et al. (2014). Waiving the age degeneracy of HD113337 remains problematic so far as we know. Even combining our interferometric radius with a technique such as asteroseismology (Creevey et al., 2007)   On both plots, evolutionary tracks for Z in the range 0.0169 to 0.0185 are displayed for three cases: best age solution (solid orange curves), lower age limit (dotted red curves), and upper age limit (dashed green curves).  the stellar mass would prove fruitless, as the two mass values corresponding to our two age solutions are already consistent together. There are hints that the rate of debris disks detected through the presence of an IR excess decreases with the stellar age (Montesinos et al., 2016), yet this does not allow us to rule out the old age solution in the specific case of HD113337.

Outer disk geometry
We adopted a simple approach to estimate the disk extent at 70 µm (Fig. 3, left panel) by assuming that the disk emission can be described by an axisymmetric model, like a Gaussian ring defined by the peak (R p ) and the width (FWHM) of the ring (R w ). The disk has a total flux, F tot , at 70 µm, and its midplane is assumed to be inclined by an angle of i from face-on (i.e., i = 0 • ), with the major axis along a position angle (P.A.). We generated a series of high-resolution model images and convolved them with the observed point spread function (PSF) derived from the calibration stars. We then determined the best-fit parameters (five free parameters) by comparing the convolved model images with the observation using a χ 2 statistic. The best-fit parameters are: R p = 85 ± 20 au, R w = 30 ± 20 au, i = 25 •+5 • −15 • , P.A. = 128 • ±5 • , and F tot = 175 ± 12 mJy. The right panel of Figure 3 shows the image residuals, all within ±1σ after the subtraction of the best-fit model.
To make sure that the best-fit disk parameters are consistent with the observed SED, we computed the model SED using the derived geometric parameters (i.e., R p and R w ). Since the disk geometric parameters are derived using the PACS 70 µm data, we only try to reproduce the SED longward of 70 µm. Assuming the disk is optically and geometrically thin and has a surface density distribution best described by the derived Gaussian ring, we are able to reproduce the bulk part of the SED using icy silicate grains (the icy grain model from Su et al., 2015, for the HD95086 system). The grain size distribution is assumed to be a power law form, ∼ a −3.5 , where a is the grain radius with a minimum a min and maximum a max cutoffs. We found that a min of Fig. 6. Spectral energy distribution (SED) of the debris around HD113337, composed of broad-band photometry and midinfrared spectrum after the removal of the stellar photosphere. Blue diamonds are the Spitzer/MIPS photometry from Moór et al. (2011), green dots (with uncertainties shown in grey area) are the Spitzer/IRS spectrum from Chen et al. (2014), purple diamonds are the Herschel/PACS fluxes from this study. Also plotted are the JCMT/SCUBA2 upper limits (Holland et al., 2017) and IRAM/MAMBO2 1.2 mm from Moór et al. (2011). The mid-and far-IR broad-band photometry can be described by a simple blackbody emission of 60 K (thin grey line); however, it is slightly too high compared to the IRS spectrum. The Gaussian ring (GR) SED is shown as the blue dashed line. The disk SED from ∼20 µm to 1.2 mm is best described by the combination of a cold GR plus a 80 K blackbody emission (see text for details).
∼2 µm and a max of 1 mm can fit the far-IR SED well (Fig. 6). The minimum grain size (2 µm) is roughly the radiation blowout size assuming a bulk density of 1.7 g cm −3 , a typical minimum grain size in a collisional cascade debris disks. A total dust mass for this cold disk is 7.3×10 −3 M Earth (up to 1 mm grains). We note that the model cold-disk SED does not fit the MIPS 24 µm and IRS data, which might be related to the warm component reported by Chen et al. (2014). To explore this possibility, we tried to fit the mid-IR part of the SED with a simple blackbody function, and found that a blackbody emission with a temperature of 80 K represents the mid-IR SED well (Fig. 6). The 80 K emission is too cold compared to the warm component derived by Chen et al. (2014). It might be an intermediate separate component in the disk structure, which will not be the first time that a complex disk structure is inferred (e.g. Eri, Su et al., 2017). Alternatively, this mid-IR emission can also arise from a small amount of dragged-in grains from the cold disk under the influence of Poynting-Robertson (P-R) drag (e.g. van Lieshout et al., 2014). Since our Herschel data do not have enough resolution to spatially resolve the inner edge of cold disk, both scenarios are possible. for the candidate companion HD113337 c we obtain m p sin i = 7.2 ± 0.5 M Jup (without significant differences between the two stellar mass solutions). Understandably, the m p sin i do not change significantly with respect to the values we derived in BO19+ (i.e. 3 ± 0.3 M Jup and 6.9 ± 0.6 M Jup for HD113337 b and c, respectively), due to the little difference on the adopted stellar mass and parallax values. HD113337 b and candidate c orbit too close to the primary to be detected or even mass-constrained with an imager such as the LBTI. However, we can use the debris disk inclination value that we derived from the disk modeling (i.e., i = 25 •+5 • −15 • ) as a possible starting point. If we assume that the GP(s) and the partially resolved outer disk orbit within the same plane, we can then estimate their true mass(es). Assuming such a system inclination, the true mass of HD113337 b would then be 7.2 +4.2 −1.5 M Jup . The true mass of the tentative HD113337 c would be 16.4 +9.6 −3.4 M Jup . These true masses would be more than twice the value of the corresponding minimal masses from BO19+, in agreement with the small inclination considered here. We emphasize that this hypothesis remains widely speculative at this stage. Yet, if confirmed, this would make HD113337 GP(s) very massive planetary companions. This would be in agreement with the trend of higher GP masses with increasing stellar masses predicted by the core-accretion theory (Kennedy & Kenyon, 2008. . The four plots correspond to our four (age; i) cases (∼20 Myr or ∼1 Gyr, and i within 10 to 30 • or uniform i distribution, respectively). On each plot, the contour colors (from white to black) correspond to a higher or lower (respectively) detection probability of additional companions, as indicated by the numbers: 1 indicates a detection probability over 99%, 0.9 over 90%, etc. We display the assumed true masses of HD113337 b and c (full and empty red dots, respectively) on the top plots (inclination assumed to be 25 • ), and their minimal masses on the bottom plots. On the top of each plot, we indicate the outer debris disk extension resolved in this study with a red band. We also indicate the disk position previously assumed from SED fits from the literature: from Moór et al. (2011) (purple triangle), Rhee et al. (2007) (yellow triangle), and the inner disk component from Chen et al. (2014) (green triangle), respectively.

Constraints on actual and possible companions
Combined mass detection probabilities -From the LBTI images, we estimated the flux-losses at each separation associated to the ADI process (Bonnefoy et al., 2014b) and derived detection limits (1-D) and detection maps (2-D) in contrast. We classically converted our detection limits in contrast into masses using the star distance, the WISE W1 magnitude (Cutri & et al., 2014) as a proxy of the L'-band star magnitude, and the COND tracks (Baraffe et al., 2003). We considered two respective ages of 20 Myr and 1 Gyr, roughly corresponding to the young and old age solutions derived from our stellar fundamental parameter analysis (Sect. 4.1). We display the derived mass detection limits in Fig. 7. The difference in terms of achieved companion sensitivity between our two age solutions (i.e. roughly one order of magnitude) highlights the importance of an accurate age determination.
We brought additional constraints on the mass of possible additional companions by combining our contrast detection (2-D) maps with radial velocities (RV). For this purpose, we used the Multi-epoch multi-purpose Exoplanet Simulation System (MESS2, Lannier et al., 2017) tool. MESS2 generates populations of synthetic planets with masses and orbital parameters within pre-defined ranges through a Monte Carlo simulation. For each of the synthetic planets, the synthetic RV signal generated at the RV observation epochs, and the simulated planet projected separation at the image's observation epoch are simultaneously compared to the RV and imaging data, respectively. With respect to "classical" mass detection limits derived from contrast (as in Fig. 7), the advantage of this approach is twofold: (1) it allows to explore different hypotheses on the companion orbital properties; and (2) it allows to assess the companion mass detection probabilities in the combined separation range covered by RV and imaging (i.e. from the star's close proximity out to the imager's field of view). MESS2 was successfully applied to famous systems such as AU Mic (Lannier et al., 2017), HD95086 (Chauvin et al., 2018), β Pictoris , and GJ504 .
We applied MESS2 to the HD113337 SOPHIE RV data set detailed in BO19+ (Sect. 2). We mainly used the RV corrected from the 2-planet Keplerian fit performed by BO19+, assuming that HD113337 c is an actual planetary companion. The RV analysis within MESS2 relies on the periodogram-based Local Power Analysis method (LPA, Meunier et al., 2012). Regarding the LBTI data, we used the contrast 2-D maps converted into masses considering the two possible age solutions of 20 Myr and 1 Gyr derived from the evolutionary track analysis (Sect. 4.1). In addition to the two age solutions, we also considered two synthetic planet inclination distributions: first, a uniform distribution between 0 • and 90 • (i.e., no assumption on i); and second, a narrow inclination range between 10 • and 30 • (assuming that the synthetic planets orbit within the same plane as the resolved outer disk). Thus, we mainly tested four MESS2 simulations. We show the deduced detection probability curves in Fig. 8.
The SOPHIE RV data set allows us to rule out any additional companion to the known GP(s) with a mass above ∼5 M Jup up to 9 au if considering a system inclination around 25 • , and still 90% of them if considering a uniform inclination distribution. We remind that we used RV residuals of a 2-planet Keplerian fit, explaining this good sensitivity. We overplotted the true masses (or minimal masses) of the known companions on Fig. 8 for the 25 • and uniform inclinations, respectively (top and bottom plots). The sudden sensitivity gap at ∼8 au is explained by the current impossibility to simulate the RV signal of synthetic planets with orbital periods longer than twice the RV data set time span within MESS2 . We are not able to fully "bridge" the gap between the RV and imaging separation domains in any of the four simulations. Understandably, the detection probabilities provided by the LBTI images are best when considering the young age solution and close to pole-on inclination (top left plot). In this case, ∼100% of additional companions above ∼10 M Jup are excluded between ∼25 and ∼35 au, and ∼90% between ∼15 and ∼45 au. Increasing the age to 1 Gyr or assuming no hypothesis on the system inclination significantly reduce our sensitivity to companions within the separation domain covered by the LBTI. Compared to our "classical" (1-D) mass detection limits (Fig. 7), our MESS2 detection probabilities show a decreased companion sensitivity within the separation range covered by the LBTI images. This can be expected as we computed the former 1-D detection limits as if assuming the HD113337 system was seen face-on (i.e. the most favourable case for imaging companions). We emphasize that we used only one epoch of observation regarding the imaging data within our MESS2 analysis. A way to increase the sensitivity to companions within the imaging separation domain is to combine multiple high-contrast images acquired at different epochs (see e.g. Lannier et al., 2017.
We additionally tested the impact of assuming the presence of only one RV-detected GP in the HD113337 system (HD113337 b, confirmed planet) on our MESS2 detection probabilities. To do so, we used the SOPHIE RV data corrected from the 1-planet Keplerian model corresponding to planet b only, while keeping the longer-term RV variability (see more details in BO19+). We performed only one simulation, considering the young age solution and an inclination between 10 • and 30 • . In this case, the detection probabilities at the shortest separations are significantly degraded (Fig. 9), increasing the gap between the separation ranges covered by RV and direct imaging. This is expected as the computation of the RV detection probabilities within MESS2 is based on the analysis of the RV Lomb-Scargle periodogram (Meunier et al., 2012, Lannier et al., 2017.

Conclusion
We combined different techniques to explore and bring constraints on various aspects of the HD113337 system. New optical long-based interferometric measurements allowed us to measure the linear radius of HD113337 with a precision better than 3%. By using the new Gaia DR2 parallax and computing the star's bolometric flux, we were able to derive two very distinct isochronal age solutions for the system. The first (young) solution corresponds to an age of ∼14-21 Myr, while the second (old) one to an age of ∼0.8-1.7 Gyr. However, as often with such age degeneracies, we could not definitely settle the question of HD113337 age. For the first time, we were able to (partially) resolve HD113337 outer debris disk and to model its radius, its radial extension and, very interestingly, its inclination. We also found hints of possible inner disk components, which would be in agreement with previous SED studies. Next, we took new high-contrast images of the system's outer environment with the LBTI imager. We used both the deduced contrast limits and previous RV data to explore the complete GP-type companion (mass, separation) range up to 80-100 au from HD113337. At the same time, we took advantage of the age solutions and disk inclination value that we found to characterize the corresponding sensitivity to companions. Interestingly, this allowed us to deduce hints of the possible true masses of the HD113337 b confirmed GP and of the candidate HD113337 c. Furthermore, we were able to bring the first constraints on the presence of additional undetected companions at larger separations using the MESS2 tool.
While it was not the main topic of this study, an important issue with the HD113337 system is to determine once and for all if the candidate companion HD113337 c is a real one (see BO19+). At this time, we are carrying out an additional long-term RV monitoring of HD113337 with the SOPHIE spectrograph to increase our RV time span. This could allow us to remove the ambiguity of the RV and spectral line profile long-term signals. Another compelling possibility is to combine our RV data with astrometric data from HIPPARCOS and/or Gaia. The tentative second GP that might orbit around HD113337 is an ideal target within this context. Furthermore, the RV + astrometric data combination would bring more constraints on the system's inclination and might even allow to confirm if the GP(s) orbit within the same plane as the outer debris disk. If present, and if seen inclined, HD113337 c would be a very massive planet, most probably formed through core-accretion while close at the same time of the commonly considered mass boundary between GP and brown dwarf companions. Finally, the combined RV + imaging analysis that we carried out in this study can be extended by using multi-epoch high-contrast imaging data, which would allow to fully "close the gap" with the RV sensitivity domain. HD113337 has already been included in several high-contrast imaging surveys. To conclude, we consider that HD113337 constitutes an exciting and rich system to further explore. It could make for a useful contribution for both stellar physics (with regard to stellar age determination), GP formation and evolution as a function of stellar properties, possibly multi-component debris disk studies, and planetary-disk interactions.