© ESO 2020

1 Introduction

Providing context to our own Solar System and the diverse atmospheres found therein is one of the goals of exoplanet research. While we have detailed constraints on the gas giants in our own Solar System thanks to the Galileo (Mahaffy et al. 2000), Cassini (Lebreton & Matson 1992), and Juno (Bolton & The Juno Science Team 2010) missions, it is only now that we have access to the precision necessary for measuring exoplanetary atmospheres and the range of confirmed exoplanets to start a more in-depth discussion on atmospheric compositions (Fletcher et al. 2019). From hot giant planets like Jupiter down to Earth-sized exoplanets, the atmosphere is accessible as a thin ring via remote transit observations, among other methods. With these observations, a transmission spectrum can be built by taking the ratio of spectra taken during the transit and out of transit. The high-resolution transmission spectrum, showing the spectral lines of the planet atmosphere resolved individually, can be used to find new elements in the atmosphere and discuss different chemical compositions and dynamic patterns.

The sodium doublet is well-suited for this endeavour because it probes the atmosphere up to the thermosphere, thus allowing for a study of many orders of magnitude in pressure. The opportunity to use sodium as a high-altitude probe comes from its characteristic resonance doublet at 589 nm (Seager & Sasselov 2000).

For hot Jupiters, a plethora of sodium detections have been confirmed, starting with Charbonneau et al. (2002) at low resolution and Snellen et al. (2008), Redfield et al. (2008) up to more recent high-resolution observations (e.g. Louden & Wheatley 2015; Ridden-Harper et al. 2016; Wyttenbach et al. 2015, 2017; Casasayas-Barris et al. 2017; Jensen et al. 2018; Khalafinejad et al. 2018; Seidel et al. 2019).

These high-resolution sodium detections show significantly broader lines than a static model of the atmosphere including pressure and temperature broadening, would suggest. In the case of WASP 76b, the detected sodium doublet is ten times broader than the instrument response function (Seidel et al. 2019). Both WASP 49b and HD 189733b show similar broad features (Wyttenbach et al. 2017; Casasayas-Barris et al. 2017). In this study, we focus on the sodium feature of HD 189733b. To observe the sodium doublet with a high enough precision to resolve the line core, we used the HARPS (High-Accuracy Radial-velocity Planet Searcher) spectrograph at ESO’s 3.6 m telescope in La Silla, Chile, with a resolution of R = 115 000 (Mayor et al. 2003). A possible explanation for this broadening is Doppler broadening from strong winds in the atmosphere, as has been claimed for CO lines which probe a different altitude in the atmosphere (Brogi et al. 2016). Theoretical studies have suggested different wind patterns in the lower atmosphere, for example: super-rotational wind, where the atmosphere rotates faster than the planet itself (e.g. Showman & Polvani 2011); day-to-night side wind, where the wind on a tidally locked planet is upwelling on the dayside and downwelling on the nightside (e.g. Flowers et al. 2019); and, lastly, vertical wind patterns pushing the atmosphere away from the planet surface (e.g. Komacek et al. 2019; Debrecht et al. 2019).

Global circulation models (GCM) build a three-dimensional atmosphere and are capable of calculating model atmospheres with equatorial jets and varying temperature profiles in latitude, longitude, and altitude (Showman et al. 2018; Flowers et al. 2019; Steinrueck et al. 2019; Deibert et al. 2019). This level of sophistication means that in almost all cases, a high computational cost limits the possible number of executions. Particularly in the case of retrieval approaches, which require executions in the order of 100 000, full GCM models are challenging. In this work, we introduce a forward model for transmission spectroscopy, including wind patterns, that operates on a one dimensional grid. The simplification of a 1D calculation allows to link the forward model to a multi-nested sampling retrieval followed by a rigorous comparison of the different proposed wind patterns to the data.

In the first section, we introduce the Multi-nested η¹ retrieval code (MERC), highlighting the implementation of wind patterns, followed by a section dedicated to benchmarking the code on simulated data. We then move on to test the code on data obtained for HD 189733b and compare our findings to the current literature. Lastly, we summarise the code and our results in the conclusions in Sect. 5.

2 Multi-nested ETA Retrieval Code (MERC)

The structure of MERC is shown as a flowchart in Fig. 1 and consists of two parts, which are individually highlighted in this section: the forward model and the nested sampling. The forward model creates a model atmosphere based on input parameters and is an adaptation of the η (Ehrenreich et al. 2006) and, subsequently, the ^πη code (Pino et al. 2018), including wind broadening. This forward model is then linked to a nested sampling retrieval to explore the full parameter space.

2.1 Forward model

In the forward model (left panel of Fig. 1), a temperature pressure profile and element abundances are used as input to calculate a model transmission spectrum. Compared to the η and ^πη code, the focus of the forward model in MERC is on the execution speed. To combine a forward model with nested sampling retrieval, the forward model has to be executed more than 100 000 times for each run, making a fast forward model indispensable.

We achieve this speed-up by parallelizing and optimizing the code, restricting the wavelength range to a single or few lines (in this case the sodium doublet), and calculating the extinction coefficient in 1D instead of 3D (already implemented in ^πη).

The atmosphere is described by a 1D grid in altitude above the planet surface. For computational reasons, an evenly spaced grid in log space with 200 cells covers the pressure from the assumed surface pressure (see Sect. 2.4) to 10⁻¹² atm, far beyond the probing depth of the sodium doublet. The number of cells in pressure was set empirically as the lowest number of cells necessary to still have unchanged transmission spectra (difference < 1%) compared to runs with much finer grids. For each of these cells a temperature is calculated according to the different models. We explore an isothermal model and a temperature gradient under the assumption of hydrostatic equilibrium. Hydrodynamical models of hot Jupiter atmospheres show no better description of the data than hydrostatic models up to the lower part of the thermosphere, implying that the hydrostatic approach is a good assumption (e.g. Yelle 2004; Koskinen et al. 2013). Additionally, we are far below the sonic point of HD 189733 b. We implement the temperature gradient T′ linear in log (P) by setting a base temperature T_base at the base of our pressure scale and a top temperature T_top at the top of our pressure scale.

This temperature-pressure grid allows the calculation of the scale height, the gravity and the absolute height above the surface under the assumption of hydrostatic equilibrium. The estimation of the pressure at surface level is complicated by the degeneracy of the surface pressure and the sodium abundance, which we address in Sect. 2.4.

From these environmental parameters, together with line parameters intrinsic to the sodium doublet – the line centre, Einstein coefficient A₂₁ and oscillator strength – we calculate the cross section κ₀ of the sodium doublet in each cell of the atmospheric grid. The cross section grid is used to create an extinction coefficient grid of our atmosphere via: $$\kappa (\lambda ,z)=\chi (z)\cdot {\kappa }_0(\lambda ,z).$$(1)

The extinction coefficient κ is, therefore, the cross-section of the molecule or element in question κ₀ weighted by the relative abundance of the element by mass, also known as the mass-mixing ratio. In our case, the element in question is sodium. The relative abundance χ is a parameter that we also address in Sect. 2.4.

At this stage, we calculated the extinction coefficient on a 1D grid corresponding to different altitudes. To compute the model transmission spectrum we need the slant optical depth (τ(λ, b)) along the line of sight. The line of sight depends on the impact parameter b(z) that quantifies how far from the planet surface the line of sight cuts through the atmosphere. To save significant computing time, the 1D grid of opacities in altitude is stored and for each calculation along the line of sight, the corresponding altitude from the surface ($\tilde{z}$ in Fig. 2) is calculated via the current impact parameter b(z) and position along the line of sight (LOS) x. The altitude above the surface, $\tilde{z}$, is then used to select the corresponding extinction coefficient from the stored grid. This creates a two-dimensional extinction coefficient grid κ(λ, x, b) from the one-dimensional κ(λ, z), which we refer to as an atmospheric slice in the following. The geometry of this procedure is shown in Fig. 2 and further elaborated on in Pino et al. (2018). This two-dimensional grid containing opacities in both z and x can then be modified by winds (see blue box in forward model in Fig. 1), which we discuss in Sect. 2.1.

With the corresponding extinction coefficient in each bin, the slant optical depth is calculated by summing over all bins along the line of sight: $$\tau (\lambda ,b)={\sum }_x\kappa (\lambda ,x,b)\cdot \Delta x.$$(2)

The contribution of the atmosphere to the transmission spectrum, called the atmospheric equivalent surface of absorption (∑ (λ)), is created from the slant optical depths by iterating through the atmosphere with the line of sight in altitude (z). The altitude grid was defined earlier as a grid of 200 cells from z(i = 0) = R₀ to z(i = 200) being the top of the atmosphere linked to the pressure grid. The nomenclature for the atmospheric equivalent surface was selected to be in agreement with Ehrenreich et al. (2006).

The atmospheric equivalent surface of absorption from Eq. (3) of Ehrenreich et al. (2006) is: $$\sum (\lambda )={\sum }_{\text{i=0}}^{200}2\pi (b(i)+R_{\text{p}})\tilde{z}(i)(1-\exp [-\tau (\lambda ,b(i))-{\tau }_{\text{rayleigh}}(\lambda ,i)]),$$(3)

where the line of sight contributions are summed up for the full atmospheric contribution utilising the symmetry of the problem. τ_rayleigh(λ, i) is the optical depth generated by Rayleigh scattering of H₂.

The transmission spectrum ℜ′(λ) is calculatedas the contribution from the atmospheric equivalent surface and the opaque disk in LOS of the planet (Eq. (4) from Ehrenreich et al. 2006): $${\Re }^{\prime }(\lambda )=-\frac{\sum (\lambda )+\pi R_{\text{P}}^2}{\pi R_{*}^2},$$(4)

where ∑ (λ) denotes the atmospheric equivalent surface of absorption and $\pi R_{\text{P}/*}^2$, respectively, the opaque disk in LOS of the planet and the star.

In the next section, we introduce an additional step when calculating the 2D atmospheric slice, a shift in the wavelength grid, which has, so far, remained constant. The wavelength grid is Doppler-shifted by wind patterns in the bins where the extinction coefficient is calculated, adding an overall broadening in the transmission spectrum.

Fig. 1 Flowchart showing the information flow from forward model to the nested sampling module and vice versa in MERC. Green arrows are executed during the continuum retrieval, blue arrows are only executed during the line retrieval. Starting with the retrieval windows in green from Fig. 6 and the parameters that can be retrieved onthe continuum (see Table 2), priors are set in the nested sampling module and handed to the forward model. One possible combination from the prior range is then used to calculate the extinction coefficient τ in a 1D grid in altitude above the planet surface. The planet surface is defined as the opaque part of the planet at the wavelength of the green region and is set at the optical white light radius. The 1D grid can be seen as a slice of the atmosphere, which is then used to create the full atmosphere (see Sect. 2.1) and calculate a model transmission spectrum. The model spectrum is handed back to the nested sampling module, where it is compared to the continuum of the data. The multi nested sampling executes this loop over the prior parameter space until the parameters converge to a value within the prior ranges. It then produces a posterior distribution of the prior and the Bayesian evidence of the model. The posterior of the continuum retrieval is then used to restrict the prior of the retrieval on the line retrieval window (shown in blue). The same process as for the continuum is started again, but now with additional wind broadening in each cell of the atmospheric slice. At the end, the Bayesian evidence of the line retrieval for each model is used to rank them. This flowchart was loosely inspired on a similar layout in Brogi & Line (2019).

Fig. 2 Schematics of the extinction coefficient calculation in the line of sight (LOS) exemplified for the outermost bin in x for a vertical wind profile. The atmospheric wind pattern is assumed to be a constant radial velocity originating from the planet’s surface. γ is the angle running from the polar position towards the observer, R the assumed radius of the planet, b the impact factor and x the sum overall bins in the line of sight x = ∑Δx.

Wind broadening

When traversing the atmosphere towards the observer, the wavelength of the light is Doppler-shifted using the amplitude of the velocity in the line of sight v_LOS in each bin in x. This velocity stems from winds in the atmosphere since the transmission spectrum is calculated in the planet rest frame. To account for these winds only, a correction for the different velocities originating from the different rest frames of observer and observed object is required. The three wind patterns most interesting for hot Jupiters are a day-to-nightside (dtn) wind (see Fig. 3), a super-rotational wind (see Fig. 4), and a vertical upward wind (see Fig. 5), which we study separately and in combination in the following.

First, we demonstrate how the velocity shift is calculated in a specific bin for the three different wind patterns. We exemplify the full calculation for a simple upward pointing constant wind velocity, see also Fig. 2 for a visualization.

Assuming a constant vertical wind originating from the planet atmosphere the line of sight velocity depends on the constant amplitude of the wind (|v_ver |) and on the angle of the wind vector and the line of sight v_LOS = |v_ver|⋅ sinγ.

With $\tilde{z}=\sqrt{{(R+b)}^2+x^2}$, the distanceof the current bin from the planet centre (see Fig. 2), the angle γ is calculated from the similarity of the two triangles in Fig. 2 as $$\sin \gamma =\frac{x}{\tilde{z}}=\frac{x}{\sqrt{{(R+b)}^2+x^2}},$$(5)

leading to a velocity in line of sight depending only on the x and y position on the grid, which is to be iterated through in the code. For a vertical constant wind the wind strength in the line of sight is calculated as follows: $$\left|v_{\text{LOS}}\right|=\left|v_{\text{ver}}\right|\cdot \frac{x}{\sqrt{{(R+b)}^2+x^2}}.$$(6)

A super-rotational wind pattern throughout the atmosphere and a day-to-night side wind are calculated the same way and the sign is adjusted depending on the respective wind direction in each hemisphere: $$\left|v_{\text{LOS}}\right|=\pm \left|v_{\text{rot/dtn}}\right|\cdot \frac{(R+b)}{\sqrt{{(R+b)}^2+x^2}}.$$(7)

The Doppler shift of the wavelength grid is calculated for each bin by $${\lambda }_{\text{shifted}}={\lambda }_{\text{grid}}\cdot \left(1+\frac{\left|v_{\text{LOS}}\right|}{c}\right),$$(8)

where c is the speed oflight.

Depending on which kind of velocity profile is implemented, the line of sight velocity changes for each bin and via its Doppler broadening the entire transmission spectrum. The three Figs. 3–5 show how the velocity shift in the line of sight (dark green bins, left side of each figure) is translated into a pseudo-3D wind pattern (right side of each figure).

In the specific case of a radial velocity field, the symmetry of the field can be utilised by calculating the shift for the side towards the observer with + |v_LOS| and the side facing away from the observer with −|v_LOS| (see Fig. 5). The cases of a super-rotational wind and a day-to-night side wind are less intuitive. For the day-to-nightside wind (see Fig. 3) both the southern and the northern hemisphere, and the morning and evening limb, behave in the same way and we rotate our line of sight calculation to compute the entire atmosphere. This implies the generalisation of a constant day-to-night side wind throughout the atmosphere without the possibility for stronger winds at the equator and weaker winds at the poles. As a consequence, the day-to-night side wind points towards the centre of the night side and does not go parallel to the equator. The same simplification applies to super-rotational winds (see Fig. 4). This is a direct consequence of our computational approach, which speeds up the execution of the code by orders of magnitude. However, we would like to highlight that the unrealistic approach of constant winds indicates that, in the case of a super-rotational wind, the actual wind speed within the jet is likely higher than our retrieved value.

However, more realistic models of atmospheres rarely show only one dominant wind pattern. We study, therefore, combined wind patterns, where a wind pattern is assumed in the lower atmosphere (either super-rotation of day-to-night side) and another one in the upper atmosphere (in this case vertical upward winds). We express the switch in orders of magnitude in altitude above the surface of the planet: $$P_{\text{switch}}=P_0\cdot {10}^h,$$(9)

where h is the variable for the retrieval (h < 0) and P₀ is the surface pressure deduced from NaX (see Sect. 2.4).

As a consequence of the Doppler shift of the profiles, a different wavelength grid is created in each bin. For the subsequent summation over all bins to produce the transmission spectrum, we have to interpolate the profile in each bin on a reference wavelength grid. We choose a reference grid at the same resolution as the HARPS wavelength calibration. The simplification of the wind profiles allows an extinction coefficient calculation in 1D, while retaining most of the 3D information of the wind patterns. This, in turn, increases computational efficiency to the degree where the forward model, including the broadening from winds, can be combined with a nested sampling retrieval.

Fig. 3 Illustration of the implementation of the day-to-night side wind pattern with a 1D calculation of the atmosphere. A polar view is shown on the left and an equatorial view on the right. The wind is indicated by the dotted blue arrows. Left: extinction coefficient is calculated in altitude along the z axis and then transposed in the x direction along the LOS (in dark green). The LOS is then iterated upward in z until the top of the atmosphere is reached and all values saved in a 2D grid (here visualised as a slice in light green). In each bin of the 2D grid the velocity and the broadened profile are calculated and stored. This is only necessary on one side of the atmosphere due to the symmetry of the problem. Right: calculated slice is then rotated to create the full atmosphere. The line of sight is shown as a dark green box in the main illustration, where the reader is in the position of the observer. Points indicate a flow towards the reader. In the day-to-night side case, the morning and evening limb winds pointtowards the observer and only one atmospheric slice has to be calculated. The slice is then rotated by 2π to create the full atmosphere. In this simplification of the atmosphere, the wind will not go parallel to the equator at all times, but point towards the centre of the night side. This reduces calculation time significantly, given that it reduces the problem from 3D to 2D, with the extinction coefficient only calculated in 1D.

Fig. 4 Illustration of the implementation of the super-rotational pattern with a 1D calculation of the atmosphere. See Fig. 3 for a detailed description of the figure and geometrical consequences on the computation. Left: since the wind is coming towards the observer for the morning terminator and away form the observer for the evening terminator, a slice of the atmosphere has to be calculated for each. Right: calculated slices are then rotated by π to get the full atmosphere (one half each). In the current model the wind does not change with latitude, creating a super-rotational stream that is constant tangentially over each atmospheric shell.

Fig. 5 Illustration of the implementation of the vertical pattern with a 1D calculation of the atmosphere. See Fig. 3 for a detailed description of the figure and geometrical consequences on the computation. Left: calculation of the LOS extinction coefficients is only necessary on one side of the atmosphere due to the symmetry of the problem. Right: calculated slice is then rotated by 2π to get the full atmosphere. The constant wind points outwards for all latitudes.

2.2 Nested sampling retrieval

Once the forward model creates a transmission spectrum for the sodium doublet, we use a nested sampling retrieval to find the best-fitting model from the selection.

We couple a nested sampling algorithm to the described forward model to explore the full parameter space of all variables (see Table 2). The main advantage of nested sampling lies in the direct calculation of the Bayesian evidence. This calculation, taking into account different parametrisations with changing number of parameters, allows for the comparison of these models equivalently (Skilling 2006; Feroz et al. 2009; Benneke & Seager 2013; Waldmann et al. 2015; Line et al. 2016), in a direct application of Occam’s Razor. The basics of these concepts are introduced in Trotta (2008); Skilling (2006). A full description of the approach taken here can be found in Lavie et al. (2017). For similar applications and examples, see Benneke & Seager (2013); Waldmann et al. (2015); Line et al. (2016).

To begin, N_live points are randomly selected from the parameter space under the constraint of our prior, where N parameters form an N-dimensional parameter space. For this first set of points, the likelihood values are calculated and in each step, the algorithm discards the point with the lowest likelihood and adds a new one until convergence is reached (see Skilling 2006; Lavie et al. 2017). Taking into account the independent Gaussian errors of our spectral observations, we calculate the likelihood as a Gaussian function (Eq. (3) in Benneke & Seager 2013). As in the implementation of Lavie et al. (2017), we use the open-source software PyMultiNest² (Buchner 2016), a wrapper for the open source MultiNest³ module (Feroz & Hobson 2008; Feroz et al. 2009, 2013).

Typical numbers of N_live range from N_live= 50–10 000 (Benneke & Seager 2013) up to 40 000 (accumulated, 400 runs with 100 live points) in Lavie et al. (2017), in this paper we use 10 000 (accumulated; 100 runs with 100 live points) each.

The Bayesian evidence for each of these models is calculated via Eq. (4) in Lavie et al. (2017) and allows us to compare the different models (here model M₀ and M₁ with Bayesian evidence Z₀ and Z₁) to each other via the logarithm of their Bayes factor (Trotta 2008): $$\ln {\mathcal{B}}_{01}=\log (Z_0/Z_1)=\log Z_0-\log Z_1.$$(10)

The Bayes factor is then judged via the so called “Jeffrey’s scale”, an empirical scale shown in Table 1 (Trotta 2008). It highlights when $\ln {\mathcal{B}}_{01}$ presents weak, moderate or strong evidence in favour of M₀ over M₁. The best-fit model allows us to rule out certain atmospheric scenarios, while taking into account the complexity of the different models and thus make quantitative statements about winds and temperature profiles in specific exoplanet atmospheres.

2.3 Combination in MERC

In MERC, we combine the forward model with the multi nested sampling module by using the priors of the nested sampling to create the temperature pressure (TP) profile and abundances needed as input for the forward model. The code then takes the result of the forward model, the model transmission spectrum, and hands it back to the nested sampling where the Bayesian evidence of the model is computed.

For the Bayesian evidence, MERC analyses the difference between the value of the model and the value of the data in each wavelength bin, meaning that each bin carries the same weight for all parameters. This comparison is performed on the unbinned HARPS data to preserve accuracy. Yet the temperature and wind strength have a strong impact only on the line core and the visible line wings, a small subset of the total number of bins in the sample, whereas the pressure at the base of the atmosphere and the constant sodium abundance are driven by the data continuum with a negligible impact from the core. To address this issue, we introduced a two-tier approach to the standard nested sampling retrieval, where we first retrieve the parameters dominated by the location of the continuum and then retrieve the additional parameters dominated by the line core.

In the first step, the data points of the sodium doublet core are blinded by leaving out the wavelength range twice the full width half maximum (FWHM), which includes the line cores. In the second step only this area containing the doublet is used for the retrieval. The two areas are highlighted in Fig. 6, where blue shows the line core retrieval range and green the continuum retrieval range. From hereon we call these two steps: continuum retrieval and line retrieval. The retrieval windows are chosen on an empirical basis so as to include all data points from the line core until the inclusion of additional points does not significantly change the retrieval results. We benchmark this approach by running the retrieval on line retrieval windows of three and four FWHM of the lines, with no changes to the result, but a significant increase in runtime.

In the flowchart of MERC (Fig. 1), the split of the data into continuum and line core is highlighted with the first, green cycle of the algorithm, the continuum retrieval, followed by the blue cycle, the line retrieval. The posterior of the continuum retrieval is used to set the prior for the line retrieval and only the line retrieval is affected by winds. The Bayesian evidence of the line retrieval is then used to rank the models and the posterior of the line retrieval gives the best-fit parameters and their uncertainties. The uncertainties are calculated as 1σ of the distribution. Given that not all posteriors follow a normal distribution, the two limits are calculated via the 0.16 and 0.84 quantile to include 68% of the distribution. The best-fit is calculated as the bin in the posterior distribution that had overall the highest Bayesian evidence. All corner plots for the posteriors with their indicated best-fits can be found in the Appendix B.

MERC is particularly well-suited for the wavelength range of the sodium doublet given its intrinsic large absorption depth, thus, capable of probing layers up to the thermosphere of the exoplanet. The different models employed to describe the atmosphere up to the thermosphere are listed in Table 2, where all retrieved parameters are shown. The continuum parameters are solely used to restrict the prior of the line retrieval and are still retrieved on the line.

Fig. 6 Spectral ratio of the hot Jupiter HD 189733b for the sodium doublet. See in grey the data obtained with the HARPS spectrograph, in black the binned data for visibility (for a detailed description of the data see Sect. 4.1). The dark blue areas highlight the retrieval zones at line centre and the light green areas the retrieval on the continuum.

2.4 Impact of the abundance pressure degeneracy

Based on the transmission spectrum geometry of Ehrenreich et al. (2006) and the theoretical work of Fortney (2005), Lecavelier Des Etangs et al. (2008) derived the effective altitude depending on wavelength as a function of extinction coefficient and partial pressure of the main absorber. Taking into account the variation of the cross section as a function of wavelength, a relation we also make use of in MERC to create a model transmission spectrum from atmospheric input parameters, they were able to relate the abundance of the absorber species to the temperature grid in altitude. However, this relation depends on the effective planetary radius at z = 0, R₀ and, subsequently, on the pressure at z = 0, P₀. Additionally, Heng & Kitzmann (2017) inferred a three-way degeneracy between P₀, R₀ and the element abundances χ, all of which change the level of the continuum in the transmission spectrum.

In a more recent work, Fisher & Heng (2018) show that this three way degeneracy is partially due to model assumptionsin semi-analytical approaches. Welbanks & Madhusudhan (2019) show that the degeneracy with the element abundances can be broken in retrieval analysis. This would require a more rigorous treatment of parameters influencing the continuum of the transmission spectra than presented here, most importantly the inclusion of H₂ −H₂ collision induced absorption (CIA), but also cloud decks, hazes, and additional absorber species.

We follow Lecavelier Des Etangs et al. (2008) in their description on the relation between the pressure scale and the height scale, assuming hydrodynamical equilibrium, and we set the reference radius at P₀ to the white light radius. We then retrieve P₀ to break the degeneracy between P₀ and R₀. Welbanks & Madhusudhan (2019) verified this approach, showing that R₀ can be set to any measured value because the retrieval adjusts P₀ accordingly.

Taking into account the simplicity of our forward model and the new findings by Welbanks & Madhusudhan (2019), we combine the two degenerate parameters P₀ and χ to set the continuum: $${\text{X}}_{\text{NaP}}=\frac{\chi }{{\chi }_{\odot }}\cdot \frac{P_0}{10 \text{bar}}={10}^{-5+\text{NaX}}.$$(11)

The right-hand side in Eq. (11) shows how the degenerate parameter is retrieved in MERC, where only the exponent that indicates the change from default values is a parameter. If NaX is 0, X_NaP is set to the default values. The abundance is then fixed throughout the atmosphere to a constant value, which is a commonly adopted approximation (Lecavelier Des Etangs et al. 2008; Agúndez et al. 2014; Steinrueck et al. 2019). The shape of the continuum is modelled via the calculation of the Rayleigh spectrum for H₂ (Lecavelier Des Etangs et al. 2008).

We consider the simplistic approach presented, which follows Heng et al. (2015), to treat the continuum sufficient for the purpose of this work, where the focus lies on the line shape and its possible origin in atmospheric wind patterns. The line shape is not significantly impacted by any of the degenerate parameters of this section. However, in follow up studies we plan to implement a more realistic approach in describing the element abundances, and the relation between pressure and height in atmospheric retrievals, following Welbanks & Madhusudhan (2019) to show their impact on line broadening.

3 Test on simulated data

We tested MERC (see Sect. 2) on simulated transmission spectroscopy data. We constructed two different simulated transmission spectra for HD 189733b with the HARPS resolution and the mean divergence and error of a HARPS transmission spectrum. We use the mean spread of the data in the continuum to introduce noise to the simulated data. MERC then creates one transmission spectra for HD 189733b with a purely isothermal atmosphere at 3600 K and no winds, and a second transmission spectra with an isothermal temperature profile at 3600 K and an added vertical wind throughout the atmosphere with a wind speed of 30 km s⁻¹. In both datasets, NaX was setto −1.3 to set the continuum. Both simulated datasets were created for arbitrary values of the parameters, one without and one with wind to highlight the ability of MERC to distinguish between a base model and a model with additional parameters. We then ran MERC on both simulated data sets, retrieving on each a isothermal temperature profile with no winds and an isothermal profile with an added vertical wind. We compared them to see if the code favoured the correct model for the simulated data sets.

3.1 Isothermal atmosphere

An overview of the parameters used to create the simulated data without winds and the two models is shown in Table 3. The isothermal model is a better fit than the model with vertical wind, but the difference in Bayesian evidence is inconclusive. Both models provide adequate fits to the simulated data since the model with added vertical winds centres the added wind around 0, making them equal (see Fig. A.2). The more complex model with the added wind is then punished for the additional retrieval parameter and has therefore a lower Bayesian evidence. Both parameters T and NaX are retrieved correctly as their original values.

3.2 Isothermal atmosphere with vertical winds

In the second test, we created a simulated dataset from a model spectrum that includes vertical upward winds at 30 km s⁻¹. We then employed MERC to model the dataset with an isothermal model and an isothermal model with additional winds, the same model used to create the data. The parameters with the posteriors from both retrievals are presented in Table 4. It shows that the isothermal model by itself is unable to sufficiently fit the provided simulated data with a Bayesian evidence 5.5 lower than the evidence for the isothermal model with added winds. This Bayesian evidence is strong evidence to support the second model (see Table 1). The second model with the added wind is additionally also able to retrieve all parameters to the correct values (best-fit compared to sim data in Table 4).

MERC is able to correctly model both simulated datasets and rank different models according to their ability to model the data. The posteriors of the retrieved parameters and the best-fit models together with the simulated data are shown in Figs. A.1–A.3.

4 Benchmarking with HD 189733b

4.1 Data reduction

We applied MERC to HD 189733b, one of the most well studied exoplanets to date. We used data from the HARPS echelle spectrograph mounted at ESO’s 3.6 m telescope in la Silla, Chile. The calculation of the transmission spectrum and the telluric correction had been undertaken by Wyttenbach et al. (2015). However, they did not correct the data for the Rossiter–McLaughlin (RM) effect. The RM effect induces distortions in the transmission spectrum due to the rotation of the star (e.g. see Rossiter 1924; McLaughlin 1924; Cegla et al. 2016; Triaud 2018).

We removed the stellar RM signature by dividing with the sum of all out-of-transit spectra (master out) after shifting it to the local velocities of the star behind the planet. The local stellar velocity was derived from Cegla et al. (2016). We compared the result to the RM corrected sodium transmission spectrum for HD 189733b from Casasayas-Barris et al. (2017). The two datasets were overplotted (see Fig. 7) and the resulting points are within 1σ of the planet transmission spectrum of the other dataset, showing that our RM correction is sufficiently accurate.

4.2 Retrieval results

We used MERC on multiple models: an isothermal temperature profile, a temperature gradient, an isothermal profile combined with different wind patterns (a day-to-night side wind, a super-rotational wind, a vertical upward wind and a combination of lower atmosphere day-to-night or super-rotational wind with vertical upward wind in the upper atmosphere).

The Bayesian evidence from all the models can be found in Table 5. All posterior distributions are in Appendix B, with the best-fit marked in blue and the mean values of the posterior printed on top of each column.

The isothermal model, as the simplest model, is used as the base model for comparison and has a Bayesian evidence of $\left|\ln \mathcal{Z}\right|=726.92\pm 0.14$. A comparison of all different models to each other is shown in Fig. 9, where the upper row of models is compared to all other models. Red stands for no evidence for a better fit, yellow for moderate evidence and green for strong evidence. The comparison scheme can be found in Table 1.

Fig. 7 Transmission spectrum of the two sodium doublet lines binned by ×5. Both spectra were generated from HARPS data of HD 189733b. In black is the spectrum from Wyttenbach et al. (2015) with the additional correction of the RM effect. In dark blue is the RM corrected spectrum from Casasayas-Barris et al. (2017) which we use as a reference with permission. The offset between the two high-resolution transmission spectra is consistently below 1σ and the smallest in the line wings, which is the most important range for wind retrieval. Therefore, we consider our RM correction to be sufficient for our purposes.

4.2.1 Isothermal profile versus temperature gradient

A comparison of the isothermal and temperature gradient models shows that there is no evidence for a better fit with the more complex temperature gradient. The posterior distributions for both models are in the appendix, Figs. B.1 and B.3. The mean of the isothermal temperature is $T(\text{mean})={3412}_{-432}^{+347}$ K. The posterior of the temperature gradient shows a wide spread for the base temperature at the surface of the planet $T_{\text{base}}(\text{mean})={2191}_{-858}^{+606}$ K and a best-fit close to the edge of the prior (T_base = 3182 K). The temperature at the top of the atmosphere shows a similarly wide spread. The comparatively poor convergence for this model stems from the degeneracy between the continuum parameter X_NaP, the base temperature T_base, and, subsequently, the top temperature T_top. To narrow down the possible solutions for the temperature gradient and compare the findings with Wyttenbach et al. (2015) and Pino et al. (2018), we set X_NaP to 0, corresponding to solar sodium abundance and a base pressure P₀ = 10 bar. Removing the degeneracy leads to a better convergence for the temperature gradient (see Fig. B.4). Plotting our retrieved TP profiles (isothermal and temperature gradient) shows broad agreement with the work of Wyttenbach et al. (2015) and Huang et al. (2017) (see Fig. 8). Work by Huang et al. (2017) has shown that an incremental fit of isothermal profiles as in Wyttenbach et al. (2015) underestimates the temperature in the higher layers of the atmosphere for atmospheres with rapidly increasing temperatures. Our work, which uses a temperature gradient, is in broad agreement with this conclusion, as our gradient is consistently on the higher temperature end of Wyttenbach et al. (2015) (see Fig. 8) and encloses results from Huang et al. (2017) at the higher range of retrieved values.

Based on the Bayesian evidence, it is clear that there is no difference for this dataset of HD 189733b whether a isothermal profile or a temperature gradient is used, therefore the wind patterns are applied with an isothermal temperature profile as the base.

4.2.2 Isothermal profiles with added wind pattern

When applying the additional basic wind profiles of a constant day-to-night side wind, a constant super-rotational wind, or a constant vertical upward wind, Fig. 9 shows that all three are better fits than the profiles without anyadditional winds, with moderate evidence for the day-to-night and super-rotational wind and strong evidence for the vertical wind pattern. All three wind patterns converge to temperatures and NaX parameters within 1σ of each other and all three wind patterns converge clearly towards a value for the velocity, with the posterior as a Gaussian distribution. In the best-fit case of vertical upward winds, this means a mean isothermal temperature of $T(\text{mean})={3426}_{-374}^{+316}$ K, NaX fitting the continuum at $\text{NaX}(\text{mean})=-{1.44}_{-0.33}^{+0.45}$ and a vertical wind of $|v_{\text{ver}}(\text{mean})|={40.55}_{-3.59}^{+3.82} \text{km}\hspace{0.17em}{\text{s}}^{-1}$. The vertical velocity is surprisingly large (see Sect. 4.3) and close to the escape velocity of HD 189733b of |v_escape| = 42.1 km s⁻¹, which suggests that more complex wind patterns combined with a vertical upward velocity in the upper layers of the atmosphere could provide an even better fit to the data. Based on this result, we implemented combinations of wind patterns. Given the current literature favouring rotational wind patterns for high pressures (up to 10⁻⁶ bar) over verticalwind patterns, we split the atmosphere into two layers. A lower part and an upper part with different wind patterns each, where the split retrieved according to Eq. (9). When splitting the atmosphere into the lower and upper atmosphere, the Bayesian evidence suggests a significantly better fit than any of the other, more simplemodels, except the model of the vertical wind pattern throughout the atmosphere (see green cells in Fig. 9). However, the posterior for the lower atmosphere day-to-night side model (see Fig. B.12) shows that the retrieval converges towards both the solution with only a vertical velocity in the atmosphere (h towards 0, the planet surface) and the split atmosphere solution (see best fit in blue in Fig. B.12), which subsequently generates a tail in the day-to-night side velocity distribution towards higher velocities. The best-fit is nonetheless found in the local likelihood minimum for the lower day-to-night side velocity pattern with h = −6. Assuming a super-rotational pattern in the lower atmosphere shows the same behaviour in the posterior distribution (see Fig. B.14) with the best-fit still for a split into lower and higher atmosphere at h = −5 and a wide spread of possible wind velocities in the lower atmosphere $\left|v_{\text{srot}}\right|={8.89}_{-5.85}^{+7.03} \text{km}\hspace{0.17em}{\text{s}}^{-1}$. Since neither lower wind pattern conclusively converges towards a wind speed and atmospheric structure, it is likely that we cannot discriminate between the wind patterns in the lower atmosphere from the current quality of the data. A retrieval of no wind pattern in the lower atmosphere and a vertical outbound wind pattern in the upper atmosphere hardens this suspicion (see Fig. B.16). The retrieval converges to a split of the atmosphere at h = −5.5, with a tail for the lower atmosphere and similar values for the other retrieved parameters when compared to the same retrieval with the vertical wind pattern throughout the atmosphere (see Fig. B.6). The vertical wind velocity remains high $|v_{\text{ver}}(\text{mean})|={40.08}_{-3.40}^{+3.72} \text{km}\hspace{0.17em}{\text{s}}^{-1}$. Additionally, the Bayesian evidence for the three retrievals with vertical winds in the upper atmosphere and different wind patterns in the lower atmosphere are nearly the same, with a small preference for super rotation in the lower atmosphere. This means we cannot make any statements about the possible wind patterns’ strength or shape in the lower atmosphere from the here presented dataset.

A plot of the data with the isothermal base model and the two best-fit, and therefore preferred, models (no lower winds with higher vertical outbound wind and super-rotation in the lower atmosphere with vertical winds in the upper atmosphere) is shown in Fig. 10.

Fig. 8 Temperature pressure profiles for HD 189733b. The points with error bars in black are from Wyttenbach et al. (2015), the red line corresponds to the profile from Huang et al. (2017) and the green and blue lines are the isothermal (see Figs. B.2 and B.1) and temperature gradient models (see Figs. B.5 and B.4) from this work. The shaded area corresponds to the 1σ area of the posterior distribution. The lowest data point in height comes from the continuum retrieval, the rest of the curve from the line cores.

Fig. 9 Difference in Bayesian evidence from Table 5. The top horizontal models are compared to each other, where dark red means no evidence, yellow means moderate evidence, and light green means strong evidence. All wind patterns that include vertical wind are strongly favoured over all models, the horizontal wind patterns show some improvement to no winds, and we cannot distinguish between an isothermal profile or a temperature gradient.

Fig. 10 Spectral ratio of the hot Jupiter HD 189733b for the sodium D2 and D1 line. In grey, the data obtained with the HARPS spectrograph at the ESO 3.6 m telescope; in black, the binned data for visibility. The red line shows the isothermal base model, compatible with findings from Wyttenbach et al. (2015), and the black line the best-fit with super-rotation in the lower atmosphere and vertical winds in the upper atmosphere. The green line shows the equally good fit of the model with no wind for pressures higher than h = −6 and solely vertical winds farther up in the atmosphere. The “bump” for both the black and green fits stems from the parts of the atmosphere where the LOS is orthogonal to the wind direction and the profile is not Doppler shifted. The parameters used to generate the fits are indicated as the blue lines (best-fit) in Figs. B.1, B.14. and B.16, respectively.

4.3 Implications

We explored an isothermal profile and a temperature gradient to fit the data, with different wind patterns added to the base model, as described in Sect. 2.1. We found no evidence that a temperature gradient provides a better fit than the isothermal base model with the current data (see Table 5 and Sect. 4.2) and we used the isothermal profile as the base for all models with added winds. The temperature gradient for the best-fit model is ${0.10}_{-0.22}^{+0.23}$ K km⁻¹, which is in agreement with the value found in Wyttenbach et al. (2015) and predicted by Heng et al. (2015), but does not allow for a validation of these results due to the high error bars. The retrieved value for the isothermal profile is in agreement with the top of the atmosphere results from Wyttenbach et al. (2015) as expected given that we retrieve on the line centre, the same part used in Wyttenbach et al. (2015) for the upper atmosphere. Both the isothermal and the temperature gradient are shown in Fig. 8, overplotted with their 1σ envelopes, together with the values from Wyttenbach et al. (2015) in black and the best-fit from Huang et al. (2017) in red. The temperatureat zero altitude was taken from the continuum retrieval and are in agreement with results from Rayleigh scattering modelling in Lecavelier Des Etangs et al. (2008) for HD 189733b.

Adding a day-to-night side wind to the isothermal base profile generates a better fit (moderate evidence) just as a super-rotational wind profile throughout the atmosphere with a wind velocity of ${3.8}_{-1.55}^{+1.51} \text{km}\hspace{0.17em}{\text{s}}^{-1}$, compatible with results from Salz et al. (2018). The strongest evidence measured against the isothermal base model is found when employing a vertical wind profile from the surface of the atmosphere to the highest probed layer. However, the wind speed for the best-fit is $\left|v_{\text{ver}}\right|={40.50}_{-3.59}^{+3.82} \text{km}\hspace{0.17em}{\text{s}}^{-1}$, speeds unlikely to exist in the lower layers of the atmosphere and in contradiction with predictions from Salz et al. (2018); Flowers et al. (2019).

For a super-rotational wind in the lower atmosphere and a vertical wind in the upper atmosphere, the best-fit gives a lower wind speed of 3.0 km s⁻¹, which changes to a best-fit vertical outwards bound wind of 38.90 km s⁻¹ (within one sigma of the other retrieved vertical wind values), switching at an altitude of 10⁻⁶ bar (corresponding to h = −5 and NaX = −1). This is roughly the same relative altitude as the atmospheric layers studied by Salz et al. (2018); Flowers et al. (2019) and, therefore, compatible with their findings of a super-rotational wind of 0− 4.5 km s⁻¹ for altitudes lower than 10⁻⁵ bar (see Fig. 2 in Flowers et al. 2019). However, this comparison is only on an order of magnitude level, given the uncertainty of the here presented retrieval values for the lower atmosphere.

The retrieved vertical wind speed in the upper atmosphere is associated with a scenario that relies on multiple assumptions and simplifications, and it suggests that additional mechanisms determine the shape the sodium line in HD 189733b. Due to the geometrical simplifications used in our forward model, we did not implement the rotation of the tidally locked planet (v_rot = 2.7 km s⁻¹). Therefore, we present the wind speed only as an upper boundary of the wind speeds needed to explain the broadening.

A full physical modelling as to the origin of this wind is beyond the scope of this paper, yet we would like to provide a possible scenario on the order of magnitude level leading to comparable wind speeds in planets similar to HD 189733b. We propose a super-rotational jet in the lower part of the atmosphere with wind speeds of the order of magnitude of a few km s⁻¹ based on Salz et al. (2018). While some work has been done on vertical shear in the super-rotational jet of HD 189733b (Brogi et al. 2016), it is not sufficient to introduce an acceleration (a) capable of propelling atoms to the velocities proposed in this paper. However, if ions at the speeds found in the super-rotational jet are subject to a magnetic field in the equatorial region, the Lorenz force creates an upwards movement of the ions. Literature on the transition zones between different atmospheric layers in Jupiter and Jupiter-like worlds is sparse, but it can be assumed that these zones are slim compared to the overall structures in the atmosphere Robinson & Catling (2014). The ionisation fraction Na⁺∕Na is computed with the Saha equation, where we use the number density of free electrons computed following a procedure from Mollière et al. (2017). We found that only a small fraction remains as neutral atoms. Assuming the Lorenz force is the potential source of these large velocities, we calculated the necessary magnetic field to propel a sodium ion to the observed speeds at the base of the upper atmosphere via: $$m\cdot a_{\text{ver}}=q\cdot v_{\text{rot}}\times B,$$(12)

where a_ver is the necessary acceleration to propel particles to the retrieved vertical wind speed over the slim transitional region between atmospheric layers (order of magnitude ~1 km), and v_rot is the velocity of particles inside the super-rotational jet. This leads to an order of magnitude for magnetic field strength (~ 50 G) orders of magnitude higher than the predicted fiels strength for planets with rotation periods of ~ 2−4 days Cauley et al. (2019). However, the necessary magnetic field strength is comparable at the order of magnitude level with the magnetic field strength measured for HD 189733b (Cauley et al. 2019). The ions recombine rapidly to neutral species, which are then at the accelerated speed of the former ions. Furthermore, the remaining ions could also drag the neutral species upwards. In this context, it is worth mentioning that the very top of the sodium line probes high altitudes, where collisions between the sodium atoms become less and less frequent, which favours a constant speed of the particles compared to a deceleration.

Additional effects that could have impacts on the sodium line broadening and change the results from our analysis are ionisation due to the close proximity to the host star and non local thermal equilibrium (non-LTE) effects. In the case of ionisation, extra ionisation sources potentially deplete the upper layers of the atmosphere of neutral sodium. The dearth of sodium creates a more shallow line centre, while the line wings keep their shape. This can be mistaken for a full sodium line with broad wings instead, creating a degeneracy between broadened line wings or strong ionisation as possible physical explanations of the line shape (Lothringer & Barman 2019). Non-LTE effects might also influence the line shape, however these effects cannot be distinguished with the current data quality and might have a smaller contribution (Fisher & Heng 2019).

Lecavelier Des Etangs et al. (2010), Lecavelier Des Etangs et al. (2012), Bourrier & Lecavelier des Etangs (2013) studied HD 189733b in Lyman α and indicated an expanding exosphere. Vidal-Madjar et al. (2003); Lammer et al. (2003) has also shown, that the planet is evaporating and that it is thought to arise from an expanding thermosphere, reaching high velocities in the upper parts. Considering the heights we probe (see Fig. 8), it is possible we reach these layers and see sodium atoms propelled to these velocities close to the escape velocity.

5 Conclusions

We introduced the MERC code - a new tool to combine 1D atmospheric forward models with pseudo 3D wind patterns together with a rigorous nested sampling algorithm. The pseudo 3D wind is created by calculating the Doppler broadening for each cell in one or more atmospheric slices and rotating them to create the full atmosphere. From this model atmosphere, a model transmission spectrum is calculated and subsequently compared to the real data via the multi-nested sampling module. To achieve a correct fitting for wind strengths and patterns, a separation of the continuum of the data from the line cores was introduced. We showed the self consistency of MERC by running it on simulated data, highlighting its ability to not only correctly retrieve the parameter values but also to distinguish between models successfully. We then applied our technique to the well-studied exoplanet HD 189733b on the sodium doublet first published in Wyttenbach et al. (2015) and corrected for the RM effect. We tested an isothermal temperature profile and a temperature gradient, highlighting that both solutions provide satisfactory fits. We then explored three additional wind profiles: a day-to-night side wind, a super-rotational wind and a vertical upward wind. Although the vertical wind throughout the atmosphere was preferred, we expanded this result by splitting the atmosphere in a lower and upper atmospheric layer. Testing the established wind patterns in the lower atmosphere with an additional vertical wind in the upper atmosphere showed that the current data quality does not allow us to distinguish wind patterns in the lower atmosphere. However, the solution with no assumptions as to the wind profile in the lower atmosphere and a vertical wind in the upper atmosphere was clearly retrieved as the best solution. Our best solution has a mean temperature of $T(\text{mean})={3412}_{-432}^{+347}$ K and a vertical upward wind of $|v_{\text{ver}}(\text{mean})|={40.08}_{-3.40}^{+3.72} \text{km}\hspace{0.17em}{\text{s}}^{-1}$ in the higher atmosphere up to the thermosphere.

The broadened sodium line in HD 189733b can, thus, be explained by strong winds in the upper layers of the atmospheres up to the thermosphere, possibly propelled by strong magnetic fields.

With the next generation of spectrographs (e.g. ESPRESSO), we hope to distinguish between different models in the lower atmosphere and classify models by their Bayesian evidence at various altitudes above the planets surface, as well as implement other possible broadening sources in MERC.

Acknowledgements

This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (project FOUR ACES; grant agreement No. 724427). This work has been carried out within the frame of the National Centre for Competence in Research “PlanetS” supported by the Swiss National Science Foundation (SNSF). L.P. acknowledges that the research leading to these results has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement no. 679633; Exo-Atmos). We thank N. Casasayas-Barris for their assistance with the RM benchmarking and L. Dos Santos for their helpful comments. We would also like to show our gratitude to Kevin Heng for discussions and insights during the course of this research, that greatly improved the understanding of theoretical issues in atmospheric modelling. We thank P. Mollière for the use of his code to calculate the ionisation fraction and the anonymous referee for their helpful comments and time.

Appendix A Posterior distributions simulated data

Fig. A.1 Posterior distributions on the simulated data. Column a: retrieved on the isothermal simulated data with no additional winds. Top panel: retrieval with an isothermal forward model, the bottom panel with an added vertical wind. Column b: retrieved on the simulated data from an isothermal temperature profile with added vertical upward winds. The top panel was retrieved with an isothermal forward model and the bottom panel with the added wind. All results were within expectations.

Fig. A.2 Simulated data from an isothermal atmosphere in grey overplotted with the best-fit model. In red the isothermal model with no winds is shown and in blue the isothermal model with added vertical wind. Both models provide an identical fit.

Fig. A.3 Simulated data from an isothermal atmosphere with vertical upward winds added in grey overplotted with the best-fit model. In red the isothermal model with no winds is shown and in blue the isothermal model with added vertical wind.

Appendix B Posterior distributions for HD 189733b

Fig. B.1 Posterior distribution of isothermal line retrieval. This was used as the base model to compare all other scenarios to.

Fig. B.2 Transmission spectrum of HD 189733b for both lines of the sodium doublet in grey, binned in black, with the best-fit retrieved by MERC in red. The best-fit was generated for the isothermal base model with no added winds and its best-fit parameters from the posterior distribution in Fig. B.1.

Fig. B.3 Posterior distribution of temperature gradient line retrieval. T base is at the surface of the planet and T top at the top of the probed atmosphere. Due to the degeneracy between the three parameters, the convergence for this model is poor.

Fig. B.4 Posterior distribution of temperature gradient line retrieval. T base is at the surface of the planet and T top at the top of the probed atmosphere. The continuum parameter XNaP is set to 0 to retrieve a converged temperature gradient and compare it to existing literature.

Fig. B.5 Transmission spectrum of HD 189733b for both lines of the sodium doublet in grey, binned in black, with the best-fit retrieved by MERC in red. The best-fit was generated for the temperature gradient model with no added winds and its best-fit parameters from the posterior distribution in Fig. B.4.

Fig. B.6 Posterior distribution of isothermal line retrieval with an added vertical wind constant throughout the atmosphere.

Fig. B.7 Transmission spectrum of HD 189733b for both lines of the sodium doublet in grey, binned in black, with the best-fit retrieved by MERC in red. The best-fit was generated for the temperature gradient model with added vertical winds and its best-fit parameters from the posterior distribution in Fig. B.6.

Fig. B.8 Posterior distribution of isothermal line retrieval with an added day-to-night side wind constant throughout the atmosphere.

Fig. B.9 Transmission spectrum of HD 189733b for both lines of the sodium doublet in grey, binned in black, with the best-fit retrieved by MERC in red. The best-fit was generated for the temperature gradient model with added day to night side winds and its best-fit parameters from the posterior distribution in Fig. B.8.

Fig. B.10 Posterior distribution of isothermal line retrieval with an added super-rotational wind constant throughout the atmosphere.

Fig. B.11 Transmission spectrum of HD 189733b for both lines of the sodium doublet in grey, binned in black, with the best-fit retrieved by MERC in red. The best-fit was generated for the temperature gradient model with added super-rotational winds and its best-fit parameters from the posterior distribution in Fig. B.10.

Fig. B.12 Posterior distribution of isothermal line retrieval with an added day-to-night side wind in the loweratmosphere and a vertical wind in the upper atmosphere. The parameter on the right h is then used to calculate the pressure where the velocity models are switched via Eq. (9). Most notably, the retrieval is able to produce a solution with a day-to-night side wind in the upper layers of the atmosphere, but due to the data quality the solution with solely a vertical wind (h closer to 0.) is degenerate (we would like to point out the two local maxima for h).

Fig. B.13 Transmission spectrum of HD 189733b for both lines of the sodium doublet in grey, binned in black, with the best-fit retrieved by MERC in red. The best-fit was generated for the temperature gradient model with added day to night side winds in the lower atmosphere and vertical winds in the upper atmosphere with its best-fit parameters from the posterior distribution in Fig. B.12.

Fig. B.14 Posterior distribution of isothermal line retrieval with an added super-rotational wind in the loweratmosphere and a vertical wind in the upper atmosphere. The parameter on the right h is then used to calculate the pressure where the velocity models are switched via Eq. (9). Just like the day-to-night side wind, the preferred solution of a super-rotational wind in the upper atmosphere (h = −5) is degenerate with a solely vertical wind (h converging towards 0).

Fig. B.15 Transmission spectrum of HD 189733b for both lines of the sodium doublet in grey, binned in black, with the best-fit retrieved by MERC in red. The best-fit was generated for the temperature gradient model with added super-rotational winds in the lower atmosphere and vertical winds in the upper atmosphere with its best-fit parameters from the posterior distribution in Fig. B.14.

Fig. B.16 Posterior distribution of isothermal line retrieval with an added vertical wind in the upper atmosphere. The parameter on the right h is then used to calculate the pressure where the vertical velocity is added via Eq. (9). The fraction of the atmosphere below that pressure has no winds. The distribution of h shows that the retrieval shows a clear preference for high vertical winds for pressures below h = −6, but is unable to discriminate between the wind patterns in the lower parts of the atmosphere.

Fig. B.17 Transmission spectrum of HD 189733b for both lines of the sodium doublet in grey, binned in black, with the best-fit retrieved by MERC in red. The best-fit was generated for the temperature gradient model with no added winds in the lower atmosphere and vertical winds in the upper atmosphere with its best-fit parameters from the posterior distribution in Fig. B.16.

References

All Tables

All Figures

Fig. 1 Flowchart showing the information flow from forward model to the nested sampling module and vice versa in MERC. Green arrows are executed during the continuum retrieval, blue arrows are only executed during the line retrieval. Starting with the retrieval windows in green from Fig. 6 and the parameters that can be retrieved onthe continuum (see Table 2), priors are set in the nested sampling module and handed to the forward model. One possible combination from the prior range is then used to calculate the extinction coefficient τ in a 1D grid in altitude above the planet surface. The planet surface is defined as the opaque part of the planet at the wavelength of the green region and is set at the optical white light radius. The 1D grid can be seen as a slice of the atmosphere, which is then used to create the full atmosphere (see Sect. 2.1) and calculate a model transmission spectrum. The model spectrum is handed back to the nested sampling module, where it is compared to the continuum of the data. The multi nested sampling executes this loop over the prior parameter space until the parameters converge to a value within the prior ranges. It then produces a posterior distribution of the prior and the Bayesian evidence of the model. The posterior of the continuum retrieval is then used to restrict the prior of the retrieval on the line retrieval window (shown in blue). The same process as for the continuum is started again, but now with additional wind broadening in each cell of the atmospheric slice. At the end, the Bayesian evidence of the line retrieval for each model is used to rank them. This flowchart was loosely inspired on a similar layout in Brogi & Line (2019).

In the text

Fig. 2 Schematics of the extinction coefficient calculation in the line of sight (LOS) exemplified for the outermost bin in x for a vertical wind profile. The atmospheric wind pattern is assumed to be a constant radial velocity originating from the planet’s surface. γ is the angle running from the polar position towards the observer, R the assumed radius of the planet, b the impact factor and x the sum overall bins in the line of sight x = ∑Δx.

In the text

Fig. 3 Illustration of the implementation of the day-to-night side wind pattern with a 1D calculation of the atmosphere. A polar view is shown on the left and an equatorial view on the right. The wind is indicated by the dotted blue arrows. Left: extinction coefficient is calculated in altitude along the z axis and then transposed in the x direction along the LOS (in dark green). The LOS is then iterated upward in z until the top of the atmosphere is reached and all values saved in a 2D grid (here visualised as a slice in light green). In each bin of the 2D grid the velocity and the broadened profile are calculated and stored. This is only necessary on one side of the atmosphere due to the symmetry of the problem. Right: calculated slice is then rotated to create the full atmosphere. The line of sight is shown as a dark green box in the main illustration, where the reader is in the position of the observer. Points indicate a flow towards the reader. In the day-to-night side case, the morning and evening limb winds pointtowards the observer and only one atmospheric slice has to be calculated. The slice is then rotated by 2π to create the full atmosphere. In this simplification of the atmosphere, the wind will not go parallel to the equator at all times, but point towards the centre of the night side. This reduces calculation time significantly, given that it reduces the problem from 3D to 2D, with the extinction coefficient only calculated in 1D.

In the text

Fig. 4 Illustration of the implementation of the super-rotational pattern with a 1D calculation of the atmosphere. See Fig. 3 for a detailed description of the figure and geometrical consequences on the computation. Left: since the wind is coming towards the observer for the morning terminator and away form the observer for the evening terminator, a slice of the atmosphere has to be calculated for each. Right: calculated slices are then rotated by π to get the full atmosphere (one half each). In the current model the wind does not change with latitude, creating a super-rotational stream that is constant tangentially over each atmospheric shell.

In the text

Fig. 5 Illustration of the implementation of the vertical pattern with a 1D calculation of the atmosphere. See Fig. 3 for a detailed description of the figure and geometrical consequences on the computation. Left: calculation of the LOS extinction coefficients is only necessary on one side of the atmosphere due to the symmetry of the problem. Right: calculated slice is then rotated by 2π to get the full atmosphere. The constant wind points outwards for all latitudes.

In the text

	Fig. 6 Spectral ratio of the hot Jupiter HD 189733b for the sodium doublet. See in grey the data obtained with the HARPS spectrograph, in black the binned data for visibility (for a detailed description of the data see Sect. 4.1). The dark blue areas highlight the retrieval zones at line centre and the light green areas the retrieval on the continuum.
In the text

Fig. 7 Transmission spectrum of the two sodium doublet lines binned by ×5. Both spectra were generated from HARPS data of HD 189733b. In black is the spectrum from Wyttenbach et al. (2015) with the additional correction of the RM effect. In dark blue is the RM corrected spectrum from Casasayas-Barris et al. (2017) which we use as a reference with permission. The offset between the two high-resolution transmission spectra is consistently below 1σ and the smallest in the line wings, which is the most important range for wind retrieval. Therefore, we consider our RM correction to be sufficient for our purposes.

In the text

Fig. 8 Temperature pressure profiles for HD 189733b. The points with error bars in black are from Wyttenbach et al. (2015), the red line corresponds to the profile from Huang et al. (2017) and the green and blue lines are the isothermal (see Figs. B.2 and B.1) and temperature gradient models (see Figs. B.5 and B.4) from this work. The shaded area corresponds to the 1σ area of the posterior distribution. The lowest data point in height comes from the continuum retrieval, the rest of the curve from the line cores.

In the text

Fig. 9 Difference in Bayesian evidence from Table 5. The top horizontal models are compared to each other, where dark red means no evidence, yellow means moderate evidence, and light green means strong evidence. All wind patterns that include vertical wind are strongly favoured over all models, the horizontal wind patterns show some improvement to no winds, and we cannot distinguish between an isothermal profile or a temperature gradient.

In the text

Fig. 10 Spectral ratio of the hot Jupiter HD 189733b for the sodium D2 and D1 line. In grey, the data obtained with the HARPS spectrograph at the ESO 3.6 m telescope; in black, the binned data for visibility. The red line shows the isothermal base model, compatible with findings from Wyttenbach et al. (2015), and the black line the best-fit with super-rotation in the lower atmosphere and vertical winds in the upper atmosphere. The green line shows the equally good fit of the model with no wind for pressures higher than h = −6 and solely vertical winds farther up in the atmosphere. The “bump” for both the black and green fits stems from the parts of the atmosphere where the LOS is orthogonal to the wind direction and the profile is not Doppler shifted. The parameters used to generate the fits are indicated as the blue lines (best-fit) in Figs. B.1, B.14. and B.16, respectively.

In the text

Fig. A.1 Posterior distributions on the simulated data. Column a: retrieved on the isothermal simulated data with no additional winds. Top panel: retrieval with an isothermal forward model, the bottom panel with an added vertical wind. Column b: retrieved on the simulated data from an isothermal temperature profile with added vertical upward winds. The top panel was retrieved with an isothermal forward model and the bottom panel with the added wind. All results were within expectations.

In the text

	Fig. A.2 Simulated data from an isothermal atmosphere in grey overplotted with the best-fit model. In red the isothermal model with no winds is shown and in blue the isothermal model with added vertical wind. Both models provide an identical fit.
In the text

	Fig. A.3 Simulated data from an isothermal atmosphere with vertical upward winds added in grey overplotted with the best-fit model. In red the isothermal model with no winds is shown and in blue the isothermal model with added vertical wind.
In the text

	Fig. B.1 Posterior distribution of isothermal line retrieval. This was used as the base model to compare all other scenarios to.
In the text

	Fig. B.2 Transmission spectrum of HD 189733b for both lines of the sodium doublet in grey, binned in black, with the best-fit retrieved by MERC in red. The best-fit was generated for the isothermal base model with no added winds and its best-fit parameters from the posterior distribution in Fig. B.1.
In the text

	Fig. B.3 Posterior distribution of temperature gradient line retrieval. T base is at the surface of the planet and T top at the top of the probed atmosphere. Due to the degeneracy between the three parameters, the convergence for this model is poor.
In the text

	Fig. B.4 Posterior distribution of temperature gradient line retrieval. T base is at the surface of the planet and T top at the top of the probed atmosphere. The continuum parameter XNaP is set to 0 to retrieve a converged temperature gradient and compare it to existing literature.
In the text

	Fig. B.5 Transmission spectrum of HD 189733b for both lines of the sodium doublet in grey, binned in black, with the best-fit retrieved by MERC in red. The best-fit was generated for the temperature gradient model with no added winds and its best-fit parameters from the posterior distribution in Fig. B.4.
In the text

	Fig. B.6 Posterior distribution of isothermal line retrieval with an added vertical wind constant throughout the atmosphere.
In the text

	Fig. B.7 Transmission spectrum of HD 189733b for both lines of the sodium doublet in grey, binned in black, with the best-fit retrieved by MERC in red. The best-fit was generated for the temperature gradient model with added vertical winds and its best-fit parameters from the posterior distribution in Fig. B.6.
In the text

	Fig. B.8 Posterior distribution of isothermal line retrieval with an added day-to-night side wind constant throughout the atmosphere.
In the text

	Fig. B.9 Transmission spectrum of HD 189733b for both lines of the sodium doublet in grey, binned in black, with the best-fit retrieved by MERC in red. The best-fit was generated for the temperature gradient model with added day to night side winds and its best-fit parameters from the posterior distribution in Fig. B.8.
In the text

	Fig. B.10 Posterior distribution of isothermal line retrieval with an added super-rotational wind constant throughout the atmosphere.
In the text

	Fig. B.11 Transmission spectrum of HD 189733b for both lines of the sodium doublet in grey, binned in black, with the best-fit retrieved by MERC in red. The best-fit was generated for the temperature gradient model with added super-rotational winds and its best-fit parameters from the posterior distribution in Fig. B.10.
In the text

Fig. B.12 Posterior distribution of isothermal line retrieval with an added day-to-night side wind in the loweratmosphere and a vertical wind in the upper atmosphere. The parameter on the right h is then used to calculate the pressure where the velocity models are switched via Eq. (9). Most notably, the retrieval is able to produce a solution with a day-to-night side wind in the upper layers of the atmosphere, but due to the data quality the solution with solely a vertical wind (h closer to 0.) is degenerate (we would like to point out the two local maxima for h).

In the text

	Fig. B.13 Transmission spectrum of HD 189733b for both lines of the sodium doublet in grey, binned in black, with the best-fit retrieved by MERC in red. The best-fit was generated for the temperature gradient model with added day to night side winds in the lower atmosphere and vertical winds in the upper atmosphere with its best-fit parameters from the posterior distribution in Fig. B.12.
In the text

Fig. B.14 Posterior distribution of isothermal line retrieval with an added super-rotational wind in the loweratmosphere and a vertical wind in the upper atmosphere. The parameter on the right h is then used to calculate the pressure where the velocity models are switched via Eq. (9). Just like the day-to-night side wind, the preferred solution of a super-rotational wind in the upper atmosphere (h = −5) is degenerate with a solely vertical wind (h converging towards 0).

In the text

	Fig. B.15 Transmission spectrum of HD 189733b for both lines of the sodium doublet in grey, binned in black, with the best-fit retrieved by MERC in red. The best-fit was generated for the temperature gradient model with added super-rotational winds in the lower atmosphere and vertical winds in the upper atmosphere with its best-fit parameters from the posterior distribution in Fig. B.14.
In the text

Fig. B.16 Posterior distribution of isothermal line retrieval with an added vertical wind in the upper atmosphere. The parameter on the right h is then used to calculate the pressure where the vertical velocity is added via Eq. (9). The fraction of the atmosphere below that pressure has no winds. The distribution of h shows that the retrieval shows a clear preference for high vertical winds for pressures below h = −6, but is unable to discriminate between the wind patterns in the lower parts of the atmosphere.

In the text

	Fig. B.17 Transmission spectrum of HD 189733b for both lines of the sodium doublet in grey, binned in black, with the best-fit retrieved by MERC in red. The best-fit was generated for the temperature gradient model with no added winds in the lower atmosphere and vertical winds in the upper atmosphere with its best-fit parameters from the posterior distribution in Fig. B.16.
In the text