Effect of improved atmospheric opacities in modelling sub-Neptunes-A case study for GJ 1214 b

Aims. We investigate the impact of updated atmospheric mean opacity input values on modelled transit radius and the distribution of interior layer mass fractions. Methods. We developed and applied a coupled interior–atmosphere model. Our straightforward semi-grey calculation of atmospheric temperature enables us to perform thousands of model realisations in a Monte Carlo approach to address potential degeneracies in interior and atmospheric mass fraction. Our main constraints are planetary mass and radius from which our model infers distributions of the internal structure of exoplanetary classes ranging from Super-Earth to Mini-Neptune. We varied the relative masses of gas, envelope, mantle, and core layers subject to constraints on the bulk density from observations, and investigated the effect of updating atmospheric mean opacities. Results. First, we validate our model output with observed temperature profiles for modern Neptune. We can reproduce the basic features in the middle atmosphere but not the temperature inversion in the upper layers, which is likely because our model lacks aerosol heating. Calculated interiors are generally consistent with modern Neptune. Second, we compare with the well-studied object GJ 1214 b and obtain results that are broadly consistent with previous findings; they suggest correlations between modelled gas, water, and core mass fractions, although these are generally weak. Updating the opacities leads to a change on the order of a few percent in the modelled transit radius. This is comparable in magnitude to the planned accuracy of the PLATO data for planetary radius, suggesting that the opacity update in the model is important to implement.


Introduction
The Kepler mission has suggested that sub-Neptunes are likely to be very common (Batalha 2014) with a wide range of masses and radii.Interior structures however are challenging to constrain since numerous, degenerate bulk compositions can fit the current data.Interior-atmosphere models that are both accurate and flexible are needed in order to explore these degeneracies and begin to quantify them.These degeneracies have now been investigated for over a decade; see for example Rogers & Seager (2010), Nettelmann et al. (2011), andValencia et al. (2013) (who all focused on GJ 1214 b).Recently, several new studies have been published that use such models to investigate this issue (e.g.Dorn & Heng 2018;Madhusudhan et al. 2020;Kite & Barnett 2020).We address this topic by applying an interior-atmosphere coupled model that compares reasonably well with more complex schemes but can be applied over a wider parameter range.Our goal is to investigate the effect of improved atmospheric parameters in order to extend our knowledge of the interior structure of sub-Neptunes.A&A 671, A65 (2023) With this isothermal profile, the radius of the atmosphere of their modelled planet grid is ≈0.1 R P .This approach neglects the effects of greenhouse heating for example, which increases the scale height and therefore the thickness of the atmosphere.In their relation, the top of the atmosphere is kept constant at an arbitrary pressure, whereas in reality the observed transit radius is a function of the atmospheric bulk properties.
Regarding atmosphere-interior models with semi-grey atmospheres, Dorn & Heng (2018) used the Guillot (2010) approximation to calculate temperature.Rogers & Seager (2010) used a broadly comparable model described in Hansen (2008).However, as their model input, these latter authors used a ratio of shortwave to longwave mean opacities (γ) based on earlier works (Freedman et al. 2008, hereafter F08) that is approximately one order of magnitude larger than the updated value in Freedman et al. (2014, hereafter F14).In general, γ values vary as a function of composition (e.g.steam, primary atmospheres) and pressure.Valencia et al. (2013), who modelled the bulk composition of sub-Neptunes, chose a mean opacity (κ) input value that yields 1000 K at 1 bar based on a line-by-line approach, and then employed an extrapolated analytical formula to vary the mean opacity with pressure, temperature, and metallicty.
Regarding complex climate-chemistry atmosphere-only models, numerous exoplanetary studies to date have applied such models to investigate for example potential habitability and atmospheric biosignatures (e.g.Segura et al. 2003Segura et al. , 2005;;Hu & Seager 2014;Meadows et al. 2018;Madhusudhan et al. 2020;Scheucher et al. 2020;Wunderlich et al. 2020).There exist degeneracies in the atmosphere-interior composition based on mass-radius measurements alone.We note that transmission spectra, for example from the James Webb Space Telescope (JWST), will likely help to address such degeneracies.
GJ 1214 b (Charbonneau et al. 2009) is one of the most wellinvestigated sub-Neptunes in terms of both planetary interior and atmosphere.Observations generally suggest a degeneracy between either a higher metallicity atmosphere or the presence of equilibrium condensate such as clouds or photochemical hazes in order to account for the relatively flat transmission spectra (Kreidberg et al. 2014;Rackham et al. 2017;Crossfield & Kreidberg 2017;Gao & Benneke 2018;Kawashima et al. 2019;Lavvas et al. 2019).The non-detection of helium (Kasper et al. 2020) suggests the possible loss of the primary atmosphere.Various works have chosen to focus on the interior, assuming a relatively simple gas layer (Nettelmann et al. 2011) or modelling the atmosphere using a climate-only model (Charnay et al. 2015;Lavvas et al. 2019) and neglecting the interior altogether.Rogers & Seager (2010) and Valencia et al. (2013) applied atmosphereinterior models to GJ 1214 b but did not explore the effect of updating input parameters such as opacity.
A main motivation of the present work is to investigate the influence of input data updates in the mean opacities from F14 updated from F08 on the simulated atmosphere-interior over a wide parameter range for GJ 1214 b.These opacities are used as input in our modelled gas layer.We compare with output from other models in the literature that use the original opacity data.As GJ 1214 b is likely to have a substantial gas layer due to its low bulk density, employing updated opacities is especially important when accurately modelling the planetary structure.
This paper is structured as follows.Section 2 provides the model description, Sects.3 and 4 present the validation and case study, Sect. 5 presents results for the effect of the opacity update, and Sect.6 presents a summary of our conclusions.Supplementary figures and tables can be found in Appendices A and B, respectively.

Model description
Our interior-atmosphere model TATOOINE (Tool for ATmospheres, Outgassing, and Optimal INteriors of Exoplanets) splits the planet into multiple layers of distinct composition: an iron core, a silicate mantle, a water-ice shell, and a hydrogen-helium gas layer.These different layers are described individually below.The interior model was discussed in detail in Baumeister et al. (2020) and here we provide just a brief description.The gas layer is a new feature and is described and validated in this paper for the first time.The different layers are described individually below followed by a brief description of the interior model.
Our main assumptions are as follows.We calculate a stratified four-layer structure of distinct compositions in hydrostatic equilibrium, with Earth-like values for interior parameters such as the magnesium number (Mg#) and the Fe/Si ratio, with only one layer in the iron core and silicate mantle, and a water layer composed of pure H 2 O (as opposed to a mixture of for example water, methane, and ammonia).Our gas layer is grain-free (no clouds or hazes), assuming a straightforward climate (see below) with no interactive chemistry, and no convection in the semi-grey portion of the atmosphere.We also do not make use of any a priori information -from for example formation models -in choosing our various model input values (Fe/Si, Mg#, atmospheric composition, or mean molecular weight, etc.).
As input, the model takes the planetary mass, layer composition, equation of state (EoS), and the mass fraction of each layer.From there, the planet is split into 2000 spherical shells.Mass, pressure, and temperature are integrated from the top downward to the planet's centre.The model is iterated until the excess mass at the planet centre is less than the preset tolerance (1 × 10 −6 ) between consecutive iterations.In order to fully characterise the inherently degenerate solution space, we perform a parameter study over a range of layer mass fractions.We adopt a Monte Carlo (MC) approach as described in the modelling method Sect.2.3.
We analyse the results of the parameter study by selecting MC model outcomes where the calculated transit radius lies within the error range for the observed radius of the planet R obs ± σ R .We then plot the posterior distribution of the mass fraction in each layer and calculate the median value where the errors are calculated at the 16% and 84% quantiles (see Figs. 1,4,A.1 and Sect. 2.4).These median values are then used as input for the core, water, and gas mass fraction in order to construct individual profiles for each model realisation (e.g.

Gas layer
The modelled gas layer has two regions: an upper layer consisting of the irradiated semi-grey atmosphere whose radiative budget is dominated by instellation, and a lower, high-pressure solar composition envelope.The composition of the atmosphere is determined by equilibrium chemistry which is calculated in the present study by the NASA Chemical Equilibrium with Applications (CEA) model (McBride et al. 1992).The CEA model uses a descent Newton Raphson method to minimise Gibbs free energy and then calculate the composition at thermodynamic equilibrium (see Table B.1).Initial element abundances are calculated based on user-input metallicity.Data are calculated for 50 reference elements including species condensation.Further details of the method and references for the thermodynamical input A65, page 2 of 15 tables are provided by McBride et al. (1992).We note that the mean opacities calculated by the Freedman studies mentioned above assume condensation of only primary species, unlike the CEA model which includes more complex condensation of secondary species.However.for the current work, we assume that the effects of this caveat are small for the rather modest range of metallicities considered.The energy budget of the radiative regime is determined on the one hand by photons emitted from the planetary interior below and on the other hand by heat release into the atmospheric region above.Within the envelope, there are two regimes: radiative and convective, the boundary of which is defined by the Schwarzschild criterion.

Irradiated atmosphere
In the irradiated atmosphere, we apply the radiative scheme of Guillot (2010) to calculate the temperature as a function of pressure for hydrogen-helium-dominated atmospheres using their Eq.( 49).As input to the Guillot model, we provide the equilibrium temperature T eq (using an albedo of zero if not otherwise known) and the intrinsic temperature T int , and γ = κ vis /κ th , where κ vis and κ th are the visible and thermal (Rosseland) mean opacities from F14, respectively.For each metallicity, the Rosseland mean opacities are provided as a table of look-up values along a grid of for example P local , T local , T eff .The Rosseland mean opacity is defined as: where the wavelength-dependent opacity (κ λ ) is calculated for each point in the pressure-temperature grid, and the Planck function (B λ ) is weighted using either the local temperature for κ th , or stellar effective temperature for κ vis .
In this work, we employ Rosseland rather than Planck mean opacities in order that we may compare with earlier modelling studies such as that of Valencia et al. (2013).The Rosseland approach assigns a higher weighting at wavelengths contributing low opacities.It is therefore particularly suited to deep atmospheric regions where photon transport is approximately diffusional.However, the Guillot method used in our work employs a single opacity that is constant throughout the atmosphere.To our knowledge, there is no theory that suggests how to apply the Planck or the Rosseland mean over the entire profile in such a case.However, we note that, in reality, changing chemical composition and pressure over the atmospheric profile could lead to variations in the value of the Rosseland opacity.These caveats are a point of investigation for future work.
The data available for a given exoplanet is generally minimal, and so as a starting point we used T eq for T local (with the associated T eff for κ vis ) because T eq can be derived observationally.We note that an aim of our work is to determine the effect of updating the original opacity data without the challenge of interpolating in p-T , which could introduce biases.For our model calculations, we therefore chose the closest available values in the input (pressure-temperature) opacity grid to our model scenarios (see below).
In the case of GJ 1214 b, we assume T eq = 555 K (assuming zero bond albedo) and a T eff = 3026 K (Charbonneau et al. 2009).The closest values in F14 are 575 K and 3000 K.For Neptune, we followed a similar procedure, but instead of using the observed T eq = 46 K from (Guillot 2005, bond albedo of 0.290), we used T 1bar = 72 K, with T local = 75 K and T eff = 6000 K being the closest value in F14.We adopted this approach because Neptune's radius is frequently defined at 1 bar (which is used here instead of the transit radius in our analysis; see Sect.3).
While this approach lends itself to estimating a value for P local in the case of Neptune, for GJ 1214 b we compared the temperature profiles using mean opacity values of different P local to profiles from different literature sources (e.g.Miller-Ricci & Fortney 2010;Valencia et al. 2013) and selected the most consistent profile, namely P local = 0.1 bar.In cases where such comparison data are not available, P local could be allowed to vary as a free input parameter.
In our model, planetary radius is taken to be the transit radius following the methodology described in Guillot (2010) and defined as the pressure level where the chord optical depth in the visible range τ ch is equal to 0.56 (Lecavelier des Etangs et al. 2008).This quantity is related to the optical depth by: where H is the scale height and R solid is the total sum of the planetary thickness of the water, mantle, and core layers.We note that the temperature and scale height dependence can lead to a 30% decrease in the optical depth at the effective planetary radius for hot Jupiters compared with warm super-Earths such as GJ 1214 b.

High-pressure envelope
Sub-Neptune atmospheres can extend over a large range of temperature and pressure such that equations of state based on the ideal gas law can begin to break down in the deeper layers.
In our work, in the high-pressure envelope, we use the EoS first described by Saumon et al. (1995) and expanded upon by Chabrier et al. (2019), hereafter collectively referred to as (SCvH).
The boundary condition between the irradiated gas layer (described above) and the underlying (radiative) envelope is determined by a pressure criterium similar to that applied by Jin et al. (2014).The switch occurs at a pressure level where incoming energy from the star has been absorbed and the main energy source is the thermal contribution of the planet; in other words, where τ ≫ 1.Using the relation in Rogers et al. (2011) τ ≫ 1/( √ 3 γ), we set the specific value of τ = 100/( √ 3 γ), which corresponds to a pressure level of approximately 10 bar.At greater pressures, the Guillot temperature and pressure are used as inputs for the non-ideal EoS density calculation as performed by SCvH.The temperature profile is calculated along the radiative gradient (∇ rad ) defined by Jin et al. (2014), as: where T int is the planetary intrinsic temperature.This procedure continues until the convective regime is reached, as defined by the Schwarzschild criterium: whereby the advective gradient (∇ ad ) in the SCvH table is compared with the radiative gradient stated above.Once this condition is reached, the entropy is fixed and the temperature A65, page 3 of 15 A&A 671, A65 (2023) and density are calculated isentropically.Appendix C discusses caveats associated with applying Eqs. ( 3) and (4) in the way described above.

Interior
We consider the following compositionally distinct layers in modelling the solid interior of each MC realisation: an iron core, a silicate mantle, and a water layer.We provide here a brief overview of the model assumptions.For more details on the model, see Baumeister et al. (2020).
The density of Earth's core suggests that it is not composed of pure iron, but that lighter alloying elements are present (Rubie et al. 2015).Even for Earth, and especially for exoplanets, the exact amount of these lighter alloys is not known well.Therefore, we follow an approach by Valencia et al. (2007) assuming a pure iron core and using a Mie-Grüneisen-Debye equation of state (Sotin et al. 2007) to calculate the density profile.
For the silicate mantle, we assume a pure composition of MgSiO 3 perovskite (bridgmanite) for the entire mantle.We model the mantle with a Vinet EoS (Vinet et al. 1989) using the suggested parameters from Oganov et al. (2001).We allow all MC realisations to have an extended water layer above the silicate mantle to cover the possibility of 'water worlds' (Léger et al. 2004;Zeng et al. 2019).We employ the AQUA EoS (Haldemann et al. 2020), which spans a pressure-temperature range suitable for interior conditions and which calculates the density and adiabatic gradients based on a suite of H 2 O EoSs.
We do not consider temperature dependence for the silicate mantle or the iron core, because the associated change in radius contribution is negligibly small for our purposes here (Seager et al. 2007).Moreover, the planetary radius in our MC realisations is mainly determined by the extended hydrogen atmosphere and envelope, whereas changes in the solid interior contribute only as second-order effects.However, it should be noted that thermal effects can have a significant impact on the thickness of the water layer (Thomas & Madhusudhan 2016).As such, we use the following adiabatic gradient to calculate the temperature profile for the water layer: dT dr = T P dP dr ln T ln P ad . (5)

Modelling method
The model takes as a boundary parameter the total mass of the planet along with starting values for mass fractions of the core, water, and gas layers, with the remaining mass in the mantle.Using a Monte Carlo approach, we distribute 5000 models over this parameter space, where the solver varies the radius randomly within user-defined minimum and maximum values.
For each mass fraction combination of core, water, and gas (with the mantle making up the rest), the model begins by calculating the temperature profile.From this, the pressure and gravity profiles are calculated, which serve as input for the EoS of each layer when calculating the density profile.With the densities and planetary radii, the solver integrates the total mass for each layer until either the pressure (for the atmosphere-envelope boundary) or mass (for the core, water, and gas layers) boundary condition is reached.Once the solver has reached the centre of the planet, it adjusts the radius for the next iteration.If mass is left over, the model increases the radius in the next iteration, and vice versa.Mean opacity input parameters Notes.Mass and radius for GJ 1214 b are taken from Charbonneau et al. (2009) to be consistent when comparing results with Valencia et al. (2013) and Rogers & Seager (2010).We assume 5% error in mass and radius for Neptune as representative of current sub-Neptune data precision. References.
This continues until the variation between consecutive iterations is less than the tolerance.
On completion of the parameter study, the results are filtered based on their final radii such that only solutions whose transit radius lies within the assumed error range of the observations (typically 1000 from the 5000 model realisations in total) are used in the analysis.These solutions are then plotted in a 'corner plot' that shows the posterior distribution of each mass fraction.The median values of these distributions are calculated along with the values at the 84th and 16th quantiles (equivalent to ±1σ for a Gaussian distribution).These median values are used to generate individual profiles for each MC run and they serve as the basis for our comparison to investigate how the input atmospheric mean opacity values influence the interior structure.Table 1 shows the key input parameters for the planets simulated.
It should be noted in general that the median values and their error bars are calculated independently from one another, which is why the sum of the median mass fractions (core, mantle, water, and gas) is not equal to one.For this reason, choosing such combinations of median values does not represent a self-consistent interior solution found by the parameter study.Alternatively, one would choose an interior structure by selecting a specific value for either core, water, or gas in order to find the mass fraction of the other layers (see Sect. 4 for an example).

Interpreting corner plots
Corner plots are useful for visualising how the distribution of solutions in a multi-variable parameter space vary with respect to one another, along with the independent distribution of solutions of each variable; see Fig. 1 for an example.In this work, we use the mass fractions (w) of each layer as free parameters.The histogram plots along the diagonal show the distribution of layer mass fractions that result in a valid solution.The contour plots show how each layer varies with respect to the other layers; that is, for a specific gas mass fraction, they show the allowed range of core, mantle, and water mass fractions.The A65, page 4 of 15 contour plots have lines drawn at each σ level (with solid, dashed, and dotted lines corresponding to 1σ, 2σ, and 3σ, respectively) and feature multiple local minima throughout the parameter space.
We also denote the median value of each distribution as a solid vertical blue line, with the upper and lower bounds shown as vertical dashed blue lines calculated at the 84% and 16% quantiles, respectively.The median and its associated errors are written above each histogram.

Validation: Neptune
We consider modern Neptune to validate our atmosphereinterior model.While there are a number of caveats involved with this comparison (see below), our main aim here is to compare observed temperature profiles with those calculated by our semi-grey scheme for Neptune conditions in the upperatmosphere irradiated regime.For this purpose, the validation is assumed to be appropriate.The main caveats are as follows: 1. Neptune's bulk mass lies just above the intended range of our (sub-Neptune) model.Furthermore, current Neptune interior models in the literature (e.g.Scheibe et al. 2019) generally apply a continuous interior modelling approach rather than a stratified planet with distinct layers as is the case in our model.As such, our model shows more variation in the final interior structures, with the "true" Neptune interior included within the range of our model output.2. Our model water layer is composed of pure water, similarly to that of Scheibe et al. (2019), rather than a mixture of water, methane, and ammonia for example, as in Podolak et al. (1995).3. The metallicity of our atmosphere is limited by the range available in our input opacity tables (F14).We use the highest available, [M/H] = 1.7 or 50× solar metallicity, although current estimates for Neptune suggest a value of ∼80 ± 20× (Mousis et al. 2020).4. Results are analysed using the radius at 1 bar as opposed to the transit radius.We assumed error bars of 5% for both the mass and radius, which we take to be representative of current sub-Neptune data.With these caveats in mind, results from our model suggest a median gas mass fraction of 2.05 M ⊕ , with a range of 1.49-2.73M ⊕ (Fig. 1 bottom right hand plot: log w gas ≈ −0.92 ⇒ 17.147 × 10 −0.92 = 2.05 M ⊕ ).Our results correspond to the gas mass fraction calculated for Neptune models in the literature, that is, from 2.3 M ⊕ to 3.45 M ⊕ (Podolak et al. 1995;Scheibe et al. 2019), which is within our ±1σ range.We used the value of 2.05 M ⊕ as an input for Neptune when modelled with the Chabrier et al. ( 2019) EoS based on previously published models that satisfy the gravity data for the second (J 2 ) and fourth (J 4 ) gravitational moments (Helled et al. 2020).The latter value was obtained based on the modelled output of J 2 and J 4 which compared best with Voyager observations of Neptune (renormalised from Owen et al. 1991).Although Neptune lies just above the mass range of exoplanets intended for our model, our results are nevertheless broadly consistent when compared to those obtained from more dedicated models.
Figure 1 presents corner plot results for the modern Neptune validation, which suggest left-skewed mass fraction distributions (top panels, blue lines) with reasonably constrained relations between the gas fraction and the other layers: w water , mantle, and core (Fig. 1, bottom row) but with less constrained (more degenerate) relations for the remaining layers (Fig. 1, second and  third rows).The tendency for the solutions to skew left towards smaller mass fractions is due to the increased number of possible layer combinations where w layer,i = 1 that result in a valid solution.
Figure 2  To test the robustness of our results, we doubled the number of MC cases from 5000 to 10 000 and compared the resulting change in the layer mass fractions output.None of the mass fractions changed by more than a few percent, which suggests that the number of MC cases used is sufficient.

Case study: GJ 1214 b
For the investigation of the atmosphere of GJ 1214 b, we compare our calculated pressure-temperature profile assuming 1× solar metallicity to the grain-free solar composition atmosphere in Valencia et al. (2013) andMiller-Ricci &Fortney (2010).For the interior validation, we compare with the model study of Rogers & Seager (2010).

Solar metallicity
Table B.2 lists the thermal mean opacity, κ th , and visible mean opacity, κ vis , from F14 as well as the associated γ (κ vis /κ th ) value assumed in our work.Valencia et al. (2013) used the F08 mean opacities for their κ th , and selected a γ value of 0.032 to give a temperature of 1000 K at 1 bar based on Miller-Ricci & Fortney (2010), which was kept fixed in their various scenarios.To select our κ th and κ vis we chose a scenario that was closest to that of Valencia et al. (2013) with γ = 0.0277, or a P local = 0.1 bar for a T local = 575 K and a metallicity of 1× solar.
Figure 3 suggests overall reasonable agreement between the modelled temperature profile from our work compared with that of the Valencia study, although there are some key differences.There is increased warming in our TATOOINE results in the upper and lower vertical sections due to the increased amounts of greenhouse gases (see Table B.1) included in our calculation of the mean opacities in F14 versus F08.We also include a highpressure envelope in our atmosphere, which contributes to the divergence from Valencia et al. (2013) at the bottom of the atmosphere, because their work varied κ with pressure, temperature, and metallicity using extrapolated values.
Figure 3 also includes profiles generated in Miller-Ricci & Fortney (2010).Rather than a semi-grey approach as in Valencia et al. (2013) and Rogers & Seager (2010), Miller-Ricci & Fortney (2010) employed a full radiative-convective model.Our atmospheric pressure-temperature profile presented in Fig. 3 shows a better fit to the profile of Miller-Ricci & Fortney (2010) than that of Valencia et al. (2013) for pressures <1 bar, which suggests that our use of the updated mean opacity values leads to an improved comparison with more sophisticated models, which do not rely on such approximations.In order to further reduce the deviation at higher pressures, one likely solution is to use a combination of a pressure-and temperature-dependent κ, such as in Valencia et al. (2013), along with the inclusion of a high-pressure envelope.

Interior structure
For our interior structure output, we compare with the model study of Rogers & Seager (2010).In their work, these authors investigated a variation of γ in the range 0.1-10.We note that they used the Planck mean opacity, while our work (and that of Valencia et al. 2013) uses the Rosseland mean (see discussion above).
In Fig. 4, the standard deviation of the solutions around the median w gas spans a range of over three orders of magnitude (0.18×10 −3 − 0.62%).While this range is about an order of magnitude smaller than those suggested in Valencia et al. (2013) and Rogers & Seager (2010), we note that our work only makes use of the median values for our comparisons.These represent the most common points in the solution space or, in other words, the layer mass fractions for which there exists the largest number of degenerate solutions that fit the bulk planetary parameters.These solutions however do not necessarily represent the most physically likely.
If we refer to Fig. 4, and look at solutions where w water → 0 in order to be comparable with the Valencia et al. (2013) study, the w gas spans a range of 10 −2 -10 −1 , approximately.For Table 2, we selected the parameter study solution with the smallest w water .The associated w gas = 1.57% is in agreement with Valencia et al. (2013) for results where w water ≈ 0.
A65, page 6 of 15 Of the three cases discussed in Rogers & Seager (2010), their 'case I' is the most similar to our 1× solar composition atmosphere case and is chosen here for further investigation.For case I, their work shows gas mass fractions with a minimum range of 0.9 × 10 −2 -0.2% and a maximum range of 3.2-6.8%,depending on the atmospheric thermal profile when the mass and radius are allowed to vary within 1σ.While the minimum range is overall consistent with this work (0.18 × 10 −3 -0.62%), the difference in both the mean opacities and treatment of γ likely contributes to the difference in the upper range of solutions.While the Planck and Rosseland mean opacities converge at pressures higher than several tens of bar, the Planck mean opacity can be several orders of magnitude higher than the Rosseland value for pressures lower than ∼1 bar depending on, for example, gas composition and temperature.Additionally, in their scenarios, γ is allowed to vary from 0.1 to 10.0 as a free parameter, while our γ value is kept fixed at 0.0277.These two factors likely account for their higher maximum gas mass fraction, which for our model is outside the ±1σ range of our results.To illustrate this point, the right hand column of Table 2 shows the parameter study solution with the greatest gas mass fraction, w gas = 3.96%, which is over 6.5σ from our median value.

Variation of metallicity
The atmospheric metallicity of GJ 1214 b is not well constrained.Miller-Ricci & Fortney (2010) assumed values ranging from 30× and 50× solar in addition to the 1× case, whereas the model study by Drummond et al. (2018) assumed up to 100× solar.The observed [Fe/H] metallicity of the host star GJ 1214 is +0.39 ± 0.15 (Harpsøe et al. 2012), which corresponds approximately to 1.74×-3.46×solar values.While this does not necessarily correlate with the planetary atmospheric composition (in our own Solar System, Neptune has a composition equivalent to greater than 50× solar metallicity), it offers a starting value when no other information is available.
Other model studies, such as Dorn et al. (2015), used the stellar metallicity as an upper limit when modelling the Fe/Si and Mg/Si ratios as input for the core and mantle.However, the main focus of our work is the planetary gas layer, and therefore changes due to variations in the core and mantle were treated in a straightforward way: the Fe/Si and Mg# were kept fixed for all runs using Earth-like values (0.19 and 0.9, respectively) and their temperature dependence was not considered; this assumption has little effect (<100 km) on the final radius.Another possibility is to allow the atmospheric metallicity to vary as a free parameter, such as in Dorn & Heng (2018), although in our work we are limited by the range of metallicities available in F14.We computed the atmospheric pressure-temperature profiles for all metallicities available in F14, and compared them to the profiles from Valencia et al. (2013) and Miller-Ricci & Fortney (2010; see Table B.2). Figure A.1 presents the corner plots for all metallicities, along with an aggregate corner plot that combines all the parameter study solutions for each metallicity.The median values for each metallicity are listed in Table B.3.A general trend can be seen in the corner plot distribution in Figs. 4 and A.1, where an increase in metallicity is associated with a decrease in gas mass fraction (e.g.Table B.4).As the metallicity increases, so too does the amount of constituent gases such as H 2 O, CO 2 , and CH 4 (see Table B.1).This leads to additional warming in the atmosphere due to increased absorption, which increases the scale height and less gas mass is needed to fit a given planetary radius.
As can be seen in Fig. 5, while our model output agrees with the results of both Valencia et al. (2013) andMiller-Ricci &Fortney (2010) for the 1× case, the profiles deviate from Miller-Ricci & Fortney (2010) for the 30× and 50× metallicities.These differences likely arose for example because Miller-Ricci & Fortney (2010) employed a local-thermal-equilibrium model developed for irradiated exoplanets, which includes various atmospheric processes -such as convection -not included in our model for our upper atmosphere, only in the high-pressure envelope.
Employing the semi-grey approximation leads to a deviation of calculated atmospheric temperature from the equilibrium temperature of the planet.This arises due to, for example, greenhouse warming and the contribution from upwelling thermal radiation from the planetary interior.For Earth, the greenhouse component is of the order ∼33 K at the surface.For our GJ 1214 b scenarios, the estimated equilibrium temperature is in the range of 393 to 555 K (Charbonneau et al. 2009).Our modelled T Tran values by comparison range from ∼490 to 760 K (see upper and lower dashed lines in Fig. 9 respectively).This difference of up to A65, page 7 of 15 A&A 671, A65 (2023) around 450 K between semi-grey and equilibrium temperatures reflects the importance of including accurate semi-grey schemes for sub-Neptunes.
To validate our interior results at higher metallicities, we compare with the so-called "Super-Earth" (SE) case in Rogers & Seager (2010).The most notable differences are the lack of He in their atmosphere, and the absence of a water layer in the planetary interior for their SE case.They calculated a gas mass fraction of 0.05%-1.2%.This is consistent with the gas mass fractions for higher metallicities (10×, 30×, and 50×) in Table B.4.The table suggests a trend with mostly a reduction in core mass fraction going from 3× to 50× metallicity.As the metallicity increases, the interior structure changes from a 'Mini-Neptune' with a large primordial atmosphere to a 'Super Earth' with a more modest gas mass fraction.
This trend is less apparent with the addition of a water layer, as can be seen in the corner plots in Fig. A.1 and Table B.3 for the median values of each mass fraction, along with their associated errors.While the same inverse relation between metallicity and gas fraction is still present, the resulting changes in the distributions for core, mantle, and water mass fractions are weaker.This suggests that it is challenging to constrain the interior based on our modelling results.Additional input constraints, such as from more sophisticated models or from other observed parameters (e.g.atmospheric spectral data or the Tidal love number k2), combined with more precise data (e.g.PLATO) are needed in order to be able to better resolve these degeneracies.

Result, comparison of 2008 and 2014 mean opacities
We investigate here the influence of the updated F14 compared with the F08 mean opacity values for different γ values.We used both the γ quoted in Valencia et al. (2013), along with the γ from F14 used previously throughout this work (Table B.2), or γ F08 = 0.032 and γ F14 ≈ 0.0277, respectively.Comparing Fig. 4 with Fig. 6 suggests that while the changes due to the updated opacities are modest, they are more significant than the differences between Fig. 6a (γ F08 ) and Fig. 6b (γ F14 ).This suggests that variations in mean opacity values as individual input parameters have a greater impact on the final model output than changing γ = (κ vis /κ th ) alone.This can also be seen in Table 3, which shows the change in each layer mass fraction and transit radius as a percentage difference.
Table 3 suggests that updating the opacities leads to an increase in the gas and water mass fractions, which is compensated by a decrease in the core and mantle mass fractions.This is different from the response discussed above (see Sect. 4.3) where increasing the metallicity led to a decrease in gas fraction due to heating through the atmosphere.The changes in Col. 1 compared to Col. 2, which demonstrate the effect of changing the mean opacities for each γ value, feature similar values of a few tens of percent.The values of Col. 3, which shows the effect of only changing γ when using the 2008 opacities, are approximately one order of magnitude smaller.This suggests that the effect of changing γ = κ vis /κ th is not large when compared to changes in κ th and κ vis .The transit radius, R P,Tran , changes by several tenths of a percent when using the median values and by up to a few percent for the specific (w core , w water , w gas ) test case (marked with an asterisk in Table 3).
Figure 7 compares temperature profiles from our model (using different combinations of 2008 and 2014 κ and γ values) with other works in the literature.Our results using the updated Fig. 6.Comparison of GJ 1214 b parameter study solutions for a 1× solar composition atmosphere using mean opacities from κ F08 (F08), and either (a) γ F08 (Valencia et al. 2013) opacities (solid purple line) lie closest to the profile from Miller-Ricci & Fortney (2010), compared to the profiles using the 2008 mean opacities, both from this model and that of Valencia et al. (2013).Updating the opacities leads to cooling in the middle atmospheric region, which leads to atmospheric contraction, which is why more gas is required to fit the observed planetary radius.This results in the increased gas mass fractions shown in Table 3.
A65, page 8 of 15 Table 3. Percentage difference due to changes in opacity and γ values of the median value and transit radius of each layer for the runs in Figs. 4, 6a, and 6b.
GJ 1214 b -F14 (1) vs. F08 (2) percentage difference w core -12.40%-13.96% -1.39% w mantle -6.04% -8.48% -2.30% w water +5.26% +6.88% +1.71% w gas +13.09% +11.18% +2.19% R P,Tran (R ⊕ ) -0.23% -0.49% -0.26% Notes.See Table B. Figure 7 suggests that the effect on P Tran of using the 2014 values (solid circle symbol) compared with the 2008 values (cross symbol) is not large.However, we note that we are comparing median values derived from two large distributions.Many of the MC realisations using the 2008 opacities lead to solutions lying within the observed radius interval.This is also the case for the MC realisations using the 2014 opacities.Comparing the resulting two median values favours smaller differences due to averaging.It is more instructive to compare MC realisations away from the median, in order to gain insight into how large individual effects can be.In addition, the median values compared have different mass fractions, which are calculated independently for each of the layers.We therefore selected a case with a combination of (w core , w water , w gas ) mass fraction which counted as a valid solution for the (κ F08 , γ F08 ) case, and used this as input to generate the comparison profiles for (κ F14 , γ F14 ) and (κ F08 , γ F14 ).This case is marked with an asterisk in Table 3.
The percentage radius difference for (κ F14 , γ F14 ) vs. (κ F08 , γ F08 ) taken from Fig. 9 is equal to where R P is the observed radius of the planet (Table 1).We note that this difference results in the transit radius (solid circle) in Fig. 9, which was calculated using the updated opacity and γ values lying outside the quoted error bars, and as such would not be considered a valid solution in our analysis.See Table 3 for the remaining percentage values.
Figure 8 shows the corner plot solutions that are unique for the 2014 input values taken from Fig. 4 (shown in red) and for the 2008 input values taken from Fig. 6 (shown in blue), with the median (solid) and quantiles (dashed) listed separately in Table 4, along with the percentage changes.A similar comparison with (κ F08 , γ F14 ) (not shown) did not differ significantly.
Figure 8 (bottom row) suggests a clear separation between a red (blue) region for lower (higher) core and mantle mass fractions and higher (lower) water and gas mass fractions.This behaviour is confirmed in the panels on the far right showing the mass fraction distributions.
difference is greater than the previous example of ∼1.7% (see Fig. 9), which was calculated for only one value, it is still not significant when the measured radius precision is σ R P ≈ 5% (Charbonneau et al. 2009).However, compared to future missions, such as the upcoming PLATO mission, differences in planetary radii on the order of a few percent will be achievable (Rauer et al. 2014).Nevertheless, this result depends on numerous other factors, such as uncertain metallicity.However, we note that PLATO does not plan to observe GJ 1214 b.We also note that observations (Kreidberg et al. 2014)   (2) ,γ F08 (3) w core 14.56 (3) Valencia et al. (2013).
spectrum for GJ 1214 b, which is likely due to the presence of clouds, which are not included in our model.Even so, our results suggest that forthcoming data will nevertheless necessitate more accurate models with improved input parameters in order to be able to accurately interpret the results.

Conclusions
The parameter space of sub-Neptunian exoplanets is large, covering a diverse range of interior-atmosphere solutions that are highly degenerate.Interior-atmosphere models must be flexible enough to investigate the large parameter space by balancing the accuracy of the results with the required computational time.
Our model, which uses the semi-grey approximation, produced results that are comparable with those of more sophisticated atmospheric models, while only requiring a fraction of the computational cost.This makes it possible to conduct thousands of model realisations, investigating the effect of a range of uncertain input parameters.By varying the metallicity as a proxy for the atmospheric composition, increased absorption of longwave radiation from the interior can raise the middle atmosphere temperature of GJ 1214 b by several hundred degrees compared with the equilibrium temperature.We note that there exist additional factors other than the input opacities, which can affect the calculated planetary transit radius.For example, the modelled atmospheric temperature profile is based on the planetary dayside of the planet, not on the terminator region where the transit radius is typically measured.Another factor is scattering in the visible and UV, which can lead to an increase in the determined transit radius.The use of additional data, such as observations over extended wavelength ranges, as well as constraints from spectral features, including the absorption bands of individual species and the Rayleigh slope (which constrains the atmospheric bulk molecular weight), can be used to constrain the atmospheric composition, potentially providing information about the planetary interior.
We find that varying the mean opacity had a larger effect on the distribution of interior solutions than varying γ.Degeneracies in the distribution of interior-atmosphere solutions cannot be constrained by the choice of atmosphere input parameters alone.Modest changes in the atmospheric gas fraction can lead to changes in transit radius that are of similar magnitude to the changes resulting from variations in the core and mantle mass fractions.The use of updated mean opacity and γ values in a semi-grey atmosphere results in improved calculations of temperature and transit radius when compared to studies that used older data.Such updates are needed for preparing the models for comparison with planned missions in exoplanet science.

Appendix C: Caveats to temperature profile calculation
Here we discuss caveats relevant to calculating the temperature profile through the atmospheric regions from top to bottom: Upper radiative regime.The central equation is taken from Guillot (2010), their equation 49.This is derived by imposing energy balance between incoming (stellar) and outgoing (planetary) radiation and assumes a mean opacity in the visible (κ vis ) and mean opacity from intrinsic heating in the infrared (κ th ).Emission from the atmosphere is assumed to have a negligible contribution to visible wavelengths.The methodic employs the Eddington approximation, assuming a static plane-parallel atmosphere in local thermodynamic equilibrium.The Guillot equation was originally derived in the context of hot Jupiters (HJs) and was applied at the time to HD 209458 b, for example.However, in our work, we simulate somewhat smaller, cooler objects such as GJ 1214 b.Table C.1 therefore compares caveat-relevant planetary and stellar data for HD 209458 b and GJ 1214 b.
Table C.1 suggests higher T eq and T int by factors of 3 and 4 for the HJ case.Planetary radius is about a factor 6 smaller for the SE case which degrades somewhat the plane-parallel assumption mentioned above.Gravity is broadly similar for the two cases in Table C.1.Kappa values are also broadly similar for the two objects considered in Table C.1.In general atmospheric hazes could change κ vis depending on particle radius whereas enhanced water vapour on super-Earths would likely lead to increased κ th values.
Switch-over from ideal to non-ideal density calculation.In our work we define a switchover boundary at t = 100/( √ 3g) which corresponds to about p>10 bar.At greater pressures, it is necessary to employ a non-ideal calculation of density.For this calculation, we use the SCvH non-ideal equation of state which yields density from a lookup table and requires temperature and pressure as input.Knowing for example the density and the volume of the gas layer enables calculation of the mass fraction, which is an output diagnostic in our cornerplots.
Radiative versus diffusional regime.A caveat of our work is that, at the switchover point of p>10 bar, we not only switch to a non-ideal density calculation but also apply equation (3) (hereafter 'the diffusional approximation') with increasing pressure for the temperature calculation until the convective regime is reached as defined by the Schwarzschild criterium.However, the diffusional approximation implicitly assumes that the radiative transfer of photons at high pressure arising from the planetary interior is limited by diffusion.However, a preferable approach which we note for future work would be to continue using the Guillot equation for temperature because this equation converges to the same solution as the diffusional approximation at high pressure.The convergence pressure depends for example on the internal temperature whereby higher internal temperatures lead to convergence at lower pressures.We now investigate the effect of this and other caveats.
The two panels (Fig. C1a and C1b) below show a range of temperature profiles for GJ 1214 b for the control run (for more details of this run see our Figure 3 above and the associated text).The left panel shows the effect of varying T eq and T int over their uncertainty range.This panel shows cooling of about 140K in the upper atmosphere radiative regime upon reducing T eq from 555 K (control) to 393 K and an increase in lower atmosphere temperatures by up to several hundred degrees while increasing T int from 24 K (control) to 50 K.The right panel shows the effect of switching from the Guillot temperature to the diffusional approximation temperature at pressures greater than 10 bar where considerably larger temperatures are calculated by the diffusional approximation in the deep atmosphere.

Fig. 1 .
Fig. 1.Neptune corner plot showing the distribution of layer mass fractions (w layer ) for an assumed 50× metallicity atmosphere.Input values are shown in the light blue panel in the upper-right corner.Solid blue circles show median values.Error bars in the uppermost panels denote the 16th and 84th percentiles based on the distribution from the MC output.

Fig. 2 .
Fig. 2. Validation of the atmospheric pressure-temperature profile for Neptune with 50× metallicity using median values from Fig. 1 as our model input (solid green line).Observations taken from Fig. 1 of Marten et al. (2005) are plotted in red.

Fig. 4 .
Fig. 4. Corner plot for GJ 1214 b showing the distribution of solutions for the parameter study with a 1× solar composition atmosphere (Table B.2).

Fig. 5 .
Fig. 5. Atmospheric pressure-temperature profiles and transit pressures (see Table B.3 for R Tran values) for GJ 1214 b for all metallicities available in F14 (solid lines), compared with profiles from Valencia et al. (2013) and Miller-Ricci & Fortney (2010; dashed lines).See Fig. A.1 for the associated corner plots.

Fig. 7 .
Fig. 7. Atmospheric temperature profiles for GJ 1214 b calculated using the median values from the parameter study runs in Fig. 4 (solid), Fig. 6a (thick dashed), and Fig. 6b (thick dash-dotted) compared to the profiles from Valencia et al. (2013; thin dashed) and Miller-Ricci & Fortney (2010; thin dotted).The transit pressures (P Tran ) are marked as solid symbols and the associated R Tran values can be found in Table B.3.

Fig. 8 .
Fig. 8. Corner plot for GJ 1214 b showing the distribution of parameter study solutions that are unique to either (κ F14 , γ F14 ) (shown in red) or (κ F08 , γ F08 ; shown in blue) refer to Figs. 4 and 6, respectively.The associated median values and error bars for each layer are listed in Table4.

Fig. 9 .
Fig. 9. Temperature profiles using both mean opacity databases and γ values for the selected (w core , w water , w gas ) (see text) and the associated transit radius (see Table B.3 for R P,Tran values).The solid grey line denotes the observed planetary transit radius of GJ 1214 b, whereas the dashed lines represent the radius error bars in Charbonneau et al. (2009).

J
Fig .Cla: Guillot Temperature Profiles

Table 1 .
Model input parameters used in this work (refer to Sect.2.1.1 regarding the mean opacity parameters).

Table 2 .
GJ 1214 b layer mass-fraction solutions associated with the minimum w water and maximum w gas from our parameter study for a 1× composition atmosphere.
Notes.Similar results for the other metallicities studied can be found in TableB.4.

Table B .
4: Comparison of solutions with either (a) the minimum w water or (b) the maximum w gas for all metallicities.

Table C .
1: Planetary and stellar parameters for HD 209458 b and GJ 1214 b. *Upper limit.Values are taken from references in this work (section 4) and from Guillot (2010) and references therein.