#  tlusty and synspec Users’s Guide IV: Upgraded Versions 208 and 54

I. Hubeny111University of Arizona, Tucson; USA<EMAIL_ADDRESS>C.
Allende Prieto222IAC, La Laguna, Tenerife, Spain<EMAIL_ADDRESS>Y.
Osorio333IAC, La Laguna, Tenerife, Spain<EMAIL_ADDRESS>and T.
Lanz444Observatoire de Cote d’Azur, Nice, France<EMAIL_ADDRESS>
###### Abstract

We present a brief description of the newly upgraded versions of the programs
tlusty, version 208, and synspec, version 54. tlusty is used to generate model
stellar atmospheres or accretion disks, and synspec produces detailed
synthetic spectra and/or opacity tables. This paper complements published
guides that provide a detailed description of earlier versions of the codes,
tlusty205, and synspec51. The main upgrades include the flexible construction
of opacity tables in synspec, and their use in producing hybrid models with
tlusty, with important species treated in NLTE, while the bulk of opacity of
atomic and molecular lines and continua are considered in LTE using a pre-
calculated opacity table.

There is also a number of additional changes and upgrades that increase the
versatility and flexibility of these codes.

###### Contents

  1. 1 Introduction
  2. 2 TLUSTY
    1. 2.1 Outline of opacities and chemical abundances
    2. 2.2 Working with opacity tables
      1. 2.2.1 Basic concepts
      2. 2.2.2 Reading and initialization of the opacity table
      3. 2.2.3 Fine points of computing NLTE models with opacity tables
    3. 2.3 New format of the opacity table
    4. 2.4 Setting frequency points
    5. 2.5 Line broadening
      1. 2.5.1 Hydrogen lines
      2. 2.5.2 Metal lines
    6. 2.6 Equation of state
    7. 2.7 Convection
    8. 2.8 Collisional rates
    9. 2.9 Fixed-density models
    10. 2.10 Changes of the output
      1. 2.10.1 Standard output
      2. 2.10.2 New output file fort.22
      3. 2.10.3 Specific intensities
  3. 3 SYNSPEC
    1. 3.1 Changes in the line list
      1. 3.1.1 Atomic line lists
      2. 3.1.2 Molecular line list
      3. 3.1.3 Specification of the format of the line lists
      4. 3.1.4 Utility programs
    2. 3.2 Constructing opacity tables
      1. 3.2.1 Input file fort.2
      2. 3.2.2 Input files fort.55 and fort.3
      3. 3.2.3 Output files
      4. 3.2.4 The utility program ngrid
    3. 3.3 Equation of state
    4. 3.4 Additional opacities
    5. 3.5 Rejection of weak lines from the line list
    6. 3.6 Updated treatment of NLTE levels
      1. 3.6.1 Modification of the atomic line list
      2. 3.6.2 Modification of the atomic data file
      3. 3.6.3 Additional input parameters
      4. 3.6.4 Input NLTE level populations
    7. 3.7 Subtle points of computing line and continuum opacity
      1. 3.7.1 Semi-explicit species
      2. 3.7.2 Quasi-explicit species
      3. 3.7.3 Summary of treating chemical species
    8. 3.8 Summary of treating line and continuum opacities by synspec
      1. 3.8.1 Input LTE model atmosphere; with no additional input
      2. 3.8.2 Input NLTE model atmosphere
      3. 3.8.3 LTE model atmosphere, with additional (NLTE) input
  4. 4 Graphical user interface programs
    1. 4.1 Plotting output from tlusty
      1. 4.1.1 IDL
      2. 4.1.2 Python
    2. 4.2 Plotting output from synspec
      1. 4.2.1 Synplot – an IDL wrapper
      2. 4.2.2 Synple – a Python wrapper
  5. 5 Downloading and installing the package
    1. 5.1 Main file
    2. 5.2 Additional files
      1. 5.2.1 Line lists
      2. 5.2.2 Iron data
      3. 5.2.3 Opacity tables
    3. 5.3 Summary of implementation
  6. 6 Compiling, testing, and running TLUSTY
    1. 6.1 Source code files and compilation
    2. 6.2 Test cases for tlusty
      1. 6.2.1 Standard model atmospheres and disks
      2. 6.2.2 Tests involving an opacity table
  7. 7 Compiling, testing, and running SYNSPEC
    1. 7.1 Source code files and compilation
    2. 7.2 Test cases for synspec
      1. 7.2.1 Tests of constructing an opacity table
      2. 7.2.2 Ordinary stellar synthetic spectra
      3. 7.2.3 Synthetic spectrum for a hybrid Opacity Table/NLTE model atmosphere
  8. 8 Summary and outlook
  9. Acknowledgements
  10. Appendix A - Adopted sets of the standard solar abundances
  11. Appendix B- Mean molecular weight and total particle density
  12. Appendix C - Thermodynamic parameters and convection
  13. Appendix D - Improved upper boundary condition for the radiative transfer equation
  14. Appendix E - Treatment of hydrostatic equilibrium within the Rybicki scheme
  15. Appendix F - Working with EXOMOL line broadening parameters
  16. Appendix G - Reducing the number of lines in a molecular line list
  17. References

##  1 Introduction

The programs tlusty, version 205, and synspec, version 51, were described
comprehensively in a series of three papers by Hubeny & Lanz (2017a, b, c;
hereafter referred to as Paper I, II, and III, respectively). In this paper,
we document recent changes in these codes. The user is strongly encouraged to
consult the previous parts of the Guide, at least Paper I, since without them
the present description cannot be fully understood.

##  2 TLUSTY

There are several basic differences between the newly developed version 208
and the previously released version 205:

  * •

A modified and upgraded treatment of opacities and chemical abundances.

  * •

Upgrades for working with opacity tables. and related changes in the treatment
of explicit atoms, levels, and transitions.

  * •

An update of the molecular equation of state, and related changes in the
format of the model output.

  * •

An improved treatment of convection.

  * •

A new and upgraded treatment of collisional rates.

  * •

An update of the Rybicki scheme for solving the linearized structural
equations, which now works for accretion disks, and even works well for NLTE
models.

  * •

A “fixed-density” option, useful for solving NLTE line formation and
ionization balance for snapshots of independent hydrodynamical calculations.

These changes are most relevant for computing model atmospheres of cool stars,
while the last change may be applicable to all types of stellar atmospheres or
accretion disks.

In addition, there are various other small changes of the layout of the
standard output, which we will briefly describe in § 2.10.

###  2.1 Outline of opacities and chemical abundances

There are three classes of opacity sources considered by tlusty, namely:

Explicit opacities

– bound-free and bound-bound transitions between explicit levels, and free-
free transitions of explicit ions. They are computed on the fly, and these are
the only ones that can be considered in NLTE.

Additional opacities

– any other opacity source that is computed on the fly. This option existed
from the very early versions of tlusty, but has been upgraded and made less
error-prone in the present version. This is described below.

Opacity from a pre-calculated table

– This option was introduced in version 205, and has been upgraded in the
present one, as described in this document.

Historically, the additional opacities, in the above meaning, were the
following: H- bound-free and free-free opacity, H${}_{2}^{+}$ bound-free and
free-free opacity, hydrogen Rayleigh scattering opacity, plus an outdated
option for including bound-free opacity of additional, higher, levels of He I
and He II. The corresponding switches. activated in the non-standard (keyword)
parameter file, are called IOPHMI, IOPH2P, IRSCT, IOPHE1, and IOPHE2. The
default values of these switches were set to 0, so these sources were not
considered unless specifically required by setting the corresponding switch to
a non-zero value, typically 1.

Table 1: Additional opacities - keywords and their values

Opacity source | Ref | Keyword | Default | Set automatically to 0 if
---|---|---|---|---
H- b-f and f-f | a | IOPHMI | 1 | H- is already an explicit ion
H${}_{2}^{+}$ b-f and f-f | a | IOPH2P | 1 | IFMOL=0 or $T>$ TMOLIM
H${}_{2}^{-}$ f-f | b | IOPH2M | 1 | IFMOL=0 or $T>$ TMOLIM
He- b-f and f-f | a | IOPHEM | 1 |
CIA H2-H2 | c | IOH2H2 | 1 | IFMOL=0 or $T>$ TMOLIM
CIA H2-He | e | IOH2HE | 1 | IFMOL=0 or $T>$ TMOLIM
CIA H2-H | f | IOH2H | 1 | IFMOL=0 or $T>$ TMOLIM
CIA H-He | d | IOHHE | 1 | IFMOL=0 or $T>$ TMOLIM
CH opacity | a | IOPCH | 1 | IFMOL=0 or $T>$ TMOLIM
OH opacity | a | IOPOH | 1 | IFMOL=0 or $T>$ TMOLIM
H Rayleigh scat. | a | IRSCT | 1 |
He Rayleigh scat. | a | IRSCHE | 1 |
H2 Rayleigh scat. | a | IRSCH2 | 1 | IFMOL=0 or $T>$ TMOLIM

References: (a) Kurucz (1970); (b) Bell (1980); (c) Borysow et al. (2001); (d)
Gustafsson & Fromhold (2001); (e) Jorgensen et al. (2000); (f) Gustafsson &
Fromhold (2003).

We have now deleted the outdated switches IOPHE1 and IOPHE2, and added more
additional opacity sources. The new opacity sources include the H${}_{2}^{-}$
free-free opacity, and four different collision-induced opacity (CIA) sources,
namely those for the H2–H2, H2–He, H2–H, and H–He collisions – see Table 1.
The default values of all of them are set to 1, that is, all sources are
considered unless specifically excluded by setting the corresponding switch to
0. In this way, one avoids accidentally omitting a potentially important
opacity source, e.g. the H${}_{2}^{+}$ or Rayleigh scattering opacity for cool
stars.555This was even done in some test cases presented in Paper III, where
Rayleigh scattering was omitted despite it provides a non-negligible
contribution to the total opacity.

Another potential source of error was the treatment of the H- opacity. If one
sets IOPHMI=1, and at the same time H- is specified as an explicit ion, its
opacity is included twice. It was entirely up to the user to check the
internal consistency of the input parameters on this regrad. The situation
when using opacity tables is even more delicate, in the sense that some of
these additional sources may have already been considered there.

Therefore, we have decided to make the overall treatment of additional
opacities much less error-prone by letting tlusty check and automatically
correct possible inconsistencies. Table 1 summarizes the relevant switches and
corrective actions.

Important note: In all cases, a particular keyword is reset to 0 if the given
opacity source was already considered in the opacity table, unless a new
keyword parameter KEEPOP is set to a nonzero value. This option is included
for pedagogical purposes, or for investigating what would be the error
introduced in a model by considering a particular opacity source twice.

We also stress that when synspec works in the table-generating mode, the
additional Rayleigh scattering opacity sources are disabled, because the
resulting opacity table should only consider the true absorptive opacity.

We have also increased a flexibility in specifying the chemical abundances of
elements. As explained in Paper III, §4̇.2, the chemical abundances are
specified in the second block of the standard input to tlusty or synspec,
namely through the parameter ABN. The second block of the standard input has
now the following structure:

The first line contains one or two numbers:

NATOMS

– as before, the highest atomic number of an element that is considered
(explicitly or non-explicitly).

IABSET

– if set, it specifies the adopted set of the standard solar abundances to be
used. If it is not set, the program will assign it to 0. The individual sets
are the following:
$\bullet\,=0$ – the adopted set of the standard solar abundances is that of
Asplund et al. (2005), which is essentially the same as that of Grevesse et
al. (2007) used in the recent large MARCS (Gustafsson et al. 2008) and Kurucz
(Meszaros et al. 2012) model atmosphere grids.
$\bullet\,=1$ – the adopted set of the standard solar abundances is a
combination of meteoritic/photospheric abundances from Asplund et al. (2009)
adopted by Husser et al. (2013) in their extensive grid of PHOENIX model
atmospheres and synthetic spectra.
$\bullet\,=2$ – the adopted set of the standard solar abundances is that of
Grevesse & Sauval (1998), which was the default in the previous versions of
tlusty and synspec.
$\bullet\,$ – any other number: IABSET is reset to 0.

Next, there are NATOMS records, each containing three items: MODE, ABN, and
MODPF. The parameter MODE will be described in § 2.2.1, and MODPF is now
outdated, so it can be set to 0. Here we describe the parameter ABN. Its
meaning is essentially the same as in the previous versions, the only
difference is that instead of dealing with a hardwired set of solar
abundances, the parameter ABN now relates the intended abundances to the
adopted standard solar abundances, specified through the parameter IABSET.

ABN

a specification of the abundance of the given species:
$\bullet\,=0$ – the standard solar abundances is assumed;
$\bullet\,<$ 0 – a non-standard abundance is assumed, abs(ABN) has now the
meaning of the abundance expressed as a multiple of the standard solar
abundance (i.e. $-0.1$ means 1/10 of the standard solar abundance, $-5$ means
5 times the standard solar abundance, etc.);
$\bullet\,>$ 0 – a non-standard abundance is assumed, expressed as $N({\rm
elem})/N({\rm ref})$, i.e. relative by number to the reference species. The
reference atom is H by default, but the reference element can be changed by
means of the optional parameter IATREF (see Paper III, § 7.4.4).
$\bullet\,>10^{6}$ – non-homogeneous (depth-dependent) abundance is assumed.
In this case, the immediately following $N\\!D$ lines should be added that
contain the individual values of the abundance (relative to hydrogen by
number), for all depth points $d=1,\ldots,N\\!D$.

For convenience, we list the adopted values of the chemical abundances for all
the three standard solar abundances sets in Appendix A. Although using more
recent sets 0 and 1 of Asplund et al. (2005, 2009), and Grevesse et al. (2007)
is preferable, we include the old set of Grevesse & Sauval (1998) for
comparison with earlier results.

###  2.2 Working with opacity tables

####  2.2.1 Basic concepts

In version 205, one could either use an opacity table describing the total
opacity, or not use an opacity table at all. The present version introduces a
third option, already mentioned in Papers I and III, of working with partial
opacity tables. This option consists in using an opacity table that contains
opacities for a subset of chemical species, while the opacity of other species
is considered explicitly, typically in NLTE (they can be considered in LTE,
but this possibility does not offer any advantages, since in such a case all
species could have been included in the opacity table). An example of using a
partial opacity table is, for instance, the construction of a solar-type model
atmosphere with H, Mg, and Ca in NLTE, while including the opacity from all
other species, atomic as well as molecular, in the opacity table. Such a table
will be referred to as a partial opacity table. We call this approach to
compute model atmospheres a hybrid opacity-table/NLTE (OT/NLTE) approach.

The basic control keyword parameter is IOPTAB, already described in Paper III
§ 7.4.1. For convenience, we repeat its description here:

IOPTAB

– a switch that specifies the use of pre-calculated opacity tables:
$\bullet\,=0$ – classical option – no pre-calculated opacity table is used;
$\bullet\,<0$ – all opacities are included by means of a pre-calculated
opacity table. In this case, one does not specify any explicit atoms, ions,
and levels. Only LTE models can be computed with this option. There are two
possibilities here:

  * •

=-1 – a full opacity table is used, but the equation of state and
thermodynamic parameters are computed on the fly.

  * •

=-2 – in addition to the opacity, the equation of state and thermodynamic
parameters are given by pre-calculated tables.

$\bullet\,>$ 0 – hybrid option, using a partial opacity table – one still
selects explicit atoms, ions, and levels, for which the opacity is computed on
the fly, while for the remaining species one uses an appropriate pre-
calculated opacity table. This option was not operational in version 205, but
it is in the present version 208.
DEFAULT: IOPTAB=0

IFRAYL

– this parameter takes effect only when using opacity tables in the full
opacity table approach, that is if IOPTAB $<0$. Since all the (thermal)
opacity is included in the table, and opacity is not calculated on the fly,
one needs to make sure that the scattering opacity is included as well.
Possible values of IFRAYL are:
$\bullet\,>0$ – the Rayleigh scattering opacity is computed on the fly, based
on the values of the Rayleigh scattering parameters IRSCT, IRSCHE, IRCSH2 for
the scattering on H, He, and H2, respectively. Since these parameters have
default values of 1 (the source is included), the Rayleigh scattering opacity
is included by default unless it is explicitly rejected.
$\bullet\,<0$ – the Rayleigh scattering opacity is read from a special table
rayleigh.tab. [This option is implemented, but the table rayleigh.tab is not
set up, so this option cannot be used at the present time.]
$\bullet\,=0$ – Rayleigh scattering is neglected completely (which is wrong,
but the option may be used to assess the importance of this scattering
mechanism).
DEFAULT: IFRAYL=1

There are some subtle points here. In the case of a full opacity table, the
frequency points selected for the run of tlusty are either the same frequency
points as those used in the table, or a different set but chosen to be within
the range covered by the table. Let $\nu_{\rm min}$ and $\nu_{\rm max}$ be the
lowest and highest frequencies in the table. In the hybrid option, one may
consider frequencies that are higher than $\nu_{\rm max}$. For example, an
opacity table used for solar-type stars has typically the maximum frequency
around the Lyman limit, $\nu_{\rm max}\approx 3.28\times 10^{15}$s-1, while
the ionization frequency of some ions treated explicitly may be higher. For
instance, when magnesium is treated explicitly,, one would include Mg I, Mg
II, and the ground state of Mg III. The ionization frequency for the ground
state of Mg II is $\approx 3.6\times 10^{15}\,$s-1, so in order to evaluate
the radiative rate for the Mg II ionization, needed for NLTE calculations of
the Mg level populations, one needs to specify the background opacity going at
least up to $\approx 7\times 10^{15}$s-1.

To this end, we have introduced two new options666Plus two more options,
MODE=4 and 5, specific to synspec  – see § 3.7.1 and 3.7.2 for treating the
chemical species, controlled by the parameter MODE, which is contained in the
block of data for chemical elements of the standard input, described in Paper
III, § 4.2. Here is its description:

MODE

– a specification of the mode of treatment of a given species:
$\bullet\,=0$ – the element is not considered at all (is completely ignored).
$\bullet\,\not=0$ – regardless of the actual value, the element is allowed to
participate in the equation of state, namely to contribute to the total number
of particles, to the total charge, and to the formation of molecules.

  * •

$=1$ – the element is treated implicitly. In this case, its opacity is not
computed on the fly, i.e. it is either neglected (in the classical approach
without an opacity table; i.e., with setting IOPTAB=0), or is contained in the
opacity table (IOPTAB$\not=0$).

  * •

$=2$ – the element is treated explicitly, i.e., selected energy levels of the
included ionization states are considered explicitly; i.e., their populations
are determined by solving the corresponding kinetic equilibrium equations. The
opacity of this element is computed on the fly, using the current level
populations.

  * •

$=3$ – the element is viewed as explicit, but its opacity is allowed to
contribute only for frequencies $\nu>\nu_{\rm max}$, i.e., beyond the range of
the opacity table. If the opacity table is not used (i.e., IOPTAB=0), then
$\nu_{\rm max}$ is set to automatically to $0$, so that in this case there is
no difference between modes 2 and 3.

  * •

$=4$ – the element is viewed as semi-explicit – see § 3.7.1. This option takes
effect only for synspec; for tlusty it is equivalent to MODE=2, i.e., the
normal explicit mode.

  * •

$=5$ – the element is viewed as quasi-explicit – see § 3.7.2. This option
takes effect only for synspec; for tlusty it is equivalent to MODE=2, i.e.,
the normal explicit mode.

  * •

$=-2$ – the element is viewed as explicit, in the sense that its opacity is
computed on the fly, but the populations of all explicit levels of all its
explicit ions are held fixed to the values from the input model atmosphere.

  * •

$=-3$ – analogous to option with MODE=3, plus the populations of all its
levels are held fixed.

The last two options allow a correct evaluation of the total opacity beyond
the frequency range of the opacity table. For instance, if a current run
intends to consider Mg in NLTE, and the hydrogen opacity is part of the
opacity table, then one sets MODE=2 for Mg, and MODE=3 for H. In this case,
the H opacity will be considered only for $\nu>\nu_{\rm max}$, which will be
satisfactory because it is the dominant source of the background opacity in
the Lyman continuum. Since H is being set as an explicit element, its level
populations will also be treated in NLTE. One can also set MODE=$-3$ for H, in
which case the opacity is evaluated as before, but the H level populations are
being held fixed to the input values (LTE or NLTE, depending on the input
model).

Finally, we remind the reader that the only way tlusty can consider molecular
line opacity is through the opacity table. There is no mechanism in tlusty to
compute molecular line opacity on the fly, or treating molecular level
populations in NLTE.

####  2.2.2 Reading and initialization of the opacity table

This is controlled by the specific input file fort.15. The file contains one
or two lines:
First line (mandatory): OPTABLE, IBINOP, where:

OPTABLE

– a character string that contains the filename of the opacity table.

IBINOP

– a switch indicating the global format of the table:
$\bullet=0\,$ – the table is stored in text format;
$\bullet\not=0\,$ – the table is stored in binary format.

Second line (optional): ISTEPT, ISTEPR, ISTEPF:

ISTEPT

– the step for taking into account the tabulated values of temperature. If
ISTEPT=1, all the temperature values are considered; if ISTEPT=2, every second
temperature value is taken into account, etc.

ISTEPR

– an analogous step for density;

ISTEPF

– an analogous step for the wavelength (frequency) grid.

If these step parameters are not specified, they are set to their default
values of ISTEPT=ISTEPR=ISTEPF=1.

####  2.2.3 Fine points of computing NLTE models with opacity tables

It should be stressed that there are two basic types of OT/NLTE hybrid models.
A simple one is to take the structure from a previously computed model
atmosphere, and to solve for NLTE line formation of selected species, while
keeping the atmospheric structure ($T$ and $N$) fixed – the so-called
restricted NLTE problem. Another possibility is to allow the atmospheric
structure to change as well, thus studying the influence of NLTE effects in
one or several species on the global atmospheric structure. In the latter
case, it is recommended to use the Rybicki scheme, set by the keyword
parameter IFRYB=1 or 2. In the present version of the code, the default
scheme, the hybrid CL/ALI method, does not work properly for this purpose.

The numerical performance of each method is demonstrated in § 6.2.2 on two
simple examples; both use an LTE model generated with a full opacity table as
a starting one, namely:
– g55nres – computes a restricted NLTE model; that is a model in which the
structure (temperature, density, electron density) is held fixed, and one
computes NLTE level populations of H, Mg, and Ca levels, using a reduced
opacity table that excludes these elements. The model is computed using the
standard Complete Linearization/Accelerated Lambda Iteration (CL/ALI) method
(see Paper II, Paper III, and Hubeny & Lanz 1995);
– g55nryb – a similar model that computes a full NLTE model in which
departures form LTE are allowed for H, Mg, and Ca, but where also the model
structure is allowed to change as a response to NLTE effects in these species.
The Rybicki scheme (see Paper II) is used here.

The former model represents a traditional procedure to examine NLTE effects in
individual chemical elements, while the latter model is more rigorous. From a
numerical point of view, both schemes have their advantages and drawbacks. The
g55nres model linearizes the kinetic equilibrium equations for individual
level populations, alongside with the radiative transfer equation, so the
changes of the NLTE level populations are determined self-consistently with
the changes of the radiation field in the iteration process. However, since in
this case one deals with over 40,000 discretized frequency points, most of
them are treated via ALI in the CL/ALI scheme [for details, refer to Hubeny &
Lanz (1995), or Hubeny & Mihalas (2014, § 18.4)], while only some 50 points
are being linearized (frequencies close to the edges of the ground-state H I,
Mg I, Mg II, and Ca II continua, and in the centers of the strongest lines
(L$\alpha$, L$\beta$, L$\gamma$ and H$\alpha$ of hydrogen, and first resonance
lines of Mg I, Mg II and Ca II). The convergence of this model is not fast,
but acceptable. requiring about 30 iterations of the CL/ALI scheme to reach a
maximum relative change in all quantities below $10^{-3}$. This relatively
slow convergence is caused by the fact that the large majority of frequency
points are not linearized. Due to a similar reason, the Cl/ALI scheme has
presently serious convergence problems when uses to compute a fully self-
consistent model where all structural equations are being solved (that is, in
addition to the radiative transfer and kinetic equilibrium, also the
hydrostatic. charge, and radiative/convective equilibrium equations are
solved). In this particular case, the interplay between convection and the
atmospheric structure is the main source of convergence problems.

In contrast, the Rybicki scheme, which essentially linearizes the radiation
field in all frequencies, and solves the radiative transfer equation together
with the radiative/convective equilibrium equation, avoids such convergence
problems. However, the kinetic + charge equilibrium equation and the
hydrostatic equilibrium equations are not a part of the linearization process,
so that they are essentially treated by means of a Lambda iteration (that is,
these equations are being solved only after a completed iteration of the
global linearization scheme). This is not a big problem for the hydrostatic
equilibrium equation because it involves the radiation intensities only if
radiation pressure is important, which is not the case for cool star model
atmospheres. However, this may be a problem for the kinetic equilibrium, that
is for determining NLTE level populations, because it is well known that the
Lambda iteration converges very slowly, or fails to converge at all, for NLTE
level populations. For these reasons, the Rybicki scheme was never applied to
compute NLTE level populations or NLTE model atmospheres.

However, tlusty does not use a traditional Lambda iteration in the formal
solution step (that is, just alternating between solving for NLTE populations
with the radiation field fixed, and solving the radiative transfer equation
with the level population fixed), but instead uses the ALI-based
preconditioning scheme of Rybicki & Hummer (1991, 1992); see also Paper II and
Paper III for additional numerical details. To make the action of this scheme
efficient, one needs to set the number of the “Lambda” iterations, the keyword
parameter NLAMBD, to a relatively large value, such as 8 as in the present
test case. This leads only to a modest increase in computing time, while the
accuracy of the NLTE level populations improves dramatically. If this is not
done, and the parameter NLAMBD is left at its default value, namely 2, the
model formally converges very fast (in 5-6 iterations), but the resulting NLTE
level populations may still be quite different from the exact solution (0.5,
up to even 0.7 dex in some cases). In contrast, when using NLAMBD$=8$, and
also to be on the safe side by decreasing the convergence criterion CHMAX to
$10^{-4}$ or $10^{-5}$, the NLTE level populations are indistinguishable from
those obtained by the CL/ALI scheme (obviously, for the temperature/density
structure being held at the structure determined by the Rybicki scheme).

The Rybicki scheme was thus shown to work well even for NLTE model
atmospheres. However, slight modifications were necessary, which are mostly
transparent to the user. The only visible modification is a somewhat extended
meaning of the keyword parameter IFRYB, which we give below:

IFRYB

– a switch controlling the basic iteration scheme:
$\bullet\,=0$ – a standard CL/ALI scheme is used;
$\bullet\,=1$ – a simple Rybicki scheme is used; to be used for LTE models;
$\bullet\,=2$ – an extended Rybicki scheme, to be used for NLTE models.
DEFAULT: IFRYB=0 (standard CL/ALI scheme)

###  2.3 New format of the opacity table

To make the computation of model atmospheres with opacity tables less error-
prone, we have introduced a new format for the opacity table. The old format
can still be used, but only in the full opacity table mode (i.e., with
IOPTAB$<0$), and if a new keyword parameter IOPOLD is set to a non-zero value.
In any case, the bulk of the table (that is, the actual opacity as a function
of frequency for all temperature-density pairs) is analogous in both cases;
the only difference is in the header of the table. Here we describe the new
header in detail. Any opacity table constructed with synspec54 will have the
new format; the possibility of using the old format is kept only for downward
compatibility with previous versions and existing opacity tables.

The purpose of adding new information to the opacity table is to avoid
inconsistencies between the parameters used in constructing the opacity table,
and those used for running tlusty. For instance, we can now avoid including
accidentally an opacity source twice, or forgetting an opacity source that may
be important, as outlined in the previous section.

The new format of the opacity table still comes in two flavors. The
traditional way is using NTEMP temperature values (for details, refer to §
3.2.1), and for each temperature NDENS values of density (or some other
density-like parameter, see § 3.2), with their numbers and values being the
same for all temperatures. Alternatively, one may choose different numbers and
different values of density for each temperature. We will call such tables
irregular tables. This will be explained in detail in § 3.2.1. We stress that
the format of the table, while being generated by synspec] by setting
appropriate input parameters, is completely transparent to the tlusty users.
The program recognizes which format of the table is being dealt with. However,
if the user intends to generate a graphical interface program to visualize the
data from the opacity table, these details are important, so that we will
point them out below. For better visibility, they are given in italics.

Here is the outline of the structure of the new opacity table:

  1. 1.

A table of abundances of the first 92 elements, used when generating the table
by synspec54. There are two sets of abundances listed there, namely those used
in the equation of state, and those used when computing detailed opacities. In
the case of using a full opacity table, that is with IOPTAB $\\!<\\!0$, these
two abundances are identical, while for a partial opacity table (IOPTAB $>0$)
the abundances of some species in the second column are set to 0, i.e. their
opacities are removed from the opacity table.

  2. 2.

A list of values for the keyword parameters IFMOL and TMOLIM used when
generating the opacity table.

  3. 3.

A list of additional opacity switches, referring to the opacities of H-,
H${}_{2}^{+}$, H${}_{2}^{-}$, He-, the four CIA opacities, CH, and OH. A value
of 1 means that the source was included when constructing the opacity table.

  4. 4.

As in the old format, the numbers of frequencies, temperatures, and densities
for the table.
For irregular tables, the number of densities (now formal) is negative. Then
there is an additional input record that for all temperatures specifies the
number of density values.

  5. 5.

As before, a list of temperatures and densities in logarithmic scale.
For irregular tables, there is a list of temperatures as usual, and then NTEMP
records with the individual values of the densities.

  6. 6.

A new list of electron densities (again, in logarithmic scale) corresponding
to the individual temperature-density pairs listed above.
For irregular tables, the format is analogous to that of listing the density
values.

  7. 7.

Exactly as before, for each frequency point, the value of the frequency, and
the logarithms of the opacity (per gram) for all temperature-density pairs.
The inner loop is over temperatures, and the outer loop over densities.
For irregular tables, the order of loops is reversed.

For completeness, we summarize below the new keyword parameters used to
describe the treatment of additional opacities and opacity tables:

IOPHMI, IOPH2P, IOPH2M, IOPHEM, IOPCH, IOPOH

– switches for additional opacity sources – see Table 1

IOH2H2, IOH2HE, IOH2H, IOHHE

– switches for additional, CIA, opacity sources – see Table 1

IRSCT, IRSCHE, IRSCH2

– switches for additional scattering sources – see Table 1

IOPOLD

– a switch for considering the old format for the opacity table (only for the
full-table mode, IOPTAB $<0$), by setting its value to 1.
DEFAULT: IOPOLD=0 (new format)

KEEPOP

– a switch for disabling the automatic correction of the additional opacity
switches to avoid considering opacity sources that may have already been
included in the opacity table. If this parameter is set to a non-zero value,
the input additional opacity switches IOPHMI, IOPH2P, IOPH2M, IOPHEM, IOH2H2,
IOH2HE, IOH2H, IOHHE, IOPCH, and IOPOH are not being checked for compatibility
with the values used when constructing the opacity table.
DEFAULT: KEEPOP=0 (an automatic correction is performed)

###  2.4 Setting frequency points

In previous versions of tlusty, choosing the frequency points was done in two
possible ways:

  * •

For classical model atmosphere construction, not using an opacity table (that
is, setting IOPTAB=0), the frequency points are set by the standard procedure,
described in Paper III, § 4.1 and 7.4.2 (for continuum frequencies), and §
11.3 (for line frequencies). The actual frequencies are set by the program;
the corresponding input data are used to determine their number and values.
The number of continuum frequencies is determined based on the input parameter
NFREAD, and their actual setting is given through the keyword parameters
FRCMIN, FRCMAX, CFRMAX, NFTAIL, and DFTAIL (Paper III, § 7.4.2). Basically,
the program sets two frequency points close to the edge of all continua
(bound-free transitions) arising from explicit levels, and places some points
in between the edges and in the region of frequencies higher than the highest-
frequency edge (determined by FRCMAX, NFTAIL, and DFTAIL).

The line frequencies are selected for each line between explicit levels, based
on the specific input in the corresponding atomic data file, as described in
Paper III, § 11.3. Once this is done for all lines, the complete set of
frequency points is checked for points corresponding to different lines that
lie very close and are thus redundant, and such frequency points are
eliminated.

  * •

For a pure opacity-table run (IOPTAB$<0$), there are no explicit transitions,
and the frequency points are set in three possible ways:
(i) to be logarithmically equidistant between FRCMIN and FRCMAX (set if the
input parameter NFREAD is negative); or
(ii) to be exactly equal to the frequencies of the opacity table (set by the
keyword parameter IFRSET=0, independently of the value of NFREAD).
(iii) to be logarithmically spaced between the minimum and maximum frequencies
in the table, set by a value of IFRSET$>0$777Obviously, IFRSET should not be
smaller than the number of frequencies of the opacity table by large factors,
in order not to deteriorate the accuracy of resulting models considerably..

However, this procedure is not satisfactory for the hybrid opacity-table mode
(IOPTAB=1). We have therefore upgraded the procedure in the following way:

  * •

In this case (IOPTAB=1), one sets the input parameter NFREAD to a positive
value, so that the frequency points corresponding to both bound-free and
bound-bound transitions between explicit levels are determined as in the
classical mode (IOPTAB=0).

  * •

One then adds the frequency points corresponding to the opacity table (and
therefore describing the opacity of all non-explicit species); these
frequencies are either exactly the same as in the opacity table, or are
logarithmically equidistant between two limiting frequencies. In the latter
case, the number of table-related frequencies and their range are optional,
and are given by the input parameters.

Here is a detailed description of the corresponding parameters:

NFREAD

(4th line of the standard input) – has essentially the same meaning as before,
namely:
$\bullet>0$ – the standard mode, i.e., the program sets the frequency points
for all explicit continua and lines;
$\bullet\leq 0$ – the subset of frequency points that are set independently of
the actual explicit continua and lines is constructed. This subset represents
either a full set of frequencies – in the case of the full opacity-table mode
(IOPTAB $<0$); or it represents a subset of frequencies (typically, the set of
table-related frequencies) to be merged with the previously determined
frequency points.The actual number of such frequency points depends on the
keyword parameter IFRSET.

IFRSET

(one of the keyword parameters) – a table-related frequency setting parameter:
$\bullet=0$ – the table-related frequency points subset is exactly the same as
the frequency points used in the table;
$\bullet>0$ – the table-related frequency set is composed of IFRSET points,
logarithmically equidistant between the minimum and maximum frequency of the
opacity table, and the values of opacities in the new set of frequencies are
interpolated from the original table. In other words, this option allows one
to reduce the number of frequencies treated through the opacity table, without
severely limiting the accuracy of the overall procedure. This option is
similar, but more flexible, than a reduction of the number of frequencies of
the table described earlier in § 2.2.2, because the present option enables
consider an arbitrary number of frequency points for the table-related subset,
while the former option only allows to consider one half, one third, etc.
frequencies from the original table. On the other hand, the latter option may
be advantageous from the point of view of reducing the total memory
requirement of the code, given that one can use a very large opacity table and
immediately remove a portion of frequencies before their corresponding
opacities are stored in memory.
$\bullet<0$ – This option is not very useful, but opens the possibility of
testing the sensitivity of a resulting model atmosphere to the frequency range
considered in the calculations. It is invoked only if NFREAD $<0$, and sets
the table-related frequency points to contain abs(NFREAD) points,
logarithmically equidistant between FRCMIN and FRCMAX. It is similar to the
option IFRSET$>0$, only that the minimum and maximum frequencies do not have
to be equal to the minimum and maximum frequency of the opacity table.
DEFAULT: IFRSET=0

This might seem complicated, but in fact it is not. To help the user, we give
some simple examples below. In all examples, let us assume that the opacity
table contains 100,000 frequencies.

(1) To compute a model atmosphere in a full opacity-table mode, and using the
exactly same frequencies as in the opacity table. In this case, one needs to
set:
IOPTAB $=-1$ in the keyword (non-standard) parameters file;
NFREAD $=0$ (or, in fact, any number NFREAD $<0$) – in the standard input;
IFRSET$=0$ – in the keyword parameter file (but it need not be specified
because 0 is the default value).

(2) Analogous as before, but using only 30,000 frequencies from the table,
logarithmically spaced between the minimum and maximum tabular values:
IOPTAB $=-1$
NFREAD $=0$
IFRSET $=30000$

(3) To compute a model atmosphere in the hybrid opacity-table mode, including
all frequencies from the table, in addition to those associated to explicit
transitions, one sets:
IOPTAB $=1$
NFREAD $=50$ (for instance, or some other non-zero value – see below)
IFRSET $=0$ (again, this does not have been specified because it is a
default).

(4) Analogous as above, but again using only a subset of 30,000 frequencies
from the opacity table:
IOPTAB $=1$
NFREAD $=50$
IFRSET $=30000$

These examples represent some of the most useful cases of actual model
construction. A note about the value of NFREAD in the case of hybrid opacity-
table mode (IOPTAB=1): For classical models, this value is usually chosen to
be larger because one needs to reproduce the continuum opacity accurately.
However, when using an opacity table, one has plenty of frequency points
included anyway, so a moderate value of NFREAD is sufficient, mostly to set up
frequencies in the tail regions beyond the span of the opacity table.

###  2.5 Line broadening

####  2.5.1 Hydrogen lines

The adopted treatment of hydrogen lines is described in detail in Paper III, §
12.1. Here we only provide a short description, highlighting differences from
the previous versions.

In both tlusty and synspec, there are several options available. A control of
these options is somewhat different in tlusty and synspec.888The reasons for
that are both historical and functional. With tlusty, one can compute very
simple models, with only a relatively small number of frequency points
considered. Moreover, for simple model atmospheres one may even use a Doppler
profile for hydrogen lines. However, synspec is used primarily for computing
detailed spectra, where many frequency points are set up, and where
considering for instance a Doppler profile for hydrogen lines does not make
much sense. In synspec, the choice is fully controlled by the keyword (non-
standard) parameter IHYDPR – see below. In tlusty, this choice is controlled,
as for any other element, by the corresponding atomic data file, namely by the
parameter IFANCY, that sets the mode of treating the line profile (see Paper
III, § 11.3.), namely:

IFANCY

– a mode of treatment of the absorption profile :
$\bullet\,=0$ – Doppler profile;
$\bullet\,=1$ or $-1$ – Voigt profile;
$\bullet\,=2$ or $-2$ – approximate Stark (+ Doppler) profile for hydrogenic
lines after Hubeny et al. (1994);
$\bullet\,=3$ or $-3$ – hydrogen line profiles given by Lemke’s (1997) tables.
$\bullet\,=4$ or $-4$ – hydrogen line profiles given by Tremblay & Bergeron
(2009) tables.

However, if the keyword parameter IHYDPR is set to 1 or 2, it overwrites the
setup specified by IFANCY. This is useful to avoid the need of modifying
manually the atomic data file for hydrogen.

IHYDPR

– a mode of treatment of the hydrogen line broadening.
$\bullet\,=0$ – in tlusty, the specification of the mode of line profile
evaluation remains to be defined by the values of IFANCY from the hydrogen
atomic data file.
– In synspec, hydrogen line broadening is treated using an approximate Stark
broadening after Hubeny et al. (1994).
$\bullet\,=1$ – hydrogen line broadening is computed using the Stark
broadening tables constructed by Lemke (1997), for all Lyman and Balmer lines
with a principal quantum number of the upper level $\leq 10$. The table
contains data for higher lines as well, but these are not recommended because
they turned out to be incorrect.
$\bullet\,=2$ – hydrogen line broadening is described using the Tremblay &
Bergeron (2009) data. The available data file, kindly suppled by
P.-E.Tremblay, is extended with respect to the one used in the previous
versions. The current file contains data for the first 20 lines of the Lyman
and Balmer series, the first 19 lines of the Paschen series, and the first 10
lines of the Bracket series. (The previous file contains data for only the
first 10 members of the Lyman and the Balmer series.)
DEFAULT: IHYDPR=0

####  2.5.2 Metal lines

Traditionally, line broadening was treated differently in tlusty and synspec.
The rationale for this approach is that in tlusty one generally does not
require detailed knowledge of the spectrum, and line broadening only enters
through its influence on evaluating radiative rates in the bound-bound
transitions, and for evaluating integrals needed to handle the radiative
equilibrium. In both cases, the quantity of interest is an integral of the
absorption profile coefficient times the mean intensity of radiation, and such
integrals are dominated by the line core, where the profile coefficient is
well approximated by the Doppler profile. This is the reason why the large
majority of explicit lines are treated with a Doppler profile.

tlusty has an option to treat selected lines with a Voigt profile, but in this
case the corresponding damping parameters should be included in the atomic
data file - see Paper I, § 11.3. The atomic data files provided in the
standard tlusty distribution contain explicit damping parameters for only very
few lines, typically strong resonance lines of the most important atoms/ions.

Such a treatment, however, leads to a somewhat peculiar situation in the case
of hybrid Opacity Table/NLTE (OT/NLTE) approach. The background opacity,
treated through the opacity table generated by synspec, is computed assuming
Voigt profiles for all lines, while the lines treated explicitly by tlusty,
which are supposedly more important, would be treated with a Doppler profile.
Moreover, the Doppler profiles are usually assumed to be depth-independent –
see Paper III, § 11.3.

It would be impractical to modify all the available atomic data files.
Instead, tlusty reads the original atomic data files, and in the OT/NLTE
approach (IOPTAB=1), modifies the setup of the basic broadening parameters for
explicit lines in the following way:

  * •

If a given line was already set for a treatment with a Voigt profile, that is
the input parameter IFANCY (later renamed to IPROF) is set to 1, then

    * –

the parameter LCOMP is set to .TRUE. \- i.e., assuming a depth-dependent
profile (which was very likely already done the original data), and the
frequency-point setting parameter INTMOD is set to 3, which invokes a
frequency grid appropriate for the modified Simpson integration rule999Defined
here as a set of three-point integrations, starting at the line center; each
consecutive one having a doubled interval between points, that is with
frequency displacements from the line center expressed in units of Doppler
width as
$0,x_{0},2x_{0},4x_{0},6x_{0},10x_{0},14x_{0},22x_{0},30x_{0},\ldots$.
Generally, for $2n+1$ points, the total extent of integration is
$(2^{n+1}-2)x_{0}$. Conversely, $x_{0}={\rm X\\!M\\!A\\!X}/[2^{({\rm
N\\!F-1})/2+1}-2]$..

    * –

the number of frequency points per half-profile (NF) is set to:
– $\max({\rm NF},25)$

    * –

the extent of coverage of frequency points from the line center expressed in
units of fiducial Doppler width, XMAX, is set to
– $\max({\rm X\\!M\\!A\\!X},4000.)$ for resonance lines

The damping parameters are left at the values supplied in the original atomic
data file.

  * •

If a given line was set to a Doppler (Gaussian) profile in the original file
(IFANCY=0), then it is automatically reset to IFANCY=1, i.e. for a treatment
with a Voigt profile, again LCOMP is set to .TRUE. and INTMOD=3, and with NF
and XMAX reset as follows: If $i$ is the index of the lower level of the
transition in the numbering of explicit levels of the given ion, that is,
$i=1$ is a ground state, then

    * –

for $i=1$, ${\rm NF}\rightarrow\max({\rm NF},23)$; and ${\rm
X\\!M\\!A\\!X}\rightarrow\max({\rm X\\!M\\!A\\!X},4000.)$,

    * –

for $i\leq 3$, ${\rm NF}\rightarrow\max({\rm NF},21)$; and ${\rm
X\\!M\\!A\\!X}\rightarrow\max({\rm X\\!M\\!A\\!X},2000.)$,

    * –

for $i\leq 5$, ${\rm NF}\rightarrow\max({\rm NF},17)$; and ${\rm
X\\!M\\!A\\!X}\rightarrow\max({\rm X\\!M\\!A\\!X},600.)$,

    * –

for $i\leq 10$, ${\rm NF}\rightarrow\max({\rm NF},13)$; and ${\rm
X\\!M\\!A\\!X}\rightarrow\max({\rm X\\!M\\!A\\!X},200.)$,

    * –

otherwise, ${\rm NF}\rightarrow\max({\rm NF},9)$; and ${\rm
X\\!M\\!A\\!X}\rightarrow\max({\rm X\\!M\\!A\\!X},60.)$.101010Such a choice of
parameters means that in all cases the setting of frequency points is very
similar. Resonance lines are allowed to contribute to a larger frequency range
than other lines, but in all cases the closest points to the line center are
separated by about $1/2$ of the fiducial Doppler width, so that the radiative
rates are evaluated with sufficient accuracy.

  * •

Damping parameters are set to (for an explanation, refer to Paper III, §
11.3): GAMAR=0, STARK1=0, VDWH=1, i.e. the damping parameters are assumed to
be given by classical expressions – see Paper I, Appendix A. For convenience,
we give the classical expressions for the natural. Stark, and Van der Waals
broadening here as well.

    * –

Natural broadening:

$\Gamma_{\rm rad}=2.67\times 10^{-22}\nu_{0}^{2},$ (1)

    * –

Stark broadening

$\Gamma_{\rm Stark}=10^{-8}n_{\rm eff}^{5/2}\,n_{\rm e},$ (2)

where $\nu_{0}$ is the frequency of the line center, $n_{\rm eff}\equiv
Z_{I}^{2}[E_{H}/(E_{I}-E_{j})]^{1/2}$ is the effective quantum number of the
upper level, $j$, of the transition in ion $I$, with the excitation energy
$E_{j}$ and the ionization energy $E_{I}$, $Z_{I}$ is the effective charge
($Z_{I}=1$ for neutrals) of the ion $I$, $E_{H}$ is the ionization energy of
hydrogen, and $n_{\rm e}$ is the electron density.

    * –

Van der Waals broadening

$\Gamma_{\rm vdW}=4.5\times 10^{-9}(2.5\,n_{\rm eff}^{4}/Z^{2})^{0.4}\,c_{w},$
(3)

with

$c_{w}=(N_{\rm H}+0.42N_{\rm He}+0.85N_{{\rm H}_{2}})(T/10^{4})^{0.45},$ (4)

where $N_{\rm H}$, $N_{\rm He}$, and $N_{{\rm H}_{2}}$ are the total number
densities of neutral hydrogen, neutral helium, and molecular hydrogen,
respectively. The contribution of the hydrogen molecule may be quite important
for cool stars.

We recall that we use here a classical formalism, where for each broadening
mechanism $\Gamma$ represents a half-intensity width corresponding to the
Lorentz profile, expressed in circular frequencies, in the atom’s rest frame
[see, e.g. Hubeny & Mihalas (2014, § 8.1)]. In the ordinary frequency units,
this width becomes $\Gamma/2\pi$. Specifically, the normalized profiles are
given by

$\phi(\omega)=\frac{\Gamma/2\pi}{[(\omega-\omega_{0})^{2}+(\Gamma/2)^{2}]},\quad{\rm
and}\quad\phi(\nu)=\frac{\Gamma/4\pi^{2}}{[(\nu-\nu_{0})^{2}+(\Gamma/4\pi)^{2}]},$
(5)

where $\omega$ is the circular frequency, $\omega_{0}$ the circular frequency
of the line center, and analogously $\nu$ and $\nu_{0}$ are ordinary frequency
and the line-center frequency; their relation being $\nu=\omega/2\pi$. We also
note that the usual Voigt parameter, that represents the damping parameter
expressed in the units of Doppler widths, is given by
$a=\Gamma/(4\pi\Delta\nu_{D})$, where $\Delta\nu_{D}$ is the Doppler width.

###  2.6 Equation of state

The treatment of the equation of state was described in Paper II, § 2.7.
However, the present versions, tlusty208, as well as synspec54, contain a
number of improvements in the treatment of molecules, both regarding the
solution of the chemical equilibrium, and the evaluation of the partition
functions.

Chemical equilibrium is treated by the Tsuji (1973) formalism. Consider, for
instance, a molecule $X_{l}Y_{m}Z_{n}$, formed by $l$ nuclei of atom $X$, $m$
nuclei of atom $Y$, and $n$ nuclei of atom $Z$. The partial pressures of the
individual constituents satisfy

$K_{P}=\frac{P_{X}^{l}P_{Y}^{m}P_{Z}^{n}}{P_{X_{l}Y_{m}Z_{n}}},$ (6)

where $P_{X}$, $P_{Y}$, $P_{Z}$, and $P_{X_{l}Y_{m}Z_{n}}$ are the partial
pressures of $X$, $Y$, $Z$, and $X_{l}Y_{m}Z_{n}$, respectively, and $K_{P}$
is the equilibrium constant. A generalization for a different number of
constituents is straightforward. For negative ions, electrons are considered
as an additional component. For positive ions, electrons are again taken as
additional component, but Eq. (6) is modified to read

$K_{P}=\frac{P_{X}^{l}P_{Y}^{m}P_{Z}^{n}}{P_{(X_{l}Y_{m}Z_{n}\\!)^{\\!+}}P_{e}},$
(7)

where $P_{e}$ is the partial pressure of electrons.

For the purpose of numerical evaluation, the equiiibrium constant is
approximated by the expression

$\log K_{P}=\sum_{i=0}^{4}a_{i}\theta^{i},$ (8)

where $\theta=5040/T$, and the fitting coefficients $a_{i}$ are given by Tsuji
(1973). This approach was used in the previous version of tlusty. For the
present version one of us (YO) used calculations from Barklem & Collet (2016)
to compute new fitting coefficients for most of molecular species considered
by Tsuji, and added new molecular species, including a number of positive and
negative molecular ions. In total, the chemical equilibrium is solved for 503
molecular species, in contrast to the original 356 species considered by
Tsuji, and 38 atomic species.

Similarly, the partition functions for the molecular species were evaluated
using the Irwin (1981) tables, wherever available, or using the Tsuji (1973)
data for molecules not considered in the Irwin tables. Again, a new Irwin-like
table was created (by YO) by fitting Barklem & Collet’s results. Many more
molecular species are considered in the new table (324) in contrast to the
original Irwin’s table (66). However, we stress that tlusty uses them only for
an evaluation of the entropy (see § 2.7 and Appendix C), but not for a
computation of the molecular line opacity which is is taken from a pre-
calculated opacity table, generated by synspec (or by another code with a
table-generating capability). The molecular partition functions are therefore
much more important in synspec than in tlusty. In synspec, we have also
implemented partition functions from the EXOMOL
project111111http://www.exomol.com/data – see § 3.3., but not in tlusty
because the available Barklem & Collet data are sufficiently accurate for this
purpose [as is seen from Eq. (36), partition functions enter only as
logarithms].

To make a model atmosphere calculation more accurate and flexible, as well as
to allow for a downward compatibility, we have introduced several new keyword
parameters:

IFMOL

– its meaning is exactly the same as in version 205, i.e., it represents a
switch for including molecular formation in the equation of state. If it is
set to a non-zero value, then the molecules are taken into account in the
equation of state, but not necessarily at all depths. This is determined by
the following, newly introduced keyword parameter TMOLIM.
DEFAULT: IFMOL=0 (but in synspec IFMOL=1 is a default!)

TMOLIM

– sets the division temperature (K) for including molecules in the equation of
state. For $T<$ TMOLIM, the equation of state includes molecules, and 38
neutral and singly ionized atomic species, as described in Paper II, § 2.7.3
and above; while for $T\geq$ TMOLIM molecular formation is not considered, and
the equation of state is treated in the standard way, as described in Paper
II, § 2.7.1.
DEFAULT: TMOLIM = 9000.

MOLTAB

– a switch to select the molecular equilibrium (Tsuji-like) table of fitting
coefficients to evaluate the equilibrium constants
$\bullet\,=0$ – the original Tsuji (1973) is used (kept for comparison
purposes);
$\bullet\,\not=0$ – an improved table based on Barklem & Collet results is
used
DEFAULT: MOLTAB=1 (i.e., the new, improved table)

IIRWIN

– a switch to consider Irwin (1981) partition function for atoms and ions as a
default for $T<16,000$ K.
$\bullet\,=0$ – the default for atoms/ions are the Hamburg (Traving et al.
1966) partition functions.
$\bullet\,=1$ – the partition functions for atoms/ions for $T<16,000$ K for
neutrals, once, and twice ionized atoms are evaluated using Irwin (1981) data.
DEFAULT: IIRWIN=1

IRWTAB

\- analogous to MOLTAB; a switch to select the table of fitting coefficients
to evaluate the molecular partition functions.
$\bullet\,=0$ – the original Irwin’s table is used. Since the improved table
is obviously preferable, we keep this option for possible comparisons with old
approaches.
$\bullet\,\not=0$ – an improved table, based on Barklem & Collet (2016) data
for diatomic molecules, is used.
DEFAULT: IRWTAB=1

IPFEXO

\- a switch to set up an evaluation of the molecular partition functions using
the EXOMOL data. It is relevant only to synspec – see § 3.3.

We stress that in version 205 when setting IFMOL=1 the molecular equation of
state was considered for all depths in the atmosphere, including the deepest,
and usually hottest, layers, where some higher ions may be formed, which were
however disregarded in the molecular state equation routine. This problem is
now relieved by introducing the keyword parameter TMOLIM. We have still
upgraded our equation of state solver, which now takes into account molecules,
neutral, once and twice ionized atoms, negative ions of some atoms and
molecules, and electrons. In the previous versions, twice ionized atoms were
not considered, which could cause some inaccuracies at temperatures close to
that specified by TMOLIM.

When NLTE is considered in a model where molecular formation is allowed for,
the total number of atoms of a given element treated in NLTE (a quantity that
enters the right-hand-side of the closing relation of the set of kinetic
equilibrium equations) is updated consistently, i.e. the nuclei that have
already been sequestered in molecules have been subtracted. This was not true
in version 205.

Another consequence of introducing an option to make the equation of state
with molecules more consistent is that the previous basic output parameters of
the model, namely the temperature, $T$, the electron density, $n_{\rm e}$, and
the mass density, $\rho$, are no longer sufficient to describe the atmospheric
state. One needs to add an additional state parameter, the total particle
number density, $N$. This is explained in more detail in Appendix B.

In order to make this change painless for users, and particularly for various
graphical tools that use the output form previous versions of tlusty, the
output on unit 8 (described in Paper III, Chap. 8) has now the following
structure:

$\bullet$ 1st line: two numbers, ND and NUMPAR, where

ND

– number of depth points

NUMPAR

\- number of input model parameters for each depth
$\bullet$ $>0$ – has the same meaning as in the previous versions, that is the
total particle number density $N$ is not part of the set of output parameters.
The first three parameters are $T$, $n_{\rm e}$, $\rho$, and there may be
additional parameters for disks (vertical coordinate $z$), and/or for NLTE
models (the set of populations of all explicit levels).
$\bullet$ $<0$ – $N$, the total particle number density, is part of the set of
output parameters, with their total number being $-$NUMPAR. The first four
parameters are now $T$, $n_{\rm e}$, $\rho$, and $N$; the rest are the same as
above.
This option is switched on automatically when the treatment of molecules is
set, i.e., when IFMOL=1.

$\bullet$ subsequent lines: analogous to the previous versions of tlusty, that
is with a line or more for each atmospheric layer, namely:
– for LTE models – just $T$, $n_{\rm e}$, $\rho$, and for models with
molecules also $N$;
– for NLTE models – these 3 or 4 parameters, and the populations of all
explicit levels (in cm-3).

###  2.7 Convection

The new version offers a somewhat modified treatment of convection, which is
however essentially transparent for the user. The main difference from version
205 is a change in the default evaluation of the adiabatic gradient,
controlled by the keyword parameter IFENTR – see below. The adiabatic gradient
is now evaluated by default using entropy instead of internal energy, which
was the default in the previous versions. We have made sure that both
approaches yield very similar results. The reason of preferring an evaluation
through the entropy is that one needs to compute numerically fewer partial
derivatives, which leads to better accuracy. For more details, and expressions
for the entropy of the material, see Appendix C.

We have also introduced two new keyword parameters for controlling the
evaluation of partial derivatives of the corresponding thermodynamic
parameters. We describe them below, and we repeat for convenience the
description of parameter IFENTR.

IFENTR

– a switch for evaluating the adiabatic gradient:
$\bullet$ $=0$ – the adiabatic gradient $\nabla_{\rm ad}$ is evaluated through
the internal energy – see Paper I, Eqs. (346) - (348).
$\bullet$ $>0$ – the adiabatic gradient $\nabla_{\rm ad}$ is evaluated through
the entropy – see Appendix C.
DEFAULT: IFENTR=1

DIFT

– numerical value of $\Delta T/T$ used for evaluating a partial derivative of
any thermodynamic quantity ($f$) with respect to $T$, namely

$\frac{\partial f}{\partial T}=\frac{f[T(1+\Delta T/T)]-f[T(1-\Delta
T/T)]}{2T\,\Delta T}$

DEFAULT: DIFT=0.01

DIFP

– analogous quantity for the derivatives with respect to the gas pressure
DEFAULT: DIFP=0.01

###  2.8 Collisional rates

This modification offers two significant improvements in the treatment of the
collisional rates of specific ions and transitions; details will be published
elsewhere (Osorio et al., in prep.). These are:

(i) Including collisions not only with electrons, but also with protons and
neutral hydrogen atoms.

(ii) Collisional rates may be supplied as tabular values instead of analytic
expressions. In all cases, the rates tabulated are upward rates, and downward
rates are evaluated internally via the detailed balance relations.

There is a subtle point here. In order to properly formulate the condition of
detailed balance, and thus to evaluate the downward rates, we need first to
consider several issues.

Atomic physics supplies a cross-section, $\sigma_{ij}(v)$, for producing the
transition $i\\!\rightarrow\\!j$ by collisions with a certain colliding
particle, or projectile (electron, proton, or neutral hydrogen), moving with
velocity $v$ relative to the atom. Then the total number of transitions is
equal to the level population times the collisional rate $C_{ij}$,

$n_{i}C_{ij}=n_{i}n_{P}\int_{v_{0}}^{\infty}\\!\sigma_{ij}(v)f(v)dv\equiv
n_{i}n_{P}\,q_{ij},$ (9)

where $v_{0}$ is the velocity corresponding to the energy of the transition,
$E_{ij}$, $n_{P}$ is the number density of projectiles, and $f(v)$ their
velocity distribution, assumed here to be Maxwellian. The quantity $q_{ij}$,
the velocity-averaged cross-sections, is also called the collisional rate
coefficient. This quantity is the actual input to tlusty. In order to evaluate
the rate of the inverse transition $j\rightarrow i$, one invokes the detailed
balance relation, which stipulates that in LTE the number of transitions
$i\rightarrow j$ is equal to number of transitions $j\rightarrow i$.

#### Collisional excitation

Regardless of the nature of the projectile, the detailed balance relation
reads:

$n_{i}^{\ast}n_{P}^{\ast}C_{ij}=n_{j}^{\ast}n_{P}^{\ast}C_{ji},\quad\Rightarrow\quad
C_{ji}/C_{ij}=(n_{i}/n_{j})^{\ast}=(g_{i}/g_{j})\exp(E_{ij}/kT),$ (10)

where the asterisk indicates an LTE population, $n_{i}$ and $n_{j}$ are the
populations of the states $i$ and $j$ of the atom under study, the $g$’s are
the statistical weights, and $E_{ij}$ is the energy of the transition
$i\rightarrow j$. Since the nature of projectiles remains unchanged during the
collision, their LTE number density cancels out. This is the standard
expression widely used in astrophysical work..

#### Collisional ionization

We consider the three possible projectiles in turn:

$\bullet$ Ionization by electrons
Again, this is a standard case, analogous to Eq. (10),

$C_{ji}/C_{ij}=(n_{i}/n_{j})^{\ast}=n_{\rm
e}\Phi_{0}T^{-3/2}(g_{i}/g_{j})\exp(E_{ij}/kT),$ (11)

where the level $j$ now corresponds to the next higher ion, usually, but not
necessarily, its ground state. The ratio of LTE populations is determined
through the Saha-Boltzmann relation, where $\Phi_{0}=(1/2)(2\pi m_{\rm
e}k/h^{2})^{-3/2}=2.0706\times 10^{-16}$, and $T$ is the temperature. For
details, refer e.g. to Hubeny & Mihalas (2014, § 4.2 and 9.3).

$\bullet$ Ionization by neutral hydrogen
From the physical point of view, the ionization processes can be either a true
collisional ionization, that is a reaction

${\rm A}_{i}+{\rm H}(1s)\leftrightarrow{\rm A}^{\\!+}_{j}+{\rm H}(1s)+e,$ (12)

with an inverse reaction being a three-body recombination; or a charge
exchange reaction

${\rm A}_{i}+{\rm H}(1s)\leftrightarrow{\rm A}^{\\!+}_{j}+{\rm H}^{-},$ (13)

where the inverse process is called mutual neutralization. Here, Ai denotes an
atom A at state $i$, and A${}^{\\!+}_{j}$ the corresponding ion at state $j$
(again, usually, but not necessarily, its ground state). One can in principle
consider collisions with hydrogen in excited states, but the cross-sections
for such collisions are poorly known, and in any case such a process is much
less important, so it is neglected here.

As shown e.g. by Barklem (2016), the later reaction, namely a charge exchange
with neutral hydrogen, is much more probable than the true ionization process,
and therefore in tlusty we only consider this process. The corresponding
detailed balance relation is therefore

$n_{i}^{\ast}\,n_{{\rm
H}}^{\ast}\,q_{ij}^{H}=n_{j}^{\ast}\,n_{\rm{H}^{\\!-}}^{\ast}\,q_{ji}^{H},$
(14)

where $n_{\rm H}$ and $n_{\rm{H}^{\\!-}}$ are the populations of the ground
state of neutral hydrogen and of H-, respectively. The superscript $H$
indicates that the rate corresponds to interaction with neutral hydrogen.
Consequently

$\displaystyle q_{ji}^{H}/q_{ij}^{H}$ $\displaystyle=$
$\displaystyle(n_{i}/n_{j})^{\ast}(n_{\rm H}/n_{\rm{H}^{\\!-}})^{\ast}$ (15)
$\displaystyle=$ $\displaystyle[n_{\rm
e}\Phi_{0}T^{-3/2}(g_{i}/g_{j})\exp(E_{ij}/kT)]\,[n_{\rm
e}\Phi_{0}g_{1\\!s}^{-1}T^{-3/2}\exp(E_{{\rm H}^{\\!-}}/kT)]^{-1}$
$\displaystyle=$ $\displaystyle(g_{i}/g_{j})g_{1\\!s}\exp[(E_{ij}-E_{{\rm
H}^{\\!-}})/kT].$

where $E_{{\rm H}^{-}}$ is the dissociation energy of H-, $E_{{\rm
H}^{-}}=0.7552$ eV.

The collisional rates for ionization by neutral hydrogen are therefore

$C_{ij}^{H}=n_{\rm H}\,q_{ij}^{H},\quad\quad
C_{ji}^{H}=n_{\rm{H}^{\\!-}}\,q_{ji}^{H},$ (16)

where $n_{\rm H}$ and $n_{\rm{H}^{\\!-}}$ are now actual, i.e. generally NLTE,
populations of the ground state of hydrogen and of H-.

Equation (16) is used in tlusty for the numerical evaluation of the rates of
collisional ionization by neutral hydrogen. We also mention that a relation
for the downward collisional rate, analogous to a relation for the downward
velocity-averaged cross-section,, Eq. (15), can be written as

$\displaystyle C_{ji}^{H}$ $\displaystyle=$ $\displaystyle n_{{\rm
H}^{\\!-}}\,q_{ji}^{H}=n_{\rm H}\,(n_{{\rm H}^{\\!-}}\\!/n_{\rm
H})\,q_{ji}^{H}$ (17) $\displaystyle=$ $\displaystyle n_{\rm H}(n_{{\rm
H}^{\\!-}}\\!/n_{\rm H})(n_{\rm H}/n_{{\rm
H}^{\\!-}})^{\ast}(n_{i}/n_{j})^{\ast}q_{ij}^{H}$ $\displaystyle=$
$\displaystyle C_{ij}^{H}(n_{i}/n_{j})^{\ast}(n_{{\rm H}^{\\!-}}\\!/n_{{\rm
H}^{\\!-}}^{\ast})(n_{\rm H}^{\ast}/n_{\rm H})$ $\displaystyle\equiv$
$\displaystyle C_{ij}^{H}(n_{i}/n_{j})^{\ast}(b_{{\rm H}^{\\!-}}\\!/b_{\rm
H}),$

where $b_{{\rm H}^{-}}$ and $b_{\rm H}$ are the NLTE departure coefficients
($b$-factors) for H-, and the ground state of hydrogen, respectively. In the
case when H- and the ground state of hydrogen are in LTE, $b=1$, the inverse
rate is given by the same simple expression as for collisional ionization with
electrons, Eq. (11), but this no longer applies in a general situation.

$\bullet$ Ionization by protons
Analogously, the collisions with protons are dominated by the charge exchange
reaction

${\rm A}_{i}+{\rm H}^{\\!+}\leftrightarrow{\rm A}^{\\!+}_{j}+{\rm H}(1\\!s).$
(18)

The corresponding detailed balance relation is

$\displaystyle q_{ji}^{p}/q_{ij}^{p}$ $\displaystyle=$
$\displaystyle[n_{i}/n_{j}]^{\ast}[n_{{\rm H}^{\\!+}}/n_{\rm{H}}]^{\ast}$ (19)
$\displaystyle=$ $\displaystyle[n_{\rm
e}\Phi_{0}T^{-3/2}(g_{i}/g_{j})\exp(E_{ij}/kT)]\,[n_{\rm
e}\Phi_{0}g_{1\\!s}T^{-3/2}\exp(E_{\rm H}/kT)]^{-1}$ $\displaystyle=$
$\displaystyle g_{i}/(g_{j}g_{1\\!s})\exp[(E_{ij}-E_{{\rm H}})/kT].$

where $E_{{\rm H}}$ is the ionization energy of hydrogen, $E_{{\rm H}}=13.595$
eV. The superscript $p$ indicates collisional rates with protons.

The collisional rates are given, analogously to Eq. (16), by

$C_{ij}^{p}=n_{{\rm H}^{\\!+}}\,q_{ij}^{p},\quad\quad
C_{ji}^{p}=n_{\rm{H}}\,q_{ji}^{p},$ (20)

Finally, a relation for the inverse collisional rate is, analogously to Eq.
(17),

$C_{ji}^{p}=C_{ij}^{p}(n_{i}/n_{j})^{\ast}(n_{\rm H}/n_{\rm H}^{\ast})(n_{{\rm
H}^{\\!+}}^{\ast}/n_{{\rm H}^{\\!+}}),$ (21)

which again differs from the standard relation, Eq. (11).

From the user’s point of view, this upgrade involves a simple modification of
the input atomic data files, described in detail in Paper III, Chap. 11.

The standard input record for each atomic transition, both bound-free and
bound-bound, may contain, in addition to the previous 9 parameters (II, JJ,
MODE, IFANCY, ICOL, IFRQ0, IFRQ1, OSC0, CPARAM), another one, NCOL. tlusty
detects121212This is detected and properly handled by synspec as well,
although synspec does not need this information, but to keep the input of data
for bound-free transitions consistent and free from confusion. whether this
parameter is present, so the code works smoothly for both the traditional
atomic data files, as well as the upgraded ones. If the parameter NCOL is
present, and has a non-zero value, then a number of subsequent parameters are
read. The meaning of the new parameters is explained below.

NCOL

– number of tables of collisional rates (corresponding to different
collisional processes); NCOL $\leq 3$. Subsequently, there are NCOL blocks of
input values, each containing:

$\bullet$ 1st record:

ITYPE

– the type of the collision process, namely for bound-bound transitions:
$\bullet=1$ – collisions with electrons
$\bullet=2$ – collisions with protons
$\bullet=3$ – collisions with neutral hydrogen atoms
while for bound-free transitions,
$\bullet=1$ – collisions with electrons – electron collisional ionization
$\bullet=2$ – collisions with protons – charge exchange with protons
$\bullet=3$ – collisions with hydrogen atoms – charge exchange with hydrogen

NCTEMP

– number of temperature values for the subsequent table

$\bullet$ following records

CTEMP(I),I=1,NCTEMP

– values of the temperature [K]

CRATE(I),I=1,NCTEMP

– values of the collision rate coefficient [cm3 s-1]

###  2.9 Fixed-density models

In version 205 and all earlier versions of tlusty, the state parameters that
could be held fixed, and the corresponding structural equation not being
solved, were:

$\bullet$ total particle number density ($N$) – not solving the hydrostatic
equilibrium

$\bullet$ temperature ($T$) – not solving the radiative/convective equilibrium

$\bullet$ electron density ($n_{\rm e})$ – not solving the charge
conservation.

Another option was also to keep fixed all three parameters, in which case only
the radiative transfer, together with kinetic equilibrium equation, is solved
(the so-called restricted NLTE problem).

In version 208, we have added another option, namely keeping the mass density,
$\rho$, fixed, while $N$ and $n_{\rm e}$ can change. This option is useful for
producing NLTE level populations and electron density for snapshots of
hydrodynamical simulations that supply temperature and mass density.
Therefore, in this option one usually considers the temperature to be fixed as
well.

The option is controlled by the following two keyword parameters:

IFIXDE

– If set to a non-zero value, the option to keep fixed density is switched on.
In this case, the fixed temperature switch, INRE, is set to zero.
DEFAULT: IFIXDE=0

TFLOOR

– represents the floor value of the input temperature. It is included because
in some hydro simulations, for instance for accretion shocks in young cool
stars, the ambient atmosphere can have very low $T$, say around 3,000 K or so,
while the accretion shock can have temperatures of the order of $10^{5}$ or
$10^{6}$ K. Since one is typically interested in these hot layers, the ambient
temperature can be artificially increased to avoid numerical problems
associated with dealing simultaneously with very high and very low
temperatures.
DEFAULT: TFLOOR = 6000.

###  2.10 Changes of the output

There are three output files that are changed or added with respect to the
previous versions, namely (i) changes of the standard output file; (ii)
changes of the output file fort.7 – a condensed model atmospheres, which were
described in § 2.6, and (iii) an additional output file fort.22 that contains
the absolute $b$-factors – see below in § 2.10.2.

Another change is an introduction of several keyword parameters that control
the amount optional tables produced in the standard output. In the previous
versions, there was a single parameter IPRINT that controlled a number of
output options, so its meaning and an actual use were quite confusing. We list
the new keyword parameters below. In all cases, if a parameter is set to zero,
the corresponding output is not produced. All parameters have a six-character
name, IPxxxx, where xxxx is an abbreviation of what kind of output it
controls. The default values of all these parameters are set to 0 (i.e., not
producing the corresponding output). These pieces of output are interesting or
useful essentially only if there are convergence or some other problems.

The first parameter controls an output of an additional table in the
initialization step:

IPTRAN

\- input parameters for explicit TRANsitions. The table is self-explanatory.
It should be invoked judiciously because it may be large.

The next two parameters control an output of tables at each iteration step of
the global iteration scheme:

IPCONF

\- output of CONvective Flux, for all depths, at each iteration. The
individual columns are: depth index; Rosseland optical depth $\tau_{Ross}$;
temperature $T$ (K); logarithmic temperature gradient $\nabla\equiv
d\ln(T)/d\ln(P)$; adiabatic gradient $\nabla_{\rm ad}$; relative convective
flux, $F_{\rm conv}/F_{\rm tot}$; relative radiative flux, $F_{\rm rad}/F_{\rm
tot}$; and finally $F_{\rm rad}+F_{\rm conv}$, which should become, at later
iterations, to be equal or at least close to 1.

IPRYBH

\- an information about a recalculation of the hydrostatic equilibrium
equation and a determination of an updated gas pressure and density in the
case of Rybicki scheme in the formal solution step, done by subroutine RYBHeq,
described in Appendix E. The individual columns are: depth index, column mass
$m$ (g cm-2, current mass density $\rho$ (g cm-3, radiation pressure $P_{\rm
rad}$ (cgs units), the ration of the current radiation pressure to the
equilibrium one, $P_{\rm rad}/(3\sigma T^{4}/4c)=P_{\rm rad}/(2.5213\times
10^{-15}T^{4})$, where $\sigma$ is the Stefan-Boltzmann constant; gradient of
the radiation pressure $dP_{\rm r}/dm$; gravity acceleration, In the case of
disks, there is an additional column, the vertical distance form the central
plane $z$ (cm).

The last two parameters control a production of two auxiliary tables at the
final iteration:

IPELCH

\- a table of the ELectron density CHeck in the case of opacity tables – see
below

IPELDO

\- a table of the contribution of the individual ELectron DOnors – see below.

####  2.10.1 Standard output

Here we describe the structure of the standard output, highlighting the
changes made in the present version. We label the individual items with “old”
or “new”, referring to whether they existed in the previous versions, or are
newly introduced here. The individual items are listed consecutively; the old
items were described in some detail in Paper III, Chap.9.

  * •

old – the basic header; namely the values of $T_{\rm eff}$ and $\log g$;

  * •

new – a verbatim copy of the content of the keyword parameter file (whose name
being specified as FINSTD in the standard input file);

  * •

old – a list of chemical elements and their adopted abundances;

  * •

old – a table of explicit ions and their basic parameters;

  * •

old – a table of explicit energy levels and their basic parameters;

  * •

new, produced only in the case when an opacity table is used (IOPTAB $\not=0$)
– a block containing a detailed information about the adopted opacity table:

    * –

the name of the file containing the opacity table;

    * –

the numbers of frequencies, temperatures, and densities in the table;

    * –

a list of chemical abundances. There are three columns. The first one lists
the abundances used by tlusty in the current run (the same ones as listed
above). The next two columns list the abundances used by synspec to generate
the opacity table; one for the equation of state, and the other one for
evaluating the opacity, as described in § 2.3. These two abundances differ
only in the case of a partial opacity table where some chemical species are
not included in evaluating the opacity;

    * –

the values of the keyword parameters IFMOL and TMOLIM used by synspec for
constructing the opacity table, and by tlusty in the current run;

    * –

information about the individual additional opacity sources, namely whether
they are included, or possibly rejected because the opacity table already
contains them;

  * •

old – several auxiliary quantities (number of overlapping transitions; an
assessment of the accuracy of the frequency integrations; number of
frequencies; a list of explicit frequencies);

  * •

old – a list of the most important keyword parameters used in the current run;

  * •

new and old – several auxiliary tables of intermediate results, not
particularly important unless something goes wrong. The appearance of these
tables is controlled by the keyword parameters IPCONF and IPRYBH – see above.

  * •

new and old – final section, entitled FINAL RESULTS, which contains three
tables:

    * –

new – if the keyword parameter IPELCH is $\geq 1$, then a table entitled CHECK
OF ELECTRON DENSITY is produced, which contains, for each depth point, the
following quantities:

      * *

depth index,

      * *

temperature,

      * *

actual electron density produced by tlusty in the current run,

      * *

electron density evaluated by the LTE equation of state for the actual values
of temperature and density at the given depth, $n_{\rm e}(T,\rho)$. Its value
may be generally different from the actual $n_{\rm e}$ for NLTE models. In the
case of partial opacity table (IOPTAB$>0$), the difference between these two
values indicates a degree of consistency of the overall calculation because
synspec uses the same procedure for determining electron density for the
tabular values of $T$ and $\rho$. In the case these two values are
significantly different, the opacities following from the opacity table are
not consistent because they were evaluated using an inconsistent electron
density;

      * *

electron density evaluated by an interpolation from the values listed in the
opacity table. A comparison of this value to the previous one indicates an
accuracy of interpolating in the opacity table, because the same interpolation
procedure is being used to evaluate opacity for actual values of $T$ and
$\rho$ from the tabulated opacity at the set of discrete values of $T$ and
$\rho$ used in the table;

    * –

new – if the keyword parameter IPELDO $\geq 1$, then a table entitled
RELATIVE CONTRIBUTIONS OF THE INDIVIDUAL ELECTRON DONORS
is produced, which lists, for each depth, relative contribution to the total
number of electrons from the most important species (H, He, C, N, O, Na, Mg,
Al, Si, S, Ca, and Fe). The first three columns of the table are the depth
index, temperature, and the value of $n({\rm H-})/n_{\rm e}$, representing a
sink of electrons due to H- (a negative contribution);

    * –

old – a table entitled
FINAL MODEL ATMOSPHERE
that presents a list of the basic state parameters of the resulting model. The
table is self-explanatory.

####  2.10.2 New output file fort.22

As explained in more detail is § 3.6.4, in some applications of synspec it is
advantageous to use as input when dealing with NLTE level populations not the
values of the populations themselves, but rather the “absolute” $b$-factors,
defined as a ratio of NLTE and (absolute) LTE populations,

$b_{i}^{\ast}\equiv n_{i}/n_{i}^{\ast\ast},$ (22)

where $n_{i}$ is a NLTE population of level $i$, and $n_{i}^{\ast\ast}$ the
absolute LTE population, obtained by solving the Saha-Boltzmann equations for
all levels of all considered ionization stages of the element.

In contrast, the traditional $b$-factors (stored as before in output file
fort.12) are defined by

$b_{i}\equiv n_{i}/n_{i}^{\ast}=n_{i}/[n_{\rm e}n_{1}^{+}\Phi_{i}(T)],$ (23)

where $n_{\rm e}$ is the electron density, $n_{1}^{+}$ the actual population
of the ground state of the next ion, and $\Phi_{i}(T)$ is the Saha-Boltzmann
factor – see Paper II, § 2.3, or Hubeny & Mihalas (2014, Chap. 9). This
definition is useful from the point of view of understanding NLTE line
formation from the theoretical point of view, while the definition we adopt
here, Eq. (22) is useful from the point of view of abundance analysis, which
in fact is one of the principal uses of synspec.

The traditional $b$-factors are stored in the output file fort.12, as in the
previous versions, while the absolute $b^{\ast}$-factors are stored in the new
output file fort.22, which has exactly the same format as fort.12, and in fact
also as fort.7 – see § 2.6. To repeat, all the three files have the same
header: the 1st line contains the number of depths and the number of physical
parameters per depth; then several lines that contain the mass-depth
coordinate (the column mass in g cm-2), and then for each depth the
temperature, electron density, mass density, total particle density (for
models with molecules), and then for all explicit levels either the level
populations (fort.7). or traditional $b$-factors (fort.12), or absolute
$b^{\ast}$-factors (fort.22). Finally, we stress that the two latter files,
fort.12 and fort.22, are generated only for NLTE models, because for LTE
models all the $b$-factors are identically equal to unity.

####  2.10.3 Specific intensities

There is a new option in tlusty to generate and store not only the emergent
flux, but also the direction-dependent specific intensities131313Before, the
specific intensities were only produced by synspec, if specifically required..
In the past, the number of frequency points considered by tlusty was often
relatively low, and in any case synspec was usually used as a post-processor
to produce a detailed spectrum. With the introduction of opacity tables in
tlusty, and the corresponding ease to produce models with large number of
frequencies considered, it became useful to obtain angle-dependent specific
intensities at the end of a run of tlusty. This is controlled by a new keyword
parameter:

INTENS

\- a switch for producing emergent specific intensities:
$\bullet\,=0$ – no specific intensities are generated;
$\bullet\,>0$ – specific intensities are generated for INTENS number of
directional cosines $\mu$, linearly spaced between 0.1 and 1.
DEFAULT: INTENS=0

The specific intensities are stored in the output file fort.18, which
contains, for each frequency point considered by tlusty, the wavelength (Å),
and INTENS values of the specific intensity $I_{\nu}(\mu)$ (erg cm-2s-1Hz-1),
starting with $\mu=0.1$.

##  3 SYNSPEC

There are several main changes between versions 51 and 54:

  * •

Changes in the line list.

  * •

An option to construct opacity tables.

  * •

A modified and upgraded treatment of additional opacities.

  * •

A modified algorithm of selecting lines from the line list.

  * •

Changes in the molecular equation of state analogous to those made in tlusty.

  * •

An improved methodology for the association of NLTE energy levels between
tlusty and synspec.

  * •

Subtle changes in the computation of the line opacity.

  * •

Allowance for spectral synthesis for LTE model atmospheres with the additional
input of NLTE level populations (more accurately, $b$-factors) taken from a
different model.

###  3.1 Changes in the line list

We stress that the structure of a line list was described in some detail in
Paper I, § 5.3.3, but only for a list providing data for lines of atoms and
ions, while the molecular line lists, which have a somewhat different format,
were not described. Since several more options were introduced recently, we
will describe the structure of the line lists in detail. We also mention that,
unlike previous versions, there is now the option of having any of the lists
in binary format, which may save a considerable amount of computer time, in
particular for evaluating opacity tables – see below.

####  3.1.1 Atomic line lists

While molecular lines can be input by means of several separate files (for
instance for some individual molecular species), synspec is designed to work
with a single atomic line list. Originally, the atomic line list had to be
transmitted to synspec with a fixed filename, namely fort.19, but now it can
have any name, specified in the input file fort.3 (see below).

Currently, there are two variants of an atomic line list, the standard one,
used in all previous versions of synspec, and an extended one, which is
preferable for NLTE model atmospheres because it provides a better association
between atomic energy levels used in tlusty and synspec – see § 3.6.

(i) The standard line list contains one input record for each spectral line.
In fact, it could contain more records, but this option is cumbersome and is
rarely used. Each input record contains the following variables:

         alam anum gf excl ql excu qu agam gs gw inext

and if the continuation indicator inext is non-zero, then the next record is

         wgr1 wgr2 wgr3 wgr4 ilwn iun iprf

These quantities have the following meaning:
$\bullet$ alam – wavelength [nm]. One follows the convention that for
$\lambda<200$ nm the wavelengths are taken in vacuum, while for $\lambda\geq
200$ nm the wavelengths are for standard air (they can however be modified to
vacuum wavelengths by setting the parameter alast in fort.55—see below—to a
negative value);
$\bullet$ anum – numerical code of the element and ion (using Kurucz’s
convention, namely atomic number, decimal point, and degree of ionization: for
instance 2.00 = He I, 6.03 = C IV, 26.01 = Fe II, etc.)
$\bullet$ gf – $\log g\\!f$ ($g$ is the statistical weight of the lower level,
and $f$ the oscillator strength)
$\bullet$ exclu – excitation potential of the lower level [cm-1]
$\bullet$ ql – the $J$ quantum number of the lower level
$\bullet$ excu, qu – analogous quantities for the upper level
$\bullet$ agam – if non-zero, $\log\Gamma_{\rm rad}$ for radiation damping
[s-1]
$\bullet$ gs – if non-zero, $\log\Gamma_{\rm stark}$ for Stark broadening
[s-1]
$\bullet$ gw – if non-zero, $\log\Gamma_{\rm vdW}$ for Van der Waals
broadening [s-1]
$\bullet$ inext – if =1, next record is needed, where:
$\bullet$ wgr1 wgr2 wgr3 wgr4 – Stark broadening values from Griem (1974)
tables, namely the values for $T=(5,10,20,40)\times 10^{3}$ K, respectively
$\bullet$ ilwn – manual setting of the index of the lower level as used by
tlusty
$\bullet$ iun – manual setting of the index of the upper level as used by
tlusty
$\bullet$ iprf – if non-zero, special procedure for evaluating the Stark
broadening (at present ony for He I – see Paper I, Appendix A).

In any of the parameters agam, gs, gw is set to zero, approximate formulae are
adopted for evaluating the corresponding broadening parameter. For details
refer to Paper I, Appendix A; see also here, § 2.5.

(ii) The extended line list. The entry for each spectral line contains, in
addition to the standard 11 parameters (ALAM, ANUM, GF, EXCL, QL, EXCU, QU,
AGAM, GS, GW, INEXT – see above, the set of 6 following new parameters:

         isql ilql ipql isqu ilqu ipqu

where
$\bullet$ isql – the quantum number $2S+1$ of the lower level
$\bullet$ ilql – the quantum number $L$ of the lower level
$\bullet$ ipql – the parity of the lower level
$\bullet$ isqu, ilqu, ipqu – analogous values for the upper level.
These parameters are explained in detail in § 3.6.

We stress that synspec automatically detects whether the line list is in the
standard, old, format or whether it has the extended format, with these 6 new
parameters. If they are present, one can still force synspec to ignore the new
information for all ions by setting the parameter INLIST (in the input file
fort.55) to a value $\geq 10$ – see also § 3.1.3.. One can also force synspec
to ignore this information for a particular ion by setting the parameter
NONSTD in the standard input file to 10, or, more precisely, if a non-zero
value of NONSTD was intended, setting NONSTD to the desired NONSTD+10.

####  3.1.2 Molecular line list

There are three variants of the molecular line list; the standard (classical-
broadening) list; its variant called a minimalistic list; and the EXOMOL-type
list, differing in the extent of the set of the supplied pressure (Van der
Waals) line broadening parameters, and consequently, in the subsequent
treatment of the pressure broadening in synspec. The latter addition was
motivated by recent progress in evaluating the broadening parameters, and an
ever increasing availability of data (e.g., Barton et al., 2017a,b). We
describe them below.

(i) Standard list. Its structure is analogous, albeit simplified, to a
standard atomic list. Each input record contains the following parameters for
each line:

         alam anum gf excl agam gs gw

The meaning of the variables is exactly equivalent to those for the atomic
line list. Since the molecular lines one always assumed to be formed in LTE,
the parameters for the upper level (EXCU and QU) are not needed. The largest
source of uncertainty in treating lines whose data are provided by such a list
is in the evaluation of pressure (Van der Waals) broadening. We stress that
Stark broadening is typically negligible under the conditions where molecules
are formed in significant amounts to make molecular lines contribute
appreciably to the total opacity, and natural broadening is usually well
approximated by the classical formula. The Van der Walls broadening parameter
is taken by a generalization of Eq. (4) in Appendix A of Paper I; see also Eq.
(4) in this paper,

$\Gamma_{\rm vdW}=10^{\tt gw}\,(N_{\rm H}+0.42N_{\rm He}+0.85N_{{\rm
H}_{2}})\,(T/10^{4})^{0.45},$ (24)

where $N_{\rm H}$, $N_{\rm He}$, and $N_{{\rm H}_{2}}$ are the number
densities of neutral hydrogen, helium, and the hydrogen molecule,
respectively, and $T$ is the temperature.

(ii) A minimalstic list. Its introduction was motivated by the observation
that in a vast majority of molecular line lists the parameters gs and gw have
the same value for all lines, and moreover for most molecular species.
Therefore, for such molecules it is a waste of computer memory to store them
in the line list; instead they are assigned by synspec when initializing line
data from the line list. Similarly, the natural broadening in such lists has
its classical value, which again does not have be stored in the line list, but
is computed by synspec instead. Consequently, the minimalistic line lists
contain only four values for each line:

         alam anum gf excl

which have the same meaning as above. The parameter agam is computed as ${\tt
agam}=2.2e13/\lambda^{2}$, with the wavelength $\lambda$ expressed in nm..
Parameters gs and gw can be specified through the keyword (non-standard)
parameters GSSTD and GWSTD, which have default values 3.1e-5 and 1.0e-7,
respectively. Therefore, if they are not specified in the input, the
parameters gs and gw will have the same values as in the usual standard lists.

(iii) The EXOMOL-type list contains for each line the following parameters:

         alam anum gf excl agam  gvdwh2 gexph2 gvdwhe gexphe

where the last 4 parameters describe the Van der Waals broadening following
Barton et al. (2017a,b) as

$\Gamma_{\rm vdW}=10^{-6}ck\,[g_{{\rm H}2}(T/296)^{n_{{\rm H}2}}N_{{\rm
H}_{2}}+g_{\rm He}(T/296)^{n_{\rm He}}N_{\rm He}]$ (25)

where the parameters $g_{{\rm H}2}\equiv{\tt gvdwh2}$, $n_{{\rm H}2}\equiv{\tt
gexph2}$, and analogously for the helium broadening. Broadening by H atoms is
ignored. More details are given in Appendix F.

We stress that synspec recognizes which type of the list is used, so that no
additional input parameter has to be specified.

####  3.1.3 Specification of the format of the line lists

In previous versions of synspec, the file name of the atomic line list was
hardwired as fort.19, and the number and the names of the molecular line lists
was set in the input file fort.55 – see Paper I, § 5.3.1. The relevant lines
in the file fort.55 are the 4th and 7th line. The fourth line contains

           ifreq  inlte  icontl  inlist  ifhe2

Here, only the parameter inlist is relevant for the line lists, in that it
specifies the format of the atomic line list:
$\bullet$ inlist = 0 or 10, all the line lists are in text format, while if
$\bullet$ inlist = 1 or 11, all the line lists are in binary format.

Notice that inlist has a dual meaning, because at the same time
$\bullet$ inlist $=0$ or 1 – the additional quantum numbers (see § 3.1.1), if
present, are being taken into accout;
$\bullet$ inlist $=10$ or $11$, the additional quantum numbers, if present,
are disregarded.

The 7th line, if it exists, contains parameters for the molecular line lists,
namely

           nmlist, (iunitm(i),i=1,nmlist)

where;
$\bullet$ nmlist – the number of molecular line lists,
$\bullet$ iunitm(i),i=1,nmlist – unit numbers of the additional molecular line
lists.
In this case, if the unit number is nn, then the explicit filename is fort.nn.

In the present version of synspec, this approach to specify the filenames of
the molecular line lists, and the mode of the atomic line list, can still be
used. However, the user may choose to use a newly introduced input file fort.3
that, if exists, overwrites the information obtained from fort.55.

The structure of the file fort.3 is as follows. For each line list, including
the atomic one, there is one input record containing 2 (in the case of an
atomic list) or 3 (in the case of a molecular line list) parameters:
AMLIST, IBIN, TMLIM (that is, only AMLIIST and IBIN for atomic lists), where

AMLIST

– a character$\ast$40 variable that specifies the name of the list

IBIN

– an indicator of the type of file:
$=0$ – a text file.
$\not=0$ – a binary file.

TMLIM

– an upper limit of temperature for which is the list being read and used.
This parameter only takes effect when evaluating the opacity tables141414The
reason why this parameter only takes effect for evaluating the opacity tables
is that in this case one reads all lists repeatedly, namely for each
temperature/density pair, while for evaluating a synthetic spectrum for a
given model atmosphere the lists are read only once, but one encounters many
different temperatures at the individual depths. If the atmosphere is
everywhere hot enough that, for instance, TiO is not formed, then one simply
does not take the TiO line list into account at all. \- see next section.

This approach is more flexible than the one using the file fort.55.

Important notes:
(1) We would like to stress that an option of using an input from the file
fort.55 to specify the line lists and their mode is kept for downward
compatibility. It is useful for users who have been using synspec in the past,
and doing spectrum calculations either for atomic lines only, or for a single
molecular line list. In that case, one can use the original file fort.55. In
all other cases, in particular when dealing with multiple line lists, it is
safer and more flexible to set up the specifications for the line lists using
information in fort.3.

(2) When the file fort.3 is used, the parameter INLIST in the file fort.55,
discussed above, does no longer have a function to set the atomic line list in
text or binary format, but it still influences the association of level
indices between tlusty and synspec; that is, if INLIST $\geq$ 10, the
additional quantum numbers, if present in the atomic line list, are
disregarded.

(3) When the file fort.3 is used, its first line has to correspond to the
atomic line list, followed by several lines corresponding to the individual
molecular line lists. Their order is arbitrary, and their number is
independent of the value of the parameters nmlist specified in the file
fort.55.

(4) In the previous versions of synspec, the basic mode IMODE (the first
number in fort.55) had to be set to a value $\geq 10$ in order molecules were
considered at all. This is now obsolete, because of an introduction of a
general keyword parameter IFMOL, transmitted to both tlusty and synspec
through the keyword parameter file. Specifically

IFMOL

– as in tlusty, it switches on the treatment of molecules in the equation of
state, independently of whether molecular lines are being considered in the
opacity or not. As for tlusty, a non-zero value allows for considering
molecules in the equation of state, for local temperatures $T<$ TMOLIM; the
latter being another keyword parameter, with exactly the same meaning as for
tlusty.

At present, synspec typically uses up to three molecular line lists:
(i) for basic diatomic molecules except TiO;
(ii) for TiO;
(iii) for H2O.

At present, the list for basic diatomics includes line data for the following
molecules: H2, CH, NH, OH, NaH, MgH, SiH, CaH, CrH, FeH, C2, CN, CO, MgO, AlO,
SiO, CaO, and VO. The line data are so far taken from the Kurucz
website151515http://kurucz.harvard.edu/molecules, although we are working on
replacing them by the EXOMOL data whenever available.

The number of lines in a full molecular line list may be very large, typically
several tens to hundreds of millions of lines, but in fact most of them do not
contribute appreciably to the total opacity. We have therefore developed a
procedure to eliminate weak lines from the line list from the outset. We
stress that if weak lines are present in the line list, they are eliminated
during the run of synspec based on the input parameter RELOP. However, working
with extended line lists is cumbersome, because (i) they take a lot of memory,
(ii) they are harder to download, and (iii) even reading an extended line
lists by synspec takes a lot of unnecessary computer time. It is therefore
desirable to devise procedures to limit the number of lines in the line list,
while keeping the accuracy of the resulting total opacity reasonably high.
This is described in some detail in Appendix G.

####  3.1.4 Utility programs

The standard distribution also contains two additional, very simple, programs
for handling the line lists, called list2bin and redulis.

$\bullet$ The program list2bin transforms a text line list to binary format.
It is quite trivial, but is included for convenience. Since there are
altogether 5 possible formats of a line list (2 for atomic, and 3 for
molecular lists), the program automatically recognizes the format of the line
list. The standard input is the text line list, and the output binary list is
on file fort.12. It is advantageous to rename this file to have the same core
name, with extension .bin. For instance, assuming that the code was compiled
as

gfortran -o list2bin.exe list2bin.f
then to rewriting the contents of the file gfATO.dat in binary format is
simply accomplished by executing

list2bin.exe <gfATO.dat

mv fort.12 gfATO.bin

$\bullet$ The program redulis performs a reduction of the line list, as
mentioned above, and is described in detail in Appendix G. The code works for
all three types of molecular line lists. The standard input contains either
just one number, $s_{\rm reject}$ (see Appendix G) (which is mandatory), or
three more quantities, $T_{\rm char}$, $\lambda_{\rm char}$, and
$\log(N/U)_{\rm char}$ – see again Appendix G. If these values are not
specified, the default values for these quantities, namely 3000, 10,000, and
10, respectively, are adopted.

The original list should be copied/linked to fort.20, and the output is on
fort.12. The standard output contains information about how many lines of the
individual molecular species are in the original list, and in the reduced one.

###  3.2 Constructing opacity tables

The option for calculating a grid of monochromatic opacities is switched on by
the parameter IMODE (1st number in the first line of the file fort.55 – see
Paper I, § 5.3.1), namely IMODE $=-3$ or $=-4$. Since there are some new
possibilities, we shall summarize the options for the parameter IMODE below:

IMODE

\- the basic mode of operation of synspec:
$=0$ – normal synthetic spectrum;
$=1$ – spectrum with only a few lines (obsolete);
$=2$ – theoretical continuum, plus H (and possibly He II) lines only;
$=-1$ – the so-called “iron-curtain option”. This is now obsolete since the
opacity table option (see below) is much more flexible;
$=-2$ – a continuum-only variant of the iron curtain option; also obsolete;
$=-3$ – a standard opacity-table option (total opacity);
$=-4$ – evaluation of an opacity table with continua + H and He II lines
(analogous to IMODE=2)

When the opacity-table mode is set up (IMODE=$-3$ or $-4$), the input values
of temperature and electron or mass density (thereafter called simply
“density”) are specified in unit 2 (fort.2). The mass density $\rho$ is
typically used as a density parameter. The reason is that, for cool models,
the electron density is less suitable as a density parameter because it is
influenced by details in the ionization equilibrium of various metals (Mg, Si,
Fe, etc.) and, moreover, for very cool models it decreases essentially to
zero.

Numerically, each temperature/density pair is considered as an “atmosphere”
represented by just 1 depth point, ND=1, with corresponding grid values of $T$
and $\rho$. synspec then computes the opacity for all these “atmospheres”,
i.e., $T/\rho$ pairs, consequently and independently of each other.

For each temperature/density pair, synspec thus computes a detailed table of
opacity as a function of wavelength. The internal wavelength grids for the
calculations are set by the program, and are generally different for each
temperature/density pair. In fact, the individual wavelength points are set by
the same procedure as for evaluating a synthetic spectrum – see Paper I. In
order to construct a universal opacity table for many temperature/density
pairs, the program transforms the computed opacities to a universal grid of
wavelength points that is also specified in fort.2. The wavelengths are
logarithmically equidistant between the starting and ending wavelengths
specified as input parameters (and therefore represent a scale with constant
velocity steps). After finishing the calculation of opacities for all
temperature/density pairs, the program outputs the resulting opacity table, in
the format required by tlusty. The name of this file is also specified in
fort.2.

The transformation of opacities from the internal to the universal wavelength
grid can be done in two different ways:

  1. 1.

By a simple interpolation in wavelength. This procedure resembles the idea of
opacity sampling used in constructing metal line-blanketed model atmospheres
(see, e.g., Hubeny & Mihalas 2014; § 18.5). A drawback of this procedure is
that if the number of wavelengths in the universal grid is too low, some cores
of the strong lines, as well as some continuum windows, may be missed.

  2. 2.

By preserving (at least approximately) the integral of the opacity in the
original and universal grids. Specifically, the opacity at the universal grid
point $i$ with wavelength $\lambda_{i}$, denoted by $\kappa^{\rm
tab}(\lambda_{i})$, is given by

$\kappa^{\rm
tab}(\lambda_{i})=\frac{2}{\lambda_{i+1}-\lambda_{i-1}}\int_{(\lambda_{i}+\lambda_{i-1})/2}^{(\lambda_{i}+\lambda_{i+1})/2}\kappa(\lambda)\,d\lambda,$
(26)

where $\kappa(\lambda)$ is the (detailed) opacity computed for the internal
wavelength points. Since the spacing of wavelengths in the internal grid is
roughly equidistant, the above equation is replaced by a mean value of
opacities in the internal wavelength points that lie in the wavelength
interval
$[(\lambda_{i}+\lambda_{i-1})/2,(\lambda_{i}+\lambda_{i+1})/2]$.

There are two basic types of an opacity table, namely:

  * •

a classical table, which also can come in three flavors, depending on how the
values of the density are selected for each temperature – see below. The
opacity is constructed for NTEMP values of temperature, and for each
temperature for a several (generally different) values of density.

  * •

a model-tailored table. Here, one takes a given model atmosphere, with ${\rm
ND}$ depth points, and constructs an opacity table just for ${\rm ND}$
temperature/density pairs, corresponding to exact values of temperature and
density of the input model, $(T_{d},\rho_{d}),\ d=1,\ldots,{\rm ND}$. Besides
the temperature and density, synspec also takes the electron density from the
input model, instead of evaluating it using the LTE equation of state. This
table may be useful for assessing the accuracy of the numerical procedure of
using pre-tabulated opacities, and for increasing the accuracy of computed
models.

To explain a potential utility of a model-tailored table, let us take an
example of computing a model atmosphere for a late-type star with H and Mg
treated explicitly, so that the opacity of these elements is computed on the
fly, and generally in NLTE161616Their opacity may be treated in LTE as well,
but such a treatment would not offer any advantage over using directly a full
opacity table with all species included in it., while the opacity of remaining
elements being taken into account through an opacity table. The table was
constructed for a discrete set of $T,\rho$ pairs. For each pair, the electron
density was determined by an LTE equation of state, $n_{\rm e}=n^{\ast}_{\rm
e}(T,\rho)$ (where the asterisk indicates an LTE value), and the opacity was
evaluated as $\chi(\nu,T,\rho)=\kappa[\nu,T,n^{\ast}_{\rm e}(T,\rho)]/\rho$,
where $\kappa$ is the absorption coefficient per length (in cm-1), while
$\chi$ is the opacity per mass (g cm2); the latter is the quantity stored in
the opacity table. The electron density for the resulting model may differ
from its LTE value because of NLTE effects in Mg (and to a smaller extent in
H). Consequently, the values of the opacity contained in the opacity table may
no longer be sufficiently accurate. To assess the effects of such errors, one
can compute a model-tailored opacity table including H and Mg in NLTE, and
rerun the model atmosphere with this table. If necessary, the overall
procedure may be iterated, but we believe that such effects should be
relatively small, so that just one recalculation of the atmospheric structure
should be entirely adequate.

Another practical use is simply an evaluation of the interpolation errors for
choosing the steps for the opacity table.

####  3.2.1 Input file fort.2

The input file fort.2 has at least 5 lines, but in some cases (see below)
more:

$\bullet$ 1st line: NTEMP, TEMP1, TEMP2:

NTEMP

– number of temperatures, and the mode of setting the opacity table:
$\bullet$ $>0$ – the opacity table is computed in the classical mode; NTEMP
represents the number of temperatures;
$\bullet$ $=0$ – the opacity table is computed in the model-tailored mode. In
this case the temperatures, densities, and electron densities are read from
input model, stored in the file fort.8  as in the standard run of synspec. The
remaining parameters on this line and the two next lines do not take effect.

TEMP1

– lowest temperature [K];

TEMP2

– highest temperature [K].

In the classical mode, the temperatures are set to be logarithmically
equidistant.

$\bullet$ 2nd line: IDENS:

IDENS

– A switch that specifies the treatment of density, and the nature of the
density parameter. It has a dual meaning:
(i) the value of mod(IDENS, 10) determines the nature of the density
parameter:
$\bullet$ mod(IDENS, 10)$=0$ – (i.e. IDENS = 0, 10, 20) – the density
parameter has the meaning of electron density;
$\bullet$ mod(IDENS, 10)$=1$ – (i.e. IDENS = 1, 11, 21) – the density
parameter has the meaning of mass density;
$\bullet$ mod(IDENS, 10)$=2$ – (i.e. IDENS = 2, 12, 22) – the density
parameter has the meaning of gas pressure.

At the same time, the value of IDENS determines
(ii) the actual values of the density parameter to be used in the opacity
table. The next block of input depends on this value:

  * •

IDENS $<10$ – (i.e. IDENS = 0, 1, 2) – the values of the density parameters
are the same for all temperatures. In other words, the pairs of
temperature/density form a regular rectangle in the $\log T/\log\rho$ space
(for simplicity, we take here the mass density for the density parameter, but
the situation is the same for all possible choices of the density parameter).
In this case the block contains just one record:
— 3rd block – one line: NDENS, DENS1, DENS2, where:

NDENS

– number of densities;

DENS1

– lowest density parameter [cgs units];

DENS2

– highest density parameter [cgs units].

  * •

IDENS $\geq 10$ and IDENS $<20$ – similar to above, but the $T/\rho$ pairs
occupy a slanted rectangle in the $\log T/\log\rho$ space. One specifies the
limiting values separately for the lowest and for the highest temperature. For
other temperatures, the limiting vales are obtained by a logarithmic
interpolation. In this case, the next block also contains one record, namely:
—3rd block – one line:
NDENS, DENSL1, DENSL2, DENSU1, DENSU2, where:

NDENS

– number of densities;

DENSL1

– low density parameter limit for the lowest temperature [cgs units];

DENSL2

– high density parameter limit or the lowest temperature [cgs units].

DENSU1

– low density parameter limit for the highest temperature [cgs units];

DENSU2

– high density parameter limit or the highest temperature [cgs units].

  * •

IDENS $\geq 20$ and IDENS $<30$ – the density limits and the number of
densities are generally different for each temperature. In other words, the
$T/\rho$ pairs occupy a completely irregular area in the $\log T/\log\rho$
space. In this case, the next block contains NTEMP input records, each
containing the limiting values as before:
— 3rd block – NTEMP lines: for each $T$, NDENS, DENS1, DENS2, where:

NDENS

– number of densities for the given temperature;

DENS1

– lowest density parameter [cgs units] for this $T$;

DENS2

– highest density parameter [cgs units] for this $T$.

Again, the densities for each temperature are spaced logarithmically
equidistantly.

We stress that for cool stars, it is most advantageous to use mass density
(e.g., IDENS=1, 11, or 21) as the density parameter.

$\bullet$ 4th block – one line: NLAMB, INTTAB, WLAM1, WLAM2:

NLAMB

– number of wavelength points in the final (universal) opacity table;

INTTAB

\- a switch for determining the mode of transformation of opacities from the
internal to the universal wavelength grid points:
$=1$ – the opacities in the universal grid are determined by an interpolation
in wavelength;
$\neq 1$ – the opacities in the universal grid are determined by approximately
preserving the integral over wavelength, i.e., by Eq. (26);

WLAM1

– shortest wavelength in the opacity table [Å];

WLAM2

– longest wavelength in the opacity table [Å].

$\bullet$ 5th block – one line: TABNAM, IBINGR:

TABNAM

– a character variable specifying the name of the file with the opacity table;

IBINGR

– the format of the opacity table:
$\bullet$ $=0$ – the table is stored in text format;
$\bullet$ $>0$ – the table is stored in binary format.

For example, a simple file fort.2 that can be used as a test, is

    -----------------------------------------------------------------
    11  3.0e3    1.5e4
    1
    7  1.0e-12  1.0e-6
    100000  0  900. 110000.
    ’optab11_7f.dat’    0
    -----------------------------------------------------------------

which constructs a table for 11 temperatures between 3000 and 15000 K, 7
values of mass density between $10^{-12}$ and $10^{-6}$g cm-3, for 100,000
wavelengths between 900 and 110,000 Å. Tabular values of the opacity are
determined by an approximate version of Eq. (26). The table is stored in the
file `optab11_7f.dat`, in text format.

If one intends to construct a table (as an example) with the same density
limits as before for the lowest $T$ (3000 K), but with limits $10^{-10}$ and
$10^{-4}$ at the highest $T$, and with interpolated values of the density
limits in between, the first three lines of fort.2 will look like:

    -----------------------------------------------------------------
    11  3.0e3    1.5e4
    11
    7  1.0e-12  1.0e-6  1.e-10   1.e-4
     -----------------------------------------------------------------

For completely irregular area covered by $T/\rho$ pairs, one can code for the
first three blocks, for instance:

    -----------------------------------------------------------------
    11  3.0e3    1.5e4
    21
    7  1.0e-12  1.0e-6
    8  1.0e-12  1.0e-6
    9  3.0e-13  3.0e-7
    9  1.0e-13  1.0e-7
    9  1.0e-13  1.0e-7
    8  3.0e-13  3.0e-7
    8  1.0e-12  1.0e-6
    7  3.0e-13  1.0e-5
    7  1.0e-11  1.0e-4
    6  3.0e-10  1.0e-4
    5  1.0e-10  1.0e-4
    -------------------------------------------------------------------

The next portion of the file fort.2 is optional. If it exists, it contains
information about changing chemical abundances or removing selected species.
It should be stressed that such changes concern only the evaluation of
opacities, and are meant as a tool for producing an opacity table for use in
the hybrid NLTE/opacity-table mode of tlusty, where some species are treated
explicitly, while the opacity of the remaining species enters through a
specific opacity table that does not contain opacity from the species to be
treated explicitly by tlusty.

The equation of state, which essentially specifies a relation between mass
density and electron density and/or the total number density, takes into
account all the species.

We stress that the elemental abundances are specified in the standard input
file for synspec. Recall that if nothing else is specified there, the
abundances are taken as the values from the chosen standard solar abundances
set – see Sec. 2.1.

If this portion of the file fort.2 is missing, all the species considered by
synspec (i.e., all species unless specifically removed by setting MODE=0 in
the standard input file – unit 5; see §2̇.2.1) are taken into account in
constructing the opacity table. If the additional part of fort.2 is present,
it contains a number of input records, one for each species that is going to
be removed (or its abundance changed). Each line contains two numbers:

IATOM RELABN

where

IATOM

– atomic number of the species to be removed/modified;

RELABN

– a factor that multiplies the adopted standard solar abundance to that used
for computing the opacity. The most usual option is RELABN=0. which removes
the element completely.

For example, the file fort.2 may look like

    -----------------------------------------------------------------
    11  3.0e3    1.5e4
    1
    7  1.0e-10  1.0e-6
    100000  0  900. 110000.
    ’optab11_7hmc.dat’    0
      1   0.
     12   0.
     20   0.
    -----------------------------------------------------------------

for constructing an opacity table analogous to that above, but removing
hydrogen, magnesium, and calcium from the resulting opacity table.

####  3.2.2 Input files fort.55 and fort.3

Let us look at an example of evaluating an opacity table, for the internal
wavelength grid that extends from 900 to 110,000 Å, with atomic as well as
molecular lines. The molecular line data data are transmitted through two line
lists, for instance, one list for several important diatomic molecules, and
the other one for TiO. One can use either an original fort.55, or the new one,
together with additional input file fort.3. We shall give both examples below.

(i) Original approach: The file fort.55 may look like this

    -----------------------------------------------------------------
       -3     1     0                        ! imode,idstd,iprint
        0     0     0     0
        0     0     0     0     0
        1     0     0     0     0
        0     0     0
     900. -110000. 50.  2000   1.0e-3   0.15
        2    20    21                        ! nmlist,iunit(ilist)
    -----------------------------------------------------------------

In this example, a negative value of the longest wavelength ($-110,000$)
signifies, as explained in Paper I, that the wavelengths of all lines are
taken in vacuum (unlike the default approach adopted in synspec that considers
air wavelengths for all lines with wavelengths larger than 2000 Å). This
setting is mandatory for constructing the opacity tables.

The 4th number on the 4th line, INLIST=0, signifies that all the line lists
are in ttext format.

Since this approach specifies only unit numbers for the molecular line lists,
one should also make symbolic links. For example, assuming that an atomic line
list gfATO.dat and two molecular line lists gfMOL.dat and gfTiO.dat171717As is
described later on, the file gfTiO.dat is a line list for TiO lines taken from
EXOMOL data, and reduced with the rejection parameter $s_{\rm reject}=-8.5$,
which contains about $8.3\times 10^{6}$ strongest lines out of the original
number of $99.9\times 10^{6}$ lines. More details are given in § 5.2.1. are
either copied or linked to the current directory, then one has to make the
following symbolic links

    ln -s -f gfATO.dat fort.19
    ln -s -f gfMOL.dat fort.20
    ln -s -f gfTiO.dat fort.21

In this case, all line lists are repeatedly read for all temperature/density
pairs, even if, e.g., TiO lines do not contribute at higher temperatures, and
thus the TiO line list is read (and computer time is being spent) in vain.

(ii) New approach. In this case file fort.55 may remain unchanged, or its last
line could be removed. As stressed above, if the file fort.3 is present, it
overwrites the line-lists-related information from fort.55. File fort.3, can
then be, for instance

    -----------------------------------------------------------------
    ’gfATO.dat’     0
    ’gfMOL.dat’     0    8000
    ’gfTiO.dat’     0    5000
    -----------------------------------------------------------------

In this case, the line list for TiO is read only for $T<5000$. Assuming now
that one also has binary versions of the molecular line lists, with a suffix
bin, then one may set up fort.3 as follows:

    -----------------------------------------------------------------
    ’gfATO.dat’     0
    ’gfMOL.bin’     1    8000
    ’gfTiO.bin’     1    5000
    -----------------------------------------------------------------

Another advantage of using input from fort.3 is that, as mentioned above, the
filenames of the molecular line lists are printed at the standard output file,
so one can have an explicit record of what line list was used, in case there
are several variants available.

.

####  3.2.3 Output files

The basic output is the opacity table, in the format described above – see §
2.3. As described there, the file contains, besides the opacity values, a
header that specifies the chemical abundances and the values of the additional
opacity switches used for constructing the table. This information is used by
tlusty to avoid possible inconsistencies in treating opacities.

The name of the file containing the opacity table is specified in the file
fort.2, as explained in § 3.2.1.

Besides the tlusty-compatible opacity table, the program produces two
additional output files that may be used for further inspection of the
results, or in the case there are problems, but these files do not have to be
stored.

fort.26

– continuum opacity, namely a table of continuum wavelengths, versus the
opacity in the continuum, for each temperature/density pair.

fort.27

– all values of opacity on the internal wavelength grid for each
temperature/density pair (i.e., before their interpolation to the universal
wavelength grid). This is a rather big file, and it is produced only if the
parameter IPRIN (specified in fort.55) is larger than 2.

fort.29

– an overview of the number of lines considered (atomic and molecular), and
the computer time for each temperature density pair. The table is more or less
self-explanatory. Briefly, the individual columns are: temperature index,
density index, value of $T$ [K], values of $\rho$ [g cm-3], electron density,
number of atomic lines selected, number of molecular lines from the individual
line lists (all in thousands), and the computing time spent to produce the
opacities on the internal grid [s].

fort.63

– a binary version of the opacity table. This can be used by tlusty if the
input parameter IBINOP is set to a non-zero value – see also § 2.2.2.

####  3.2.4 The utility program ngrid

Since computing the monochromatic opacities on the internal wavelength points
is the most time-consuming part of the construction of an opacity table, it is
sometimes advantageous to generate and store the auxiliary output file
fort.27, and subsequently use it for generating another opacity table with a
different universal wavelength grid and/or a different mode of transformation
of opacities from the internal to the universal grid. This can be done with
the program ngrid.

The program requires the following three input files:

  1. 1.

Standard input file – that contains three input records, the first two being
analogous to the last two records of the input file fort.2 used to generate
the original table, namely

$\bullet$ 1st line: NLAMB, INTTAB, WLAM1, WLAM2,
where

NLAMB

– number of wavelength points in the new opacity table;

INTTAB

\- a switch for determining the mode of transformation of opacities from the
internal to the new universal wavelength grid:
$=1$ – the opacities in the new table are determined by an interpolation in
the wavelengths;
$\neq 1$ – the opacities in the new table are determined by an approximate
version of Eq. (26);

WLAM1

– shortest wavelength in the opacity table [Å];

WLAM2

– longest wavelength in the opacity table [Å].

$\bullet$ 2nd line: TABNAM, IBINGR,
where

TABNAM

– a character variable specifying the name of the file with the new opacity
table;

IBINGR

– the format of the new opacity table:
$\bullet$ $=0$ – the table is stored in text format;
$\bullet$ $>0$ – the table is stored in binary format.

$\bullet$ 3rd line: OPTABLE, IBINOP,
where

OPTABLE

– a character variable specifying the name of the previously generated opacity
table.

IBINOP

– the format of the original (old) opacity table:
$\bullet$ $=0$ – the table is stored in text format;
$\bullet$ $>0$ – the table is stored in binary format.

  2. 2.

File fort.27 – as described above, this is the file that contains the set of
opacities in the internal wavelength points, generated by a previous run of
synspec.

  3. 3.

An old opacity table, whose name is specified by the input parameter OPTABLE.
This file is needed just for providing a header for the new opacity table,
tracking the information from the original calculation.

###  3.3 Equation of state

All the improvements in the equation of state (EOS), and the corresponding
keyword parameters controlling them, used in tlusty208 and described in § 2.6,
apply to synspec54 as well. In fact, the improvements in the equation of state
are more important in synspec than in tlusty, because the most important
effect of molecules, namely the molecular line opacity, is evaluated by
synspec and enters tlusty only as an input, pre-calculated, opacity table.

To evaluate the molecular line opacity, not only the molecular number
densities are needed, but also accurate values of the partition functions. We
will now describe the treatment of molecular partition functions in detail.

In previous versions, the partition functions for the molecular species were
evaluated using the Irwin (1981) tables, whenever available, or using the
Tsuji (1973) data for molecules not included in the Irwin tables. In this
version. a new Irwin-like table was created (by YO) by fitting Barklem &
Collet’s (2016) results. Many more molecular species are considered in the new
table (324) compared to the original Irwin’s table (66).

We have also implemented partition functions provided by the EXOMOL
project.181818http://www.exomol.com/data In the current version, we include
data for the partition functions of 29 diatomic molecules – AlO, C2, CH, CN,
CO, CS, CaH, CaO, CrH, FeH, H2, HCl, HF, MgH, MgO, N2, NH, NO, NS, NaH, OH,
PH, SH, SiH, SiO, SiS, TiH, TiO, and VO; and 3 triatomic molecules – H2O, H2S,
and CO2. We stress that our current molecular line lists do not yet contain
data for some of these molecules, namely CS, HCl, HF, N2, NS, PH, SH, SiS,
H2S, and CO2.

The EXOMOL tables contain values of the partition functions for temperatures
between $T=1$ and $T=T_{\rm max}$ with a step of 1 K. The value of $T_{\rm
max}$ varies for the individual species; ranging from 10,000 (for H2O, H2,,
C2, and FeH) down to 3,000 (for CN, CS, CaH, CrH, MgH, and H2S), and being
typically 9,000, 8,000 or 5,000 for the rest. In order to provide values for
$T>T{\rm max}$, we use the following expression

$U(T)=U^{\rm Irwin}(T)\frac{U^{\rm EXO}(T_{\rm max})}{U^{\rm Irwin}(T_{\rm
max})},\quad{\rm for}\quad T>T{\rm max},$ (27)

which scales the (available) Irwin values at $T>T_{\rm max}$ to match the
EXOMOL data at $T_{\rm max}$.