© The Authors 2024

1 Introduction

We have developed a model atmosphere code consisting of an updated version of the classical Model Atmospheres in Radiative and Convective Scheme (MARCS) code (Gustafsson et al. 1975; Gustafsson et al. 2008) combined with the cloud formation code StaticWeather (Helling & Woitke 2006; Helling et al. 2008b) and the chemical equilibrium code GGchem (Woitke et al. 2018). We name the new code and the corresponding model atmospheres MSG, an abbreviation of its three major components (MARCS-StaticWeather-GGchem). For warmer models, the atmospheric temperature versus gas pressure structures of the MSG models are quite similar to the corresponding models of the MARCS grid from 2008 (Gustafsson et al. 2008), which we refer to in this work as ‘the classical MARCS models’. Minor differences are found when one approaches the lower end of effective temperatures of T_eff = 2500 K in the 2008 classical MARCS grid. There are no published classical MARCS models below T_eff = 2500 K, but we show what they would have looked like in this work and illustrate what was missing in the classical models to describe these lower temperatures and what the differences to the new grid would, approximately, have been. The new grid is at present computed down to effective temperatures of T_eff = 300 K, which will allow for studies of the Earth-like exoplanetary temperature regime, and in models still to be computed, it will enable study of the atmospheric effect of rocky surface development and biological evolution.

A basic philosophy behind the MARCS models, as well as the StaticWeather and GGchem computations, has been to keep all computations fully self-consistent, with no free or adjustable parameters to improve fits to observed spectra. This philosophy is inherited into the MSG grid. We believe that this approach maximises the visibility of unknown (or non-included) physics and missing data, with a long-term goal of revealing, for example, a potential variety of known or unknown biological effects on exoplanetary atmospheric structures. The existing grid of classical models available at the MARCS website¹ contains a few tens of thousands of models covering a broad range of dwarf and giant stars of many different metallicities in the effective temperature range from 8000 K to 2500 K, as described in detail in Gustafsson et al. (2008). This huge number of easily available models has made it possible to test the model predictions against a large variety of stellar objects, and the more than 3000 citations in the scientific literature to the two basic model grid papers (Gustafsson et al. 1975, 2008) serves as evidence of a broadly applied and well-tested use of the classical MARCS models.

The MARCS code was originally written (Gustafsson et al. 1975) with ‘solar-type’ (i.e. F, G, K type) stars in mind. Such stars are well represented by plane parallel, local thermodynamic equilibrium (LTE) models in chemical equilibrium, with relatively few opacity sources, and they converge well in a scheme where variations in the opacity are calculated as a function of temperature and electron pressure, and the adiabatic index is computed based on only a few molecular and other opacities. During the almost five decades that have passed since the initial 1975 grid, the models have been modified, updated, expanded, and augmented several times by the introduction of new physics and basic input data suitable for a much larger parameter space, including the introduction of more complete atomic and molecular opacities (e.g. Jørgensen et al. 1985), as well as by updating of the opacity distribution function (ODF) method to the opacity sampling (OS) method in handling of the absorption coefficients (Saxner & Gustafsson 1984; Jørgensen 1992), which all made it possible to compute stars of a wider range of effective temperatures (for an overview see Gustafsson & Jørgensen 1994). Inclusion of spherical geometry has made it possible to reach larger stars of lower surface gravity (Plez et al. 1992; Jørgensen et al. 1992), and the first attempts of a description of cloud formation in the models was introduced in Juncher et al. (2017).

To assure the reliability of the classical models, it was decided to limit the 2008 grid to temperatures above T_eff = 2500 K. On top of this, models usually had convergence problems at the lowest temperatures, which could further support the suspicion that something was physically unrealistic at the lowest effective temperatures. In the present paper, we show that this suspicion was correct to a large extent, and we introduce new data and methods that are consistent with the classical models for T_eff = 2500 K and above, while at the same time extends the grid from its present lower limit of T_eff = 2500 K down to T_eff = 300 K. This extension completes the grid to include the coolest known stars and brown dwarfs at temperatures overlapping with the bulk of the hot-Jupiters and Earth-like exoplanets too, including all types of M, L, and T-dwarfs and probably reaching the bulk of the disputed Y-types, too (see, e.g., Burningham et al. 2008; Cushing et al. 2021 and our discussion of NH₃ in Section 6). A corresponding extension of the grid to warmer models than the present upper limit of T_eff = 8000 K is in progress by Edwardson et al. (in preparation).

The extension of cloud-free models down to T_eff = 300 K has been made possible by the inclusion of a considerably more complete set of molecular opacities assembled from the Exo- Mol database (Tennyson & Yurchenko 2012; Tennyson et al. 2016; Tennyson et al. 2020; Chubb et al. 2021), by inclu- sione of the extensive new chemical equilibrium routines of GGchem (Woitke et al. 2018) that allow for computation of the relevant species down to T = 100 K, and by including a new way of handling the electron pressure in the iterative convergence scheme, and by other numerical improvements. The corresponding cloudy models have been made possible by the integration with the self-consistent cloud formation program StaticWeather (Helling & Woitke 2006; Helling et al. 2008b). The MSG computations are able to compute, in a self-consistent and homogenous way, both the host stars and their irradiated exoplanets based on an improvement of the more simple radiation schemes that began with the series of papers by Nordlund and Vaz for the description of binary stars already in the 1980s (Vaz & Nordlund 1985; Nordlund & Vaz 1990).

The present updates as well as the planed extensions have been made possible due to two major grants (mentioned in the acknowledgments) that form the backbone of the Centre for ExoLife Sciences². Under the auspices of this project, the effects of irradiation (Amadio et al., in prep.), details of the atmospheres of cloudy exoplanets (Campos Estrada et al. 2024), and application to rocky exoplanets and its consequences for the effects of biological evolution on the atmospheric structures (D’Alessandro et al., in preparation) are all part of the ongoing extension of the grid.

In the present paper, we describe the overall results, the basic methods, and the grid of cloud-free models and corresponding synthetic spectra, and we introduce the ongoing extensions to cloudy, irradiated, and non-equilibrium exoplanet models. The paper is structured with Section 2 giving a short summary of the basic assumptions and methods of the three merged codes; Section 3 outlining the degree of consistency with the 2008 classical MARCS-grid; Section 4 focusing on some challenges, describing the addition of the new molecular line lists, partition functions, and chemical equilibrium routines; Section 5 discussing the importance of the correct treatment of the line profiles; Section 6 focusing on aspects of the total opacities; Section 7 describing effects of irradiation; Section 8 describing effects of cloud formation; Section 9 providing considerations about chemical non-equilibrium; Section 10 describing the boundary between M-type and L-type as well as the transition from L-type to T-type; Section 11 demonstrating the difference between oxygen-rich and carbon-rich models. Finally, the conclusions and outlook are outlined in Section 12.

2 Summary of the three merged model codes

2.1 MARCS

The classical MARCS models are cloud-free, 1D, stratified, flux-constant, LTE models, computed in radiative-convective equilibrium in an optical depth scale, and solved iteratively by the Newton-Raphson method. With some modifications, they form the basic structure for the solution of the energy balance in the MSG models too. The models are characterised by the effective temperature, T_eff, the surface gravity, g, and the relative abundances of the elements. Typically, the gravity is given as log10(g) in cgs-units such that the log(g) of the Sun, Jupiter, and the Earth are 4.44, 3.39, and 2.99, respectively, and therefore the log(g) of the respectively cool stellar and brown dwarfs, exo-Jupiters, and exo-Earths can be considered as being ~4.5, 3.5, and 3.0. The code can compute the models in plane parallel as well as in spherical geometry. In plane parallel approximation planetary (or stellar) radius is an irrelevant parameter, and only gravity needs to be specified. If sphericity is chosen, mass is specified in the input too such that the radius can be computed from mass and gravity. The full spherical atmospheric structure is then solved by use of the numerical method developed by Nordlund (1984) and further described in, for example, Jørgensen et al. (1992), Heiter & Eriksson (2006). It is a computationally fast, iterative method, giving better than 1% accuracy in the source function after only two to three iterations, and therefore it requires typically only a factor of five to ten times more CPU time than a corresponding plane parallel computation. The atmospheric structure of extended stars of low gravity, such as red giants, are substantially affected by sphericity, but for the type of objects discussed here, sphericity effects are in general negligible.

The input variables for the calculation of the thermodynamic quantities and the opacities are temperature, T (τ), and electron pressure, P_e (τ). The use of P_e as input variable in the MARCS-based part of the MSG scheme is somewhat modified to compensate for the extremely low abundance of free electrons in the coolest MSG models.

The dependent variables include the gas pressure, P_g, and radiative pressure, P_rad, as well as the integrated radiative and convective fluxes, F_rad and F_conv, and the convective velocity v_conv. Together these variables can be expressed in a highly nonlinear closed system of equations that includes the hydrostatic equilibrium, the radiative-convective flux-constancy, expressions for the gas pressure and radiative pressure, calculation of the convective flux from mixing-length theory, and calculation of the radiative flux in LTE, based on the Feautrier-method (Feautrier 1964), with boundary conditions as described in Nordlund (1984). The adopted formalism of the mixing length follows Henyey (1965) who besides the mixing length parameter α itself (in units of the scale height) also includes a parameter y which relates to the temperature distribution within the convective elements, and ν which relates to dissipation by the turbulent viscosity, and which in all our models are fixed to (α, y, ν) = (1.5, 0.076, 8), following the recommendation in Henyey (1965), but easily can be varied in the input file of the models. Additionally, also turbulent pressure P_t can be included as a dependent variable in the system of equations, but following the arguments in Gustafsson et al. (2008) all the models in the present MSG grid are computed with P_t =0. More details about the specific equations, the choices of methods, and the introduced approximations can be found in Gustafsson et al. (2008).

The code solves for the radiative energy balance (weighted sum over heating and cooling for each wavelength point) above a height where the convective flux is negligble, while below it solves for constancy of total flux, in both cases using linearisation of the relevant equations. The point of transition between the two methods is not critical, except it should occur where the convective flux is zero or insignificant. With the second order Feautrier method, constant radiative flux and radiative balance are numerically equivalent, given suitable definitions and centering of optical depth increments, but numerical precision is better for radiative equilibrium at small optical depths, and for constant flux at large optical depths. Due to the linearisation of the coupled equations (which while not including everything does include the most important dependencies), the temperature corrections decrease rapidly towards the end of iterations. We typically stop the iterations when the temperature corrections from one iteration to the next are smaller than 1 K in all layers, and tests have shown that this corresponds closely to remaining temperature errors. When clouds form in the model structure, it is the sum of the increased cloud opacities plus the reduced gas opacities that defines the energy balance (as described below in Section 2.3), and this usually requires additional iterations compared to cloud free models and often demands a stepwise introduction of the cloud particles in order to avoid too abrupt changes in the opacities from one iteration to the next.

In the overlap regions between the classical MARCS grid and the corresponding parts of the PHOENIX grid (Hauschildt et al. (1999), and in particular the following later versions) and the grids by Kurucz and collaborators (Kurucz 1979; Castelli & Kurucz 2003), there is a good general agreement between models from those three grids. This is obviously very reassuring given the large degree of independence between the three codes, concerning adopted methods, approximations, and (to some extent) adopted input data. We will comment in detail below about the level of agreement between the classical MARCS models and our MSG models in the overlap region between these two grids.

The extent of the atmosphere is defined by the input choice of the Rosseland optical depth, τ_ross, in the top (where P_g = 0) and bottom of the atmosphere. It should be noted that although τ_ross is a very poorly defined physical quantity in gasses where the opacity is dominated by line absorption, the adoption of τ_ross is only used as a means for the discretisation of the equation system, and it does not have any effect on the physical solution. As such, it can be thought of mainly as a convenient labelling of the individual discrete layers of the stratified atmosphere. The step length between the layers in any physical meaningful variable (e.g. gas pressure, temperature, height) can be adjusted in the input by adjusting the step length in units of τ_ross between any of the individual layers. In gaseous objects of our grid, we typically define the bottom of the atmosphere as the layer where log(τ_ross) = 2.5, while for solid objects (e.g. Earth-like exoplanets), we obviously define the bottom of the atmosphere as the planetary surface. In the bulk of the models presented in the following, we define the ‘top’ of the atmosphere as the layer where log(τr_oss) = –5, and we use a step length of Δlog(τ_ross) ranging between 0.1 and 0.15 in various regions of the atmosphere, resulting in typically 53 layers in our ‘standard’ models. For stellar atmospheres, it is a customary practice to call the ‘top’ of the atmosphere for the ‘surface’ of the star. In order to avoid confusion with the obvious use in planetary practice, with the ‘bottom’ of the atmosphere also being the ‘surface’ of the (solid) planets, we recommend, and try to be consistent with, using only the terms ‘top’ and ‘bottom’ of the atmosphere, avoiding the term ‘surface’ in the stellar-community meaning of the word as the ‘top’ of the atmosphere.

2.2 GGchem

We used the publicly available thermo-chemical equilibrium code GGchem of Woitke et al. (2018) to calculate the concentrations of all neutral and single ionised atoms, electrons, molecules and molecular ions. In all models of the presented MSG grid, higher than single ionised atoms are irrelevant, but in order to be consistent with the classical MARCS grid, and in order to compute exoplanets and their irradiating host star self-consistently, there is an input option to include higher ionisation from the classical MARCS code (subroutine JON) if desired. GGchem is based on the principle of minimisation of the total Gibbs free energy, including condensates, and applicable in a wide temperature range between 100 K and 6000 K. However, in this paper we use GGchem solely for the purpose to calculate the gas phase particle densities based on temperature, pressure and the gas phase elemental abundances. In MSG the solid state condensates are obtained from subsequent calls to the more advanced StaticWeather computations of the dust grains (see next sub-section). Hence, the gas phase elemental abundances change layer by layer though the MSG iterations and convergence, keeping everything self-consistent.

The thermo-chemical data has been collected from different sources, in particular from the JANAF thermochemical tables (Chase 1996, 1998) and from Barklem & Collet (2016). The data have been carefully compared to other data sources by Worters et al. (2017). We consider 22 elements (H, He, Li, Be, C, N, O, F, Na, Mg, Al, Si, P, S, Cl, K, Ca, Ti, V, Cr, Fe, Zr). Each element can be present in the form of neutral atoms, or singly charged positive or negative ions. According to this selection of elements, GGchem finds additional 543 molecules, molecular ions and cations in its database. GGchem has been benchmarked down to 400 K for a pure neutral gas phase chemistry against TEA (Blecic et al. 2016), and down to 500 K for chemistry, ionisation and condensation against Sharp & Huebner (1990). Below these temperatures, the other codes start to have problems to converge. GGchem uses a hierarchical approach of pre-iterations to stepwise increase the number of elements, combined with quadruple-precision arithmetics, to solve the thermo-chemical equilibrium down to 100 K. GGchem is written in ultra-fast Fortran-90 and needs about 0.003 CPU-sec per call in this setup.

2.3 StaticWeather

The kinetic cloud formation model StaticWeather (also known as DRIFT) was developed by Woitke & Helling (2003, 2004). The model calculates the rate of nucleation of seed particles from the gas phase by modified classical nucleation theory (Gail 1984), and considers an explicit set of surface chemical reactions to calculate the growth rates of all cloud particles by moleculesurface collisions, and the reverse rates, the sublimation rates. Each pair of growth and sublimation rates is related to each other via the supersaturation ratio S of the considered solid material in the atmosphere. In order to calculate the molecular densities required for the effective growth rates and the supersaturation ratios, the GGchem code is used internally (see Sect. 2.2). The StaticWeather code uses the dust moment method of Gail & Sedlmayr (1988) to consider a cloud particle size distribution function that changes with height in the atmosphere. The size-dependent gravitational settling of the cloud particles is part of the modelling, which leads to a downwards transport of condensable elements.

The element abundances in the gas phase are hence consumed by nucleation and growth, and replenished by sublimation and convective mixing. For the replenishment by convective mixing, the StaticWeather code uses a simple approach based on the mass exchange timescales derived from the 3D hydrodynamical models for M-star atmospheres by Ludwig et al. (2002). At each atmospheric layer, the gas with cloud particles is exchanged by the cloud-free gas from the deep atmosphere, on a timescale that depends on atmospheric height.

In praxis the inclusion of clouds into the self-consistent modelling of the MSG atmospheric structure is done by starting with computing (or reading in) a cloud-free model structure of given elemental abundances (that are the same in all layers) as specified in the input file. Then StaticWeather computes the amount of gas that will condense into solids for each layer of this structure. In next iteration the opacities from the computed abundances of the solids and from the corresponding (reduced) gas species are computed and added up, and the radiative transfer is calculated. A new structure is then computed, and so on until a self-consistent convergence is reached. After full convergence, the input relative elemental abundances of the gas phase are only conserved in the regions outside the cloud formation area. Inside the cloudy region, the gas phase elemental abundances of those elements that participate in the mineral cloud particles are reduced, sometimes by large factors. The total elemental abundances present in the solid and the gaseous states in each cloudy layer do not in general sum up to the input elemental abundances, but depends on the local condensation, sublimation, and mixing time scales. The emerging spectrum of such fully self-consistent models can be substantially different from non-self-consistent modelling.

Helling & Woitke (2006) have generalised the StaticWeather model to mixed condensates, assuming that all stable solid materials grow on the same cloud particle surfaces. The model is anchored in experimental laboratory data, all the way from the chemical reactions it uses to the optical properties of the condensates it considers. A hand-picked set of about 20 condensates is selected, including several magnesium-silicates, solid iron, a few high-temperature Al-Ca-Ti oxides, and some simple metal oxides and sulphides. In atmospheres that are cool at the top and hot at the bottom, the nucleation kick-starts the formation of “dirty” cloud particles materials high in the atmospheres, but as these particles settle down into warmer regions, they stepwise “purify”, until only the most stable condensates remain – mineral clouds – before even these materials sublimate at the cloud base. Helling et al. (2008a) have compared the StaticWeather model to other cloud formation models. Helling et al. (2008b) and collaborators have included the StaticWeather into PHOENIX stellar atmosphere models of Hauschildt et al. (1999) for brown dwarfs, as well as into gas-giant atmospheres (e.g. Helling et al. 2021; Lee et al. 2015; Lee et al. 2016; Juncher et al. 2017; Samra et al. 2022), including 3D GCM models (Helling et al. 2016).

2.4 Computational speed

The computational speed of the code scales to a first approximation with the number of adopted frequencies, n_λ, in the radiative transfer computation, times the number of layers in the atmosphere, nτ, times the number nµ of individual directions in the radiation field adopted in the computation of the mean intensity. In order to self-consistently include clouds into the model computation, the MSG code allows for varying relative elemental abundances in the gas phase chemistry from layer to layer and from iteration to iteration, which obviously adds considerable time to the computation of the opacities (as does the fact that the MSG models include many more individual opacity sources than the classical MARCS models), but a combination of increased core memory and computer capacity and improved coding has allowed the total time per iteration to remain approximately the same as 20 years ago (for a corresponding classical MARCS model). This is very reassuring because it allows one to use the MSG models not only for comparison of models over an extensive range of parameter space but also as a test base for self-consistent inclusion of advanced physics and chemistry into dynamically advanced models such as modern 3D GCM models. With typically 100 000 λ-points, 50 depth-layers and 6 directions (in the plane parallel approximation, and about 5 times more in the spherical case), and 10 iterations for full convergence (depending strongly on the applied input model for the first iteration), a typical cloud-free model convergence time when executing at a single core of our present hpc computing system (Intel(R) Xenon(R) CPU E5-2680 v4 @2.40 GHz) is of the order of 15 minutes.

3 Consistency between non-irradiated, cloud-free MSG models and the classical MARCS models

As described above, the classical MARCS grid was limited to models of T_eff = 2500 K and warmer, among other things due to the lack of sufficient input data for cooler models. Such data have now become available due to the huge work on molecular line data as computed and collected in the ExoMol database (Tennyson et al. 2020), due to the extensive development in cloud formation data and theory (Helling & Woitke 2006), and the improved chemical equilibrium computations (Woitke et al. 2018), as well as due to extensively updated atomic databases (Grimm et al. 2021) and updated data in general. In this section we will compare the cloud-free, non-irradiated, plane parallel, chemical equilibrium, solar metallicity MSG models with the corresponding classical MARCS models available at the Uppsala-MARCS-webpage³ and described in Gustafsson et al. (2008) and the DRIFT-MARCS models described in Juncher et al. (2017).

The classical MARCS models include the line opacity from 16 molecules (HCN, H₂O C₂, CH, CN, CO, CaH, FeH, MgH, NH, OH, SiH, SiO, TiO, VO, ZrO) from various sources that were available at that time, as listed in the 2008 grid paper, plus atomic lines from the VALD database. The DRIFT-MARCS models from Juncher et al. (2017) include the same 16 molecules (but not in all cases from the same sources) plus 6 additional molecules (CO₂, TiH, CrH, NO, LiH, H₂). Cloud formation and cloud opacities were included in the Juncher et al. (2017) models, but in general atomic line data were not included. The standard version of the present MSG model grid includes the same 22 molecules as in Juncher et al. (2017) plus additionally 28 molecules (SO₂, AlCl, PO, NaCl, SiS, LiF, LiCl, PH₃, KCl, HCl, CP, NaH, AlF, BeH, MgF, CaF, CH₄, H₂CO, NH₃, AlO, PN, CS, KF, PS, SN, AlH, NaF, HS). All included molecules are now from the ExoMol database. References to the original work behind the individual line lists adopted from the ExoMol database are given in Table 1 below together with the original sources for the individual atomic line list date adopted from the DACE database. The data for continuum and collision induced opacities are the same as used in the classical MARCS models, and listed in Table 1 of Gustafsson et al. (2008). Some of the classical MARCS and the DRIFT-MARCS models included two more carbon rich molecules (C₂ H₂ and C₃ ), that are not generally included here. They only affect carbon rich atmospheres and not the solar metallicity models that are the focus of the present paper. They and other molecules will be included in forthcoming versions of the MSG model grid where relevant.

The upper panel of Fig. 1 shows that for models of T_eff = 2500 K there is a good general agreement between the MSG, the classical MARCS, and the 2017 models, computed based on their adopted input as mentioned above. For higher values of T_eff the agreement can be expected to be even closer, since the additional data added in the MSG grid computations will have a smaller importance than for the T_eff = 2500 K models. This is very reassuring, because it shows that the new grid of cloud- free and non-irradiated MSG models and the classical MARCS models are consistent and in general agreement with one another at the overlap region between the two grids.

The lower panel of Fig. 1 zooms in on the uppermost layers of the same models as in the upper panel of the figure, to quantify clearly the region where the various models differ the most. It shows that the models in the most extreme cases differ from one another with up to 70 K. The classical MARCS models fall close to the middle of the range.

We consider the presently best MSG models to be those that are based on line data from the 50 molecules of the ExoMol database mentioned above plus from 20 different neutral atoms available in the DACE database (Grimm et al. 2021) as listed in Table 1, plus of course updated chemical equilibrium computations from GGchem and other updates as well. For T_eff = 2500 K, Fig. 1 shows that this model is only ~20 K cooler in the upper layers than the shown classical MARCS model (for a given value of P_gas). We note that most of the classical MARCS grid is computed for models with microturbulence of ξ=2 km/s as well as ξ=5 km/s, while the T_eff = 2500 K model was computed for ξ=2 km/s only, and the present MSG grid is computed for a fixed value of ξ=3 km/s. Typically, cool models get ~3 K hotter per increase of 1 km/s in ξ, so we conclude that the preferred T_eff =2500 K MSG model is approximately 25 K cooler than the corresponding classical MARCS model (in the upper layers and for given value of P_gas). The choice of microturbulence determines the broadness of the Gaussian line profiles used in the computation of the molecular opacities, and it is this effect that accounts for the difference between the ξ=2 km/s and ξ=5 km/s models. We will discuss in detail in a separate section below the effect of the important choice of line profile in general for the resulting structure of cool-star and exoplanetary models.

In order to understand what effect various updates have played for the models on the boundary to the classical MARCS grid, we show in Fig. 1 also a number of T_eff =2500 K MSG models computed with a variety of chosen molecular and atomic input data. This will hopefully add some hints about where future improvements in exoplanetary modelling in general could be most beneficial.

The model labeled ‘G2008MARCS 2km/s’ is the classical MARCS model (full drawn red line). The group of three models in Fig. 1 that are a few degrees cooler than the classical MARCS model are based on line opacities from respectively the set of 22 molecules that were also adopted in the 2017 grid (but here all take from the ExoMol database), all the 50 molecules taken from the ExoMol database, and all these 50 plus the 20 atoms from the DACE database (full drawn black line). This latter is our preferred combination, used as “standard” for the MSG grid. All these three models differ less than ∼10 K from one another in the uppermost layers of the atmosphere, and it is therefore concluded that the 28 newest molecules added from the ExoMol database only have a marginal effect (of slightly more than 5 K) on the model structure (of a T_eff =2500 K model) and the addition of atomic lines has an even smaller effect. This latter conclusion (about the atomic lines) has always been implicitly assumed in our older cool MARCS models (although not in the 2008 grid), but has never been thoroughly and quantitatively tested before the results shown in Fig. 1. We will return to it in a separate section below.

Of the group of three models in Fig. 1 that are all approximately 10 to 15 K warmer than the classical model for a given P_g, the two coolest are computed with respectively the 16 and 22 molecules from (approximately) the line lists adopted in the 2008 grid and in the 2017 grid, while the warmest of the three is computed based on the 16 molecules from the classical grid augmented with the 20 atoms from the DACE database that we use in the MSG grid (labeled ’MSG 16mol+20atoms’). The difference between these three models is very small, and we conclude that the effect of the difference in the opacities included in the 2008 and 2017 grids were very small for these values of T_eff.

Finally, the group of the three warmest models shown in Fig. 1 are almost 50 K warmer than the published classical MARCS model. The two (almost) identical of these models are based on the 16 molecules in the classical MARCS grid or the 22 from the 2017 grid, respectively, and both with the water line opacity from the BT2 ExoMol list (Barber et al. 2006), while the third of these models are with the Prokazatel H₂O line list from ExoMol (Polyansky et al. 2018). We conclude from these three models that H₂O is likely to be the most important individual molecule for the opacity at these temperatures (and compositions), and that the difference between choosing the BT2 or the Prokazatel ExoMol water line list give rise to only minor changes.

We summarize from the comparisons in Fig. 1 described above that the coolest end of the classical 2008 MARCS grid models and the cloud-free, non-irradiated MSG models in the present grid are in very good agreement with one another. Our T_eff = 2500 K, log(g)=4.5, solar composition, cloud-free, nonirradiated MSG models with our “standard” input choices are ~25 K cooler in the upper layers than the corresponding classical MARCS models at a given value of the gas pressure, and has a ∼0.1dex lower gas pressure at the formal top layer of log(τ_ross) = –5.0. Further we conclude that the effect of the additional 28 molecules that were added in the MSG computations, but not included in the DRIFT-MARCS models, is marginal (at T_eff = 2500 K and solar composition), and that the effect of atomic line opacities only have a negligible effect on the model structure at T_eff = 2500 K (and cooler). We therefore conclude, to our great satisfaction, that the conclusions of the many tests and analyses of observational data that have been done based on the classical MARCS grid will remain unchanged by the introduction of the MSG grid; even for the coolest models of the published MARCS grid. It is only for the models of lower effective temperatures there are big differences between the old (unpublished) MARCS models and the new MSG models. We do, however, caution that the large amount of different molecular line opacity sources individually gives rise to heating as well as to cooling of the atmospheric layers, and that the effect of different combinations of the available line computations demonstrated in Fig. 1 combined give rise to what might or might not be seen as an uncertainty in the results of ∆T~70 K and ∆log(𝑔)~0.2 dex in the upper layers of the models, and corresponding spectral interpretations. We have chosen to hopefully lower this uncertainty by consistently choosing all molecular line data from the Exo- Mol database and all atomic line data from the DACE database as listed in Table 1. The imposed uncertainty introduced by the choice of line profile applied to these data will be discussed in a separate section below.

Fig. 1 Comparison of the classical 2008 MARCS model at the lower boundary of the 2008 grid (labeled ‘G2008-MARCS’) with models computed with the MSG code and based on various levels of completeness and choices of input data as described in the text. Models labeled ‘MSG 22…mol..’ are based on molecular opacities from the 22 molecules included in the DRIFT-MARCS computations, while models labeled ‘MSG 16mol…’ are based on opacities from the 16 molecules included in the classical MARCS grid. Upper panel show the T–Pg structure of the full models from log(τross) = –5.0 to +2.5 while the lower panel zooms in on the upper part of the atmosphere.

4 The cooler MSG models

To a first order approximation, the physical temperature at a given atmospheric optical depth decreases almost linearly with the effective temperature. With decreasing physical temperature, an increasing number of atoms will combine into molecules. Since molecules can have electronic transitions (as atoms) as well as a cascade of associated vibrational and rotational transitions, the absorption coefficient per atom will generally increase substantially when neutral atoms combine into diatomic and polyatomic molecules. Cooling the atmosphere will therefore usually result in an increasing opacity per gram of stellar material. As a response to the increased opacity, the atmosphere will usually expand, and the gas pressure at any given optical depth will decrease. We note that τ_ross only has a strict physical meaning when the gas opacity is dominated by continuum sources, and therefore the “surface” or “top” of the atmosphere, defined here as where τ_ross = 10^–5, is mainly to be considered as a computational label of the atmospheric top, and the more meaningful comparison of individual models is for example what the physical temperature is at a place of a given gas pressure.

Since the input data to the classical MARCS models were known to lack opacity sources from species that are observed in stars cooler than T_eff = 2500 K, it could be suspected that the models would be too compact at the lowest effective temperatures, due to the missing opacities, and the result would be not only that the computed spectra would be absent of the missing opacity sources, but the computed spectra would also be correspondingly too strong in the other spectral features because of too high gas pressure (and hence partial pressures) in the models. This would translate into the model spectra predicting lower than real abundances when used for interpreting observed spectra (at least to a first approximation), or alternatively giving erroneous estimates of the effective temperatures. A now classical and extreme example of this effect was the analysis of the hydrogen abundance in giant stars (Jørgensen 1989), and it would have been in effect for the classical MARCS models too, if they had been squeezed below T_eff ~ 2500 K.

Fig. 2 illustrates both how the physical temperature scales with the effective temperature as well as quantify the effect of insufficient inclusion of molecular opacities in the classical grid for temperatures below T_eff =2500 K. For T_eff =2500 K the T – P_g and T – τ structures of the classical MARCS models (with 22 molecules in the opacity) and the present models (with 50 molecules in the opacity) are seen to be virtually identical. This is of course reassuring, since it shows that the classical grid as published (i.e. down to T_eff =2500 K) is not affected by lacking (now known) molecular opacities. For T_eff=2000 K the difference between the models with respectively 22 and 50 molecules in the opacity are noticeable, and for T_eff=1500 K the difference amounts to more than 100 K all the way from τ_ross=1 and upwards in the atmosphere, i.e. basically throughout the whole spectrum forming region.

Also derived structures such as the physical size of the atmosphere are for T_eff = 2500 K almost identical in the classical MARCS and the MSG grid. For T_eff = 2000 K, however, although the T – τ structure of the old classical MARCS and the new MSG models are quite similar, the gas pressure at the top of the atmosphere (defined as τ_ross = 10^–5) of the MSG models are ~30% lower than in the ones simulating the classical MARCS models (for the reasons described above, due to the larger molecular opacities), giving rise to a quite different abundance determination when compared to observed data. At T_eff = 1500 K the MSG models are completely different from the classical MARCS ones, with the gas pressure in the top layers (τ_ross ~ 10^–5) reduced by two orders of magnitude in the MSG models compared to models simulating the classical MARCS ones, and the physical size of the atmosphere (from log(τ_ross) = 2.5 to –5) increased with 40% (from 140 to 190 km for log(g)= 4.5 models). Below T_eff = 1500 K, it was no longer possible (nor meaningful) to attempt to converge the models with input data and computational methods corresponding to what was used in the classical MARCS models.

In general, one sees from Fig. 2 that the temperature (and from Fig. 3 the electron pressure as well) of the MSG models decrease with a larger and larger amount compared to the classical MARCS models with decreasing effective temperature, resulting in a slight increase in the physical temperature for deeper atmospheric layers with values of log(τ_ross) larger than 0 (best understandable from Fig. 2 lower panel for the T_eff =1500 K models). Taken isolated, these changes would result in a decreased continuum opacity and an increased molecular line (and pseudo-continuum) opacity, which isolated seen would result in an increased intensity in the predicted spectral lines and hence a lower abundance estimated when compared to observed spectra. As usual in stellar atmosphere and self-consistent planetary atmosphere theory, it is, however, not possible to give a general scaling, because the problem is non-linear.

Fig. 2 Upper panel shows the T–P𝑔 structure and lower panel the T- τ structure from log(τross) = +2.5 to –5.0 for models with Teff from 2500 K to 500 K. For each Teff, we show a classical MARCS model with 22 molecules in the opacity calculation as well as the models from the present grid (i.e. 50 molecules plus atoms included in the opacity), except for the Teff =1000 K and Teff =500 K models where it has not been possible to converge the models with only 22 molecules.

Fig. 3 Electron pressure versus optical depth for models with Teff from 2500 K (upper curves) to 500 K (lower curve), computed from log(τross) = +2.5 to –5.0. For the warmer of the models are shown both the results based on including only the original 22 molecular opacities in the equilibrium computation as well as the full set of the 50 molecules, as discussed in the text. For the coolest models only those with opacities from all the 50 molecular opacities adopted in the MSG models could converge.

5 Three “standard challenges” in self-consistent cool star and planet atmosphere modelling

We will discuss here three “standard challenges” in deriving the final gas-pressure versus temperature structure in atmospheric computations for a sufficiently wide range of input parameters to make a general-purpose grid. The first is somewhat specific to the way classical MARCS models were computed, while the two others are of more general character.

5.1 The choice of independent variables

The first challenge is to choose “the right” independent variables in the convergence scheme to make the models converge efficient and accurately for a large temperature range. The original MARCS code (Gustafsson et al. 1975) introduced temperature and electron pressure as the independent variables in the iteration scheme towards flux-constant atmospheric model structures. The same choice is done in for example Gray’s text book on stellar atmospheres (Gray 2022). This choice is very convenient and valuable for solar type and warmer stars and even for giant stars in globular clusters, which was an early focus of the application of the MARCS models. When the stars become considerably colder, or we enter the sub-stellar and planetary regime, however, electrons are sparse and the iteration on electron pressure becomes a challenge.

At high temperatures the opacity is dominated by continuum sources, including electron scattering and H^– b-f and f-f transitions. At somewhat lower temperatures, more and more ions recombine with electrons and become neutral atoms, and the electron pressure (and therefore also H^-) decreases. At even lower temperatures the atoms combine to form molecules, and in chemical equilibrium the electron pressure rapidly approaches zero. Fig. 3 shows how the electron pressure in the uppermost layers drop 20 orders of magnitude from 10^–7 dyn/cm² to 10^–27 dyn/cm² when lowering the effective temperature from T_eff = 2500K to T_eff = 1000K (for models of log(g) = 4.5 and solar metallicity). While this is not a physical problem, it creates a numerical challenge in the computation of the derivative of the opacity as function of electron pressure, δO/δPe, and in derived thermodynamical quantities. Since the convergence scheme in MARCS builds on δO/δT and δO/δPe, the MARCS code traditionally has problems converging at low temperatures.

However, for the low temperatures, one can take advantage of the fact that almost all non-noble gas elements are bound in molecules and the small number of atoms that are not bound in molecules are basically all neutral. Hence the gas pressure is to a very good approximation the sum of the partial pressures of molecules and noble-gas atoms only (ie no electrons and no ions). The small number of free electrons can therefore be computed from a simplified Saha equation that sums over just a few atoms. Adopting the nomenclature from Gray (2022), one can write the electron pressure as $$P_e={\displaystyle \sum _{j}k}T\cdot N_{ej}={\displaystyle \sum k}T\left(N_j\frac{{\text{Φ}}_j(T)/P_e}{1+{\text{Φ}}_j(T)/P_e}\right),$$(1)

with $${\text{Φ}}_j(T)=1.20\cdot {10}^9\left(u_{1j}/u_{0j}\right){\theta }^{-5/2}{10}^{-\theta I_j},$$(2)

where N_j is the number density of element j, N_ej is the number density of electrons coming from element j, u_0j and u_1j are the partition functions of respectively the neutral and first ionised state of element j, I_j is the ionisation potential from the neutral to the first ionised state of element j (in units of electron volt, eV). θ = 5040/T is the “temperature” in units of per eV when T is measured in K.

For temperatures where the electron pressure is very small compared to the gas pressure, the correct temperature versus gas pressure model structure is not dependent on the actual value of P_e used for the convergence, as long as the correct electron pressure is used for the opacity computations. In order of convenience and to make the MSG code efficient and accurate for high as well as low temperatures, we therefore kept the convergence method from the classical MARCS models, but we allowed the MSG convergence scheme at low temperatures to be a function of P_e, P_ae = f (P_e) instead of the actual electron pressure itself, with f (Pe) computed from Eq. (1) as described above.

Once the gas pressure has been related to the artificial electron pressure, P_ae, we used GGchem to calculate the true electron pressure, P_e, which is then subsequently used for all calculations of physical quantities (such as the continuum opacities, the molecular opacities, etc.) to estimate δO(Pe)/δP_ae. Since Φ_j(T) (and hence also P_ae) is a linear function of ${10}^{-\theta I_j}$ (Eq. (2)), we have essentially changed the code to iterate over gas pressure instead of electron pressure, where the real electron pressure P_e has now become a dependent variable instead of the artificial electron pressure P_ae.

At the very lowest temperatures, this same reasoning can be applied to allow for further efficiency in the convergence by artificially lowering the ionisation potential, I_j, of the relevant atoms (which for solar composition is potasium only), and in this way, we secured efficient convergence of the model structure for all models where no physical layer reaches a temperature close to the convergence criterium of GGchem at around 100 K. For non-irradiated models, this corresponds to T_eff≈ 300 K, and for irradiated models, it corresponds to a lower value of T_eff. Once the low-temperature models have converged, the thermodynamic quantities (such as heat capacities) are now self-consistently calculated using the NASA polynomials from data/nasa9.dat (McBride et al. 2002). This assures that the correct adiabatic index (and hence boundary between layers in radiative and convective equilibrium) and the correct gas pressure etc., are computed at low temperatures while still allowing for high accuracy by including both single and double ionised atoms (and only few molecules) at the highest temperatures, as in the original MARCS code. In this way the flexibility, computational speed, and high accuracy in the convergence to a final T–P_g structure is assured for a sufficiently wide range of T_eff to compute the exo- planetary system and its host star together with the same code and input data.

Several tests have assured us that for the models that could converge with both options in handling of P_e the resulting structures were identical, as is illustrated in Fig. 4 where the computed electron pressure as function of τ_Ross is seen to be 7 orders of magnitude higher for P_ae than P_e (upper panel) in the “old convergence scheme” (models labeled “…p1…”) than the new (which is further seen to be independent of the assumed value of the ionisation potential of potasium, as it should), but yet, as expected and explained in the text above, leading to exactly the same T – P_g model structure (as is seen in the lower panel).

5.2 The number of transitions to include in the opacity computation

A second “standard assumption” one needs to take when computing a self-consistent atmospheric model, is how many transitions to include in the opacity computation (unless one just adopt an opacity table accessibly on the web, but then somebody else have taken the decision instead). Choosing too few lines will make the outgoing radiation escape between the lines, with the effect that the atmospheric structure becomes too compact. A too compact model will not have too much effect on the observed spectrum of the dominating opacity source(s) that is (are) computed wrong. This is because the computed atmospheric structure will balance itself such that the hydrostatic equilibrium, flux constancy, and energy balance is preserved and hence also the emerging flux from the top of the atmosphere (the spectrum). However, the opacity from the sources (mainly molecules because they for cool objects carry most of the opacity) that are not dominant will be wrong, often by a substantial factor, and the sensitivity to for example abundance analysis will therefore often be sensitive to the choice of number of lines included in the opacity computation of the model structure used for the analysis. On the other hand choosing too many lines can have the opposite effect. The problem is the same concerning the choice of line profiles that we will return to below. One could be tempted to think that one can never choose too many lines, because one could just choose “all of the lines” and be on the safe side. The problem here is that there is no such concept as “all the lines”. The following will illuminate why such a concept does not exist.

In modern quantum mechanics, transition moments (“intensities of electronic, vibrational or rotational transitions”) are computed as integrals over the corresponding wavefunctions, but the factors that define how many lines a given electronic or vibrational transition ends up producing relies on the quantum numbers associated to the eigenstates as well as the definition of the energy of the highest energy level in the energy potential, both of which relies on quantum numbers from classical mechanics derived from the harmonic oscillator theory or some more or less empirical modifications of it. The same is true if one attempts to evaluate the completeness of an observed line list. The highest eigenstates are most affected by the incorrectness of the harmonic oscillator theory, and one will often argue that since the highest eigenstates are least populated in a Boltzmann distribution, then it is probably not too important to compute them accurately, or one could just disregard them. In the pure harmonic oscillator theory, there are infinitely many lines originating from the uppermost states in the potential, each approaching infinitely small intensity. The “unimportance-argument” is probably true concerning the total integrated absorption coefficient, but not necessarily concerning how much energy escapes between the included lines, which affect the model structure in a complex and non-linear way which is not straight forwards to predict. The only realistic way forwards may be a rather cumbersome comparison between high-resolution observed spectra (which has to include transitions from high and low excitation energies in opacity-dominant as well as opacity un-important sources) of high-temperature objects (stars and planets) of various types compared to theoretical results based on different attempts to do something “reasonable”. We tested and described in detail in Jørgensen et al. (2001) the numerical effect on the model structure of successively including from a few million water lines to a few billion water lines, and we described in detail in Popovas & Jørgensen (2016) various choices one can make for including or excluding lines from excitation energies above a certain level, and calculated as an example what could be a reasonable choice for H₂, which was later tested and found in good agreement with laboratory experiments for H₂ (Retter et al. 2020).

Fig. 4 Models with Teff = 1000 K computed based on the ‘old’ way (labeled ‘…Pe1K4’) of estimating the electron pressure as well as the new way with Ij equal to 2.339 eV (labeled ‘…Pe2K2’) and 4.339 eV (labeled ‘…Pe2K4’), respectively. The upper and lower panels respectively show, as a function of log(τross), Pe and Pg.

5.3 The choice of line profiles in the opacity computation

A third “standard assumption” one needs to take when computing a self-consistent atmospheric model, is the selection of line profiles in the opacity (and spectrum) computation. This challenge is quite analogous to the “second standard assumption” described above since it also concerns how much flux one chooses to allow the object to escape between the lines (in the opacity).

In all the computed MARCS atmospheric models so far, the molecular line profiles were treated as Gaussians; that is, the wings of the full Voigt profile were ignored. If one assumes that a line list from a given molecule (or the sum of line lists of all the molecules included in the computation of the model structure) includes, for example, a billion lines from somewhere in the visual region to λ equal to say 10 microns, then we have 10 000 lines per angstrom. With a microturbulence-dominated Gaussian line width of say ξ = 3 km/s, then the number of lines per Gaussian half-width that we have is of the order of 1000 lines. Intuitively, one could guess that the sum of the opacity of such a huge number of even, relatively weak line cores would always dominate over the wings of any even relatively strong single line. This would intuitively justify the exclusion of the line wings, but it has never really been possible to test this assumption.

One reason for excluding the wings is that they are temperature and pressure dependent, while the Gauss profile is temperature dependent only. One can therefore save substantial computing time (and core memory) by limiting the computation to Gaussian profiles. Another reason for doing it was of more practical character, namely that the necessary data for computing the full Voigt profile for the huge number of lines were unavailable. This problem was overcome with the creation of the ExoMol database, and the test of the ‘correct’ line profile is therefore now in principle possible. Opacity computations based on the full list of relevant transitions in the ExoMol database has been done by us for the present work assuming Gaussian line profiles with ξ = 3 km/s, and by Chubb et al. (2021) for Voigt profiles with ξ = 0 km/s. The use of the two set of line lists in the opacity computation can lead to differences in the model structures larger than the differences in the model structures based on Gaussian profiles with respectively 2km/s and 5km/s that in the classical MARCS models were implicitly assumed to encompass the uncertainty due to line profiles. It is not obvious which of the two existing approaches (excluding the Lorenz broadening or excluding the microturbulence) are most severe (most likely none of them are really “correct”) and which values of the microturbulence are most correct for the sub-stellar objects, only that the difference between the two approaches seems to give rise to larger structural differences than the choice of ξ=2 and 5km/s do. A proper solution may require an observational estimate of the line widths as a function of sub-stellar parameters, such as was fruitfully done through decades of observational studies for stars.

Fig. 5 illuminates the complexity of the problem. For low temperatures (here 1000K) and low gas pressures (here 10⁻⁸ bar), the Lorenz broadening gives rise to a larger blanketing at all wavelengths, whereas for high temperature and high pressure (here 10⁻³ bar) the effect of the microturbulence broadening of the Gauss profile on the blanketing strongly exceeds the effect of the Voigt wings at wavelengths below 1 µm. For higher temperatures (here 2900K) the effect is opposite, but also much smaller. Qualitatively, this shows us that we will have blanketing from the microturbulence line broadening mainly in the deeper layers (short-wavelength Gaussian blanketing at high pressures) and therefore cooling of most of the atmospheric structure, while the effect of the Lorentz wings in the Voigt profile will dominate at lower pressures and therefore block more light in the upper layers and backwarm the bulk of the lower atmosphere. Generally, we therefore expect that excluding the microturbulence cools the models (the lines or line cores get more narrow, and more radiation therefore escapes between the lines), while substituting the Gauss profiles with Voigt profiles heats the models (there are fewer regions with ultra-low opacities, so more radiation is blocked). Experimenting with the two sets of opacities showed us that the effect is larger for the cooler models than for the warmer models. One has, however, to be cautious to generalize this conclusion because the effect is complex and non-linear, but it could be because fewer levels of the individual molecules are populated at lower temperatures, but on the other hand there are also more different molecules present. We also noted that the bulk of the effect generally moves downwards in the models when one goes from warmer to colder models and that for the cooler models introducing Voigt profiles heats the inner structure and removing microturbulence cools the inner structure, which most likely reflects a temperature dependence of the degree of backwarming due to the pseudo-continuum of the veil of weaker lines. It could therefore be that the effect of lacking modelling of the Lorenz wings might have had a somewhat larger effect in the upper layers of existing models than we have usually assumed.

Fig. 5 Absorption coefficient of the OS lines as a function of wave number, computed from the ExoMol line list under the assumption of Gaussian molecular line profiles (blue) and Voigt profiles (light brown), respectively, and for two different values of temperature (T = 1000 K and 2900 K) and two different values of Pg (10−8 and 10−3 dyn/cm2).

6 The sampling method and the changing partial pressures

In order to transform the line lists into molecular line absorption coefficients and further into a total monochromatic molecular absorption coefficient (e.g. in units of cm²/molecule = cm/molecule/cm⁻¹) that can be used for the radiative transfer calculations, some kind of statistical approach to the resulting opacity (e.g. in units of cm²/g_*, where g_* is grams of stellar material) has to be adopted, which here (as in all MARCS models computed after ∼1990) is chosen to be the Opacity Sampling scheme (as described for example in Jørgensen 1992). We typically use a sampling density of R = λ/∆λ = 15 000 which is sufficient for medium resolution spectrum calculations, and then sample each 10^th point for the radiative transfer calculation. In the OS approximation the line frequencies for the transfer computation has to be chosen randomly, but this is well satisfied with a fixed step length as long as the line positions can be regarded as random relative to the chosen OS frequencies.

The chosen sampling is a good compromise between speed of computation and accuracy of the final model structure, and it corresponds to solving the radiative transfer problem in ∼10 000 frequency points. The computing time scales approximately linearly with the number of frequency points. Such sampling gives essentially the same result as the ODF method (which has some similarity to the often used correlated K-method; Goody et al. 1989) that was used and described in the oldest versions of the MARCS code (Gustafsson et al. 1975). The OS method does not have the problem of lacking frequency correlation between individual species that the ODF method has (Saxner & Gustafsson 1984) and is therefore more flexible and considerably faster when many opacity species are included. Moreover, it is also relatively easy to estimate how many frequency points are needed to reach a correct solution (Helling & Jørgensen 1998). It further has the advantage that the same OS’s can be used consistently in the model iteration and the spectrum computation.

Once the individual monochromatic absorption coefficients are calculated and tabulated, it is then simple to compute the individual molecular opacities in the sampling points by multiplying with the individual partial pressures (for example in units of number density of species per number density unit of stellar mixture, ρ_molecule/ρ_*, which when applying the ideal gas equation as the equation of state is identical to the molecular partial pressures divided with the stellar gas pressure) such that the opacities of the individual species becomes in units of per stellar mixture (for example cm² per gram of stellar material). These opacities can then be added linearly since they are all in units of per gram stellar material, and the computing time of the stellar opacity therefore scales only approximately linearly with the number of opacity sources included.

Figure 6 shows the total opacity as function of gas pressure in units of cm² per gram of stellar material. The figure illustrates well the reason behind the changes in the P_g − τ structure described above and the changing T − τ structure shown in Fig. 2. It is seen that for the lower limit (T_eff = 2500 K) of the classical 2008 grid, the inclusion of the many new molecules from Exo- Mol compared to the only 16 molecules included in the classical grid, does not change the stellar opacity at any depth of the atmosphere, and hence the model structure does not change either. For T_eff = 2000 K the opacity changes with approximately a factor of 2 (0.3 dex), while for T_eff = 1500 K the new opacities are 2 orders of magnitude larger than the old ones. For T_eff = 1000 K we must expect the differences to be even larger, and it is somehow reassuring that it was not possible, as described in connection with Fig. 2, to converge models with such unrealistic input data.

It should be stressed that the opacities shown in Fig. 6 are not physical entities that as such could be measured at any place in the atmosphere. They are the Rosseland-mean opacities, arising from the sum of the individual monochromatic absorption coefficients multiplied with the respective partial pressures for the relevant molecules, which only has a physical meaning for non-monochromatic absorption coefficients. They are therefore only illustrations of the approximate reason for the changing model structure as function of effective temperature. While the real opacity that goes into the model structure computation is monochromatic (in the OS statistical sense), the Rosseland opacity (sometimes a bit misleading called the Rosseland mean absorption coefficient, κ_ross), is one over the integral of one over the monochromatic opacity times the planck function. It is a bit like the way we talk about thin and thick cloud cover on Earth, meaning how the opacity is qualitatively sensed with our eyes, through the sensitivity function of our eyes multiplied with the intensity distribution of the sun light. It is therefore also to be understood that the two orders of magnitude decrease in the gas pressure in the top of the atmosphere mentioned above, because of the 100 times increased value of the opacity seen in Fig. 6 for the new line lists, is influenced by the definition of “top of the atmosphere”, but the figure serves well its illustrative purpose.

As described above, the opacity consist of a product between the absorption coefficients and the corresponding partial pressures. Fig. 7 shows the sum of the partial pressures of the 22 molecules that went into the computation of our classical DRIFT-MARCS models, compared to the sum of the partial pressures of the full set of 50 molecules introduced in the present grid, for the new models of T_eff =2500 K, 2000 K, 1500 K and 1000 K, respectively. The effect on the final atmospheric structure depends both on how large the partial pressures are and how strong the absorption coefficients are, and sensitively on how the absorption coefficient is distributed as function of frequency, but the relative change in the partial pressures between the models gives a good feeling for the reasons behind the changing model structures.

H₂ is the dominant contributor to the gas pressure at the temperatures and compositions discussed here, but H₂ in itself has a very small absorption coefficient concentrated on a few narrow frequency intervals. In order to illustrate most clearly the effects of the 28 new molecules added to the present grid, the sum of the partial pressures of the 22 opacity bearing molecules in the old grid and the 50 in the new grid are therefore shown as the sum of these minus the partial pressure of H₂ (normalised to the gas pressure of the respective models). Figure 7 indicates, as expected, that the sum of the 22 partial pressures (black dotted line) and 50 partial pressures (black dashed line), defined as described above, are almost identical for the (new) model of T_eff = 2500 K, while for the (new) models of 1500 K (green lines) and 1000 K (blue lines) the sum of the partial pressures of the 50 opacity bearing species (dashed lines) is substantially higher than the corresponding sum of the 22 molecules (dotted lines) in the classical models, in particular for the upper and coolest atmospheric layers at temperatures around and below 1000 K.

Figures 8 and 9 show how the distribution of the gas pressure onto various individual species develops from T_eff =2500 K (Fig. 8) to T_eff =1000 K (Fig. 9). From top to bottom in the two figures, the long-dashed black line and the short-dashed green line show the sum of the partial pressures of all the neutral molecules together and the partial pressure of H₂ alone, respectively. For the T_eff = 1000 K model these two lines are identical and practically equal to the total gas pressure (log(pp/p_gas) ≈0), while for the T_eff = 2500 K model they are also identical but deviate a bit from the total gas pressure in the deeper layers of the model because of the relative increase in the sum of all neutral atoms (dash-triple-dot purple line) which mainly consist of hydrogen (full drawn black line) and helium (yellow dotted line) in solar metallicity models.

With three orders of magnitude lower partial pressures we find (in Figures 8 and 9) the sum of the partial pressures of all the 22 molecules that contribute to the opacity of the classical models (full drawn red line) and the corresponding 50 in the new models (full drawn green line). In the T_eff = 1500 K (not shown) and T_eff = 1000 K (Fig. 9) models the partial pressures of the 50 molecules (minus H₂) are of the order of twice that of the original 22 opacity carrying molecules in the old models. It is worth noticing how small a fraction of the total number of molecules per cm³ of the stellar gas that actually carries all of the opacity and hence determines the structure. For solar metal- licity objects it is approximately 0.1% of the gas pressure that carries 100% of the opacity and for metal poor stars this percentage is correspondingly lower. This is in fact approximately the same in the Earth’s (cloud-free) atmosphere where N₂ and O₂ carries approximately 99% of the gas pressure but almost no opacity (because of their homonuclear structure). It is the reason for why so relatively small changes (in temperature or chemical composition) have so dramatic consequences for the structure (for stars as well as for Earth).

Right below the two lines with the sum of the 22 or 50 opacity bearing molecules (which are identical to one another in the T_eff = 2500 K case) we find the partial pressures of CO and H₂O. In the T_eff =2500 K model the relative partial pressures of CO and H₂O are an almost constant fraction of the gas pressure at all atmospheric depths, with CO dominating, while the ratio of the two form a more complex pattern in the cooler models, as is well illustrated at T_eff =1000 K.

It is seen by comparing Figs. 8 and 9 that while the partial pressure of CH₄ (methane) is orders of magnitudes below those of CO and H₂O in the T_eff = 2500 K model, even at such low temperatures as 1500 K, the partial pressure of methane approaches that of CO and water already at 2000 K in the T_eff = 1000 K model, and supersedes the one of CO at T ∼ 1500 K. The reason for this difference has to do with the gas pressure, which is 5 orders of magnitude higher at 1500 K in the T_eff = 1000 K model than in the T_eff = 2500 K model, and the general tendency of the chemical equilibrium to shift towards larger molecules at higher pressures. This illustrates in an extreme form why the transition between L-type spectral class and T-type spectral class (defined as the visibility of methane in the spectrum) must depend not only on T_eff , but also on log(g). If the effective temperature is well determined, the ratio between the intensity of CH₄ and CO bands can therefore in principle be used to determine log(g) if the object is a star or a brown dwarf and the mass of the object if it is a transiting exoplanet.

It is the same phenomenon that explains the huge difference between the H₂O (light blue dashed-triple-dot line) to OH (dotted green line) ratio at similar temperatures in the T_eff =2500 K and T_eff =1000 K models. We will further see (in a later section) that TiO will be much stronger in the spectrum of the warmer of the models in the present grid than in the colder ones, and we see here that it is not because the partial pressure of TiO is larger in the warmer models (in fact it is very low in both of the models presented in Figs. 8 and 9), but because the partial pressure of TiO increases outwards in the warmer models and inwards in the cooler models, while the opposite effect is present for NH₃ (ammonia). This difference in the behaviour in the partial pressures of the two molecules makes them not only a good indicator of T_eff, but the difference in the reduced observed band intensities of the two molecules compared to what can be computed from cloud free models, as those presented here, becomes a direct indicator of the height in the atmosphere of the presence of clouds.

Titanium oxide, TiO, has a relatively low partial pressure in both models, although several orders of magnitude lower in the uppermost layers of the T_eff = 1000 K model than in the T_eff= 2500 K model. However, its high absorption coefficient per molecule and its many electronic transitions in the visual part of the spectrum, makes it an important contributor to the structure and spectrum for the T_eff = 2500 K models, and it will also play a role even in the T_eff = 1000 K models. In addition, it plays a crucial role in the temperature inversion for irradiated models, as discussed below, and therefore is also a good diagnostic in this. Hydroxide, OH, is the most important species for initiating atmospheric oxidation on Earth, and (in chemical equilibrium) it is considerable more abundant in high temperature models than in low temperature models.

Finally, the blue dashed lines show the partial pressure of the sum of all positively charged atoms and the dash-dot light blue lines show the sum of the partial pressures of all negatively charged atoms. The ratio between the two curves is a direct illustration of the origin of the electron pressure, and is another way of quantifying the strongly decreasing electron pressure in the cooler models already discussed in connection with Fig. 3. The partial pressures of the negative ions as well as the positive have become very small in the T_eff = 1000 K model, and the partial pressures of the positive and negative ions are almost identical in the upper atmospheric layers, hence the electron pressure (not shown) goes towards zero as already discussed above.

Fig. 6 Rosseland mean opacity as function of gas pressure, for models with Teff from 2500 K (upper curves) to 1000 K (lower curves), based on the full set of 52 molecular opacities (full drawn lines) and the original 22 molecular opacities, respectively, in the radiative transfer. Colour coding correspond to the effective temperature of the models, as shown in the figure legends.

Fig. 7 Logarithm of the sum of the partial pressures (minus the partial pressure of H2 , and normalised to the gas pressure) of respectively the original 22 opacity bearing molecules (dotted lines) and the 22 plus the new 30 ones (dashed lines), for the new models presented here of Teff = 2500 K (black lines), 2000 K (red lines), 1500 K (green lines), and 1000 K (blue lines).

Fig. 8 Distribution of the total gas pressure onto its various components for a model of Teff =2500 K. As indicated with the legend, the curves show the partial pressures of the molecules H2 , CO, H2O, CH4 , NH3 , TiO, and OH, the neutral H and He atoms, the sum of all positively charges species (“pel+”), all negatively charged species (“pel-”), sum of the partial pressures of all neutral atoms together (“pnat”), all neutral molecules (“pnmo”), and finally the sum of all the opacity bearing molecules in the classical MARCS models (“p24m”) as well as in the new models presented here (“p54m”). We note that the curves for H2 and all neutral molecules (pnmo) as expected plot on top of one another, as do the curves p24m and p54m (therefore p24m/p54m show as a single full drawn green curve, and the only red curve is for the partial pressure of TiO). The upper black curve is for H and the lower for NH3 .

Fig. 9 Same as Fig. 8 but for at model of Teff =1000 K. We note that for this model the curves p24m and p54m do deviate, in particular for the upper layers, illustrating the importance of the additional molecular opacities in the MSG models compared to the classical MARCS models.

7 The effect of irradiation

The influence of irradiation of one component onto the other in eclipsing binary stellar systems was already included into the MARCS code by Vaz & Nordlund in the 1980s (Vaz & Nordlund 1985; Nordlund & Vaz 1990) in a plane parallel treatment. We use here these early initiatives as a starting point to develop an irradiation scheme in the MSG code that can realistically describe the effects of host star irradiation on an exoplanet (as well as being able to treat the atmospheric structure and spectrum formation in binary stars in a more up-to-date radiative transfer and geometry description). The MSG scheme includes angle dependent irradiation from the host star, and can in principle be computed in plane parallel as well as full spherical geometry of the host star and/or the irradiated exoplanet atmospheres.

In plane parallel geometry of the planet, the optical path of the incident irradiation converges towards infinity for 90 degrees angle towards zenith (i.e. at the “morning and evening terminator”, Eq. (5)). To compensate this, we followed in our plane parallel computations the approximations introduced by Malik et al. (2019) (see in particular their equation 27 and associated figure 2). Since the models are static, and therefore do not include rotation and dynamical effects as in a GCM 3D modelling, we term the method 1.5D atmospheric modelling, in the sense that there will be no difference in the angle dependence of the atmospheric structure between planetrary longitude and latitude. In our future irradiated models we plan to implement a parametrisation of computationally extensive GCM results (which was presented within our team in, e.g., Schneider et al. 2022), which will allow us to keep the high flexibility and the speed that 1D self-consistent studies offer for a large grid of exoplanets with complex physical, chemical, and biological effects. A route forwards for this coupling of 3D modelling results into the 1D computations was developed in detail in Hubeny (2017), in particular for irradiated models in his section 2.2. More details of our own experiments with hybrid coupling of the MSG modelling with 3D results is described in Kiefer et al. (2024). Here we describe the most basic features of our irradiated models as they are included in the present grid, with illustration of a few examples in Figs. 10 and 11.

To a first approximation, the irradiation intensity J_ν,star reaching layer k of an exoplanet from a host-star shining onto the planet from and angle Θ_star with the zenith direction at a given longitude and latitude of the planetary surface, can be expressed as $$J_{star}(\nu ,\mu ,k)=B_{\nu }\left(T_{star}\right)\text{ΔΩ}exp\left(-{\tau }_{star}\right),$$(3)

where B_ν(T_star) is the black body flux of a host star of effective temperature of T_star . ∆Ω is the solid angle subintended by the star onto the planet, $$\text{ΔΩ}=\pi {\left(\frac{R_{star}}{a_{au}}\right)}^2\frac{1}{4\pi },$$(4)

where a_au is the semi-major axis of the exoplanet’s orbit and R_star is the radius of the host star. In the plane parallel approximation, µ_star is a simple function of Θ_star, $${\mu }_{star}=1/\cos \left({\Theta }_{star}\right),$$(5)

and τ_star is the optical depth scale that the stellar light travels inside the planetary atmosphere to reach layer k of the zenith direction at the local position of the planet, τ_star = τ_kµ_star.

Since the MSG grid can be used to compute the structure (and spectrum) of the exoplanets as well as their host stars, B_ν(T_star) can be substituted with the real stellar spectrum to compute a correct and fully frequency dependent planetary transit spectrum, which can eventually be folded with a filter function to give a transit light curve for any arbitrarily chosen filter, without any assumed relation between scale height, radius and frequency.

A more detailed description of the effects of irradiation onto the models will be given in separate papers (Amadio et al, in preparation, D’Alessandro et al, in preparation). Here we give just two examples of preliminary results from this part of the grid. Figure 10 illustrates how the model structure changes as function of irradiation angle of the host star relative the planetary surface. Quantitatively the figure is a simulation of a cloud- free exoplanet with parameters as those of WASP-39b (i.e. a 0.28 M_jup mass planet orbiting a T_eff = 5400 K G8 type star at a semi-major axis of 0.0486 au), with the internal model temperature estimated from the relations between equilibrium temperature and internal temperature given by the formulas in Thorngren et al. (2019). Within the limitations of 1.5D modelling, the 3 curves show the gas-pressure temperature structure of the atmosphere “at noon” (the sub-stellar irradiation corresponding to Θ = 0°; the right-most curve), Θ = 15° (middle curve) and Θ = 30° (left most curve). Figure 11 illustrates the temperature inversion caused by the strong absorption of the irradiated star light by TiO and VO in our model of WASP-43b (a 2 Jupiter-mass planet orbiting a K7V star of T_eff = 4500 K at an orbital semi-major axis of 0.015 AU). The effect is equivalent to the temperature inversion due to absorption of solar irradiated energy by ozone in the Earth’s atmosphere.

Fig. 10 Temperature versus gas pressure structure of our irradiated 1.5D model simulation of an exoplanet with the parameters of WASP-39b at the sub-stellar point (right-most curve), at Θ = 15° and at Θ = 30° (leftmost curve).

Fig. 11 Temperature inversion caused by TiO and VO at the sub-stellar point in irradiated MSG 1.5D model simulations of WASP-43b.

8 The effects of microphysical cloud formation

To study the effects of microphysical cloud formation in sub-stellar atmospheres, the MSG code couples MARCS and StaticWeather to one another in a self-consistent manner (see Section 2). This implies the cloud radiative feedback is accounted for in the radiative transfer scheme. The cloud radiative effect is added to the radiative transfer by considering the cloud’s opacity contribution, as well as the gas-phase elemental depletion caused by the cloud formation process. Given the cloud is composed of a mix of materials, we calculate its opacity with spherical particle Mie theory (Mie 1908) combined with effective medium theory.

The effective optical constants for the material mixtures were obtained with effective medium theory. We generally used the Bruggeman method (Bruggeman 1935), with the exception of rare cases where we found non-convergence and therefore used the analytic Landau–Lifshitz–Looyenga [LLL] method (Looyenga 1965) (see Section 2.4.1. of Lee et al. (2016) for more details). The clouds extinction and scattering coefficients are computed with Mie theory using the routine developed by Wolf & Voshchinnikov (2004), which is based on the widely used Bohren & Huffman (1983) routine. However, the Wolf & Voshchinnikov (2004) routine allows for the consideration of larger dust grain sizes.

Here we consider 12 condensates: TiO₂[s], Mg₂SiO₄[s], MgSiO₃ [s], Al₂O₃ [s], Fe[s], SiO[s], SiO₂[s], FeO[s], FeS[s], Fe₂O₃ [s], MgO[s] and CaTiO₃ [s]. We use TiO₂ [s] as the nucleation species. For the dust optical constants, we use the data tables compiled by Kitzmann & Heng (2018) (see Table 1 in Kitzmann & Heng (2018) for references). We note we use the amorphous (sol-gel) data for Mg₂SiO₄ [s] and MgSiO₃ [s], and the amorphous data for SiO₂ [s]. The most relevant data are given in Table 2 below, and more details are given in Campos Estrada et al. (2024).

We ran MARCS and StaticWeather iteratively until we found a converged solution. Besides the typical convergence criteria in the radiative-transfer scheme, we required both the cloud opacity and the emission spectra between the last two iterations to be converged.

The upper panel of Figure 12 shows the effect of clouds on the atmospheric temperature versus gas pressure structure of a MSG model of a self-luminous object with T_eff = 1500 K and log(g) = 4.0, at solar abundances and C/O. The radiative feedback of the clouds has a warming effect in the entire atmosphere (solid curve labeled MSG) compared to a cloud-free atmosphere of the same T_eff and log(g) (dashed curve labeled MARCS). We note in Figure 12 that the radiative regions are plotted with a thinner line width than the convective regions which are plotted with a thicker line width. The cloudy model presents a detached convective zone between 0.01 and 0.1 bar. This detached convective zone has its origin on the cloud’s back warming effect. We discuss the emergence of detached convective zones in our cloudy models in detail in Campos Estrada et al. (2024). Such detached convective zones are also observed in other model grids (e.g. Morley et al. (2024)). The lower panel in Figure 12 shows the average cloud particle size and cloud particle number density for the model with clouds. We note where the average cloud particle size is larger, there is a drop in the number density as physically expected.

The composition of the clouds and the gas elemental depletion are shown in Fig. 13. At the very top of the atmosphere (TOA) we see that TiO₂ [s] is the major component of the cloud particles as it is the nucleation species in this specific model presented (Fig. 13 upper panel). This also explains why Ti in its gas-phase is highly depleted at the TOA (Fig. 13 lower panel). Although it is not as visible, O is also depleted at the TOA. However, O is a lot more abundant than Ti, and therefore this depletion is not as drastic. After TiO₂ [s] nucleates, other cloud species grow on top of the TiO₂[s] seeds. Throughout most of the atmosphere, the bulk of the cloud is composed of Mg₂SiO₄ [s], MgSiO₃ [s], SiO[s] and SiO₂ [s]. Deep in the atmosphere where the temperatures become too high for silicate to be thermally stable high-temperature condensates take over the bulk of the cloud (Al₂O₃ [s] and Fe[s]). Towards the bottom of the atmosphere, we can see regions where there is an enrichment of the cloud forming elements in their gas-phase. These enrichments occur when particular cloud species evaporate: for example, we can see that there is an enrichment of Mg and Si when the silicates and magnesium-silicates evaporate. This is because these elements were trapped in the cloud particles and transported down until those cloud species evaporates, at which point we see an enrichment of the elements involved in that region of the atmosphere. We note that the gas-phase equilibrium chemistry is always computed after we have computed the gas-phase element depletions from cloud formation.

Figure 14 shows the relative number density of some important molecules and atoms in sub-stellar atmospheres, including H₂ O, CO and CH₄ , along the atmosphere for the MSG and MARCS models shown in Figure 12. The gas-phase chemistry is of course different in each model, not only because of the effect of the cloud formation surface reactions but also because of the distinct temperature gas pressure structures (see CH₄ for example).

Clouds present a challenge to atmospheric characterisation as they can often hide spectral features of the gaseous components of atmospheres. Ignoring the presence of clouds can lead to a wrongful determination of various atmospheric parameters, such as metallicity.

Figure 15 shows the comparison between the emission spectra of the cloud-free and cloudy models shown in Fig. 12. The spectrum of the cloudy model shows that a number of the molecular features present in the cloud-free spectrum are diminished by the presence of clouds. This is due to a combination of the decrease of molecular abundances due to the higher atmospheric temperatures of the cloudy model, as well as the contribution of the cloud opacity, which is significant in the observable atmosphere.

In the final version of the cloudy MSG models, to be added to the present grid, we will present and analyse atmospheric models of a larger range in T_eff (Campos Estrada et al. (2024)), as well as discuss the importance of the cloud radiative feedback, the cloud composition, the nucleation, the mixing timescales, and other intricacies of the microphysics processes.

Fig. 12 Effects of clouds in a Teff =1500K model atmosphere. Upper panel: pressure–temperature profiles for MSG models at Teff =1500K and log(g)=4.0 respectively with clouds (labeled MSG) and without clouds (labeled MARCS). We note that the convective part of the atmosphere is represented with a thicker line. Lower panel: average cloud particle size (solid curve) and cloud particle number density (dashed curve) for the MSG model shown in the upper panel. We note that the spike in the number density at the bottom of the cloud is numerical and due to the rapid decrease of the average cloud particle size.

Fig. 13 Cloudy MSG model at an effective temperature of 1500 K and log(g)=4.0. Upper panel: cloud volume fraction of the 12 cloud species present in the atmosphere. Lower panel: relative elemental depletion in the gas phase of elements affected by the cloud formation.

Fig. 14 Relative number density of important molecular and atomic gas species in the cloudy MSG (solid curves) and cloud-free MSG (dashed curves labeled MARCS) models at 1500 K and log(g)=4.0 (shown in Figure 12).

Fig. 15 Cloudy (blue) and cloud-free (grey) synthetic spectra (binned to a rsolution of 1000) of a self-luminous body with an effective temperature of 1500 K and log(g)=4.0, with solar metallicity and solar C/O ratio.

9 The effect of non-equilibrium chemistry

A standard adopted assumption in classical (stellar) atmospheric computations, including most versions of the MARCS code, is the local thermo-dynamical equilibrium, LT E, approximation, where the relative abundances of the chemical species is assumed to be in equilibrium with the radiation field “locally”, i.e. layer by layer the kinetic and radiative temperature are assumed to be identical to one another. In optically very thin atmospheres, for example giant or supergiant stars, the radiation field in a given layer is often a combination of local Planck-like radiation augmented with a non-negligible radiation field from neighbouring layers of different temperatures, giving rise to what in astrophysics classically is termed NLTE (NonLocal-Thermo-dynamical Equilibrium). In a broader sense the non-local radiation field could even be a layer of a neighbouring object, such as we know it from the Earth’s ozone-layer that is affected by absorption of UV-radiation from the Sun because the uppermost layers are transparent to solar-temperature UV radiation.

In exoplanetary model atmospheres the chemical equilibrium in the higher layers are often assumed to be a combination of LTE and a gas-component from a wind-system that moves the gas fast enough upwards to add an LTE component from an atmospheric layer of a non-local temperature (and density), which in its simplest form is often approximated by an almost free (in the case of exoplanets with unknown topography and wind-systems) parameter, expressed by the so-called eddy diffusion coefficients.

In a more detailed description of the deviation from LTE one would like to include a full network of forwards and backwards reaction rates (to obtain a local steady state) and a realistic dynamical mixing of the layers. Only with such a more complex treatment of the gas chemistry would it be possible to realistically quantify further contributions from potential entropy-reducing entities (including biology), which is the final goal of the additions of more complex MSG models to the basic grid presented here. Earth’s atmosphere reflects the existence of all of these non-equilibrium processes, from classical NLTE (in the common astrophysical use of the term) represented most extreme by the effect of the solar radiation field on the ozone formation in the upper atmosphere and the photochemistry of the OH radical, over the dynamical non-equilibrium estimated from empirically confined eddy mixing, to the introduction of (biologically and non-biologically produced) oxygen and methane on dynamical timescales faster than the corresponding reaction rate timescales. Figures 16, 17, and 18 illustrate preliminary results from our first steps of including steady state computations as well as dynamical mixing into the MSG model grid towards the aim of including self-consistent biological effects in the modelling.

As a first illustration (Fig. 16) of the effect of including time-dependent reaction rates into the MSG modelling, we have constructed a limited 4-reaction Chapman mechanism chemical network capable of computing the O, O₂ and O₃ abundances. The input starting abundances (left side of the figure; time zero) are the LTE abundances from GGchem, and the network is then propagated for long enough in time to assure that steady state has been reached, as shown in the right side of Fig. 16. The figure shows the result from a single layer of an MSG model atmosphere of solar abundance, T_eff = 2500 K and log(g) = 4.5. The chosen layer had T = 1950 K and P_gas = 5.5 10⁴ dyn/cm². The photochemistry is driven by irradiation from a solar type star at 1 AU from the modelled object (for example an exoplanet), and the rate coefficients were taken from the NIST database. We note how the stellar irradiation in a timescale of a few thousand years increases the relative ozone abundance with a factor of ∼10¹⁰ compared to an LTE modelling.

As a next example (Fig. 17) is shown the O₃ partial pressure as function of total gas pressure, calculated by assuming respectively thermo-dynamic equilibrium and non-equilibrium from a 4-reaction Chapman-like mechanism (as in Fig. 16), now for a surface-less exoplanet with Earth-like atmospheric elemental abundances, log(g) = 3.0 and T_eff = 1000 K, illuminated by a Solar-type star (of T_eff = 5772 K and solar radius) at 1 AU from the planet. We note how the ozone abundance starts deviating from the LTE abundances at gas pressures below a few bars, and that the non-equilibrium ozone abundance reaches more than 25 orders of magnitude increased partial pressure compared to its LTE abundance in the uppermost layers of the model. Combining the irradiation and non-equilibrium modelling capabilities of the new MSG code with the Earth-like elemental abundances allows for modelling of the photochemical Chapman-mechanism responsible for formation of the ozone layer in Earth’s atmosphere. The photolytic dissociation of O₂ results in formation of oxygen atoms, which combine with O₂ to form O₃ . This results in a local maximum in O₃ concentration comparable to the ozone layer in Earth’s atmosphere. At an effective temperature of 1000 K, the effect of the ozone layer is too small to affect the pressure-temperature structure of the model, but it could be visible in the spectrum.

To allow for modelling of more general steady-state conditions, we have developed a (beta-)version of the MSG that is coupled to the KROME package, designed specifically to include time-dependent chemistry into astrophysical models (Grassi et al. 2014). It attempts to solve the chemistry by treating all relevant chemical processes with their respective rate coefficients, as a system of ordinary differential equations. The relevant equation that needs to be solved would therefore be $$\frac{d{n}_i}{dt}={\displaystyle \sum _{j\in F_i}\left(k_j{\displaystyle \prod _{r\in R_j}{n}_{r(j)}}\right)}-{\displaystyle \sum _{j\in D_i}\left(k_j{\displaystyle \prod _{r\in R_j}{n}_{r(j)}}\right)},$$(6)

where n_i is the i-th species in the network, and the first sum represent all the processes that would create the i-th species, with their respective rate coefficients k_j and the number densities of the reaction partners n_{r( j)} , whilst the second sum represents all processes that would destroy the i-th species, with again their respective rate coefficients k_j and the number densities of the reaction partners n_{r( j)}. For more details on how the system of differential equations is setup and solved, and which processes KROME is able to include, we refer again to the original publication Grassi et al. (2014).

Given a chemical reaction network in the KROME format, the MSG species included in the reaction network are treated kinetically in each layer starting from the GGchem equilibrium concentrations, as illustrated in Figs. 16 and 17. The MSG atmospheric structure can then be iterated until self-consistent convergence is obtained with the non-equilibrium steady state concentrations of the selected species. All remaining species retain their equilibrium abundances.

As a final example of non-equilibrium, we show in Fig. 18 effects of the often used quenching approximation in exo- planetary atmosphere computations (for example PICASO-3.0; Mukherjee et al. 2023), which assumes that vertical flow of air under some conditions can bring bubbles of air from deeper layers to higher layers fast enough that the chemical equilibrium of the deeper layer is conserved in the higher layers in spite of their different temperatures (and gas pressures). Hence, the assumption is that the convective (quencing) timescales are shorter than the chemical reaction timescales. A typical introduction of quencing could be to let the chemical abundance of the first layer (from bottom and upward) where the vertical convective timescales was found to be smaller than certain reaction timescales in that layer, persist throughout the remaining part of the atmosphere (upward).

In Earth’s atmosphere such movements could typically take place when a mountain range forces air flow to move quickly upward, expressed by eddy diffusion with the diffusion coefficients fitted to atmospheric observations. Corresponding coefficients for a given exoplanet could be treated as a free parameter that could be adjusted to obtain a desired structure or reproduce an observed spectrum. We tested quenching based in MSG computations as an option for the 8 species CO, CH₄, H₂O, NH₃, N₂, CO₂, HCN and PH₃ . The typical chemical timescales for conversion are calculated based on parametrisations (Zahnle & Marley 2014 and Visscher et al. 2006) and compared to the timescale of vertical mixing as approximated by a fixed eddy diffusion coefficient, K_ZZ . In regions above the lowest layer where the vertical mixing time scale for the chosen K_ZZ -value in this sense is lower than the estimated timescale for chemical conversions, the mixing ratios of the selected species are then kept fixed. An example is shown in Fig. 18 for a chosen fixed value of K_ZZ = 10¹⁰ cm²/s. This option will allow for testing of the quenching approximation against a fully self-consistent chemical non-equilibrium layer-by-layer computation, which will be the preferred solution, especially when no particular topology and wind system is known.

Fig. 16 Example of the temporal evolution of the O, O2 , and O3 number densities from their equilibrium abundances (left) to steady-state (right) for a 4-reaction Chapman mechanism chemical network.

Fig. 17 Local thermodynamic equilibrium and non-equilibrium steadystate partial pressures of O3 as a function of total gas pressure (shown as full-drawn and dotted lines, respectively) for a surface-less exoplanet with Earth-like atmospheric elemental abundances, log(g) = 3.0 and Teff = 1000 K, illuminated by a Solar-type star (of log(g) = 4.4, Teff = 5772 K, and solar radius) at 1 AU from the planet.

Fig. 18 Atmospheric mixing ratios of the molecules for which quenching has been implemented into MSG. Solid lines show the equilibrium values, while the dotted lines show the mixing ratios with a fixed eddy diffusion coefficient of Kzz = 1010 cm2/s. The model has Teff = 2500 K and log(g) = 4.5.

10 Partial pressures and spectral diagnostics

The presented MSG grid covers the temperature range of late- M, L, T, and the bulk of Y spectral type objects, representing the temperatures of cool dwarf stars, brown dwarfs, hot to temperate gas and rocky exoplanets, including the temperature range of the atmospheres of Earth-like planets. This diversity of objects cannot be described by a simple sequence of T_eff, log(g), and scaled solar chemical composition. The ambition of this first basic MSG grid of these complex objects therefore is also only to present some fundamental features, understanding, and methods relevant for the sequence, focusing on the basic gas phase chemistry and self-consistent radiative transfer. Ongoing work on inclusion of more complexity into the grid includes the issues illustrated in the three previous sections. In these last two short sections we will based on the existing basic grid first touch upon the fundamental reasons behind the transitions between the empirically defined late-type spectral classes (the present section), and some of the expected major changes in the self-consistent atmospheric structure when shifting from solar-like to Jupiter-like and Earth-like compositions (next section).

The radiative properties of self-consistent, cloud-free, static atmospheres are determined by the opacities, which are determined by the line transitions, line-profiles and sampling techniques as described in previous sections, which determine the absorption coefficients (Fig. 19) and hence the partial pressures (Fig. 20) and the radiative transfer leading to the spectrum (Fig. 21).

Fig. 19 shows the monochromatic absorption coefficient for the most important molecules needed for the understanding of the M-, L-, T- (and Y-)sequence. Due to its many electronic transitions, TiO stands out from most other species by the enormously strong absorption coefficient per molecule in the visual and photographic infrared region. It also has a relatively strong transition in the region from 10 to 12 µm, that per molecule is comparable to that of water, methane and ammonia. Due to the low abundance of Ti relative to C, N, and O, it is usually only the visual transitions that are discussed in observational spectra, but comparison of the visual and infrared absorption holds a rather unexplored diagnostic potential. VO has an absorption coefficient relatively similar to that of TiO, and also the absorption coefficient of FeH is similar to that of TiO at least in the photographic infrared. The solar abundance of Fe is more than two orders of magnitude higher than that of Ti, while the abundance of V is almost an order of magnitude below than that of Ti.

The integrated absorption coefficients per molecule are rather similar for H₂O, NH₃ and CH₄ , but due to the difference in the internal configuration of the three molecules, water is the most temperature dependent of the three and CH₄ the least. For water, the amount of weak lines increase rapidly with temperature, making the absorption coefficient stronger and more “continuum-like” for high temperatures, giving it an important structural role for stellar and sub-stellar objects of a wider range of effective temperatures than the two other abundant polyatomic molecules. The difference in the C-H and C-O stretching (and bending) energies happens to be so that the most pronounced H₂O and CH₄ bands alternate with one another throughout large parts of the spectrum (which is not the case for NH₃) and hence give rise to a much stronger structural effect than would have been from each one of H₂O and CH₄ alone. This is also the basic reason for why methane in even tiny amounts have a relatively large green-house effect on Earth’s atmosphere, and this will make it possible to quantify even small biological contributions to exoplanetary atmospheres by measuring trace elements in (coming) high-resolution exoplanetary spectra if the planet has water and its biology includes metabolic methane as it does on Earth.

CO also has a strong absorption per molecule. Being a diatomic molecule it has (as TiO) much fewer vibrational transitions than water, methane and ammonia and hence more discrete infrared bands, so even though the individual bands can be strong, CO has less importance for the structure in cool stellar and sub-stellar environments than the corresponding polyatomic molecules. Most pronounced is its fundamental vibrational transition around 5 µm, but also the overtones in the more accessible JHK region are strong. Due to its unusual high binding energy it will often consume most of the less abundant of either C or O, and therefore remain with a relatively constant band intensity in the spectrum independent of the atmospheric C/O ratio. The ratio of the H₂O and the CH₄ band intensities towards the CO band intensities are therefore good indicators of the atmospheric C/O ratio in cool stellar and sub-stellar environments.

CO₂ has a strong near- and mid-infrared absorption coefficient, and is well known from Earth, Venus and Mars to be able to produce a strong absorption around 5 µm, even at low abundances. If it does show up in the spectrum simultaneously with CO and CH₄ it has enormous biological diagnostic potential, and its presence together with methane, water, and ozone in Earth’s infrared spectrum is considered evidence of photosynthesis on Earth. At high temperatures the absorption coefficient becomes “continuum-like” resembling that of water a wavelengths beyond 10 µm. Its strongest transition overlap with the corresponding CO transition, and depending on temperature, gas pressure and relative oxygen abundance, it should be expected to have smaller structural effects on habitable exoplanets than methane.

While the absorption coefficients are basic for identifying the spectral features, it is the product of the absorption coefficient and the partial pressure that determines the opacities and hence the atmospheric structure and the quantitative spectroscopy. The development of the partial pressures (normalised to the gas pressures) through the M-, L-, T-sequence is shown in Fig. 20 as function of optical depth. Colour coding representing the different molecules is indicated with the legends, and the results are shown for models of T_eff = 1500 K (full drawn lines), 1000 K (dash-dot lines), and 500 K (dashed lines).

In the observational classification of the M-, L- and T-type (and Y-type) objects, the transition between M-type and L-type is seen as the gradual disappearance of the TiO bands from their maximum intensity around M7 (T_eff ∼2700 K) towards a disappearance (in the visual to near-IR) at spectral type L3 to L6 (T_eff ∼1700 K). VO plays a similar role in the classification, although appearing and disappearing over a more narrow range of spectral classes, usually not seen before in the very latest M-types. Parallel with the disappearance of the small oxides, is observed an intensified appearance of metal hydrides and hydroxides (FeH, CaH, MgH, CaOH and others) reaching maximum intensities for mid L-types and disappearing again with late L-types.

Somewhat surprising, atomic lines, that would usually have been associated with much warmer objects, show up in the observationally defined L-type sequence, and gets broader and broader towards L5 to L8. In spectral type L5 the Na-triplet diverges into what is interpreted as one broad line almost 1400 Å wide, and the K-triplet similarly reach a broadness of 1000 Å in type L8. It is hard to understand this as a temperature effect, and one suggestion has been that the decreasing electron density results in decreasing continuum opacities, which in turn result in a more transparent atmosphere, such that one sees deep into the atmosphere through a large column density of neutral atoms, which would increase the line intensities and in particular contribute to broader Van der Walls wings due to the higher pressures in the warmer layers at the bottom of this clear photosphere. This interpretation is not immediately in agreement with our MSG models. Our models do show a particularly marked decreasing electron pressure (by more than 15 orders of magnitude) in what could be interpreted as the late L-type sequence (Fig. 3), but our continuum is not very dependent on the contributions from electrons already at higher values of T_eff (Fig. 21). The interpretation is open for more detailed analysis by use of the present grid. The understanding obviously also raises the question of the general assumption of more dust-clouds influence on the cooler atmospheres, and it all underlines the complexity of the processes compared to warmer objects. At what is interpreted as even lower effective temperatures, the observed metal lines disappear from spectra again.

At lower effective temperatures than the L-type sequence, methane (CH₄ ) appears in the spectrum, and the traditional classification feature of the beginning of the T-type sequence is the appearance of the CH₄ overtone bands at 1.6 and 2.2 microns, but the fundamental band at 3.3 microns is of course even stronger, and we see it in our cloud free models already in what is traditionally classified as mid-L type objects at T_eff ∼1500 K (Fig. 21). The T-type objects are sometimes called the methane brown dwarfs, and the whole T-type sequence belongs to sub- stellar objects. For even cooler temperatures, ammonia (NH₃) becomes abundant, although these objects (some times called Y-type objects) usually are so dim that the spectra are open for several interpretations. Our cloud-free solar metallicity models predicts the NH₃ bands to be very strong in the T_eff = 500 K models, but overlapping strongly with the CH₄ bands and probably best identifiable around 6 µm (Fig. 21).

Just as for the changing depth distribution of TiO in the atmosphere (Figs. 8 and 9), we have seen in Figs. 20 and 21 that in the coolest MSG models presented also NH₃ changes depth distribution with changing values of T_eff rather than changing total atmospheric abundance, but in the opposite depth direction that TiO (and neutral atoms). Depending on the height of cloud formation, NH₃ may therefore form above a potential cloud layer in the coolest layers (just as in the much cooler atmosphere of Uranus), and the relative intensity of the atomic lines and the strongest TiO bands, on the one hand, and the ammonia bands, on the other, may therefore be valuable detailed diagnostics of the level of cloud formation in self-consistent modelling. The relative weakening of the TiO and NH₃ bands will depend on how high in the atmosphere the clouds form and the T to Y transition is therefore not solely a measure of a temperature sequence, but depend among many factors also on the height in the atmosphere that the clouds form. High level clouds will weaken or erase the NH₃ features, while low lying clouds weakens the TiO bands but not the NH₃ bands. In our preliminary cloud models (Juncher et al. 2017), the first thin clouds can already form in the upper layers of models of T_eff≈2500K and will therefore move the M-L (and L-T) transition to higher values of T_eff than in the cloud-free models.

Fig. 19 Absorption coefficient at T = 1500 K of the most important molecules (listed in the legend) for the structure and understanding of the coolest dwarfs.

Fig. 20 Partial pressures of the most important molecules in models of Teff = 1500 K, 1000 K, and 500 K, respectively.

Fig. 21 Development of the spectrum of H2O and TiO (the two left panels), CH4, and NH3 (middle panels) and FeH and VO (right panels) for models of Teff = 1500 K (red), 1000 K (blue), and 500 K (green). The dashed curves are the true continua (i.e. the spectrum computed by allowing only the continuum opacities into the spectrum computation, but all opacity sources into the model computation), while the full drawn lines likewise are the spectra of each individual molecule listed in the legend, computed by allowing the opacity from only that molecule in the spectrum computation, but all opacity sources into the model computation. Remark that for the coolest models (and the plotted resolution) the combined spectra do not at any place reach even close to the (true) continuum, which is the continuum that should be the reference in abundance analyses.

11 Solar versus Jupiter-like and Earth-like chemical compositions

Primordial planetary atmospheres of all kinds, as well as present day gas-giant exoplanetary atmospheres, may well be close to (scaled) solar composition for many elements. Jupiter’s atmosphere has close to the solar relative abundance of C, N, S (and the noble gasses Ar, Kr, and Xe), although they all are a factor of approximately 3 increased relative hydrogen, compared to their solar ratios. He and P are almost solar (relative to hydrogen) while Ne and O are of the order a factor 10 down compared to their solar values. Expressed in the standard є-scale (i.e. as the logarithm of the number abundance of an element A per 10¹² hydrogen atoms, є(A) = log(A/H)+12), then Jupiter’s abundances (Atreya et al. 2003) would translate into epsilon values of є(He, C, N, O, P, S, Ne, Ar, Kr, Xe) = (10.896, 9.021, 8.593, 7.930, 5.486, 7.607, 7.090, 6.957, 3.638, 2.640). Since these abundances are not widely different from the solar values, є(He, C, N, O) = (10.93, 8.39, 7.78, 8.66), the corresponding T-Pg structures are not widely different either, as illustrated in Fig. 22 for models of T_eff = 1000 K and 2500 K. One sees that the solar and the Jupiter composition models span approximately the same temperature and gas pressure ranges, with the (hot, cloud-free) Jupiter composition models having only slightly higher/lower gas pressures for given temperature in the warmest/coldest of the models presented in Fig. 22, caused by the only slightly different (and temperature dependent) opacities from the two different compositions. Jupiter itself has log(g) close to 3.5 rather than the solar value of (close to) log(g) = 4.5 that is illustrated in Fig. 22. The figure is meant to illustrate the effect of chemical composition by isolating the elemental parameter differences. A main difference between the temperature response of an object with a Jupiter-like (C/O>1) and a solar-like (C/O<1) elemental abundance pattern is hidden in the different temperature dependence of dominating carbon molecules and oxygen molecules, as was illustrated for TiO, H₂O, and CH₄ in Fig. 20. In Fig. 22, this effect is reflected in which one of the two models (Jupiter-like or solar-like) that is the warmest at given gas pressure for the two values of T_eff presented.

The epsilon values of Earth’s atmosphere, on the other hand, are widely different from the solar composition, mainly because of the presence of biology, and can be expressed as є (He, C, N, O, Ne, Ar, Kr) = (8.418, 10.319, 13.892, 13.219, 8.958, 11.670, 7.756) as converted from the Earth’s molecular abundances given in the so-called NASA standard Earth (dry) atmosphere (Williams 2024). The T-Pg structure is therefore also widely different from that of corresponding solar or Jupiter composition models. Figure 22 therefore also includes an illustration of the physical effect of such a chemical composition. The main difference of the Earth-composition models from the Jupiter or solar composition models, is a substantially higher gas pressure for given physical temperature. This is a direct effect of the high nitrogen and oxygen abundances which give rise to Earth-composition models being dominated by dipole in-active homonuclear molecules, making the dry atmosphere almost completely transparent. Together with the very low hydrogen and carbon abundances in Earth-composition models, this also makes the temperature gradient as function of height very sensitive to small changes in the carbon and hydrogen abundances, which is well known as the Earth’s extreme greenhouse sensitivity to tiny changes in carbon and/or hydrogen abundances. The much higher temperature in the bottom layers of Earth composition models seen in Fig. 22 is not a greenhouse effect, but just an artefact of the way Rosseland optical depth is defined, the physical important difference between the models being the changing gas pressure values at given temperature, which defines the changing transparency (spectrum) of the models.

The gravity g of the Sun in units of cm/s² is close to log(g)=4.5, Jupiter has log(g)≈3.5, while Saturn and Earth both have log(g)≈3.0. A change of one order of magnitude in the value of g will give rise to a change in P_g of a similar amount as the differences seen in Fig. 22 due to the abundance difference between Earth-like and Jupiter/solar-like compositions. However, while the changing abundances will give rise to substantial changes in the resulting spectrum, the changing value of log(g) will in most cases hardly be visible in the spectrum, as long as the opacity carrying species remain unchanged. This is a well known and well described effect for cool stars, because it makes it difficult to estimate the mass of cool stars based on their spectra (as opposed to hotter stars). The reason being that in a self-consistent energy balance, the amount of energy absorbed in a column of atmospheric material (the ‘integrated opacity’) depends on the flux (the effective temperature) and not the gravity, such that a change in the gravity will adjust the height structure until the absorbed energy reflects the flux, and hence to a first approximation, it will reach the same spectrum of the opacity dominated species at the top of the atmosphere independent of log(g). As already discussed above, the spectra of the trace elements will, however, be widely different in two such self-consistent models of different values of log(g) and therefore immediately reveal, for example, nitrogen or oxygen biological features even if oxygen itself would not be visible, provided the energy balance in the analysing model atmosphere is computed self-consistently. In exoplanetary spectra we will probably have to await large enough instruments, such as the ELT, to be able to identify the spectral features from trace molecules in the atmosphere.

Fig. 22 Temperature versus gas pressure structure of models with Teff = 1000 K (black) and 2500 K (red), with respectively solar, Jupiter- like, and Earth-like chemical atmospheric compositions. All models are for log(g) = 4.5.

12 Outlook and conclusions

We have spent a considerable amount of effort throughout this paper to investigate the quality of the assumptions that were used in the previous versions of the MARCS code, by comparing results from computations where these assumptions were included with the corresponding MSG models where it is now possible to relax these assumptions. The logic behind these comparisons was to secure consistency with the previous MARCS grids and to study the effects of various limitations that are no longer needed (for example due to the existence of more advanced input data or due to new computational methods). We hope that this detailed comparison will also help the huge community of users of previous versions of the MARCS code (judged from the more than 3000 citations to the code in the literature) as to whether some specific analyses and results may or may not deserve revision.

By comparing the models based on opacities from the limited number of molecules that were available in 2008 with models based on a more extensive number of molecules now available in the ExoMol database, we conclude that the limitation to including only 16 molecules in the classical 2008 grid (and 22 in the 2017 DRIFT-MARCS computations) did not affect the accuracy of the models, even for the coolest models of the 2008 grid of T_eff = 2500 K (Figs. 1 and 23). We further conclude that it has been valid to exclude atomic line opacities in any of the applications (that we know of) where it was ‘just intuitively’ assumed to be acceptable (i.e. typically for solar metallicity models of T_eff ≈ 3500 K and cooler). Both of these conclusions could of course be challenged for compositions deviating substantially from solar composition and not included in the tests performed in the present study. For the sake of consistency and continuity, all MSG models made available do include the full set of molecular species listed in this paper as well atomic line opacities, even when some of it is expected to be without importance for the structure.

The intuitive assumption in all MARCS models that it was sufficient to use Gaussian line profiles in the computation of the molecular opacities was, however, challenged. One might have assumed that the uncertainty in the model structure of older models due to lacking molecular line profile data could be simulated by introducing models with ξ = 2 km/s and 5 km/s in the Gauss profile, which might however underestimate the true uncertainty. We found that the difference between models based on line opacities from Gaussian profiles of ξ = 3 km/s and based on extensive Voigt profile datasets now available was larger than the difference between models based on the ξ = 2 km/s and ξ = 5 km/s in the molecular Gaussian line profiles. However, the presently available Voigt profile based opacity data for molecules in the ExoMol database that we tested here were computed only for ξ = 0, so more work is needed on this issue in order to reach a final conclusion. This uncertainty goes of course not only for MARCS and MSG models but for all other self-consistent exoplanetary models in the literature as well.

We have focused in the present paper on a description of the basic principles that have made it possible to expand the previous grids of well-cited MARCS models to lower effective temperatures, to a greater variety of chemical compositions, to irradiated models with cloud formation, and to models of time-dependent chemistry that are necessary for a self-consistent description of the effect of biological activity on exoplanetary atmospheres and spectra. In the present paper, we have given examples of these models, and we illustrate in Fig. 23 the rough changes in T − P_g structure through the new grid of the basic models. In further MSG models, to be added to the present basic grid, we will present models that include larger ranges in T_eff , log(g), chemical elemental compositions, stellar irradiation, surface structures, as well as clouds of additional nucleation and dust species.

The code in its present form is able to self-consistently compute the host star atmospheric structure and radiation field of a wide range of stars and sub-stellar objects. At the same time it can calculate the atmosphere and spectrum of the associated system of planets from Earth-like temperatures to even the hottest clear and cloudy gas-giant exoplanets. It can include structural and spectral effects of chemical non-equilibrium, including those that can be suspected to arise from biological activity. Future versions of the grid will include a wider range of atmospheric models for rocky exoplanets of various surface compositions, to facilitate studies of different kinds of formation and outgassing scenaria and biological activity patterns.

Fig. 23 Temperature versus gas pressure of basic MSG models from Teff = 400 K to 2500 K (plus the coolest classical MARCS model for comparison as well).

Acknowledgements

We acknowledge funding from the European Union H2020- MSCA-ITN-2019 under grant no. 860470 (CHAMELEON) and from the Novo Nordisk Foundation Interdisciplinary Synergy Programme grant no. NNF19OC0057374. We are thankful for many detailed and insightful comments from an anonymous referee.

References

All Tables

All Figures

Fig. 1 Comparison of the classical 2008 MARCS model at the lower boundary of the 2008 grid (labeled ‘G2008-MARCS’) with models computed with the MSG code and based on various levels of completeness and choices of input data as described in the text. Models labeled ‘MSG 22…mol..’ are based on molecular opacities from the 22 molecules included in the DRIFT-MARCS computations, while models labeled ‘MSG 16mol…’ are based on opacities from the 16 molecules included in the classical MARCS grid. Upper panel show the T–Pg structure of the full models from log(τross) = –5.0 to +2.5 while the lower panel zooms in on the upper part of the atmosphere.

In the text

Fig. 2 Upper panel shows the T–P𝑔 structure and lower panel the T- τ structure from log(τross) = +2.5 to –5.0 for models with Teff from 2500 K to 500 K. For each Teff, we show a classical MARCS model with 22 molecules in the opacity calculation as well as the models from the present grid (i.e. 50 molecules plus atoms included in the opacity), except for the Teff =1000 K and Teff =500 K models where it has not been possible to converge the models with only 22 molecules.

In the text

Fig. 3 Electron pressure versus optical depth for models with Teff from 2500 K (upper curves) to 500 K (lower curve), computed from log(τross) = +2.5 to –5.0. For the warmer of the models are shown both the results based on including only the original 22 molecular opacities in the equilibrium computation as well as the full set of the 50 molecules, as discussed in the text. For the coolest models only those with opacities from all the 50 molecular opacities adopted in the MSG models could converge.

In the text

	Fig. 4 Models with Teff = 1000 K computed based on the ‘old’ way (labeled ‘…Pe1K4’) of estimating the electron pressure as well as the new way with Ij equal to 2.339 eV (labeled ‘…Pe2K2’) and 4.339 eV (labeled ‘…Pe2K4’), respectively. The upper and lower panels respectively show, as a function of log(τross), Pe and Pg.
In the text

	Fig. 5 Absorption coefficient of the OS lines as a function of wave number, computed from the ExoMol line list under the assumption of Gaussian molecular line profiles (blue) and Voigt profiles (light brown), respectively, and for two different values of temperature (T = 1000 K and 2900 K) and two different values of Pg (10−8 and 10−3 dyn/cm2).
In the text

	Fig. 6 Rosseland mean opacity as function of gas pressure, for models with Teff from 2500 K (upper curves) to 1000 K (lower curves), based on the full set of 52 molecular opacities (full drawn lines) and the original 22 molecular opacities, respectively, in the radiative transfer. Colour coding correspond to the effective temperature of the models, as shown in the figure legends.
In the text

	Fig. 7 Logarithm of the sum of the partial pressures (minus the partial pressure of H2 , and normalised to the gas pressure) of respectively the original 22 opacity bearing molecules (dotted lines) and the 22 plus the new 30 ones (dashed lines), for the new models presented here of Teff = 2500 K (black lines), 2000 K (red lines), 1500 K (green lines), and 1000 K (blue lines).
In the text

Fig. 8 Distribution of the total gas pressure onto its various components for a model of Teff =2500 K. As indicated with the legend, the curves show the partial pressures of the molecules H2 , CO, H2O, CH4 , NH3 , TiO, and OH, the neutral H and He atoms, the sum of all positively charges species (“pel+”), all negatively charged species (“pel-”), sum of the partial pressures of all neutral atoms together (“pnat”), all neutral molecules (“pnmo”), and finally the sum of all the opacity bearing molecules in the classical MARCS models (“p24m”) as well as in the new models presented here (“p54m”). We note that the curves for H2 and all neutral molecules (pnmo) as expected plot on top of one another, as do the curves p24m and p54m (therefore p24m/p54m show as a single full drawn green curve, and the only red curve is for the partial pressure of TiO). The upper black curve is for H and the lower for NH3 .

In the text

	Fig. 9 Same as Fig. 8 but for at model of Teff =1000 K. We note that for this model the curves p24m and p54m do deviate, in particular for the upper layers, illustrating the importance of the additional molecular opacities in the MSG models compared to the classical MARCS models.
In the text

	Fig. 10 Temperature versus gas pressure structure of our irradiated 1.5D model simulation of an exoplanet with the parameters of WASP-39b at the sub-stellar point (right-most curve), at Θ = 15° and at Θ = 30° (leftmost curve).
In the text

	Fig. 11 Temperature inversion caused by TiO and VO at the sub-stellar point in irradiated MSG 1.5D model simulations of WASP-43b.
In the text

Fig. 12 Effects of clouds in a Teff =1500K model atmosphere. Upper panel: pressure–temperature profiles for MSG models at Teff =1500K and log(g)=4.0 respectively with clouds (labeled MSG) and without clouds (labeled MARCS). We note that the convective part of the atmosphere is represented with a thicker line. Lower panel: average cloud particle size (solid curve) and cloud particle number density (dashed curve) for the MSG model shown in the upper panel. We note that the spike in the number density at the bottom of the cloud is numerical and due to the rapid decrease of the average cloud particle size.

In the text

	Fig. 13 Cloudy MSG model at an effective temperature of 1500 K and log(g)=4.0. Upper panel: cloud volume fraction of the 12 cloud species present in the atmosphere. Lower panel: relative elemental depletion in the gas phase of elements affected by the cloud formation.
In the text

	Fig. 14 Relative number density of important molecular and atomic gas species in the cloudy MSG (solid curves) and cloud-free MSG (dashed curves labeled MARCS) models at 1500 K and log(g)=4.0 (shown in Figure 12).
In the text

	Fig. 15 Cloudy (blue) and cloud-free (grey) synthetic spectra (binned to a rsolution of 1000) of a self-luminous body with an effective temperature of 1500 K and log(g)=4.0, with solar metallicity and solar C/O ratio.
In the text

	Fig. 16 Example of the temporal evolution of the O, O2 , and O3 number densities from their equilibrium abundances (left) to steady-state (right) for a 4-reaction Chapman mechanism chemical network.
In the text

Fig. 17 Local thermodynamic equilibrium and non-equilibrium steadystate partial pressures of O3 as a function of total gas pressure (shown as full-drawn and dotted lines, respectively) for a surface-less exoplanet with Earth-like atmospheric elemental abundances, log(g) = 3.0 and Teff = 1000 K, illuminated by a Solar-type star (of log(g) = 4.4, Teff = 5772 K, and solar radius) at 1 AU from the planet.

In the text

	Fig. 18 Atmospheric mixing ratios of the molecules for which quenching has been implemented into MSG. Solid lines show the equilibrium values, while the dotted lines show the mixing ratios with a fixed eddy diffusion coefficient of Kzz = 1010 cm2/s. The model has Teff = 2500 K and log(g) = 4.5.
In the text

	Fig. 19 Absorption coefficient at T = 1500 K of the most important molecules (listed in the legend) for the structure and understanding of the coolest dwarfs.
In the text

	Fig. 20 Partial pressures of the most important molecules in models of Teff = 1500 K, 1000 K, and 500 K, respectively.
In the text

Fig. 21 Development of the spectrum of H2O and TiO (the two left panels), CH4, and NH3 (middle panels) and FeH and VO (right panels) for models of Teff = 1500 K (red), 1000 K (blue), and 500 K (green). The dashed curves are the true continua (i.e. the spectrum computed by allowing only the continuum opacities into the spectrum computation, but all opacity sources into the model computation), while the full drawn lines likewise are the spectra of each individual molecule listed in the legend, computed by allowing the opacity from only that molecule in the spectrum computation, but all opacity sources into the model computation. Remark that for the coolest models (and the plotted resolution) the combined spectra do not at any place reach even close to the (true) continuum, which is the continuum that should be the reference in abundance analyses.

In the text

	Fig. 22 Temperature versus gas pressure structure of models with Teff = 1000 K (black) and 2500 K (red), with respectively solar, Jupiter- like, and Earth-like chemical atmospheric compositions. All models are for log(g) = 4.5.
In the text

	Fig. 23 Temperature versus gas pressure of basic MSG models from Teff = 400 K to 2500 K (plus the coolest classical MARCS model for comparison as well).
In the text