Planet and star synergy at high-spectral resolution. A rationale for the characterization of exoplanet atmospheres The infrared

Context. Spectroscopy of exoplanet atmospheres at high-resolving powers is rapidly gaining popularity for measuring the presence of atomic and molecular species. While this technique is particularly robust against contaminant absorption in the Earth’s atmosphere, the non-stationary stellar spectrum, in the form of either Doppler shift or distortion of the line proﬁle during planetary transits, creates a non-negligible source of noise that can alter or even prevent detection. Aims. Our aim was to use state-of-the art three-dimensional stellar simulations to directly remove the signature of the star from observations prior to cross correlation with templates for the planet’s atmosphere, which are commonly used to extract the faint exoplanet signal from noisy data. Methods. We computed synthetic spectra from 3D simulations of stellar convection resolved both spatially and temporally, and we coupled them with an analytical model reproducing the correct geometry of a transiting exoplanet. We applied the method to the early K-dwarf, HD 189733, and re-analyzed transmission and emission spectroscopy of its hosted exoplanet. In addition, we also analyzed emission spectroscopy of the non transiting exoplanet 51 Pegasi b, orbiting a solar-type star. Results. We ﬁnd a signiﬁcant improvement in planet detectability when removing the stellar spectrum with our method. In all cases, we show that the method is superior to a simple parametrisation of the stellar line proﬁle or to the use of 1D stellar models. We show that this is due to the intrinsic treatment of convection in 3D simulations, which allows us to correctly reproduce asymmetric and blue-shifted spectral lines, and intrinsically model center-to-limb variation and Rossiter-McLaughlin effect potentially altering the interpretation of exoplanet transmission spectra. In the case of 51 Pegasi b, we succeed in conﬁrming a previous tentative detection of the planet’s K -band spectrum due to the improved suppression of stellar residuals. Conclusions. Future high-resolution observations will beneﬁt from the synergy with stellar spectroscopy and can be used to test the correct modeling of physical processes in stellar atmospheres. We highlight key improvements in modeling techniques and knowledge of opacity sources to extend this work to shorter wavelengths and later-type stars.


Introduction
The remote atmospheric characterisation of planets outside our solar system (exoplanets) is considered a key milestone in unravelling their physical and chemical composition (Miller-Ricci et al. 2009), their formation scenarios (Madhusudhan 2012;Piso et al. 2015;Eistrup et al. 2018), and ultimately the presence of conditions amenable to life (Schwieterman et al. 2018). Albeit challenging, these observations have reached a level of maturity where multiple observing techniques can be applied to both space-and ground-based telescopes. Among these, highresolution spectroscopy (HRS) at resolving powers R > 25 000 led to the detection of molecular (CO, H 2 O, CH 4 , HCN, TiO) and atomic (H, He, K, Na, Mg, Fe, Ti) species in a dozen exoplanets (see e.g. Birkby 2018, and references therein for a recent review). There are two unique aspects of HRS that make it particularly suitable for exoplanet characterisation. Firstly, at high resolving powers, molecular species are partially resolved into a dense forest of lines that can be robustly identified by linematching techniques such as cross-correlation. Secondly, the planet is subject to a detectable Doppler shift as it moves along the orbit, which can be used to solve for the orbital inclination of non-transiting systems similarly to spectroscopic binary stars (e.g. Brogi et al. 2012). Doppler shift and broadening also constrain atmospheric winds (Louden & Wheatley 2015;Brogi et al. 2016;Flowers et al. 2019) and overall bulk-planet rotation (Snellen et al. 2014;Schwarz et al. 2016).
Initially limited to the brightest stars, HRS has recently achieved a substantial increase in sensitivity due to the advent of new, performing infrared spectrographs mounted at large and medium-size telescope facilities, including GIANO (Origlia et al. 2014), CARMENES (Quirrenbach et al. 2014), SPIRou (Artigau et al. 2014), and soon CRIRES+ (Follert et al. 2014). Coupled with the progressive discovery of more planets orbiting bright stars, either from radial-velocity or transit surveys (e.g. Motalebi et al. 2015;Gandolfi et al. 2018), this means that the sample of exoplanets potentially detectable with HRS is now on the order of dozens, and will likely be of the order of hundreds after the end of the TESS mission (Barclay et al. 2018). Importantly for future studies, it has been shown that HRS with ground-based facilities and space-borne, low-resolution spectroscopy are highly complementary and can be combined Brogi & Line 2019) to improve contraints on the chemical make up and thermal structure of atmospheres.
Owing to the ability to disentangle sources with different Doppler signature, HRS is particularly suitable to diagnose potential sources of spurious signals, which can severely complicate the interpretation of exoplanet spectra. One of these is the non-uniformity of the planet-hosting stars. Their photosphere is covered with a complex and stochastic pattern associated with convective heat transport (i.e. granulation). Convection manifests in the surface layers as a particular pattern of downflowing cooler plasma and bright areas where hot plasma rises (Nordlund et al. 2009). Convection is a difficult process to understand, because it is non-local, and three-dimensional, and it involves nonlinear interactions over many disparate length scales. In this context, the use of numerical three-dimensional (3D) radiative hydrodynamical (RHD) simulations of stellar convection is crucial, but has only become possible in recent years with the increase of computational power, resulting in large grids of simulations covering a substantial portion of the Hertzsprung-Russell diagram (Magic et al. 2013;Trampedach et al. 2013;Beeck et al. 2013a;Ludwig et al. 2009). The use of those simulations has proven that the convection-related surface structures have different size, depth, and temporal variations, depending on the stellar type (Tremblay et al. 2013;Beeck et al. 2013b;Magic & Asplund 2014). More importantly, the related activity (in addition to other phenomena such as magnetic spots, rotation, dust, etc.) has an impact in stellar parameter determination Creevey et al. 2012;Chiavassa et al. 2012), radial velocity (Bigot & Thévenin 2008;Chiavassa et al. 2011;Allende Prieto et al. 2013), chemical abundances determinations (Asplund et al. 2005Caffau et al. 2011), photometric colours (Chiavassa et al. 2018;Bonifacio et al. 2017), and on planet detection (Magic et al. 2015;Chiavassa et al. 2017).
In the context of exoplanet studies, the use of 3D RHD simulations has already been endeavoured in several works. Chiavassa et al. (2015) showed that 3D simulations are better suited than ad-hoc limb-darkened laws for the interpretation of transit light curves in terms of ingress and egress slopes as well as the emerging flux. Chiavassa et al. (2017) continued this work evaluating the impact of granulation from optical to IR wavelengths for several planet/star systems and determined the granulation error budget on the determination of the planetary radius. Cegla et al. (2016) studied the impact of the stellar surface on the Rossiter-Mclaughlin effect, differential rotation, and the resulting star-planet alignment. Cegla et al. (2018) inspected the impact of the granulation on the resulting line profiles, and found induced center-to-limb variations in shape and net position. Dravins et al. (2017Dravins et al. ( , 2018 scrutinised the possibility of studying stellar microvariability, highlighting small surface segments, hidden behind the planet transits, with high spectral-resolution observations.
In this work, we demonstrate that 3D RHD simulations can already be applied to correct existing HRS observations of exoplanets and lead to a noticeable improvement in the detectability of their atmospheres compared to uncorrected spectra, or spectra corrected with 1D stellar models. Section 2 introduces the numerical methods used, Sect. 3 evidences the stellar intrinsic spatial and temporal variability, Sect. 4 explains the observational tools developed, and Sects. 5 and 6 report two direct applications of our approach.

Simulations
The stellar surface convection simulations used in this work are calculated using the STAGGER CODE, (Collet et al. 2011;Nordlund et al. 2009) which is a state-of-the-art (magneto) hydrodynamic code that solves the time-dependent hydrodynamic equations for mass, momentum, and energy conservation, coupled with the 3D radiative transfer equation, in order to correctly account for the interaction between the radiation field and the plasma. The equations are solved on a staggered mesh where the thermodynamical scalar variables (density, internal energy, and temperature) are cell-centered, while the fluxes are defined on the cell faces. This scheme has several numerical advantages when simulating surface convection. It is robust against shocks and ensures conservation of the thermodynamic variables. The domain of the simulation contains an entropy minimum at the surface that is sufficiently deep to ensure a flat entropy profile at the bottom. The code uses periodic boundary conditions horizontally, and open boundaries vertically. At the bottom of the simulation, the inflows have a constant entropy and pressure. The outflows are not tightly constrained and are free to pass through the boundary. The code is based on a sixth-order explicit finitedifference scheme, and a fifth-order interpolation. It employs realistic input physics: the equation of state is an updated version of the one described by Mihalas et al. (1988), and the radiative transfer is calculated for a large number over wavelength points merged into 12 opacity bins (Nordlund 1982;Skartlien 2000;Magic et al. 2013). They include continuous absorption opacities and scattering coefficients from Hayek et al. (2010), as well as line opacities described in Gustafsson et al. (2008), which in turn are based on the VALD-2 database (Vienna Atomic Line Database; Stempels et al. 2001) of atomic lines. A solar chemical composition is assumed ).

Synthetic spectra in the infrared region
We used the 3D pure-LTE radiative transfer code OPTIM3D (Chiavassa et al. 2009) to compute synthetic spectra from the snapshots of the RHD simulations reported in Table 1 and taken from the STAGGER-grid (Magic et al. 2013). The code takes into account the Doppler shifts due to convective motions. The radiative transfer equation is solved monochromatically, using pre-tabulated extinction coefficients as a function of temperature, density, and wavelength. The lookup tables were computed for the same chemical compositions as the RHD simulations using the same extensive atomic and molecular continuum and line-opacity data as the latest generation of MARCS models (Gustafsson et al. 2008). We assumed microturbulence equal to zero, since the velocity fields inherent in RHD simulations are expected to self-consistently and adequately account for nonthermal Doppler broadening of spectral lines (Asplund 2000). More details about OPTIM3D can be found in Chiavassa et al. (2009Chiavassa et al. ( , 2010. As in Chiavassa et al. (2018), we computed spectra with a constant resolving power of λ/∆λ = 300 000 from 22 850 to 23 900 Å, casting vertical rays through the computational box Table 1. 3D RHD simulations from STAGGER-grid (Magic et al. 2013)  Notes. T eff are from Chiavassa et al. (2018). In addition, a temporal average is also performed over ten snapshots adequately spaced so as to capture several convective turnovers. The final resulting spectra are: (i) a temporally averaged intensity spectrum at different µ and φ-angles, which will be used to model the changes in the disk-averaged stellar spectrum during the transit of exoplanet HD 189733 b (Sec. 5.1); (ii) a temporally-and spatially-averaged flux ( Fig. 1), which will be used to correct the stellar spectra in datasets targeting the thermal emission of exoplanets (Sects. 5.2 and 6).

Intrinsic stellar variability and its impact in the infrared
The granulation pattern is associated with heat transport by convection. The bright areas on the stellar surfaces, the granules (Fig. 2, bottom panel), are the locations of upflowing hot plasma, while the dark intergranular lanes are the locations of downflowing cooler plasma. Additionally, the horizontal scale on which radiative cooling drives the convective motions is linked with the granulation diameter (Nordlund et al. 2009). One piece of evidence of the convective-related surface structures comes from the unresolved spectral lines, in particular at high spectral resolution. In fact, they combine important properties such as velocity amplitudes and velocity-intensity correlations, which affect the line shape, shift, and asymmetries. These structures derive from the convective flows in the solar photosphere and solar oscillations (Asplund et al. 2000a;Nordlund et al. 2009).  Table 1).
The solid red (intergranular lane) and blue (granule) lines display two particular positions (colored star symbols) extracted from the intensity map. Figure 2 displays an example of the spatially-resolved intensities corresponding to different regions across the granulation pattern. Correlations of velocity and temperature cause characteristic asymmetries of spectral lines as well as net blueor red-shifts depending on the area probed (Dravins 1987; A100, page 3 of 10 A&A 631, A100 (2019) Fig. 3. Spatial and temporal average spectrum (red) of one CO line (same as in Fig. 2) for K-dwarf simulation (Table 1) together with a 1D hydrostatic MARCS model (Gustafsson et al. 2008) with the same stellar and spectral line parameters. Gray 2005). Bright and rising convective elements (granule) being blue-shifted and contributing more photons than the cool dark shrinking gas (intergranular lanes) that are red-shifted (Dravins 1982). As a consequence, the resulting spatiallyaveraged absorption line appears blue-shifted, like in Fig. 3 (red line).
An evident and important difference when using RHD simulation with respect to traditional 1D, stationary, hydrostatic model stellar atmospheres such as MARCS (Gustafsson et al. 2008), ATLAS (Kurucz 2005), or PHOENIX (Husser et al. 2013) is the treatment of convection. One-dimensional models can treat convective energy transport in approximate manners (e.g. mixing-length theory, Böhm-Vitense 1958) that are all dependent on a number of free parameters. However, since the convection zone reaches up to the optical surface in main sequence stars, convection directly influences the spectrum formation, both by modifying the mean stratification and by introducing inhomogeneities and velocity fields in the photosphere. It is therefore important to properly account for their effects in order to extract accurate and precise information from the analysis of stellar spectra. This is displayed in Fig. 3, where the shape of the CO line is entirely symmetric and centered to 0 km s −1 in the case of the 1D model (black line) and blue-shifted in the case of RHD simulation (red line). Each line has unique fingerprints in the spectrum that depend on line strength, depth, shift, width, and asymmetry across the granulation pattern depending on their height of formation and sensitivity to the atmospheric conditions. In this context, the line strength plays a major role (Asplund et al. 2000b).
One last important point of concern is the temporal variation of spatial-averaged spectral lines during planet transits. Figure 4 displays the uncorrelated temporal variations over slightly more than one hour of simulation time. This variability has already been detected on several occasions. For example, it has been observed that the Sun's total irradiance varies on all timescales relevant for transit surveys, from minutes to months (Aigrain et al. 2004). Moreover, granulation manifested on photometric variability of the SOHO quiet-Sun data (ranging from 10 to 50 ppm, Jenkins 2002;Frohlich et al. 1997), which were explained by Chiavassa et al. (2017) with the same RHD simulations used in this work. In the next sections, we apply synthetic spectra computed for 3D RHD simulations to correct existing HRS observations of exoplanets and check the improvement in the detectability of their atmospheres.  Table 1). The increasing time changes with colors from blue to red (ten snapshots and 70 min spanned in time).

Observations
We applied the models described in the previous section to correct the stellar spectrum imprinted on real spectroscopic data of bright exoplanet host stars. We devote this section to detailing the instrumental setup in common between data sets, and the basic data analysis applied to the spectra. The spectra were acquired with the Cryogenic Infrared Echelle Spectrograph (CRIRES, Kaeufl et al. 2004) mounted at the Nasmyth-A focus of the 8.2-m ESO Very Large Telescope UT1. The spectrograph was set to observe with a relatively narrow slit of 0.2 in order to achieve the maximum resolving power of R ≈ 100 000. A standard observing sequence with ABBA nodding pattern was used for these observations, with the telescope slewing by 10 along the slit between nodding positions A and B. Different images (A-B or B-A) were therefore used to subtract the contribution from sky emission lines and thermal backgrounds. One peculiar aspect of these observations is that the temporal information is preserved, meaning each couple of A-B or B-A frames is combined, and the corresponding 1D spectrum extracted, which differs from the common practice of co-adding all the spectra obtained during an observing night. This strategy allows us to have a time resolution between one and two minutes per combined exposure, enough to prevent the change in orbital radial velocity of exoplanets from broadening the planet spectral lines.
Extraction of the 1D spectra is obtained via the standard ESO pipeline. Each dataset described below has been published at different stages of the lifetime of CRIRES, and therefore different pipeline versions have been used. The data of 51 Pegasi ) and the thermal emission data of HD 189733  were processed with version 1.11.0 of the pipeline. The transit data of HD 189733 were instead processed with version 2.1.3 (Brogi et al. 2016) of the pipeline.
CRIRES images spectra on four detectors are physically separated by a few mm. This means that, although our data cover the interval between 2.27 and 2.35 µm, there are corresponding small gaps in the spectral coverage. Our data target most of the 2-0 rovibrational band of carbon monoxide (CO). In hightemperature exoplanet spectra, water vapour is also expected to produce significant spectral lines and had indeed been reported multiple times. Methane could also contribute with a dense forest of weaker lines, but so far its presence has not been confirmed at these wavelengths.   (Brogi et al. 2012), HD 209458 in red , 51 Pegasi in gold , and HD 189733 in blue (Brogi et al. 2016). The amplitude of stellar CO lines increases with decreasing effective temperature as expected, while their position varies due to the different barycentric + systemic velocity of the targets for different targets and/or nights of observations. In contrast, lines formed in the Earth's atmosphere (telluric lines) are stationary in wavelength and only slightly vary in depth due to airmass effects.
Interestingly, since the same CO lines are also present in the spectrum of the two stellar hosts examined in this paper, there is a chance of contaminating the planet signal with uncorrected stellar lines when the planet-star differential radial velocity is small. For planets on circular orbits, this happens during transit but also close to superior conjunction. A small fraction of a CRIRES spectrum is shown in Fig. 5, with telluric lines -mostly due to methane in this case -absorbing at a fixed wavelength and with approximately constant intensity, while the stellar CO lines vary in position and intensity according to the stellar systemic velocity and spectral type.

Data calibration pipeline
To enable accurate removal of the stellar spectrum, not only do we need to obtain a reliable wavelength solution, but also to measure the instrumental profile of CRIRES for each individual spectrum. In preparation for these steps, we use the ESO Sky Model calculator to precompute a coarse grid of telluric models in airmass (in the range 1.4-2.9, enough to cover the observations) and precipitable water vapour (PWV, in the range 0.5-3.5 mm).
The pixel-wavelength solution of CRIRES is known to drift by up to a pixel (approximately 0.1 Å at the wavelength of these observations) during a night. This means that arc frames taken at the end of the night are insufficient to provide a stable and accurate solution. The position and known wavelength of telluric lines imprinted on each spectrum is instead used as simultaneous reference to solve for the pixel-wavelength outcome. If spectral lines from the parent star are also prominent, these can lower the quality of the solution due to blends and to the non-negligible change in the barycentric velocity of the observer during the observations. Although stellar spectral lines have been neglected in previous literature, in this work, we revised the calibration pipeline in order to include the effects of a variable (in Doppler space) stellar spectrum. We proceeded as follows: (i) we started from a guess wavelength solution, which is the solution published in the literature for the three sets of data. As in past work, we assume a quadratic pixel-wavelength solution for each of the CRIRES detectors. This means that any wavelength solutions can be completely described by a triplet of points (x 1 , λ 1 ), (x 2 , λ 2 ), (x 3 , λ 3 ), where x is the pixel position and λ is the corresponding wavelength. We adopt (x 1 , x 2 , x 3 ) = (268, 512, 768) to evenly sample the 1024 pixels of a CRIRES detector. (ii) To assess the quality of the wavelength solution, we compared the position of spectral lines in the observed spectra to those in model spectra via cross correlation. As results are almost exclusively dependent on line position, we did not need to exactly match the amplitude of spectral lines at this stage. Therefore, we chose a telluric model as close as possible to the average airmass of each night of observations, and with PWV = 1.5 mm, the median value at Cerro Paranal. (iii) We launch a Markov chain Monte Carlo (MCMC) using the Python package emcee (Foreman-Mackey et al. 2013) where we vary the λ i (3 parameters) at each step of the Monte Carlo. We interpolate the telluric and stellar models (the latter Doppler-shifted accounting for the systemic, barycentric, and orbital radial velocity) to the tested solution, cross correlate the combined model with the data, and translate the cross correlation value into log-likelihood, as in Zucker (2003), to drive the MCMC. (iv) We re-grid the data to a common wavelength solution with constant steps in ∆λ/λ. This is necessary in order to correctly estimate the instrumental profile, as explained below.
The use of a MonteCarlo technique allows us to robustly assess the precision of the wavelength calibration. Each of the λ i can be determined with a typical error between 5 × 10 −4 and 10 −3 nm at these wavelengths, which is four-eight percent of the size of a CRIRES pixel.
We approximate the Instrumental Profile (IP) with a Gaussian profile with variable FWHM. The width of the profile depends on the changing resolving power of the instrument, due to, for example, varying slit illumination, pointing jitters, pressure and temperature changes, and any inaccuracies in combining A and B frames. Accurate IP estimation requires an unbroadened model as close as possible to the actual data. This is why, at this stage, we also fit for airmass and PWV to adjust the depth of the telluric spectrum to the observations. The stellar spectrum is further broadened with a rotational kernel corresponding to the literature values of the stellar projected rotational velocity, or v sin(i). Attempts to fit for v sin(i) as a free parameter resulted in no meaningful constraints within ±1 km s −1 .
We launch an MCMC on each observed spectrum where airmass, resolving power, and PWV are fitted at once (3 parameters). As for the wavelength calibration, the procedure uses the cross-correlation-to-likelihood mapping of Zucker (2003). Each telluric model spectrum is now obtained by bi-linearly interpolating across the pre-computed grid of models in logspace, and both telluric and stellar spectrum are broadened by the variable IP prior to cross correlation. Figure 6 shows the typical behaviour of resolving power and PWV over three different nights of observations, with error bars reporting the 1-σ confidence intervals from the MCMC.
At the end of the calibration process, the final wavelength solution and measured IP are used to remove the stellar spectrum from the data. This is done by shifting the appropriate stellar model (see details in each individual section below) based A100, page 5 of 10  on the three radial-velocity components listed above (systemic, barycentric, gravitational). The shifted model is then broadened by convolution with the IP and divided out from the spectra. Since the stellar model had its continuum previously normalised, no re-scaling is necessary before dividing it out.

Transmission spectra
The spectral sequence presented in this section include spectra taken during one transit of exoplanet HD 189733 b and thoroughly described in Brogi et al. (2016) and Flowers et al. (2019), to which we point the reader for further details. In this context, it is worth mentioning that out of the 45 spectra in the sequence, the first six are out of transit (before ingress), and the remaining 39 cover ingress, mid-transit, and most of the egress of the planet.
Already, Brogi et al. (2016), pointed out that the stellar CO lines were a significant contaminant in transmission spectra, as shown by the left panel of Fig. 7. They attempted a correction of the stellar spectrum by modeling the average line profile with micro-turbulence, macro-turbulence, and rotational broadening. This allowed them to pinpoint the CO planetary absorption at the expected planet maximum radial velocity (dashed lines, middle panel), albeit at relatively low significance, and with residual signal from the star. A more recent analysis of the same data based on the technique described in this paper (Flowers et al. 2019) achieved a much better correction and resulted in a unique and unambiguous identification of the planetary signal in CO alone (i.e. without the aid of the extra cross-correlation signal from water vapour), as shown in the right panel of Fig. 7. This improvement is reflected into a refined inference on the rotational rate and wind speed of exoplanet HD 189733 b.
As the process of stellar removal is only briefly explained in Flowers et al. (2019), we provide more details in this section. The metallicity of the star HD 189733 is calibrated on a small set of unblended stellar CO lines far from telluric lines. We find that the literature value of [Fe/H] = −0.03 is a good match to the strength of the stellar CO lines, although a solar metallicity spectrum is not visually distinguishable at the noise level of these data.
The application of the stellar model spectra to transit observation requires a numerical model capable of reproducing the correct geometry of the transit. The model is based closely on the work of Brogi et al. (2016), and it involves sub-dividing the stellar disk into a grid of pixels. Each pixel has a different stellar spectrum assigned, based on its value of (µ, φ) as defined in Sect. 2, and obtained by bi-linearly interpolating across the temporally-averaged grid of 3D simulations. Furthermore, the spectrum at each pixel is Doppler shifted based on the projected stellar rotational velocity of v sin(i) = 3.3 km s −1 (Triaud et al. 2009). Although in this work we assume rigid rotation for the star, the model makes it possible to incorporate differential rotation, which might be necessary once measurements of stellar radial velocity can constrain this quantity (see, e.g., Cegla et al. 2016). Additional radial velocity components, such as a gravitational redshift of 0.68 km s −1 , the reflex motion due to the orbiting exoplanet (with semi-amplitude 202 m s −1 ), and the combined systemic and barycentric radial velocity, are also applied to each time-resolved spectrum.
For observations outside of the planet transit, the total stellar spectrum is given by the sum of the spectra of each individual pixel. This temporally and spatially averaged spectrum will be utilised in Sect. 5.2 as well to correct dayside observations of the same system. During transit, the correct planet-star position needs to be computed in order to determine which pixels of the stellar surface are occulted by the planet disk. As explained in Brogi et al. (2016), we use the known spin-orbit angle of the system (i spin = 0.85 +0.28 −0.32 degrees, Triaud et al. 2009), transit impact parameter (b = 0.6631 ± 0.0023, Agol et al. 2010), and planet orbital phase to solve for the geometry. Once the occulted pixels of the stellar disk are zeroed, the corresponding stellar spectrum at that specific planet phase is computed again by summing all the non-zero pixels. The output stellar model, which is now time-resolved during the planet transit, naturally incorporates the distortions in the line profile due to center-to-limb variations and Rossiter-McLaughlin effect. As each 3D RHD simulation spectrum is normalised by the local continuum, the correction of the stellar spectrum is obtained by dividing the sequence of observed spectra through the sequence of modeled spectra, with no further scaling required.

Emission spectra
110 spectra of HD 189733 with planet "b" just before secondary eclipse (orbital phases 0.38-0.48) were analyzed by de . They reported a detection of the planet's thermal spectrum in CO at a S /N = 5, but also a significant contamination from stellar CO, which they solved by masking those portions of the spectra corresponding to strong stellar absorption. This was possible due to the significant difference between the stellar and planetary radial velocity far from secondary eclipse, however it still produced residual stellar noise due to the sharp edges of the masking filter applied to the data.
We use the 3D stellar spectra of HD 189733 (Table 1) Fig. 7. Total cross-correlation signal from carbon monoxide in transmission spectrum of HD 189733 b, as function of stellar removal applied to data. Left panel: uncorrected spectra show a candidate signal at expected planet maximum orbital velocity (K P ∼ 151 km s −1 , white dashed lines), but also significant stellar contamination at all values of K P from uncorrected CO stellar lines. Middle panel: application of a parametric 1D stellar model partially mitigates the problem and allows recovery of the CO signal at 5σ (Brogi et al. 2016). There is, however, still a residual stellar signal that hinders detection. Right panel: application of phase-resolved 3D stellar models accounting for the geometry of transit (this work and Flowers et al. 2019) allows us to unambiguosly detect the planet CO signal at >7σ of confidence level. observations as well. In this case, since the planet is not transiting, we do not need to accurately model the transit geometry, and a model averaged both spatially and temporally is adopted. This is once again shifted by the systemic and barycentric velocity at the time of observations, and gravitational redshift. Figure 8 shows the results of the analysis in the case of an uncorrected stellar spectrum (top panel), and after the stellar correction. In the latter case, we recover the signal of exoplanet HD 189733 b in CO at a S /N = 4.5, consistent with de , and no stellar residual above the S /N = 3. Conversely, without a correction of the stellar spectrum, the planet signal would be completely outshone by stellar residuals, preventing its unambiguous identification.

Confirming the atmospheric signature of exoplanet 51 Pegasi b
The last set of spectra analyzed in this work are described in detail in Brogi et al. (2013). It consists of three sequences of spectra taken on three different nights just before, around, and after superior conjunction of the planet. As 51 Pegasi b is not transiting, its identification in velocity space is complicated by the unknown orbital inclination and poorly-known time of inferior conjunction (or any equivalent reference time), which means a large parameter space in the rest-frame velocity (V rest ) versus semi-amplitude of the planet radial velocity (K P ) diagrams, such as those in Sects. 5.1 and 5.2, is allowed. The presence of stellar contaminating signal can thus significantly complicate or even prevent detection. In their original study, Brogi et al. (2013) attempted to scale an empirical solar spectrum to remove a large fraction of the stellar CO spectrum. This approach worked sufficiently well on two of the three nights, meaning, when the planet was far from superior conjunction, and its radial velocity departed from the systemic velocity by tens of km s −1 . However, on the third night (close to superior conjunction) no signal was detected, and no compelling evidence was found to justify the non-detection. We decided to revisit these data and remove the spectral signature of the star 51 Pegasi with application of the models described in this paper. Figure 9 shows the observed spectrum of 51 Pegasi after telluric lines have been removed with the model spectra in output of the ESO Sky Model Calculator, with the fitted values of airmass, PWV, and fitted IP. Overplotted are the 3D RHD simulations used to correct these spectra, which is clearly a good A&A 631, A100 (2019)

Wavelength [nm]
Normalised stellar flux Fig. 9. Average observed spectrum of 51 Pegasi obtained from entire sequence of CRIRES observations presented in Sect. 6, after removal of a telluric model (black line). The corresponding 3D spectrum from Table 1, broadened by v sin(i) = 1.5 km s −1 and by the measured CRIRES instrumental profile, is overplotted in red. match to most of the stellar lines, especially the CO lines clearly identifiable by the band-head near 2.323 µm. As 51 Pegasi has nearly solar physical parameters, we use solar models for this work, with the exception of stellar metallicity, which is set by visual inspection to the value of [Fe/H] = −0.15. Figure 10 should be directly compared with Fig. 2 of Brogi et al. (2013), and it shows the total cross-correlation signal from their best-fitting planetary model, which contains both CO and H 2 O when the stellar CO lines are not corrected (left panel), and when they are removed with the analysis described in this paper (right panel). Similarly to their work, the stellar correction manages to bring the residual stellar cross-correlation below the level of the noise, while the planet signal, barely visible on individual nights as a bright stripe of positive correlation around V rest = 0 and K P = 132 km s −1 , is preserved in the process. The Although the exoplanet cross-correlation signal is detected even without removing the stellar spectrum, it is confused in the strong stellar residuals (left panel). Conversely, after removal of the stellar spectrum, the planetary signal is the only unambiguous signature detected, and all three nights co-add constructively as expected.
overall performances of the stellar removal appear superior to the analysis in Brogi et al. (2013), and this is further demonstrated when the individual observing nights are co-added. Figure 11 shows the combined cross-correlation signal from the three observing nights (top row) and from the nights of October 16 and 17 only (bottom row). The stellar residuals make a striking difference in the co-added signal. Although a clear positive deviation in the total cross correlation signal appears at the known position of the planet regardless of whether the star is removed or not, in the latter case, its planetary nature cannot be unambiguously determined. The presence of strong stellar residuals near the planet position on the night of October 25, which is a consequence of the similar planet and stellar-radial velocities near-superior conjunction, makes the planet detection doubtful.
Conversely, when the stellar spectrum is removed, the planet 51 Pegasi b remains the only unambiguous source detected. Remarkably, the planet is now detected in all three nights of observation, with a noticeable gain in S/N (about 10%) when adding the night of October 25. This is a clear improvement from the non-detection originally reported by Brogi et al. (2013) Fig. 9) with 1D hydrostatic models (light blue) and 3D RHD simulation (black). From top to bottom: different CRIRES detectors, except detector number four, which is known to suffer from odd-even effect, and likely has an asymmetric instrument profile. Each bisector results from one of the 166 exposures observed on the night of October 16, 2010. on this night, which allows us to confirm the measurement of molecular absorption in K-band spectra of this planet for the first time. In Sect. 3, we showed that convection plays a crucial role in the formation of spectral lines and deeply influences the shape, shift, and asymmetries of lines in late-type stars. One way to detect the asymmetries in the line is the bisector, defined as the locus of the midpoints of the spectral line. a symmetric profile has a straight vertical bisector, and is the consequence of a complete absence of velocity-brightness correlations (i.e. the case of 1D hydrostatic models), while bisectors with C-shape are formed mostly in the upflows (granules), or reverse C-shape bisectors are generally formed in downflows (Dravins et al. 1981). In the literature, bisectors of several kinds of stars have already revealed the presence of asymmetries and wavelength shifts caused by the granulation (e.g. Ramírez et al. 2008;Gray 2009).
In this section, we aim to check whether our use of numerical stellar models is correct. Initially, we cross-correlate the stellar spectrum 51 Pegasi, resolved temporally, with both the 3D RHD simulation from Table 1 and the 1D hydrostatic model from MARCS package (Gustafsson et al. 2008) with the same stellar parameters as the 3D simulation. We run the cross-correlation at R = 10 6 , meaning 10 × the resolution of CRIRES, resulting in cross-correlation functions sampled with a step size of 0.3 km s −1 .
Subsequently, we compute the bisectors for the resulting cross-correlated functions as reported in Fig. 12 (light blue for 1D and black for 3D). A straight and zero-centered bisector would indicate that the model we used to correct the stellar signal perfectly matches the stellar properties of the star. The first point of concern is the shape of bisectors: blue and dark curves are similar across different detectors. This indicates that the CO lines are comparable, in term of shape, in 1D and 3D. The second point, more importantly, is the bisector shifts. The use of 3D simulation systematically returns values between zero and 50 m s −1 . This velocity shift with respect to zero might be due to the gravitational redshift (Lindegren & Dravins 2003). In this work, we assumed a gravitational redshift of 650 m s −1 for 51 Pegasi, and added this extra radial-velocity component when comparing the synthetic spectra to observations. However, the exact values for 51 Pegasi is unknown even though, for late-type dwarfs (log g ≈ 4.5), the gravitational shift ranges between 0.7 and 0.8 km s −1 , and it dramatically decreases with surface gravity down to 0.02-0.03 km s −1 for K giant stars with log g ≈ 1.5 (Allende Prieto et al. 2013). Another possibility for the residual shift of bisectors obtained from 3D spectra is the similar uncertainty in the systemic velocity, which is −33.218 +0.076 −0.077 km s −1 according to Birkby et al. (2017). In terms of our measurements, gravitational redshift and systemic velocity could technically be separated, because the former only influences the star, whereas the latter influences both stellar and planetary radial velocity. However, this would require a precision less than 100 m s −1 in the planetary radial velocity, which is nearly one order of magnitude below the current precision (just below 1 km s −1 for the best planet detections, see, e.g., Flowers et al. 2019). Therefore, for the purpose of this study, one can consider systemic velocity and gravitational redshift indistinguishable and producing just a net shift of the stellar lines.
Whatever the case, cross correlation functions obtained with the 1D model show a residual blueshift between −150 and −200 m s −1 , indicating that 3D models are substantially more adequate in accurately reproducing the photospheric velocity field of the star.

Conclusions
We used 3D RHD stellar convection simulations to provide synthetic stellar spectra resolved both spatially and temporally, and we coupled them with an analytical model of a transiting exoplanet accounting for the variable portion of the occulted stellar disk. We applied the method to VLT/CRIRES observations and removed the spectrum of two bright exoplanets' host stars: the early K-dwarf HD 189733 (covering transmission and emission spectroscopy of its transiting planet), and the solar-type 51 Pegasi (covering emission spectroscopy of its non-transiting planet).
Removing the stellar spectrum with our method led to a significantly improved detection of the exoplanet atmosphere. We show that our modeling is superior to a simple parametrisation of the stellar line profile or to the use of 1D stellar models. This is due to the fact that 3D RHD simulations correctly describe the stellar properties and dynamics. Three-dimensional spectra match the observed asymmetry and/or blue-shifted in spectral lines, and, intrinsically, they also model center-to-limb variation and the Rossiter-McLaughlin effect, potentially altering the interpretation of exoplanet transmission spectra. In the case of 51 Pegasi, we succeeded in confirming a previous tentative detection of the planet's K-band spectrum due to the improved modeling of the stellar residuals.
To further establish the synergy between stellar modeling and exoplanet observations, it is also fundamental to extend this work to bluer wavelengths covering the helium triplet at 10 830 Å (Oklopčić & Hirata 2018;Nortmann et al. 2018;Allart et al. 2018), and the alkali lines in the optical (Na and K doublet, see, e.g., Wyttenbach et al. 2015;Casasayas-Barris et al. 2017). However, this is also potentially complicated by the necessity of treating non-LTE effects (e.g., Amarsi et al. 2018), and by the general limitation in the knowledge of atomic and molecular linelists (Jofré et al. 2019).
Another important aspect is also the extension of the method described in this work to M-dwarf stars that are increasing in terms of number of detections, due to their abundance in the solar neighborhood, and the favorable planet/star contrast (Anglada-Escudé et al. 2016;Gillon et al. 2017). Preliminary 3D global (i.e., including the whole convective surface layers within the numerical box, star-in-a-box configuration) RHD simulations have been achieved in Allard et al. (2013) in the presence of rotation, essential to resolving the cloud-surface distribution and to induced potential variability. However, M dwarfs are magnetically active possibly with local suppression of convection by magnetic field lines emerging at the surface. Thus, global Magneto-RHD simulations are needed to identify the mixing length to use in classic 1D models to compensate for this suppression (Allard et al. 2013).
In conclusion, three-dimensional RHD simulations are now established as realistic descriptions for the convective photospheres of various classes of stars. The good and time-dependent representation of the background stellar disk is a natural and necessary step forward toward a better understanding of stellar properties. We have now proven that, in the context of exoplanet science, 3D RHD simulations are also extremely useful for a detailed and quantitative analysis of the atmospheric signatures of transiting and non-transiting planets.