A Grid of Synthetic Spectra for Subdwarfs: Non-LTE line-blanketed atmosphere models

A new grid of detailed atmosphere model spectra for hot and moderately cool subdwarf stars is presented. High-resolution spectra and synthetic photometry are calculated in the range from 1000{\AA} to 10,000{\AA} using Non-LTE fully line-blanketed atmosphere structures. Our grid covers eight temperatures within 10,000<Teff [K]<65,000, three surface gravities in the range 4.5<logg [cgs]<6.5, two helium abundances matching two extreme helium-rich and helium-poor scenarios, and two limiting metallicity boundaries regarding both solar ([Fe/H] = 0) and Galactic halo ([Fe/H] = -1.5 and [{\alpha}/Fe] = +0.4). Besides its application in the determination of fundamental parameters of subdwarfs in isolation and in binaries, the resulting database is also of interest for population synthesis procedures in a wide variety of stellar systems.


INTRODUCTION
The evolutionary scenario of subdwarfs is still emerging. Hot subdwarf stars evolve from low-to intermediate-mass progenitors and reach far beyond the main sequence, at the blue end of the horizontal branch (HB) or mixed with post-HB stars (Heber 2009). Even though they lie mostly between the main sequence and the white dwarf (WD) cooling sequence in the Hertzsprung-Russell diagram, they are contaminated with WDs.
The most accepted formation scenario is of an extreme HB star that lost its envelope after the He-burning phase, evolving directly to the WD-cooling sequence by avoiding the asymptotic giant branch (AGB) (Wade et al. 2005;Heber 2009;Fontaine et al. 2012). These objects experienced mass transfer followed by common envelope ejection in a binary system leaving the stellar core exposed (Green et al. 2008;Geier et al. 2010). More recent investigations support this scenario, such as the observational evidence of hot O-and B-type subdwarf formation (sdO and sdB, respectively) from the product of binary interactions (Pelisoli et al. 2020). Also, the expected evolutionary path suggests that sdOs evolve from sdBs (Heber 2009). An alternative evolutionary scenario suggests that the merger of two low-mass He-core WDs might form isolated H-rich B-type hot subdwarfs (Hall & Jeffery 2016;Schwab 2018).
A small fraction of sdBs in close binaries (Wade et al. 2005) fall out of the red giant branch before He ignition, perhaps because they are low-mass stars that do not support the burning of He in their cores (Heber 2009).
The hot subdwarfs can be separated into three main groups according to Saffer et al. (2001). i) Those nonbinaries that have no detectable spectral lines from a cool companion, and show only small or insignificant velocity variations (35%). ii) Those single-lined spectroscopic binaries which have significant or large velocity variations and probable orbital periods on the order of a day (45%). iii) Those showing additional spectral lines from a cool FGK main-sequence or subgiant companion that have slowly varying or nearly constant velocities, indicating periods of many months to years.
There is observational evidence indicating that subdwarfs are the most likely candidates to account for the ultraviolet (UV) upturn observed on spectroscopic and photometric analyses of globular clusters and elliptical galaxies (Yi et al. 1999;Busso et al. 2005;Green et al. 2008).
Analysis from asteroseismology (Charpinet et al. 2005), spectroscopy (Wade et al. 2005;Dorsch et al. 2018Dorsch et al. , 2019, and photometry (Johnson et al. 2014) of field subdwarf stars found nonsolar He abundances. Herich and He-poor sequences, as a function of effective temperature, were proposed by Edelmann et al. (2003) and updated by Németh et al. (2012) and Lei et al. (2019). Those sequences also support the evidence that He-enriched sdOs are more common than the Hedeficient ones (Heber 2009).
It is ideal to have homogeneous and widely available model spectra for these elusive systems, in order to improve our knowledge of their fundamental physical parameters and evolution. Such models should be computed as fully as possible regarding the radiative transfer, considering non-local thermodynamic equilibrium (NLTE), fully blanketed atmosphere structures Lanz et al. 1997).
Considerable improvements in NLTE atmosphere models have been achieved over the last five decades. The first NLTE spectra of hot atmospheres in the early 1970s (Mihalas & Hummer 1974) were built based on pure H atmosphere structures, and the metal lines were considered only in the spectral synthesis. Analysis by  took H, C, and Fe into account for the opacity of hot and high-gravity NLTE atmosphere models. They found critical differences in effective temperature determinations and line profiles arising from the inclusion of explicit atomic species (i.e. with the kinetic equilibrium equation solved) in the atmospheric structure. The subdwarfs' atmospheres have more evident NLTE effects due to their lower surface gravities (when compared to WDs; Lanz et al. 1997). They built spectra considering full line blanketing with H/He, CNO, Si, Fe, and Ni as implicit ions and concluded that this approach significantly improved the line profile analysis. Hubeny et al. (1998) showed that NLTE metal-line-blanketing effects produce a conspicuous difference (typically 20%) in the line profiles seen in the UV spectra of hot WDs. By considering explicit metal ions such as H/He, CNO, Si, and Fe in the atmosphere structure a better emergent flux match was achieved.
Overall, the subdwarfs' spectra are built for specific targets and studies. We have hundreds of hot subdwarf stars cataloged (Geier et al. 2017;Heber et al. 2017;Geier et al. 2019;Geier 2020), with their spectral anal-ysis focused on classification and/or kinematic properties (Luo et al. 2019). A lot of work was done to analyse the He-abundance (Fontaine et al. 2019) sequences (Lei et al. 2019). The sequences also show differences between field subdwarfs and extreme HB subdwarfs in globular clusters (Lei et al. 2020).
NLTE structures and subdwarf spectra were calculated by Nemeth et al. (2014), covering specific temperature and gravity ranges. The opacities of H, He, and a few metal ions were included to compute the atmosphere structure in these models. However, only a small and nonhomogeneous collection of detailed subdwarf model spectra is currently available. High-temperature spectral grids may be also calculated by other codes e.g. the Tübingen NLTE Model Atmosphere package (Rauch et al. 2018) In this work we present an extensive high-spectralresolution grid of NLTE synthetic spectra for subdwarfs in the optical and UV. The models are fully line blanketed with H, He, C, N, O, Ne, Mg, Al, Si, S, and Fe as opacity sources (see table 1). We also consider nonsolar chemical abundances to better sample the subdwarf parameter space. Low-temperature convective atmospheres were also included in the square grid. The complete grid is available on our SpecModel website 2 as well as on the SVO's Theoretical Spectra Web Server 3 and the Vizier 4 databases.
Section 2 describes the atmosphere structure models, in particular the differences between the LTE and NLTE assumptions made. Section 3 describes the synthetic spectra, comparing the results with observational data. Section 4 describes the construction of the synthetic magnitudes and their application in color-color diagrams. Section 5 presents the concluding remarks.

ATMOSPHERE STRUCTURE
The atmosphere structure models were computed by TLUSTY code v205 and v208 (for convective 10,000 K models ;Hubeny 1988;Hubeny & Lanz 2011a). It calculates a self-consistent solution of the equations describing the radiative (or radiative plus convective) transfer and physical state of the atmosphere. The geometry of the model is plane-parallel with homogeneous chemical abundances .
The most popular subdwarf models available have been built in LTE. In LTE all energy partitioning such as atomic, ionic, and molecular level populations is given by Saha-Boltzmann equations, being defined by the local temperature. The LTE conditions may differ from those derived in actual statistical equilibrium, therefore level populations in NLTE were calculated. Departures from LTE may be specially significant for high-temperature subdwarf atmospheres Lanz et al. 1997) as also found in DAtype WD atmospheres (Levenhagen et al. 2017).
The starting point was to compute a LTE gray structure model, which is then used as an input to build the NLTE line-blanketed models (Hubeny 1988;Hubeny et al. 1994). We considered H, He, C, N, O, Ne, Mg, Al, Si, S, and Fe as explicit ions (see table 1) contributing to the opacity. The set of ionization states for the selected atomic species were chosen on the basis of ionization energies, model-effective-temperatures and candidate-permitted transitions from Willians & Livio (1995). The atomic data of each ion included as an explicit species in the NLTE atmosphere models are summarized in table 1 (see details in Lanz & Hubeny 2003, Lanz & Hubeny 2007. As described by Lanz & Hubeny (2007), the LTE fluxes are about 10% higher than the NLTE predictions, and this is most noticeable in the near-UV. Lanz & Hubeny (2007) also mentioned that the lower NLTE fluxes result from the overpopulation of the H i (n = 2) level at the depth of formation of the continuum flux, hence implying a larger Balmer continuum opacity.
The adopted abundances are expressed as solar ([Fe/H] = 0.0; Asplund et al. 2009) or Galactic halo ([Fe/H] = -1.5 and [α/Fe] = +0.4) as well as Herich (Edelmann et al. 2003) and He-poor abundances (Németh et al. 2012;Lei et al. 2019). Our choice of using solar and halo metallicity is not intended to describe subdwarf typical abundances. Instead, these solar and halo points just set a broad Z range, probing its effect on the continuum and line profiles. It is well established that atmospheric subdwarf metal abundances have little or virtually no memory of their parent main sequence abundances (e.g. Moehler 2001;Németh et al. 2012;Byrne et al. 2018). Metals at the subdwarf surface vary by large amounts depending on previous binary evolution and diffusion processes, leading to diverse abundance patterns (e.g. O'Toole & Heber 2006) with loose observational constraints. He-rich and He-poor abundance variations were chosen based on linear fits as a function of temperature based on Lei et al. (2019) (see table 2). Note that TLUSTY (Hubeny 1988;Hubeny & Lanz 2011a) uses the solar chemical abundances as de-fined by Grevesse & Sauval (1998) as default values, and we rearranged all the chemical species to the solar abundance as defined by Asplund et al. (2009). The atmosphere structures (available upon request to the authors) may be interpolated in metallicity and specific modified metal abundances used in custom spectral synthesis.
The time spent on computational work increases and the convergence is more difficult when the complexity of the heavy elements is included in the TLUSTY atmosphere code (Hubeny 1988;Hubeny & Lanz 2011a). As a reference, each one of the 96 structures were computed and analyzed in a typical time scale of a day.
The broadening of lines was treated carefully considering the Lyman and Balmer line profile tables from Tremblay & Bergeron (2009). The high-frequency limit for lines that are taken into account in the opacity sampling was set to 6 × 10 11 × T eff or the highest boundfree edge frequency, if greater. Opacity-sampling steps smaller than the thermal broadening were set. The microturbulent velocity was set to 10 km s −1 . The opacitysampling approach was used to compute the superline cross sections with better accuracy, also, the iron peak lines were treated in a line-blanketed model.
The only exception to the description above was the cooler atmosphere models with 10,000 K which require a convective treatment (Fontaine et al. 1981;Bergeron et al. 1992); a mixing-length theory parameter equal to the pressure-scale height (α MLT = 1.0) was used. For the 10,000 K case, the final structure model was computed in LTE without a microturbulent velocity.
Most models were computed using the hybrid complete-linearization and accelerated lambda-iteration (CL/ALI) method Hubeny 2003), which is the default procedure for computing fully consistent, NLTE metal-line-blanketed atmosphere models in TLUSTY (Hubeny 1988). When the convergence was not achieved directly from CL/ALI we computed the atmosphere structure models using the Rybicki scheme (Hubeny & Mihalas 2014), which is also suitable for treating high-temperature and high-gravity atmospheres, with better convergence behavior in some cases.
The effective temperature as function of the mass depth is shown in figure 1. The outermost mass depth corresponds to the Rosseland optical depth τ = 10 −7 while the innermost depth has τ = 100, logarithmically sampled in 70 layers. The color-bar scale indicates the effective temperature from bottom 10,000 K in yellow to top 65,000 K in purple. The hottest structures are more extended toward the inner thick region compared to the coolest model. The structure models in figure 1 (a) have a surface gravity (log g) equal to 4.5, so they are more extended if we compare them to the structures in figure 1 (b), which represent more compact objects with log g equal to 6.5. The dotted lines represent the starting model computed using LTE assumptions, where we can see the almost isothermal behavior in the external region for models with different effective temperatures. On the other hand, the dashed lines represent the structure computed using the NLTE assumptions and line blanketing of the explicit species. In the latter, a small temperature inversion on their outer regions can be seen in some cases, which is important to the line-core profile formation. The inner region converges to the solution with minor differences if we compare the LTE and NLTE structures. Note that mass depths are shown, if we look at a specific line optical depth the differences between LTE and NLTE structure models are much more significant.
The structure models are highly sensitive to the chemical abundances via opacities, electron density, and mass density. The explicit atoms and ions included to compose the level and superlevels are the most important ingredients to compute detailed NLTE structure models. They have an impact on the cross section computation, justifying the careful choice of the most important ions for each model. The electron density and the mass density are shown in figure 2, having a similar profile as a function of the mass depth. The effects of the surface gravity in the atmosphere height are evident as already mentioned above. The coolest structure models such as the example of 10,000 K in figure 2 a) and b) have an inversion on the density profiles near log [Depth (mass)] = -1, this is due to the effect of convection that was considered for  these specific temperatures. The hot structure models as the example of 65,000 K in figure 2 c) and d) have densities with a linear dependency of the mass depth. An Inglis-Teller diagram for the grid was produced, which is a classical tool for evaluating model sequences' behavior over a gravity range. It involves the electron density and the maximum n level that originates a distinguishable Balmer absorption line, which is useful for diagnosing our model's quality. We interpolated the model electron density to a specific optical depth τ = 0.1, counting the highest visible term in the final synthetic spectra. A linear least-squares fit of the electron density, and the maximum number of the absorption lines near the Balmer's discontinuity, is shown in figure 3.  We also performed convergence analyses for all the 96 subdwarf-structure models computed. The most critical convergence criterion is the magnitude of the relative changes of the components of the state vector, which is defined as a set of all structural parameters (e.g. tem-perature, particle number densities, and the mean radiation intensities in discretized frequency points) in a given discretized depth point (Hubeny & Lanz 2011a). The necessary condition for convergence is that the maximum relative change of all state vector components in all the 70 depths is smaller than 10 −3 . However, a supplementary condition is the conservation of the total flux concerning the total theoretical flux, σT eff 4 . The output parameters such as the number of depths, column mass, temperature, and densities in each depth and relative changes between iterations are used in the convergence analysis, as well as in the emergent flux in all frequency points used by the SYNSPEC code (Hubeny & Lanz 2011b) to build synthetic spectra.

SYNTHETIC SPECTRA
We computed the grid of synthetic spectra with the SYNSPEC code designed to synthesize the emergent spectra from atmosphere model structures. We used the NLTE atmosphere structure from TLUSTY models (Hubeny 1988;Hubeny & Lanz 2011a) described in section 2. We considered NLTE assumptions for the evaluation of the level and superlevel populations. The line profiles were carefully treated in a special computation for H and He (Lanz & Hubeny 2003, 2007. The reference atomic line list is Coelho (2014) The hydrogen and hydrogenic lines are treated as a part of the continuum and their profiles are computed using tables from Tremblay & Bergeron (2009). The quasi-molecular satellites of Lα, Lβ, and Hα (λ = 1215.67, 1025.18 and 6562.79Å, respectively), are considered. In that case, additional input files containing the corresponding data were used (Allard et al. 2009).
The four He i triplet lines, (λ = 4026, 4387, 4471, and 4921Å) were treated using special line-broadening tables (Barnard et al. 1974;Shamey 1969). The He ii lines are treated as the approximate hydrogenic ion by analytical values of the Stark + Doppler profile (Hubeny et al. 1994), which improves the accuracy of the line profile for T eff > 10,000 K, and the line profiles are given by the Stark-broadening tables of Schöning & Butler (1989). We are also considering Stark broadening computed by Tremblay & Bergeron (2009).
The spectral grid coverage is between 1000 -10,000Å with steps of 0.01Å. The resulting spectrum was subsequently processed with the ROTIN code (  2011b), which resamples the original synthetic spectrum but considering that no rotational velocity and no instrumental degradation have been taken into account.
We show the normalized emergent flux (in arbitrary units) given as a function of wavelength for the spectral coverage from 1000 to 9000Å in figure 4. A forest of lines is seen at this particular sampling, where no instrumental degradation is included. They are synthetic spectra with T eff equal to 65,000,45,000,35,000,30,000,25,000,20,000,15,000 and 10,000 K (from top to bottom) and log g [cgs] = 5.5 with halo and He-rich abundances.
In figure 5 we can compare the different spectral types in the near-UV region between 1000 and 2000Å, where for visual effect it is degraded to a gaussian instrumental profile with full width at half maximum (FWHM) equal to 5Å. The hotter spectra have a bluer continuum and the He lines are dominant. The Stark broadening is more evident on the coolest spectra, where a stronger Lyman series is present, besides its quasi-molecular absorption features.
In figure 6 we can compare the different spectral types on the Balmer break region between 3500 -6750Å, where for visual effect the resolution is degraded by a gaussian instrumental profile with FWHM equal to 5Å. As expected, the H is mostly ionized and the Balmer series is weaker in the hotter spectra. The coolest models are not favorable to the formation of strong He lines. Moreover, the spectra presented in figure 6 follow the Heabundance sequence from Németh et al. (2012), which shows low He abundance even for the He-rich sequence of cool subdwarfs.
With the purpose of illustrating the grid, we present a simple comparison of model spectra with the spectra of the Lamost 442708048, Lamost 183405148, and HD 4539 subdwarfs. No interpolation in the grid was performed as the intent is not to determine exact abundances or stellar parameters, but rather to exhibit the grid potentialities. The targets were selected from the subdwarf catalog in Lei et al. (2019) and Hubble Space Telescope's (HST)/Space Telescope Imaging Spectrograph pc. The closest model spectrum adopted corresponds to T eff = 45,000 K, log g = 5.5 [cgs], [Fe/H] = 0, and log (n He /n H ) = 0.92 (see figure 7 (b)). HST/STIS Legacy Archive data on the bright subdwarf HD 4539 were used to illustrate the models in the FUV. Parameters for this sdB are as follows: T eff = 26,000±500 K, log g = 5.2±0.1, log (n He /n H ) = -2.32±0.05, E(B − V ) = 0.04 ± 0.01 (Sale et al. 2008) and d = 171.6 ± 2.1 pc (Gaia EDR3; Gaia Collaboration et al. (2021)). The models were scaled to match the continuum flux in the middle of each spectral range. Hepoor branch and halo low-Z abundances were assumed here while T eff and log g were linearly interpolated from model nodes to the literature values (see figure 8). Rota- tion, instrumental resolution, and exact chemical composition were neglected, which would explain the significant differences found in the lines and continua offsets. By using the Gaia distances and reddening values above, stellar radii compatible with typical values for subdwarfs in eclipsing binaries (e.g. Rebassa-Mansergas et al. (2019) could be found. A line and/or continuum fitting of spectra can be performed with the XTGRID facility (Nemeth 2019).

SYNTHETIC MAGNITUDES
Synthetic magnitudes have been computed for several standard photometric bands to trend the grid's behavior in the color indices space and provide a comparison with photometry data. We used the filter response functions available on the Filter Profile Service at the Spanish Virtual Observatory 7 to convolve our synthetic integrated fluxes in the AB and Vega systems. The filters were interpolated in the spectra wavelength steps within the Shannon-Whittaker scheme. Photon-counting integrated fluxes were assumed (Bessell et al. 1998). Vega's zero-points were previously evaluated from Calspec's standard spectrum (Bohlin et al. 2014). Figure 9 illustrate our models in two color-color panels: panel (a) was chosen to trace the overall continuum inclination, and panel (b) traces the behavior of the Balmer jump. In panel (a) we show F469N -F673N versus FQ757N -FQ750N colors of the HST/Wide Field Camera 3 (WFC3) photometric system. The dependence on effective temperature is clear. Those indices measure how lower-temperature synthetic models have a flatter continuum in the optical region. The point-size scale represents the surface gravity, where larger symbols relate to lower gravity (log g = 4.5), and smaller symbols stand for higher gravity (log g = 6.5), which also reveal a linear dependence in the color space. Finally, the circles represent solar abundance and downpointing triangles represent low halo metallicity, which has no clear separation on the diagram.
In figure 9 (b) the Strömgren photometric system indices were used. Zero-points were derived from Vega ubvy magnitudes (Hauck & Mermilliod 1998) and its calibrated spectrum from the STScI CALSPEC database (Bohlin et al. 2014). The color-color diagram of figure 9 (b) was constructed with u (349.6nm)v (410.3nm) and v -b (466.6nm) colors. It shows a temperature dependence of the Balmer break, which is more evident and also gravity-dependent for cooler atmospheres.
The study of the 5874 hot subdwarf stars with Gaia Data Release 2 by Geier 2020 is the most complete sample of subdwarfs. From this catalog we selected data from 3450 targets (1482 with determined effective temperature and colors) observed by the Sloan Digital Sky Survey (SDSS) photometric system to compare with the grid synthetic magnitudes. The color-color diagram in figure 10 is composed of g -r and u -g colors without any color cutoff. Observational data with a determined effective temperature are shown as dots scaled by hue, while those with an undetermined effective temperature are gray crosses. Our synthetic colors are displayed in the same hue, style, and size scales as in figure 9, matching the sdO, sdB, and sdOB previously classified by Geier 2020. The upper data sequence represents the subdwarf composite binaries with main-sequence stars companions.  Figure 10. Color-color diagram (g -r vs. u -g) composed of observational data on subdwarfs from SDSS (see text) with a hue scale representing determined effective temperatures, while the gray crosses shown an undetermined effective temperature. The synthetic colors from our grid are plotted with the same hue, style, and size scales as in figure 9.

SUMMARY
We presented a grid of NLTE, fully blanketed theoretical spectra and synthetic photometry for hot and moderately cool subdwarf stars. The atmosphere models were computed considering line blanketing of H, He, C, N, O, Ne, Mg, Al, Si, S, and Fe. The effective temperatures are T eff = 10,000, 15,000, 20,000, 25,000, 30,000, 35,000, 45,000, and 65,000 K, while the surface gravities are log g [cgs] = 4.5, 5.5, and 6.5. The two representative chemical abundances are solar and Galactic halo, each one with two extreme scenarios for He-rich and He-poor stellar atmospheres.
The main differences between LTE and NLTE atmosphere structures are shown for these objects. They have significant differences in the outermost atmospheres, leading to distinct line-core profile formation.
The complete high-resolution spectral synthesis is performed from the UV to near-IR (1000 to 10,000Å) in 0.01 A steps. We provided an illustrative analysis for the UV and the optical regions by comparing our models with observed spectra from LAMOST and HST/STIS Legacy Archive data.
The behavior of the color indices were analyzed using the HST/WFC3 and the Strömgren photometric systems. A clear separation in effective temperature can be seen, as well as gravity for lower-temperature models, as provided by the Balmer discontinuity. We also matched our synthetic magnitudes against SDSS subdwarf data with fair agreement.
These results pave the way for both spectroscopic and photometric analyses of fundamental parameters in isolated or binary objects which, in turn, may provide a more detailed insight into the atmosphere models themselves. In addition, classical stellar population synthesis analysis can make use of the homogeneous spectral grid to better understand the blue and UV properties of old stellar populations. The full spectral grid and synthetic indices are available in the IAG-USP 2 , SVO 3 , and Vizier 4 databases.

ACKNOWLEDGMENTS
We thank Ivan Hubeny who kindly released the new versions of TLUSTY 208 and SYNSPEC 54 and boosted the discussions about the convective models. We thank the anonymous referee for refereeing this paper and for the suggestions that improved the work.
This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior -Brasil (CAPES), Finance Code 001. M.P.D. acknowledges support from CNPq under grant #305033. P.C. acknowledges support from Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) under grant #310041/2018-0 and from Fundação de Amparo a Pesquisa do Estado de São Paulo (FAPESP) process #2018/05392-8. This research has made use of the SVO Filter Profile Service (http://svo2.cab.inta-csic.es/theory/ fps/) supported from the Spanish MINECO through grant AYA2017-84089. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC; https://www.cosmos.esa.int/web/gaia/ dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.