ensemble of non-Hermitian but real random matrices. Each element of real
random matrices is independent and produced by the same Gaussian distribution.
For the GinSE, non-Hermitian random matrices are defined to satisfy
$\sigma_{y}H^{*}\sigma_{y}=H$ with

$\displaystyle\sigma_{y}=\begin{pmatrix}0&-\text{i}\\\ \text{i}&0\\\
\end{pmatrix},$ (S.115)

so that $H$ has the following symplectic structure:

$\displaystyle H=\begin{pmatrix}X&Y\\\ -Y^{*}&X^{*}\\\ \end{pmatrix},$ (S.116)

where $X$, $Y$ are generic non-Hermitian random matrices with no constraint.

As shown in Figs. S.6 and S.7, the DoS of non-Hermitian disordered systems in
class AI has a sharp peak on the real axis, while the DoS in class AII has a
soft gap on the real axis. In order to understand these behavior, we study the
DoS of the GinOE and GinSE. The eigenenergy distribution and the DoS along the
imaginary axis are shown in Fig. S.9. The DoS around the real axis for non-
Hermitian disordered Hamiltonians in classes AI and AII is qualitatively
consistent with the DoS of the GinOE and GinSE, respectively.

(a) $E=0$, (2,3,0,1)

(b) $E=0.5{\rm i}$, (3,3,0,1)

Figure S.1: Normalized localization lengths $\Lambda$ as a function of the
disorder strength $W$ for the 3D O(1) model in class AI at (a) $E=0$ and (b)
$E=0.5\text{i}$. The points with the error bars are the numerical data with
the different system sizes $L$. The colored curves are the fitted curves with
the expansion order $(m_{1},n_{1},m_{2},n_{2})$.

(a) 2D, (1,3,0,1)

(b) 3D, (2,3,0,1)

Figure S.2: Normalized localization lengths $\Lambda$ as a function of the
disorder strength $W=W_{r}=W_{i}$ for the (a) 2D and (b) 3D SU(2) model in
class AII† with $E=0$. The points with the error bars are the numerical data
with the different system sizes $L$. The colored curves are the fitted curves
with the expansion order $(m_{1},n_{1},m_{2},n_{2})$.

(a) 2D, $E=0.01{\rm i}$, (3,3,0,1)

(b) 3D, $E=0$, (2,3,0,0)

(c) 3D, $E={\rm i}$, (1,3,0,1)

Figure S.3: Normalized localization lengths $\Lambda$ as a function of the
disorder strength $W=W_{r}=W_{i}$ for (a) the 2D SU(2) model in class AII with
$E=0.01\text{i}$, (b) the 3D SU(2) model in class AII with $E=0$, and (c) the
3D SU(2) model in class AII with $E=\text{i}$. The points with the error bars
are the numerical data with the different system sizes $L$. The colored curves
are the fitted curves with the expansion order $(m_{1},n_{1},m_{2},n_{2})$.

(a) $E=0$, 2D, class CII†, (1,3,0,0)

(b) $E=0$, 2D, class DIII, (1,2,0,0)

Figure S.4: Normalized localization lengths $\Lambda$ as a function of the
disorder strength $W\equiv W_{i}$ for the 2D SU(2) model in (a) class CII† and
(b) class DIII. The points with the error bars are the numerical data with the
different system sizes $L$. The colored curves are the fitted curves with the
expansion order $(m_{1},n_{1},m_{2},n_{2})$.

(a) 3D class AI

(b) 3D class AII

(c) 2D class AII

(d) 2D class AII, $E=0.01{\rm i}$

Figure S.5: Phase diagrams of the O(1) and SU(2) models for (a) 3D class AI,
(b) 3D class AII, and (c) 2D class AII in terms of the disorder strength $W$
and the imaginary part of eigenenergies $E$. The real part of $E$ is set to
$0$. The blue squares and red circles are the phase boundaries for the
Anderson transitions. The phase boundaries are determined by the localization
lengths. (d) Normalized localization length as a function of $W\equiv
W_{r}=W_{i}$ for the 2D class AII model at $E=0.01\text{i}$.

(a) $W=0$

(b) $W=0$

(c) $W=10$

(d) $W=10$

(e) $W=20$

(f) $W=20$

Figure S.6: (a), (c), (e) Eigenenergy distributions in the complex plane. Each
point corresponds to each eigenenergy in one sample. The color of the points
describes $1/I$ with the inverse participation ratios $I$ for the
corresponding eigenmodes. (b), (d), (f) Density of states (DoS) for the
imaginary part of eigenenergies. The numerical calculations are performed for
the 3D O(1) model in class AI with the cubic system size $L=20$, under the
periodic boundary conditions, and with the disorder strength $W=0,10,20$.

(a) 2D

(b) 2D

(c) 3D

(d) 3D

Figure S.7: (a), (c) Heat maps of eigenenergy density in the complex plane.
(b), (d) Density of states (DoS) for the imaginary part of eigenenergies.
Eigenenergies are calculated for the 2D and 3D SU(2) model in class AII under
the periodic boundary conditions, with $W_{r}=W_{i}=1$, and over the 640
samples with different disorder realizations. The system sizes are $70\times
70$ for 2D and $16\times 16\times 16$ for 3D.

(a) 2D

(b) 3D

Figure S.8: Eigenenergy distributions in the complex plane. Each point
corresponds to each eigenenergy in one sample. The color of the points
describes $1/I$ with the inverse participation ratios $I$ for the
corresponding eigenmodes. The eigenenergies are calculated for the SU(2)
models in class AII† with the system sizes $60\times 60$ for 2D and $16\times
16\times 16$ for 3D. The periodic boundary conditions are imposed in all the
directions for both 2D and 3D. The disorder strength is set to $W_{r}=W_{i}=4$
for 2D and $W_{r}=W_{i}=7$ for 3D.

(a) Ginibre Orthogonal ensemble

(b) Ginibre Orthogonal ensemble

(c) Ginibre Symplectic ensemble

(d) Ginibre Symplectic ensemble

Figure S.9: Heat maps of eigenenergy density in the complex plane for (a) the
Ginibre orthogonal ensemble and (c) the Ginibre symplectic ensemble. The color
describes the density of states (DoS) in the complex plane. The DoS for the
imaginary part of eigenenergies of (b) the Ginibre orthogonal ensemble and (d)
the Ginibre symplectic ensemble. The numerical data come from the 640
realizations of random matrices with the size $10^{4}\times 10^{4}$.


# Explaining the differences in massive star models from various simulations

Poojan Agrawal1,2,3, Dorottya Szécsi,4,5, Simon Stevenson1,2, Jarrod Hurley1,2
1 Centre for Astrophysics and Supercomputing, Swinburne University of
Technology, Hawthorn, VIC 3122, Australia
2 OzGrav: The ARC Centre of Excellence for Gravitational Wave Discovery,
Hawthorn, VIC 3122, Australia
3 McWilliams Center for Cosmology, Department of Physics, Carnegie Mellon
University, Pittsburgh, PA 15213, USA
4 Institute of Astronomy, Faculty of Physics, Astronomy and Informatics,
Nicolaus Copernicus University, Grudziądzka 5, 87-100
Toruń, Poland
5 I. Physikalisches Institut, Universität zu Köln, Zülpicher-Str. 77, D-50937
Cologne, Germany
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The evolution of massive stars is the basis of several astrophysical
investigations, from predicting gravitational-wave event rates to studying
star-formation and stellar populations in clusters. However, uncertainties in
massive star evolution present a significant challenge when accounting for
these models’ behaviour in stellar population studies. In this work, we
present a comparison between five published sets of stellar models from the
BPASS, BoOST, Geneva, MIST and PARSEC simulations at near-solar metallicity.
The different sets of stellar models have been computed using slightly
different physical inputs in terms of mass-loss rates and internal mixing
properties. Moreover, these models also employ various pragmatic methods to
overcome the numerical difficulties that arise due to the presence of density
inversions in the outer layers of stars more massive than 40 M⊙. These density
inversions result from the combination of inefficient convection in the low-
density envelopes of massive stars and the excess of radiative luminosity to
the Eddington luminosity. We find that the ionizing radiation released by the
stellar populations can change by up to $\sim$15 percent, the maximum radial
expansion of a star can differ between $\sim$100–2 000 R⊙, and the mass of the
stellar remnant can vary up to 20 M⊙ between the five sets of simulations. We
conclude that any attempts to explain observations that rely on the use of
models of stars more massive than 40 M⊙ should be made with caution.

###### keywords:

stars: massive – stars: evolution – gravitational waves – stars: formation –
stars: black holes – galaxies: stellar content

††pubyear: 2021††pagerange: Explaining the differences in massive star models
from various simulations–Explaining the differences in massive star models
from various simulations

##  1 Introduction

Stellar evolutionary model sequences serve as input for a broad range of
astrophysical applications; from star-formation (e.g. Gatto et al., 2017) to
galaxy evolution (e.g. Weinberger et al., 2020), from cluster dynamics (e.g.
Heggie & Hut, 2003) to gravitational-wave studies (e.g. Vigna-Gómez et al.,
2018). These sequences provide an easy and powerful way to account for both
individual stars (e.g. Schneider et al., 2014) and stellar populations (e.g.
Brott et al., 2011b) in a given astrophysical environment.

One dimensional (1D) model sequences (from now on: stellar models) can be
computed from first principles (Kippenhahn & Weigert, 1990) and have become a
household tool in astrophysical research. But when it comes to stars more
massive than $\sim$ 9 M⊙ – those that, despite being rare, provide the bulk of
the radiation, chemical pollution and the most exotic death throes in the
Universe (Woosley et al., 2002) – stellar models are still riddled with large
uncertainties.

High-mass stars are born less often than their low-mass counterparts
(Salpeter, 1955) and have comparatively shorter lives (Crowther, 2012).
Consequently, observational constraints on their evolution are more difficult
to obtain. The situation is further complicated by many massive stars being
observed to be fast rotators (Ramírez-Agudelo et al., 2013), which breaks down
perfect symmetry, and to have a close-by companion star (Sana et al., 2012),
breaking the assumption of perfect isolation. Even for isolated, non-rotating
single stars, the physical conditions both inside (Heger et al., 2000) and
around (Lamers & Cassinelli, 1999) the star are so peculiar and complex that
developing appropriate numerical simulations becomes highly challenging. This
is why the evolution of massive stars remains an actively studied field to
this day.

Much progress has been made in the last few decades concerning massive stars
and their evolution. Mass loss in the form of high-velocity winds from massive
stars is being intensively studied and accounted for in the models (Smith,
2014; Sander & Vink, 2020). Observations of massive stars from the Large and
Small Magellanic Clouds are being used to constrain the efficiency of interior
mixing processes (Brott et al., 2011a; Schootemeijer et al., 2019). 1D stellar
models have also been updated to account for the effects of rotation (Maeder,
2009; Costa et al., 2019) and magnetic fields (Maeder & Meynet, 2005;
Takahashi & Langer, 2021) which can significantly change their evolutionary
paths (Walder et al., 2012; Petit et al., 2017; Groh et al., 2020).

Despite the progress, there are still many open questions surrounding the
lives of massive ($\gtrsim$9 M⊙) and ‘very’ massive (here designated as
$\gtrsim$40 M⊙) stars, and in the absence of well-defined answers, stellar
evolution codes make use of different assumptions. Earlier studies comparing
models of massive stars from different codes (e.g. Martins & Palacios, 2013;
Jones et al., 2015) have already established that the differences in the
physical parameters such as mixing and mass-loss rates adopted by various
stellar evolution codes can affect the evolutionary outcome of these stars.

Here we highlight another major uncertainty arising due to the numerical
treatment of low-density envelopes of very massive stars. These stars have
luminosities close to the Eddington-limit, so changes in the elemental
opacities during their evolution can lead to the formation of density and
pressure inversions in the stellar envelope (Langer, 1997). The presence of
these density inversions can cause numerical instabilities for 1D stellar
evolution codes. To deal with these instabilities, the codes use different
pragmatic solutions whose interplay with mixing and mass loss can further vary
the evolution of massive stars.

The role of the Eddington limit and the associated density inversions in
massive stars is well known within the stellar evolution community but remains
relatively unknown outside the field. With the surge in the use of massive
star models, for example, in gravitational-wave event rate predictions and
supernova studies, it has become important to be aware of this issue. Our goal
is to present the broader community with a concise overview, including how it
affects the evolutionary properties such as the radial expansion and the
remnant mass of very massive stars. To this end, we compare models of massive
and very massive stars from five published sets created with different
evolutionary codes: (i) models from the PAdova and TRieste Stellar Evolution
Code (PARSEC; Bressan et al., 2012; Chen et al., 2015); (ii) the MESA
Isochrones and Stellar Tracks (MIST; Choi et al., 2016) from the Modules for
Experiments in Stellar Astrophysics (MESA; Paxton et al., 2011); (iii) models
(Ekström et al., 2012) from the Geneva code (Eggenberger et al., 2008); (iv)
models from the Binary Population and Spectral Synthesis (BPASS; Eldridge et
al., 2017) project; and (v) the Bonn Optimized Stellar Tracks (BoOST; Szécsi
et al., 2020) from the ‘Bonn’ Code.

We describe the major physical ingredients used in computing each set of
models in Section 2 and Section 3. In Section 4, we compare the predictions
from each set of models in the Hertzsprung-Russell diagram and in terms of the
emitted ionizing radiation, as well as the predictions for the maximum radial
expansion of stars, and their remnant masses. Finally we draw conclusions in
Section 5.

##  2 Physical inputs

Table 1: Summary of input parameters used in the computation of the models of
massive stars from different codes. See Section 2 for details.

Stellar Model |  Z⊙ |  Hot Wind |  Cool Wind |  Wolf-Rayet Wind |  Convective boundary |  $\alpha_{\rm MLT}$ |  Overshoot type |  $\alpha_{\rm ovs}$ |  $\alpha_{\rm semi}$
---|---|---|---|---|---|---|---|---|---
BPASS |  0.020 |  Vink et al. (2000, 2001) |  de Jager et al. (1988) |  Nugis & Lamers (2000) |  Schwarzschild (1958) |  2.0 |  Pols et al. (1998) |  0.12e |  –
BOOST |  0.008a |  Vink et al. (2000, 2001) |  Nieuwenhuijzen & de Jager (1990) |  Hamann et al. (1995)b |  Ledoux (1947) |  1.5 |  step |  0.335 |  1.0
GENEVA |  0.014 |  Vink et al. (2000, 2001) |  de Jager et al. (1988)c |  Nugis & Lamers (2000) |  Schwarzschild (1958) |  1.6d |  step |  0.1 |  –
MIST |  0.014 |  Vink et al. (2000, 2001) |  de Jager et al. (1988) |  Nugis & Lamers (2000) |  Ledoux (1947) |  1.82 |  Herwig (2000) |  0.016e |  0.1
PARSEC |  0.015 |  Vink et al. (2000, 2001) |  de Jager et al. (1988) |  Nugis & Lamers (2000) |  Schwarzschild (1958) |  1.74 |  Bressan et al. (1981) |  0.5e |  –

  * Notes.

  * a

For calculating mass-loss rates and opacities, Z⊙$=0.017$ is used.

  * b

Reduced by a factor of 10.

  * c

For ${\rm logT_{eff}/K\leq 3.7}$, mass-loss rates from Crowther (2000) are
used.

  * d

For stars with initial mass $\geq$ 40 M⊙, $\alpha_{\rm MLT}=1.0$ is used but
with a different scale height (see Section 3).

  * e

The rough equivalent in the step overshooting formalism is 0.2, 0.25 and 0.4
for the MIST, PARSEC and BPASS models respectively.

###  2.1 Chemical composition

The chemical composition of the Sun is often used as a yardstick in computing
the metal content of other stars. However, the exact value remains
inconclusive and has undergone several revisions since 2004 (see Basu, 2009;
Asplund et al., 2021, for an overview). Therefore, different stellar models
often make use of different abundance scales.

The BPASS, PARSEC, Geneva and MIST models base their chemical compositions on
the Sun, while the BoOST models use a mixture tailored to the sample of
massive stars from the FLAMES survey (Evans et al., 2005) with
$Z_{Gal}=0.0088$. For stellar winds and opacity calculations, BoOST models use
$Z=0.017$ from Grevesse et al. (1996) as the reference solar metallicity. The
BPASS models use solar abundances from Grevesse & Noels (1993) with Z⊙$=0.02$.
The PARSEC models follow Grevesse & Sauval (1998) with revisions from Caffau
et al. (2011) and Z⊙$=0.01524$. Geneva models use Asplund et al. (2005)
abundances with Ne abundance from Cunha et al. (2006) and Z⊙$=0.014$ while,
finally, the MIST models base their abundance on Asplund et al. (2009) with
Z⊙$=0.01428$.

For the purpose of comparison here, we use $Z=0.014$ models for each set
except for BoOST where we use the Galactic composition, $Z=0.0088$ as the
closest match. Iron is an important contributor to metallicity, as numerous
iron transition lines dominate both opacity and mass-loss rates, therefore
directly affecting the structure of massive stars (e.g., Puls et al., 2000).
The stellar models compared here have similar iron content, with the
normalised number density ($A_{\rm Fe}={\rm log(N_{Fe}/N_{H})+12.0})$ ranging
from $7.40$ (for BoOST models) to $7.54$ (for MIST models).

###  2.2 Mass-loss rates

Stellar mass is a key determinant of a star’s life and evolutionary outcome.
It can, however, change as stars lose their outer layers in the form of
stellar winds, and through interactions with a binary companion. Consequently,
mass loss can affect not only the structure and chemical composition of the
star, but is also important in determining its final state (Renzo et al.,
2017).

For massive stars, the effects of mass loss are even more pronounced. The mass
loss experienced by hot massive stars (O type stars and Wolf–Rayet stars) is
known to be line-driven (Lamers & Cassinelli, 1999) while that of cool massive
stars (red supergiants, Levesque, 2017) is suggested to be dust-driven. Both
types of mass loss are an intensively studied subject. However, the complexity
of the problem of atomic and molecular transitions in the wind together with
the rarity of stars at these high masses means that the model assumptions are
usually based either on a few observations (a small sample of stars) or on
what we know about the wind properties of low-mass stars.

All models in this study follow Vink et al. (2000, 2001) for hot wind-driven
mass-loss. The PARSEC, Geneva, BPASS and MIST models follow de Jager et al.
(1988) for cool dust-driven mass-loss and Nugis & Lamers (2000) for mass loss
in the naked helium star phase. The Geneva models further switch to Crowther
(2000) for hydrogen-rich stars with ${\rm logT_{eff}/K\leq 3.7}$. The BoOST
models follow Nieuwenhuijzen & de Jager (1990) for cool winds and Hamann et
al. (1995, reduced by a factor of 10) for stars with surface helium mass
fraction $>\,0.7$. For computing mass-loss rates of stars with surface helium
mass fraction between 0.4 and 0.7, BoOST models linearly interpolate between
the mass-loss rates of Vink et al. (2001) and Hamann et al. (1995, reduced by
a factor of 10).

To account for the dependence of the mass-loss rates on the chemical
composition, BoOST, MIST and BPASS models scale the mass-loss rates by factor
of $Z^{0.85}$ (Vink et al., 2001) 111Note that for some models the factor
$Z^{0.69}$ is quoted, depending on whether the dependence of the terminal
velocity on Z is explicitly considered or not. Geneva and PARSEC models also
use additional mass-loss as described in Section 3.

###  2.3 Convection and overshooting

Internal mixing processes such as convection and overshooting play an
important role in determining both the structure and evolution of massive
stars (e.g., see Sukhbold & Woosley, 2014). Similar to mass loss, these
processes represent another major source of uncertainty in massive stellar
evolution (Schootemeijer et al., 2019; Kaiser et al., 2020). In 1D stellar
evolution codes convection is modelled using the Mixing-Length Theory (Böhm-
Vitense, 1958, MLT) in terms of the mixing length parameter
${\rm\alpha_{MLT}}$. However, 3D simulations suggest that convection in
massive stars might be more sophisticated and turbulent than described by MLT
(Jiang et al., 2015).

The BoOST, Geneva, PARSEC and BPASS models used here follow standard MLT (Cox
& Giuli, 1968) for convective mixing with mixing length parameter $\alpha_{\rm
MLT}$ = $(1.5,1.6,1.74,2.0)$ respectively. MIST follows a modified version of
MLT given by Henyey et al. (1965) with $\alpha_{\rm MLT}=1.82$. Convective
boundaries in PARSEC, Geneva and BPASS models are determined using the
Schwarzschild criterion (Schwarzschild, 1958). BoOST and MIST use the Ledoux
criterion (Ledoux, 1947) for determining convective boundaries with
semiconvective mixing parameters of 1.0 and 0.1, respectively. For determining
convective core overshoot, Geneva and BoOST use step overshooting with
overshoot parameter $\alpha_{ov}$ = $(0.1,0.335)$. MIST uses exponential
overshooting following Herwig (2000) with $\alpha_{ov}$ = 0.016. PARSEC uses
overshoot from Bressan et al. (1981) with $\alpha_{ov}$ = $0.5$. BPASS uses
the overshoot prescription from Pols et al. (1998) with $\alpha_{ov}$ = 0.12.
For comparison, the rough equivalent in the step overshooting formalism would
be 0.2, 0.25 and 0.4 for the MIST, PARSEC and BPASS models respectively (See
Choi et al., 2016; Pols et al., 1998; Bressan et al., 2012, for details of
each method.)

MIST and PARSEC also include small amounts of overshoot associated with
convective regions in the envelope. However, apart from modifying surface
abundances, envelope overshoot has a negligible effect on the evolution of the
star (Bressan et al., 2012). Rotational mixing also plays an important role in
the evolution of massive stars, however, for simplicity we only compare non-
rotating models for PARSEC, MIST, Geneva and BPASS. For BoOST, although we do
use slowly rotating (100 km s-1) models, this small difference in the initial
rotation rate is again not relevant from the point of view of the overall
evolutionary behaviour (Brott et al., 2011a).

Major input parameters used in each set of models are summarized in Table 1.

Figure 1: Hertzsprung–Russell diagrams of the massive single star models
analysed in this work. All models have near solar composition. Symbols mark
every $10^{5}$ years of evolution. Only the core-hydrogen- and core-helium-
burning phases are plotted. The dashed red line marks the observational
Humphreys–Davidson limit (Humphreys & Davidson, 1979) where relevant. The
tracks become more varied with increasing initial mass. This is because the
codes apply various treatments for the numerical instabilities associated with
the Eddington-limit proximity, cf. Section 3.

##  3 Eddington luminosity and the numerical treatment of density inversions

The Eddington luminosity is the maximum luminosity that can be transported by
radiation while maintaining hydrostatic equilibrium (Eddington, 1926). In the
low-density envelopes of massive stars changes in the elemental opacities
during the evolution of stars can cause the local radiative luminosity to
exceed the Eddington luminosity (Langer, 1997; Sanyal et al., 2015). To
maintain hydrostatic equilibrium, density and pressure inversion regions form
in the stellar envelope. In the absence of efficient convection (which is also
typical for the low-density envelopes, Grassitelli et al. 2016), this can lead
to convergence problems for 1D stellar evolution codes (Paxton et al., 2013).
Owing to numerical difficulties, the time-steps become exceedingly small,
preventing further evolution of the star. While less massive stars are only
affected by this process in their late evolutionary phases (Harpaz, 1984; Lau
et al., 2012), very massive stars can exceed the Eddington-limit already
during the core-hydrogen-burning phase (Gräfener et al., 2012; Sanyal et al.,
2015) and inhibit computation of their evolution. Therefore, 1D stellar
evolution codes have to employ various solutions to compute further evolution
of very massive stars.

In PARSEC models, density inversions and the consequent numerical difficulties
are avoided by limiting the temperature gradient such that the density
gradient never becomes negative (see Sec. 2.4 of Chen et al., 2015; Alongi et
al., 1993). Limiting the temperature gradient prevents inefficient convection
and the evolution of the stars proceeds uninterrupted. Also, the models
include a mass-loss enhancement following Vink (2011) whenever the total
luminosity of the star approaches the Eddington-luminosity.

MIST models suppress density inversions through the MLT++ formalism (Paxton et
al., 2013) of MESA. In this method, the actual temperature gradient is
artificially reduced to make it closer to the adiabatic temperature gradient
whenever radiative luminosity exceeds the Eddington luminosity above a pre-
defined threshold. This approach again increases convective efficiency,
helping stars to overcome density inversions. Additionally, radiative pressure
at the surface of the star is also enhanced in the MIST models to help with
convergence (Choi et al., 2016).

In the extended envelopes of massive stars, the density scale height is much
larger compared to the pressure scale height (which is typically used for
computing the mixing length). Therefore, setting the mixing length to be
comparable with the density scale height helps avoid density inversion
(Nishida & Schindler, 1967; Maeder, 1987). The Geneva models include this
treatment when computing models with initial masses greater than 40 M⊙ with
$\alpha_{\rm MLT}=1.0$ (see Sec. 2.3 of Ekström et al., 2012). Additionally,
the mass-loss rates for the models are increased by a factor of three whenever
the local luminosity in any of the layers of the envelope is higher than five
times the local Eddington luminosity.

BoOST models do not include any artificial treatment to prevent massive stars
from encountering density inversions. Instead, their models undergo envelope
inflation when massive stars reach the Eddington limit (Sanyal et al., 2015).
On encountering the density inversions in their envelopes, the computation of
very massive stars becomes numerically difficult. Further evolution of such
stars is then computed through post-processing. It involves removing layers
from the surface of the star (which would anyway happen due to regular mass
loss) while correcting for surface properties such as effective temperature
and luminosity (Szécsi et al., 2020).

BPASS models also allow density inversions to develop in the envelope of
massive stars. However, these models are able to continue the evolution
without numerical difficulties, most likely due to the use of a non-Lagrangian
mesh (see Stancliffe, 2006, for an overview) and the resolution factors being
lower than in other models (J.J. Eldridge, private communication, also see
Eggleton, 1973).

##  4 Comparing the models

###  4.1 Differences between models in the Hertzsprung–Russell diagram

The evolution of stars can be easily represented through tracks on the
Hertzsprung-Russell (HR) diagram, depicting the evolutionary paths followed by
a series of stars. Figure 1 presents the HR diagram of stars of various
initial masses from the five simulation approaches. The observational analogue
to the Eddington-limit is the Humphreys–Davidson limit or HD-limit (Humphreys
& Davidson, 1979). Since the luminosity of a star depends on its mass, more
massive stars also are more luminous. This means they can easily exceed the
Eddington-limit, develop density inversions and require the use of numerical
solutions as discussed in Section 3.

From Figure 1 we see that the tracks of the 25 M⊙ (or 24 M⊙ in some cases)
stars agree well during most of the evolution. This is because stars of this
mass do not exceed the Eddington limit and are thus not affected by the
related numerical treatments. The minor differences in their tracks are due to
the difference in physical inputs (Section 2) between the simulations. For
example, the differences in the position of the main-sequence (MS) hook
feature in the HR diagram arise due to the varied extent of convective
overshoot used in each set of models.

A 40 M⊙ star is clearly affected by the numerical treatment employed during
the post-main sequence phase of its evolution, as evidenced by the difference
in the tracks in the HR diagram shown in Figure 1. More massive stars, i.e,
those with initial masses 80/85 M⊙ and 120/125 M⊙, can exceed the Eddington-
limit in their envelopes while on the MS and therefore their simulations
differ significantly from each other. At these masses, the mass-loss rates can
be as high as 10-3 $-$ 10-4 M⊙ yr-1, completely dominating over every other
physical ingredient in determining the evolutionary path. While all tracks
have been computed with similar prescriptions for wind mass loss (cf. Section
2), the actual rates can be strongly modified by the numerical methods adopted
by each code in response to numerical instabilities (Section 3), resulting in
vast differences in the tracks.

Table 2: Ionizing photon number flux [s-1] in the Lyman continuum emitted on average by the stellar models during their lives, cf. Section 4.2. The last column provides the amount of Lyman radiation (number of photons [s-1]) that a 107 M⊙ population (e.g. a starburst galaxy or a young massive cluster in the Milky Way) containing these massive stars would emit.  Mini [M⊙] | 24/25 | 40 | 80/85 | 120/125 | pop.
---|---|---|---|---|---
PARSEC | $3.7\times 10^{48}$ | $1.3\times 10^{49}$ | $5.5\times 10^{49}$ | $1.0\times 10^{50}$ | $1.08\times 10^{54}$
MIST | $3.3\times 10^{48}$ | $1.5\times 10^{49}$ | $5.1\times 10^{49}$ | $1.1\times 10^{50}$ | $1.06\times 10^{54}$
Geneva | $3.5\times 10^{48}$ | $1.2\times 10^{49}$ | $4.6\times 10^{49}$ | $7.8\times 10^{49}$ | $9.27\times 10^{53}$
BPASS | $3.6\times 10^{48}$ | $1.3\times 10^{49}$ | $4.5\times 10^{49}$ | $7.7\times 10^{49}$ | $9.34\times 10^{53}$
BoOST | $3.7\times 10^{48}$ | $1.2\times 10^{49}$ | $4.4\times 10^{49}$ | $7.4\times 10^{49}$ | $9.14\times 10^{53}$

###  4.2 Ionizing radiation and synthetic populations

The ionization released by a stellar population in e.g. a cluster or galaxy is
influenced by the contribution of the most massive stars (Topping & Shull,
2015). As shown above however, these are the stars for which the simulations
give the most diverse predictions.

To demonstrate this effect, we calculate the ionizing radiation emitted by a
simple stellar population, supposing a Salpeter (1955) initial mass function
and a star-forming region of 107 M⊙ total mass – which is aimed to represent
either a typical starburst galaxy, or a young massive cluster in the Milky
Way. In Table 2 we list the Lyman photon flux predicted by the individual
stellar models analysed here. In the absence of spectral synthesis models
computed for all five sets (cf. Wofford et al., 2016), we have opted to simply
use black body estimation. To correct for optically thick winds, we follow the
method explained in Chapter 4.5.1 of Szécsi (2016, which relies on ). The
results of the simple population syntheses based on all five sets are also
reported in Tab. 2.

We find that, in terms of how much Lyman flux is emitted by a given synthetic
population, the model predictions can differ as much as $\sim$15 percent
between simulations. This supports earlier findings (e.g. Topping & Shull,
2015) that relying on the ionizing properties of massive stars from
evolutionary models should be done with caution. Indeed one should keep in
mind that the behaviour of the most massive models, those that dominate the
radiation profile of any star-forming region, is weighted with large
uncertainties – the source of which is the treatment of the Eddington limit,
explained in Section 3.

###  4.3 Predictions of maximum stellar radii

The radial expansion of a star plays a significant role in determining the
nature of binary interaction as it can lead to episodes of mass transfer in
close interacting binaries. Recent studies indicate that most of the massive
stars occur in binaries (Sana et al., 2012; Moe & Di Stefano, 2017).
Therefore, predictions of stellar radii become even more important for
determining the binary properties of massive stars.

Figure 2 shows the maximum radial expansion achieved by massive stars from
each simulation. For stars with initial mass up to 30 M⊙, all simulations
predict the formation of a red supergiant. The maximum difference in the
radius predictions here is $\lesssim$1000 R⊙. For higher initial masses, the
predictions for maximum stellar radii become more divergent as proximity to
the Eddington-limit increases and numerical treatments adopted by each code
modify the mass-loss rates.

The greatest difference in the maximum radius predictions ($\gtrsim$1000 R⊙)
occurs for stars with initial masses between 40 and 100 M⊙. Above $\sim$ 100
M⊙, stars have even higher mass-loss rates which can completely strip a star
of its envelope before it can become a red-supergiant. Such stars evolve
directly towards the naked helium star phase and have much smaller radii.
Therefore, for stars with initial masses more than 100 M⊙, the difference
between the maximum radius predictions by each simulation reduces to
$\lesssim$100 R⊙. The predictions in this mass range seem to further converge
into two groups: PARSEC and MIST represent one group predicting smaller radii
compared to the second group which consists of BoOST and BPASS models. Since
the maximum initial mass of a star in the Geneva set is 120 M⊙, we cannot
comment on the behaviour of their models outside of this range.

Figure 2: Maximum stellar radii as a function of the initial mass of the star.
Similar to Figure 1, differences in the physical inputs and the numerical
methods adopted by each code can lead to a difference of more than $1\,000\,$
R⊙ in predictions in terms of the maximum radial expansion achieved by the
stars.

###  4.4 Remnant mass predictions

Stellar evolutionary models provide an easy way of estimating the properties
of stellar remnants such as black holes and neutron stars, which are needed in
many fields including supernova studies (e.g. Aguilera-Dena et al., 2018;
Raithel et al., 2018), gamma-ray bursts progenitors (e.g. Yoon et al., 2006;
Szecsi, 2017), and gravitational-wave event rate predictions (e.g. Stevenson
et al., 2019; Mapelli et al., 2020).

Following Belczynski et al. (2010) we show in Figure 3 how the uncertainties
in the models we compare here also pose a challenge for the predictions of
remnant properties. Remnant masses have been calculated from the carbon-oxygen
core mass and the total mass of the star at the end of the core helium-burning
phase using the prescription of Belczynski et al. (2008) (same as the
StarTrack prescription in Fryer et al., 2012).

The remnant masses are heavily influenced by the modelling assumptions (cf.
Section 2) and the numerical methods (cf. Section 3) especially above Mini $=$
40 M⊙ where the most massive black holes are predicted. We find that the mass
of the black holes predicted by the different sets of models can differ by
$\sim$ 20 M⊙. The maximum black hole mass varies from about 20 M⊙ for BPASS
models to about 38 M⊙ for MIST and PARSEC models and 35 M⊙ for the BoOST
models. Since all the models we study have near-solar metallicity and
therefore rather high mass-loss rates, none of them stays massive enough at
the end of their lives to undergo pair-instability.

Similar variability is found in the core mass and the final total mass of the
star (from which the remnant masses have been calculated). BPASS models
consistently predict the lowest values of final core mass and the total mass
for most of the massive stars while BoOST predictions are the highest. MIST
and PARSEC models predict values in between. The Geneva models do not provide
information on core masses, hence we omit them from the analysis of the core
mass and the remnant mass. However, we do show their predictions for the final
total mass which, for stars up to 100 M⊙, are similar to the predictions of
BPASS models.

Similar to the predictions of maximum stellar radii in Figure 2, the
predictions of final masses for stars with initial mass $\gtrsim$ 100 M⊙
separate into two groups, although the groups are slightly different this
time, BPASS and MIST represent one group predicting smaller final masses
compared to the second group consisting of BoOST and PARSEC models.

---
Figure 3: Final masses of stars as a function of their initial mass, ${\rm
M_{ZAMS}}$. The top panel shows the mass of stellar remnants as predicted by
the BPASS, BoOST, MIST and PARSEC models. The middle panel shows the carbon-
oxygen core mass and the bottom panel shows the total mass of the star at the
end of the core helium-burning phase, as used in the calculation of the
remnant masses. Differences in the evolutionary parameters for massive stars
can cause variations of about 20 M⊙ in the remnant masses between the stellar
models from various simulations.

##  5 Conclusions

We compare 1D evolutionary models of massive and very massive stars from five
independent simulations. Focusing on near-solar composition, we find that the
predictions from different codes can differ from each other by more than 1000
R⊙ in terms of maximum radial expansion achieved by the stars, by $\sim$15
percent in terms of ionizing radiation, and about 20 M⊙ in terms of the
stellar remnant mass. The differences in the evolution of massive stars can
arise due to physical inputs like chemical abundances, mass-loss rates, and
internal mixing properties. However, very massive stars, that is stars with
initial masses 40 M⊙ or more, show a larger difference in evolutionary
properties compared to lower mass stars. For these stars, the differences in
the evolution can be largely attributed to the numerical treatment of the
models when the Eddington-limit is exceeded in their low-density envelopes.

The different methods used by 1D codes to compute the evolution of massive
stars beyond density inversions (or to avoid the inversions) can modify the
radius and temperature of the star, and can therefore affect the mass-loss
rates. A phenomenological justification for the mass loss enhancement comes
from the fact that there are stars observed with extremely high, episodical
mass loss, i.e. luminous blue variables (Bestenlehner et al., 2014; Sarkisyan
et al., 2020). However, other studies, such as the recent measurement of a
approximately 20 M⊙ black hole in the Galactic black hole high-mass X-ray
binary Cyg X-1 (Miller-Jones et al., 2021), suggest that the mass-loss rates
for massive stars at near-solar metallicity may be lower than usually assumed
in the 1D stellar models (Neijssel et al., 2021). The exact nature of wind
mass loss for very massive stars remains disputed (Smith & Tombleson, 2015).
Moreover, variation in remnants masses in Figure 3 shows other uncertainties
in massive star evolution can lead to differences at least as large as
variations in mass-loss rates, which could also easily explain the formation
of a 20 M⊙ black hole in Cyg X-1 in the Galaxy.

None of the solutions that the BoOST, Geneva, MIST and PARSEC models employ
can currently be established as better than the others. In each case they have
been designed to address numerical issues in 1D stellar evolution. However,
the interplay of these solutions with mass-loss rates and convection further
adds to the uncertainties in massive stellar evolution. Therefore, a
systematic study to untangle the effect of the treatment of the Eddington-
limit from other physical assumptions will be conducted in an upcoming paper
(Agrawal et al. in prep) .

In the case of BPASS the stellar models evolve without requiring any numerical
enhancement. Whether this is a result of using a non-Lagrangian mesh (the
‘Eggletonian’ mesh, which is more adaptive to changes in stellar structure) or
if this is an artifact of bigger time steps (that helps stars skip problematic
short-lived phases of evolution) is currently not known. A separate study to
explore the effect of the ‘Lagrangian’ versus the ‘Eggletonian’ mesh structure
for massive stars (similar to Stancliffe et al. 2004 study for low and
intermediate mass stars) is highly desirable.

In conclusion, it is crucial to be aware of the uncertainties resulting from
numerical methods whenever the evolutionary model sequences of massive stars
are applied in any scientific project, such as gravitational-wave event rate
predictions or star-formation and feedback studies.

We only focus on massive stars as isolated single stars in this work. However,
there is mounting evidence that massive stars are formed as binaries or
triples, thus treating them as single stars might not be correct (e.g.,
Klencki et al., 2020; Laplace et al., 2021). Several studies have shown that
binarity can heavily influence the lives of massive stars through mass and
angular momentum transfer (de Mink et al., 2009; Marchant et al., 2016) and
can therefore help in avoiding density and pressure inversion regions in
stellar envelopes (Shenar et al., 2020).

We also limit our study to massive stars at near-solar metallicity, where due
to high opacity, the numerical instabilities related to the proximity to the
Eddington-limit are maximum. Since opacity decreases with metallicity, opacity
peaks become less prominent at lower metallicity. Nevertheless, stars with
low-metal content also reach the Eddington-limit, although at higher initial
masses (Sanyal et al., 2015, 2017). While progenitors of currently detectable
gravitational-wave (GW) sources may have been born in the early Universe where
the metal content is sub-solar (cf. Santoliquido et al., 2021), high star-
formation rates at near solar metallicities offer a fertile ground for the
formation of more GW sources, although less massive compared to sub-solar
metallicity (Neijssel et al., 2019). As such, there is good motivation for
studying the behaviour and reliability of massive star models across a wide
range of metallicities (Agrawal et al. in prep).

Collecting observational data as well as improvements in 3D and hydrodynamical
modelling will help us better constrain the models of massive stars in the
future. Until then, however, we urge the broader community to treat any set of
stellar models with caution. Ideally, one would implement all available
simulations as input into any given astrophysical study, and test the outcome
also in terms of stellar evolution related uncertainties. With tools such as
METISSE (Agrawal et al., 2020) and SEVN (Spera et al., 2019), this task is
becoming feasible.

## Data Availability

All the stellar models used in this work are publicly available:
PARSEC: people.sissa.it/$\sim$sbressan/parsec.html,
MIST: waps.cfa.harvard.edu/MIST,
Geneva: obswww.unige.ch/Research/evol/tables_grids2011,
BPASS: bpass.auckland.ac.nz/9.html,
BoOST: boost.asu.cas.cz.

## Acknowledgements

DSz has been supported by the Alexander von Humboldt Foundation. PA
acknowledges the support from the Australian Research Council Centre of
Excellence for Gravitational Wave Discovery (OzGrav). We also thank Stefanie
Walch-Gassner, Debashis Sanyal, Ross Church and Jan Eldridge for useful
comments and discussions.

## References

  * Agrawal et al. (2020) Agrawal P., Hurley J., Stevenson S., Szécsi D., Flynn C., 2020, MNRAS, 497, 4549
  * Aguilera-Dena et al. (2018) Aguilera-Dena D. R., Langer N., Moriya T. J., Schootemeijer A., 2018, ApJ, 858, 115
  * Alongi et al. (1993) Alongi M., Bertelli G., Bressan A., Chiosi C., Fagotto F., Greggio L., Nasi E., 1993, A&AS, 97, 851
  * Asplund et al. (2005) Asplund M., Grevesse N., Sauval A. J., 2005, in Barnes III T. G., Bash F. N., eds, Astronomical Society of the Pacific Conference Series Vol. 336, Cosmic Abundances as Records of Stellar Evolution and Nucleosynthesis. p. 25
  * Asplund et al. (2009) Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481
  * Asplund et al. (2021) Asplund M., Amarsi A. M., Grevesse N., 2021, arXiv e-prints, p. arXiv:2105.01661
  * Basu (2009) Basu S., 2009, in Dikpati M., Arentoft T., González Hernández I., Lindsey C., Hill F., eds, Astronomical Society of the Pacific Conference Series Vol. 416, Solar-Stellar Dynamos as Revealed by Helio- and Asteroseismology: GONG 2008/SOHO 21. p. 193
  * Belczynski et al. (2008) Belczynski K., Kalogera V., Rasio F. A., Taam R. E., Zezas A., Bulik T., Maccarone T. J., Ivanova N., 2008, ApJS, 174, 223
  * Belczynski et al. (2010) Belczynski K., Bulik T., Fryer C. L., Ruiter A., Valsecchi F., Vink J. S., Hurley J. R., 2010, ApJ, 714, 1217
  * Bestenlehner et al. (2014) Bestenlehner J. M., et al., 2014, A&A, 570, A38
  * Böhm-Vitense (1958) Böhm-Vitense E., 1958, Z. Astrophys., 46, 108
  * Bressan et al. (1981) Bressan A. G., Chiosi C., Bertelli G., 1981, A&A, 102, 25
  * Bressan et al. (2012) Bressan A., Marigo P., Girardi L., Salasnich B., Dal Cero C., Rubele S., Nanni A., 2012, MNRAS, 427, 127
  * Brott et al. (2011a) Brott I., et al., 2011a, A&A, 530, A115
  * Brott et al. (2011b) Brott I., et al., 2011b, A&A, 530, A116
  * Caffau et al. (2011) Caffau E., Ludwig H. G., Steffen M., Freytag B., Bonifacio P., 2011, Sol. Phys., 268, 255
  * Chen et al. (2015) Chen Y., Bressan A., Girardi L., Marigo P., Kong X., Lanza A., 2015, MNRAS, 452, 1068
  * Choi et al. (2016) Choi J., Dotter A., Conroy C., Cantiello M., Paxton B., Johnson B. D., 2016, ApJ, 823, 102
  * Costa et al. (2019) Costa G., Girardi L., Bressan A., Marigo P., Rodrigues T. S., Chen Y., Lanza A., Goudfrooij P., 2019, MNRAS, 485, 4641
  * Cox & Giuli (1968) Cox J., Giuli R., 1968, Principles of Stellar Structure: Physical principles. Principles of Stellar Structure, Gordon and Breach, https://books.google.com.au/books?id=TdhEAAAAIAAJ
  * Crowther (2000) Crowther P. A., 2000, A&A, 356, 191
  * Crowther (2012) Crowther P., 2012, Astronomy and Geophysics, 53, 4.30
  * Cunha et al. (2006) Cunha K., Hubeny I., Lanz T., 2006, ApJ, 647, L143
  * Eddington (1926) Eddington A. S., 1926, The Internal Constitution of the Stars. Cambridge Science Classics, Cambridge University Press, doi:10.1017/CBO9780511600005
  * Eggenberger et al. (2008) Eggenberger P., Meynet G., Maeder A., Hirschi R., Charbonnel C., Talon S., Ekström S., 2008, Ap&SS, 316, 43
  * Eggleton (1973) Eggleton P. P., 1973, MNRAS, 163, 279
  * Ekström et al. (2012) Ekström S., et al., 2012, A&A, 537, A146
  * Eldridge et al. (2017) Eldridge J. J., Stanway E. R., Xiao L., McClelland L. A. S., Taylor G., Ng M., Greis S. M. L., Bray J. C., 2017, Publications of the Astronomical Society of Australia, 34, e058
  * Evans et al. (2005) Evans C. J., et al., 2005, A&A, 437, 467
  * Fryer et al. (2012) Fryer C. L., Belczynski K., Wiktorowicz G., Dominik M., Kalogera V., Holz D. E., 2012, ApJ, 749, 91
  * Gatto et al. (2017) Gatto A., et al., 2017, MNRAS, 466, 1903
  * Gräfener et al. (2012) Gräfener G., Owocki S. P., Vink J. S., 2012, A&A, 538, A40
  * Grassitelli et al. (2016) Grassitelli L., Fossati L., Langer N., Simón-Díaz S., Castro N., Sanyal D., 2016, A&A, 593, A14
  * Grevesse & Noels (1993) Grevesse N., Noels A., 1993, in Prantzos N., Vangioni-Flam E., Casse M., eds, Origin and Evolution of the Elements. pp 15–25
  * Grevesse & Sauval (1998) Grevesse N., Sauval A. J., 1998, Space Sci. Rev., 85, 161
  * Grevesse et al. (1996) Grevesse N., Noels A., Sauval A. J., 1996, in Holt S. S., Sonneborn G., eds, Astronomical Society of the Pacific Conference Series Vol. 99, Cosmic Abundances. p. 117
  * Groh et al. (2020) Groh J. H., Farrell E. J., Meynet G., Smith N., Murphy L., Allan A. P., Georgy C., Ekstroem S., 2020, ApJ, 900, 98
  * Hamann et al. (1995) Hamann W. R., Koesterke L., Wessolowski U., 1995, A&A, 299, 151
  * Harpaz (1984) Harpaz A., 1984, MNRAS, 210, 633
  * Heger et al. (2000) Heger A., Langer N., Woosley S. E., 2000, ApJ, 528, 368
  * Heggie & Hut (2003) Heggie D., Hut P., 2003, The Gravitational Million-Body Problem: A Multidisciplinary Approach to Star Cluster Dynamics. Cambridge University Press
  * Henyey et al. (1965) Henyey L., Vardya M. S., Bodenheimer P., 1965, ApJ, 142, 841
  * Herwig (2000) Herwig F., 2000, A&A, 360, 952
  * Humphreys & Davidson (1979) Humphreys R. M., Davidson K., 1979, ApJ, 232, 409
  * Jiang et al. (2015) Jiang Y.-F., Cantiello M., Bildsten L., Quataert E., Blaes O., 2015, ApJ, 813, 74
  * Jones et al. (2015) Jones S., Hirschi R., Pignatari M., Heger A., Georgy C., Nishimura N., Fryer C., Herwig F., 2015, MNRAS, 447, 3115
  * Kaiser et al. (2020) Kaiser E. A., Hirschi R., Arnett W. D., Georgy C., Scott L. J. A., Cristini A., 2020, MNRAS, 496, 1967
  * Kippenhahn & Weigert (1990) Kippenhahn R., Weigert A., 1990, Stellar Structure and Evolution. Springer
  * Klencki et al. (2020) Klencki J., Nelemans G., Istrate A. G., Pols O., 2020, A&A, 638, A55
  * Lamers & Cassinelli (1999) Lamers H., Cassinelli J., 1999, Introduction to Stellar Winds. Cambridge University Press
  * Langer (1989) Langer N., 1989, A&A, 210, 93
  * Langer (1997) Langer N., 1997, in Nota A., Lamers H., eds, Astronomical Society of the Pacific Conference Series Vol. 120, Luminous Blue Variables: Massive Stars in Transition. p. 83
  * Laplace et al. (2021) Laplace E., Justham S., Renzo M., Götberg Y., Farmer R., Vartanyan D., de Mink S. E., 2021, arXiv e-prints, p. arXiv:2102.05036
  * Lau et al. (2012) Lau H. H. B., Gil-Pons P., Doherty C., Lattanzio J., 2012, A&A, 542, A1
  * Ledoux (1947) Ledoux P., 1947, ApJ, 105, 305
  * Levesque (2017) Levesque E. M., 2017, Astrophysics of Red Supergiants. 2514-3433, IOP Publishing, doi:10.1088/978-0-7503-1329-2
  * Maeder (1987) Maeder A., 1987, A&A, 173, 247
  * Maeder (2009) Maeder A., 2009, Physics, Formation and Evolution of Rotating Stars. Springer Berlin Heidelberg, doi:10.1007/978-3-540-76949-1
  * Maeder & Meynet (2005) Maeder A., Meynet G., 2005, A&A, 440, 1041
  * Mapelli et al. (2020) Mapelli M., Spera M., Montanari E., Limongi M., Chieffi A., Giacobbo N., Bressan A., Bouffanais Y., 2020, ApJ, 888, 76
  * Marchant et al. (2016) Marchant P., Langer N., Podsiadlowski P., Tauris T. M., Moriya T. J., 2016, A&A, 588, A50
  * Martins & Palacios (2013) Martins F., Palacios A., 2013, A&A, 560, A16
  * Miller-Jones et al. (2021) Miller-Jones J. C. A., et al., 2021, Science, 371, 1046
  * Moe & Di Stefano (2017) Moe M., Di Stefano R., 2017, ApJS, 230, 15
  * Neijssel et al. (2019) Neijssel C. J., et al., 2019, MNRAS, 490, 3740
  * Neijssel et al. (2021) Neijssel C. J., Vinciguerra S., Vigna-Gómez A., Hirai R., Miller-Jones J. C. A., Bahramian A., Maccarone T. J., Mandel I., 2021, ApJ, 908, 118
  * Nieuwenhuijzen & de Jager (1990) Nieuwenhuijzen H., de Jager C., 1990, A&A, 231, 134
  * Nishida & Schindler (1967) Nishida M., Schindler A. M., 1967, PASJ, 19, 606
  * Nugis & Lamers (2000) Nugis T., Lamers H. J. G. L. M., 2000, A&A, 360, 227
  * Paxton et al. (2011) Paxton B., Bildsten L., Dotter A., Herwig F., Lesaffre P., Timmes F., 2011, ApJS, 192, 3
  * Paxton et al. (2013) Paxton B., et al., 2013, ApJS, 208, 4
  * Petit et al. (2017) Petit V., et al., 2017, MNRAS, 466, 1052
  * Pols et al. (1998) Pols O. R., Schröder K.-P., Hurley J. R., Tout C. A., Eggleton P. P., 1998, MNRAS, 298, 525
  * Puls et al. (2000) Puls J., Springmann U., Lennon M., 2000, A&AS, 141, 23
  * Raithel et al. (2018) Raithel C. A., Sukhbold T., Özel F., 2018, ApJ, 856, 35
  * Ramírez-Agudelo et al. (2013) Ramírez-Agudelo O. H., et al., 2013, A&A, 560, A29
  * Renzo et al. (2017) Renzo M., Ott C. D., Shore S. N., de Mink S. E., 2017, A&A, 603, A118
  * Salpeter (1955) Salpeter E. E., 1955, ApJ, 121, 161
  * Sana et al. (2012) Sana H., et al., 2012, Science, 337, 444
  * Sander & Vink (2020) Sander A. A. C., Vink J. S., 2020, MNRAS, 499, 873
  * Santoliquido et al. (2021) Santoliquido F., Mapelli M., Giacobbo N., Bouffanais Y., Artale M. C., 2021, MNRAS, 502, 4877
  * Sanyal et al. (2015) Sanyal D., Grassitelli L., Langer N., Bestenlehner J. M., 2015, A&A, 580, A20
  * Sanyal et al. (2017) Sanyal D., Langer N., Szécsi D., -C Yoon S., Grassitelli L., 2017, A&A, 597, A71
  * Sarkisyan et al. (2020) Sarkisyan A., et al., 2020, MNRAS, 497, 687
  * Schneider et al. (2014) Schneider F. R. N., Langer N., de Koter A., Brott I., Izzard R. G., Lau H. H. B., 2014, A&A, 570, A66
  * Schootemeijer et al. (2019) Schootemeijer A., Langer N., Grin N. J., Wang C., 2019, A&A, 625, A132
  * Schwarzschild (1958) Schwarzschild M., 1958, Structure and evolution of the stars.. Princeton University Press
  * Shenar et al. (2020) Shenar T., Gilkis A., Vink J. S., Sana H., Sander A. A. C., 2020, A&A, 634, A79
  * Smith (2014) Smith N., 2014, ARA&A, 52, 487
  * Smith & Tombleson (2015) Smith N., Tombleson R., 2015, MNRAS, 447, 598
  * Spera et al. (2019) Spera M., Mapelli M., Giacobbo N., Trani A. A., Bressan A., Costa G., 2019, MNRAS, 485, 889
  * Stancliffe (2006) Stancliffe R. J., 2006, MNRAS, 370, 1817
  * Stancliffe et al. (2004) Stancliffe R. J., Tout C. A., Pols O. R., 2004, MNRAS, 352, 984
  * Stevenson et al. (2019) Stevenson S., Sampson M., Powell J., Vigna-Gómez A., Neijssel C. J., Szécsi D., Mandel I., 2019, ApJ, 882, 121
  * Sukhbold & Woosley (2014) Sukhbold T., Woosley S., 2014, Astrophys. J., 783, 10
  * Szécsi (2016) Szécsi D., 2016, PhD thesis, Mathematisch-Naturwissenschaftlichen Fakultät der Universität Bonn, doi:10.5281/zenodo.998070
  * Szecsi (2017) Szecsi D., 2017, in XII Multifrequency Behaviour of High Energy Cosmic Sources Workshop (MULTIF2017). p. 65 (arXiv:1710.05655)
  * Szécsi et al. (2020) Szécsi D., Wünsch R., Agrawal P., Langer N., 2020, arXiv e-prints, p. arXiv:2004.08203
  * Takahashi & Langer (2021) Takahashi K., Langer N., 2021, A&A, 646, A19
  * Topping & Shull (2015) Topping M. W., Shull J. M., 2015, ApJ, 800, 97
  * Vigna-Gómez et al. (2018) Vigna-Gómez A., et al., 2018, MNRAS, 481, 4009
  * Vink (2011) Vink J. S., 2011, Ap&SS, 336, 163
  * Vink et al. (2000) Vink J. S., de Koter A., Lamers H. J. G. L. M., 2000, A&A, 362, 295
  * Vink et al. (2001) Vink J. S., de Koter A., Lamers H. J. G. L. M., 2001, A&A, 369, 574
  * Walder et al. (2012) Walder R., Folini D., Meynet G., 2012, Space Sci. Rev., 166, 145
  * Weinberger et al. (2020) Weinberger R., Springel V., Pakmor R., 2020, ApJS, 248, 32
  * Wofford et al. (2016) Wofford A., et al., 2016, MNRAS, 457, 4296
  * Woosley et al. (2002) Woosley S. E., Heger A., Weaver T. A., 2002, Reviews of Modern Physics, 74, 1015
  * Yoon et al. (2006) Yoon S.-C., Langer N., Norman C., 2006, A&A, 460, 199
  * de Jager et al. (1988) de Jager C., Nieuwenhuijzen H., van der Hucht K. A., 1988, A&AS, 72, 259
  * de Mink et al. (2009) de Mink S. E., Cantiello M., Langer N., Pols O. R., Brott I., Yoon S.-C., 2009, A&A, 497, 243


# Artificial Intuition: Efficient Classification of Scientific Abstracts

Harsh Sakhrani, Naseela Pervez, Anirudh Ravi Kumar, Fred Morstatter
Information Sciences Institute, Viterbi School of Engineering, University of
Southern California
Alexandra Graddy Reed
Sol Price School of Public Policy, University of Southern California
Andrea Belz
Information Sciences Institute, Viterbi School of Engineering, University of
Southern California

###### Abstract

It is desirable to coarsely classify short scientific texts, such as grant or
publication abstracts, for strategic insight or research portfolio management.
These texts efficiently transmit dense information to experts possessing a
rich body of knowledge to aid interpretation. Yet this task is remarkably
difficult to automate because of brevity and the absence of context. To
address this gap, we have developed a novel approach to generate and
appropriately assign coarse domain-specific labels. We show that a Large
Language Model (LLM) can provide metadata essential to the task, in a process
akin to the augmentation of supplemental knowledge representing human
intuition, and propose a workflow. As a pilot study, we use a corpus of award
abstracts from the National Aeronautics and Space Administration (NASA). We
develop new assessment tools in concert with established performance metrics.

##  1 Introduction

Analyzing technical documents is a crucial strategic task, enabling the
management of research portfolios, tracking investment trends, and exploring
scientific advancement. On a more tactical level, it can aid the preliminary
screening of scientific abstracts in systematic reviews (Buchlak et al., 2020;
Rios and Kavuluru, 2015; Ambalavanan and Devarakonda, 2020).

Several approaches are possible. First, authors can label their own work, but
this presents several challenges: (1) authors that self-label their own texts
may make idiosyncratic decisions, (2) authors in close disciplines may use
different terms for related concepts, such as “robotics” and “autonomy”, and
(3) multidisciplinary projects may require novel or multiple labels.

A second method is to impose an external ontology. However, these schemes
often have both fine- and coarse-granularities (e.g., “networks” versus “ad-
hoc networks”). Another concern is that the scheme simply lacks appropriate
labels, especially for emerging fields.

Automated processes do exist. Those with a large number of parameters are now
customizable at lower computational cost (Hu et al., 2021; Ben Zaken et al.,
2022). Although dedicated pre-trained models can yield robust results, they
incur significant expenses in manual annotation due to reliance on supervised
learning (Beltagy et al., 2019; Chang et al., 2008; Cohan et al., 2020).

In summary, we face two distinct needs in the analysis of scientific
documents: (1) a unified, coarse-grained, non-overlapping taxonomy, tailored
to uniquely classify a set of documents; and (2) an unsupervised methodology
that circumvents the reliance on manual annotation while effectively managing
the peculiarities of scientific text. These challenges are exacerbated for
abstracts.

In manual labeling, an expert’s rapid progress often hinges on integrating
prior knowledge, crucial for effective comprehension Reid Smith and Hammond
(2021). In so doing, the expert rapidly identifies the phrases conveying the
most information and uses those for classification. Importantly, this process
is not a simple frequency or statistical analysis; indeed, the most important
phrase may appear only once. Moreover, multigrams carrying high semantic value
may not appear systematically in the same place in a sentence or paragraph.

Here we describe “artificial intuition,” a method mimicking the expert’s
process to execute two objectives: generating an optimal label space and
producing accurate predictions within this new space. We integrate tools into
a novel workflow to identify important terms, augment them with relevant
background information, then aggregate these enhanced documents into clusters
for classification purposes.

As a pilot case to evaluate our methodology, we analyze award abstracts of
federally funded projects from the National Aeronautics and Space
Administration (NASA) Small Business Innovation Research (SBIR) Program. We
obtain domain knowledge by extracting and ranking the abstract’s keywords /
keyphrases (which we will collectively refer to as “keywords”). We generate
metadata for these keywords in a zero-shot setting and derive embeddings for
the keyword-metadata concatenations using a pre-trained Sentence Transformer.

For label space generation, we implement a clustering process that represents
the task of organizing awards into funding themes. This method not only
clarifies the thematic organization of the documents but also reveals the
hierarchical relationships between different topics. We introduce a novel
evaluation scheme to assess whether the label set comprehensively spans the
document space and can serve as a set of basis vectors.

To predict labels, we reinterpret the multilabel-classification problem as a
semantic matching challenge wherein the document space is characterized by the
keyword-metadata concatenation and the label space is described by the element
closest to the centroid for each cluster. This retrieval-based perspective
allows for flexibility in adapting to new label spaces without the need for
retraining.

This framework accommodates various levels of parsimony, which we explore
extensively in our experiments. Finally, using our test sample, we demonstrate
the efficacy of our prediction methodology and quantify the performance.

##  2 Related Work

Various methodologies have been proposed for text classification. Bayesian
approaches (Tang et al., 2016) classify the text by extracting features. One
method is to first select document features with discriminative power, then
compute the semantic similarity between features and documents (Zong et al.,
2015), but this becomes more difficult as the number of features grows.
Support Vector Machines (SVMs) can be used for document classification (Cai
and Hofmann, 2004). However, these approaches are constrained by the
requirement for manual feature engineering, limiting their ability to capture
the complexity of natural language.

New deep learning techniques have advanced scientific document classification.
Neural network-based architectures (Lee and Dernoncourt, 2016), particularly
Convolutional Neural Networks (CNNs) (Sun et al., 2019) and Recurrent Neural
Networks (RNNs) (Xun et al., 2019; Liu et al., 2016), outperform some
traditional machine learning methods. These models automatically learn feature
representations from data, capturing both the semantic and syntactic nuances
of text.

These methods presume that documents are related to only one label. Newer
approaches (e.g., Liu et al. (2017); Song et al. (2022); Xiao et al. (2019);
Blanco et al. (2019); Chang et al. (2020)) classify documents with multiple
labels, and one alternative attempts to map 10,000 fine-grained labels for
scientific documents (Zhang et al., 2022a) although most methods consider
10-50 coarse labels. These models are incompletely validated because many
real-world datasets will have limited or poorly labeled data.

Weakly supervised learning and zero-shot learning (ZSL) models do not use
annotated data. Some pre-trained language models demonstrate impressive
performance in zero-shot document classification (Devlin et al., 2019; Beltagy
et al., 2019; Liu et al., 2019) and can be used to assign multiple labels to a
given document (Yin et al., 2019). On the other hand, hierarchical multi-class
methods can use just class names - without training examples - as supervision
(Shen et al., 2021; Zhang et al., 2022b). Large language models trained on
scientific data, such as Galactica (Taylor et al., 2022) and SciNCL
(Ostendorff et al., 2022), can be used to assign labels to a scientific
document. Many approaches use metadata, such as generic descriptions, as
supervision for further classification (Zhang et al., 2023). However, these
methods are still potentially subject to noise. Here we describe a method to
identify keywords and derive context-specific metadata to improve
classification accuracy, particularly for short abstracts.

##  3 Approach

###  3.1 Problem Formulation

The scientific literature tagging task can be conceptualized as a multi-label
classification (each paper can be relevant to more than one label) problem,
where all candidate tags (e.g., “Aerodynamics,” “Superconductance/Magnetics”)
constitute a label space $Y$ of arbitrary size. We seek to:

  * •

Construct a new label space $Y$ comprising coarse-grained labels and
aggregating correlated labels (e.g., merging “Optics” and “Photonics” into
“Optical technologies”).

  * •

Develop an unsupervised multi-label classifier that can effectively map an
abstract to the new label space $Y$.

A simplistic approach would utilize a pre-trained language model to encode
each document and label, generate their embeddings, and then conduct a nearest
neighbor search in the embedding space. However, this method encounters two
primary challenges: (1) the existing language models are largely trained on
general English text that does not discern technical terms, and (2) analogous
labels (e.g., “networking” and “ad-hoc networks”) confound the results. One
might augment the label embedding process with generic metadata, such as a
brief description from Wikipedia or using solutions like Positive Instance
Feature Aggregation (PIFA) Yu et al. (2022).

Instead, we seek to generate a context-specific glossary. This has the added
advantage that the labels can be fine-tuned, converting a multi-label problem
into a simpler system. For instance, a thermal protection system (TPS)
consists of materials suited to handle extremely high temperatures. In a
conventional classification scheme, this might require two labels, such as
“materials” and “temperature.” In contrast, we create a system such that
“thermal protection system” is itself sufficient to serve as the only label.
This is possible only with a label space customized to the knowledge domain.

###  3.2 Implementation Components

  * •

Yet Another Keyword Extractor (YAKE) (Campos et al., 2020) is a lightweight,
unsupervised keyword extraction algorithm that uses statistical properties and
contextual information.

  * •

Mistral 7B is a Large Language Model (LLM) with strong performance of Llama-2
13B on key benchmarks Jiang et al. (2023).

  * •

Maximal Marginal Relevance (MMR) (Carbonell and Goldstein, 1998) iteratively
selects candidate items that simultaneously maximize their relevance to the
query and their novelty compared to previously selected items.

Figure 1: Label Space Generation flowchart. The clusters are named with the
keyword closest to the cluster centroid.

  * •

Sentence Transformer (S-Transformer)111https://huggingface.co/sentence-
transformers/all-mpnet-base-v2 Reimers and Gurevych (2019) constructs dense
vector representations of sentences to enable efficient comparison of text
semantics.

###  3.3 Document Corpus

The NASA SBIR program publishes abstracts of funded projects. We used 1,230
abstracts from 2010 to 2015 extracted online from the publicly available
archive222sbir.gov. The average abstract length is about 450 words. All
abstracts were pre-processed by removing stop words, which were variations and
combinations of: “NASA”, “space”, “mission(s)”, “research”, “SBIR”,
“spacecraft”, “future”, and “science”. These words and multigrams comprising
these words appeared in a large number of the abstracts, and therefore they
provided little information to assist in classification. We randomly drew 100
abstracts (roughly 10%) for manual classification, described below.

###  3.4 Label Space Generation

We generate the label space as illustrated in Figure 1. Initially, pre-
processed abstracts are submitted to YAKE to extract keywords. One
hyperparameter of our workflow is the number of keywords $\hat{c}$. We
estimated that $\hat{c}$ should be approximately 5 as it represents 1-2% of
the abstract length. We confirmed that the F1 results, described in more
detail below, showed a general lack of sensitivity to this parameter (Figure
2), and therefore we set $\hat{c}=5$ in our main analyses.

Figure 2: Variation of F1 score with the number of keywords at the threshold
of top 1%.

We sought to supplement these keywords with contextual definitions to form
metadata. We used Mistral-7B Instruct v0.2 with hyperparameters set at default
values and submitted the following prompt:

> Given the scientific abstract and the keywords that have been extracted for
> the document, provide a concise meta data/prior information for every
> keyword in context of the document. Incorporate any extra knowledge that can
> help classify the document to relevant topics.

This combined data-keyword concatenation is processed using the S-Transformer
model to produce embeddings. A critical aspect of this process is that the
metadata generated for each keyword is tailored specifically to the context of
the related document, ensuring that the embeddings are context-specific rather
than generic. We use k-means clustering (Habibi and Popescu-Belis, 2015) to
partition these embeddings into clusters, represented by the keyword closest
to their centroids and effectively summarizing each cluster’s thematic focus.
This method approximates the scheme by which such abstracts would be sorted in
a funding portfolio.

Unlike $\hat{c}$, the number of clusters $\hat{k}$ requires closer
examination. We seek a parsimonious model that minimizes the number of labels
per document. In practice, we seek to organize approximately 1,000 documents
into approximately 10-20 classes. In addition to making this a tractable
problem, it adequately represents the portfolio management process.

###  3.5 Annotation Task Design

We conducted a manual annotation task to label the test set of the NASA SBIR
abstracts. We presented the annotator, a NASA expert, with a scientific
abstract and the generated label set. The annotator was instructed to assign a
label to the scientific abstract only if one of the presented labels was
appropriate, and to leave it unlabeled otherwise. The same documents were
labeled for each configuration for consistency.

##  4 Results

###  4.1 Label Space Orthogonality: Redundancy

Our first task is estimate the degree of overlap within the label space. To do
so, we define the redundancy, $\mathcal{R}$, as a measure of the orthogonality
between labels. This figure of merit (FOM) is intrinsic to the label space and
assessed independently of individual document projections.

The labels are transformed into normalized embeddings using the S-Transformer
model, resulting in a label matrix $\mathcal{L}$ of dimensions $\hat{k}\times
v$ (in our case, $v=768$). Each element $\mathcal{L}_{ij}$ represents the
$j$-th dimension of the $i$-th label embedding.

To measure the orthogonality, we calculate the cosine similarity between each
pair of distinct label embeddings. If the labels are orthogonal and distinct,
the cosine similarity should approach 0; on the other hand, two labels
capturing closely related ideas will give a cosine similarity that approaches
1. Formally, for normalized label vectors $\mathcal{T}_{i}$ and
$\mathcal{T}_{j}$ in $\mathcal{L}$, we define redundancy $\mathcal{R}$ as the
maximum cosine similarity among all pairs:

$\mathcal{R}=\max_{i\neq j}(\text{cosine
similarity}(\mathcal{T}_{i},\mathcal{T}_{j}))$ (1)

where

$\text{cosine
similarity}(\mathcal{T}_{i},\mathcal{T}_{j})=\frac{\mathcal{T}_{i}\cdot\mathcal{T}_{j}}{\|\mathcal{T}_{i}\|\|\mathcal{T}_{j}\|}$

A value of $\mathcal{R}$ close to 0 is desirable because orthogonal label
embeddings suggest that each label contributes unique information without
redundancy. Conversely, a value of $\mathcal{R}$ approaching 1 shows that at
least one pair of labels shares a high degree of overlap. Overlap implies that
multiple labels may be describing similar features within the documents, thus
complicating the interpretability and utility of the label space. Our goal is
to represent each key concept with a unique label.

To understand the redundancy in our basis vector set, we executed the label
space generation process but systematically varied $\hat{k}$. We then
evaluated $\mathcal{R}$ for each label space. $\mathcal{R}$ increased with
$\hat{k}$ (Figure 3), as expected. Notably, we identified three general
regimes: At low $\hat{k}$, $\mathcal{R}$ was fairly flat and low. The labels
do not overlap. At approximately $\hat{k}=8$, the redundancy increased to a
new plateau. At much higher values of $\hat{k}$ (18 and higher), this FOM
entered a regime in which the value dramatically oscillated.

We therefore conclude that at very low cluster numbers ($\hat{k}<6$), the
severely reduced $\mathcal{R}$ indicates that the labels are probably
insufficient to describe the document set. At higher values of $\hat{k}>18$,
the risk of overlapping labels increases substantially, but the likelihood
that each document is at least minimally described also increases.

Figure 3: Variation of redundancy $\mathcal{R}$ with the number of clusters
$\hat{k}$.

###  4.2 Spanning the Document Space: Coverage

We defined the redundancy $\mathcal{R}$ to characterize the orthogonality of
our proposed label space basis vectors. Next, we study how comprehensively
these labels describe the documents, essentially determining if our labels can
span the document space.

Figure 4: Analysis workflow and use of the coverage matrix $\mathcal{W}$. In
one application (final step in green), the element with the maximum value is
used to generate the Coverage. The second usage (blue final step) is to
extract those values exceeding a specific threshold $T$ for the label
prediction task.

We architected a second workflow (Figure 4). Again we begin with YAKE usage
for a single document. We submit these keywords to Mistral-7B for document-
specific contextual definitions as supplementary metadata. Both the document
itself and the keyword-metadata concatenations are subsequently processed
through the S-Transformer model to generate their individual embeddings,
refined using MMR. This forms a new keyword embedding matrix $\mathcal{C}$ of
dimensions $v\times\hat{c}$, where $v$ (768 in our case) represents the
embedding dimension, and the extracted keywords are still parameterized by
$\hat{c}$.

Likewise, we still have the label embedding $\mathcal{L}$ of dimensions
$\hat{k}\times v$. As our goal is to understand the overlap between the labels
and the corpus embeddings, we define a new matrix, termed “coverage”,
$\mathcal{W}$ with elements $w_{ij}$:

$w_{ij}=\sum_{v}{L_{iv}C_{vj}},$ (2)

The coverage matrix $\mathcal{W}$ has the resulting dimension $\hat{k}$ x
$\hat{c}$, where $\hat{k}$ represents the number of labels and $\hat{c}$
represents the number of keywords. In other words, $\mathcal{W}$ is the
projection of the keywords onto the label space. (Strictly speaking, the
S-Transformer embeddings of length $v$ can be understood as creating a
coordinate system to facilitate projections.) Each $w_{ij}$ element ranges
from -1 to 1.

A high value of any element $w_{ij}$ indicates that a label and keyword are
highly aligned. Therefore, finding the maximum value that appears in this
matrix $\mathcal{W}$ will signify how well the label space describes the
keywords of an individual document in the best case. Consequently, we define
the coverage $\mathcal{S}$ for a given document $d$ (where $d$ is a member of
the document corpus $\mathcal{D}$):

$\mathcal{S}^{d}=\textrm{max}(w^{d}_{ij})$ (3)

The coverage for the corpus $D$ is simply the average of the individual
documents’ coverage:

$\mathcal{S}^{D}=\frac{\sum_{D}S^{d}}{D}$ (4)

Figure 5: Variation of coverage $S$ with $\hat{k}$.

This proposed figure of metric, coverage, provides critical validation that
the new label space is pertinent to the knowledge domain encompassed by the
documents. One would expect for coverage to be small if the label space is not
large enough - namely, for small values of $\hat{k}$. An intermediate regime
would appear if each new label adds significant new information. Eventually, a
final regime would be reached wherein the new information provided by an
additional label is marginal as the segmentation becomes finer, such as
comparing ‘Chemical Propulsion Technologies’ and ‘Electronic Propulsion
Technologies’. In other words, a general analytical form for coverage should
start near 0, then experience rapid growth until the space is largely covered
and it tapers off. The corpus coverage is bounded by 1 because the individual
documents’ coverage is given by a cosine similarity of two normalized vectors,
thus limited to 1.

We tested this concept by varying the number of clusters $\hat{k}$ from 2 to
28 and evaluating coverage $\mathcal{S}$ for each newly developed label space.
As $\hat{k}$ increased, the labels did indeed relate well to the documents, as
represented by keywords (Figure 5). In addition, the variation revealed a
generally asymptotic form, as expected.

Table 1: Labels at $\hat{k}=15$

Advanced Optical Systems
---
Advanced Photovoltaic Systems
Aeroservoelastic Analysis and Aircraft Systems Analysis
Aeroservoelastic Analysis Tools
Electric Propulsion Systems
Electrolyzers
High Energy Density Electronics
LIDAR Remote Sensing
Multifunctional Composite Materials
Optical Communications Technology
Radiation-Hardened Electronics
Robotic Science Missions
Technologies Fault Management
Thermal Protection Systems
Unmanned Aircraft Systems

###  4.3 Label Assignment: Precision and Recall

Table 2: Labels at $k=25$

Advanced Aeroservoelastic Analysis and Rotorcraft Aeromechanics
---
Advanced Composite and Ceramic Matrix Materials
Advanced ESR Technologies for Space Exploration
Advanced Energy Storage and Power Systems
Advanced Fluid and Thermal Management Technologies
Advanced Laser and Optical Communication Technologies
Advanced Manufacturing Technologies for Aerospace
Advanced Microwave and Remote Sensing Technologies
Advanced Optical Systems for Scientific Missions
Advanced Structural Sensors and NDE Technologies
Advanced Thermal Protection Systems
Airborne Measurement and Sensing Systems
Automation and Control in Robotic Science Missions
Fault Management Technologies
High-End Computing and Data Handling
Highly Capable Propulsion Systems
Innovative Aerospace Structural Design
Innovative Fiber-Optic and Navigational Technologies
International Space Station
LIDAR Remote Sensing Technologies
Mars Sample Return Missions
Radiation-Hardened Electronics and Sensors
Regenerative Life Support Systems
Solar Power Technologies for Advanced Energy Solutions
Unmanned Aircraft Systems Operations

We seek to create a label space with high coverage, indicating relevance; and
low redundancy or overlap. However, these measures act in opposition as higher
coverage naturally can lead to greater redundancy. That is, these two measures
form a trade space in which we strive to optimize $\hat{k}$.

We revisited the workflow generating the coverage matrix (Figure 4) and
developed a prediction pipeline, mirroring the initial process through the
creation of the coverage matrix $\mathcal{W}$.

In the coverage study, we took the maximum $w_{ij}$ value to characterize the
space. Here, we seek to find all relevant values of $w_{ij}$. To
operationalize this, we analyze the distribution of all $w_{ij}$ values and
establish a threshold $T$, which defines the minimum percentile to be used as
a filter for the $w_{ij}$ values, effectively distinguishing between
significant and negligible overlaps. For instance, setting $T=1\%$ means we
retain only the top 1% of the $w_{ij}$ values, which is more restrictive than
setting $T=10\%$. In practical terms, for a system of 5 keywords and 15
labels, a 1% threshold would retain just one label (top 1% of 5x15 = 75 matrix
elements results in one). On the other hand, a 10% threshold retains seven
elements that could be distributed in various ways. For instance, all five
keywords might describe label 1, with two of those keywords linked to label 2;
or only one keyword could be associated with each of seven labels. As the
threshold $T$ gets larger, the variability in possible outcomes increases.

For each document, we select labels associated with the values of $w_{ij}$
that exceed the threshold $T$. However, to accurately evaluate the
classification, a set of ‘true’ labels is required. While NASA maintains its
own taxonomy of approximately 200 labels that could theoretically serve this
purpose, the inconsistency in this taxonomy year-to-year and the excessive
number of labels compared to our needs complicate its use. Instead, as noted
in Section 3.5, we manually aligned the abstracts with our new labels.

Using the three regimes of Figure 3 as a guide, we considered three values for
$\hat{k}$ \- 4, 15, and 25 - and estimated the usual classification measures
of precision, recall, and F1. Moreover, we varied the threshold $T$,
hypothesizing that at low restrictive values of $T$, these measures should
improve as only the most significant overlaps in the coverage matrix would be
retained.

At $\hat{k}$ = 4, the labels were: Propulsion Technologies, Remote Sensing
Technologies, Thermal Protection Systems, and Unmanned Aircraft Systems.
However, the manual classification task failed because the labels simply did
not describe the abstracts.

At $\hat{k}=15$, the labels consisted of words generally associated with space
technologies (Table 1). Similarly, the $\hat{k}=25$ generated labels related
to space (Table 2); however, in this case the word “Advanced” preceded nearly
half the technical topics, suggesting that the semantic content of that word
decreases in this context. (Notably, the word “advanced” has been linked to
other technical contexts where its semantic content is diluted (Belz et al.,
2023)).

To evaluate our method’s quality, we set aside $\hat{k}=4$ and considered
differences between $\hat{k}=15$ and $\hat{k}=25$. We focused on the F1 score
and found that $\hat{k}=15$ consistently yielded higher scores than the
overdetermined space represented by $\hat{k}=25$ (Figure 6). As a result, we
concluded that $\hat{k}=15$ represented a better set of labels to describe
this space.

Figure 6: Variation of F1 scores for assigned labels with weights $w$
exceeding the percentile threshold $T$, as defined in the text.

Our final task was to demonstrate the advantage of augmenting the abstract
with the metadata extracted from the additional analysis of the keywords.
Using the $\hat{k}=15$ label space described above, we evaluated the
performance of our model with and without the metadata generated by the LLM.
We found that the LLM consistently improved the F1 score (Figure 7) for all
tested values of the threshold $T$. This was due to improvement primarily in
the precision (Table 3).

Figure 7: Variation of F1 scores for assigned labels with weights $w$ exceeding the percentile threshold $T$, as defined in the text, for $\hat{k}=15$. Table 3: Precision, recall and F1 scores for $k=15$ at varying thresholds ($T$). $T$ (%) | Precision | Recall | F1
---|---|---|---
1 | 0.56 | 0.56 | 0.56
5 | 0.35 | 0.79 | 0.49
10 | 0.22 | 0.81 | 0.35
15 | 0.16 | 0.86 | 0.27
20 | 0.15 | 0.91 | 0.25

##  5 Discussion and Future Research

Scientific communication is designed to efficiently carry rich information
between experts. The abstract of a grant or publication is perhaps the most
striking example, wherein sophisticated concepts are conveyed in a relatively
dense, short vehicle. Years of study generate a large body of knowledge to
guide the expert in a classification process. Indeed, this additional material
and the associated judgments underpin the rapid decision-making characteristic
of human intuition.

We have sought to replicate that process in an automated methodology. Our
unsupervised approach is robust and flexible, enabling its use in various
domains. Its independence from specific label sets underscores its
adaptability and broad applicability. Our contributions range from applied
text processing tasks to economics and public policy, with several interesting
directions ahead.

First, we have tested this approach on a relatively narrow set of abstracts by
selecting a NASA corpus of documents as the first test case. This exercise
should be conducted on benchmark datasets such as Maple (Metadata-Aware Paper
colLEction) 333https://github.com/yuzhimanhua/MAPLE/tree/master/. This would
demonstrate the generalizability of our approach.

Second, a different validation would be to compare these results with those of
longer documents. For instance, one could analyze both publication abstracts
and the full text. It is not clear if the publications would contain more
noise; or perhaps the complete text would carry the metadata such that the LLM
task would be less necessary.

In addition, here the manual classification exercise assigned only one label
to each abstract as a rigorous test. We have not explored the opportunity to
generate multiple labels for a single abstract. Indeed, the $\hat{k}=25$ data
set points to this, as some of the labels (such as “Advanced Thermal
Protection Systems”) addressed the technology itself, while others described
the intended application (e.g., “Mars Sample Return missions”). In the future,
we can develop a new weighting scheme addressing this complex classification.

Finally, our method opens lines of inquiry in business or public policy, as we
could use this labeling method to generate metadata for the abstracts
themselves. In this fashion, the labels could form a variable to be used in
further assessment, such as patterns in funding, research direction, patents,
or other corpora where scientific documents are condensed in short summaries.
Using this method with public company reports could create entirely new
industry categories, updating existing schemes (Shweta and Belz, 2021; Hoberg
and Phillips, 2010, 2016). Moreover, these data could be combined with other
tags, such as principal investigator, institution, or other bibliometric
characteristics to create a complex profile. Such a data set could be used to
track a number of interesting trends.

##  6 Conclusion

For labeling short scientific documents, such as abstracts, pre-existing
domain-specific taxonomies are ambiguous. Defining a label space spanning the
set of documents is an important task that humans execute easily. In this
paper, we demonstrate that the text of the documents is insufficient to either
define the label space or predict the labels. We present evidence that an LLM
can provide critical metadata to address this gap, forming the basis for
artificial intuition. Additionally, we propose both an architecture to address
this and two novel measures to evaluate the constructed label spaces. Testing
our model with a corpus of NASA award abstracts, we demonstrate a workflow
that integrates the LLM’s supplemental data successfully.

## References

  * Ambalavanan and Devarakonda (2020) Ashwin Karthik Ambalavanan and Murthy V. Devarakonda. 2020.  Using the contextual language model bert for multi-criteria classification of scientific articles.  _Journal of Biomedical Informatics_ , 112:103578.
  * Banerjee et al. (2020) Soumya Banerjee, Debarshi Kumar Sanyal, Samiran Chattopadhyay, Plaban Kumar Bhowmick, and Partha Pratim Das. 2020.  Segmenting scientific abstracts into discourse categories: a deep learning-based approach for sparse labeled data.  In _Proceedings of the ACM/IEEE Joint Conference on Digital Libraries in 2020_ , pages 429–432.
  * Beltagy et al. (2019) Iz Beltagy, Kyle Lo, and Arman Cohan. 2019.  SciBERT: A pretrained language model for scientific text.  In _Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP)_ , pages 3615–3620, Hong Kong, China. Association for Computational Linguistics.
  * Belz et al. (2023) Andrea Belz, Alexandra Graddy-Reed, FNU Shweta, Aleksandar Giga, and Shivesh Meenakshi Murali. 2023.  Deterministic bibliometric disambiguation challenges in company names.  In _2023 IEEE 17th International Conference on Semantic Computing (ICSC)_. IEEE.
  * Ben Zaken et al. (2022) Elad Ben Zaken, Yoav Goldberg, and Shauli Ravfogel. 2022.  BitFit: Simple parameter-efficient fine-tuning for transformer-based masked language-models.  In _Proceedings of the 60th Annual Meeting of the Association for Computational Linguistics (Volume 2: Short Papers)_ , pages 1–9, Dublin, Ireland. Association for Computational Linguistics.
  * Blanco et al. (2019) Alberto Blanco, Arantza Casillas, Alicia Pérez, and Arantza Diaz de Ilarraza. 2019.  Multi-label clinical document classification: Impact of label-density.  _Expert Systems with Applications_ , 138:112835.
  * Buchlak et al. (2020) Buchlak, Quinlan, and Leveque JC. 2020.  Machine learning applications to clinical decision support in neurosurgery.  _an artificial intelligence augmented systematic review._
  * Cai and Hofmann (2004) Lijuan Cai and Thomas Hofmann. 2004.  Hierarchical document categorization with support vector machines.  In _Proceedings of the Thirteenth ACM International Conference on Information and Knowledge Management_ , CIKM ’04, page 78–87, New York, NY, USA. Association for Computing Machinery.
  * Campos et al. (2020) Ricardo Campos, Vítor Mangaravite, Arian Pasquali, Alípio Jorge, Célia Nunes, and Adam Jatowt. 2020.  Yake! keyword extraction from single documents using multiple local features.  _Information Sciences_ , 509:257–289.
  * Carbonell and Goldstein (1998) Jaime Carbonell and Jade Goldstein. 1998.  The use of mmr, diversity-based reranking for reordering documents and producing summaries.  In _Proceedings of the 21st Annual International ACM SIGIR Conference on Research and Development in Information Retrieval_ , SIGIR ’98, page 335–336, New York, NY, USA. Association for Computing Machinery.
  * Chang et al. (2008) Ming-Wei Chang, Lev-Arie Ratinov, Dan Roth, and Vivek Srikumar. 2008.  Importance of semantic representation: Dataless classification.  In _Aaai_ , volume 2, pages 830–835.
  * Chang et al. (2020) Wei-Cheng Chang, Hsiang-Fu Yu, Kai Zhong, Yiming Yang, and Inderjit Dhillon. 2020.  Taming pretrained transformers for extreme multi-label text classification.
  * Cohan et al. (2020) Arman Cohan, Sergey Feldman, Iz Beltagy, Doug Downey, and Daniel Weld. 2020.  SPECTER: Document-level representation learning using citation-informed transformers.  In _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics_ , pages 2270–2282, Online. Association for Computational Linguistics.
  * Devlin et al. (2019) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2019.  Bert: Pre-training of deep bidirectional transformers for language understanding.
  * Fortunato et al. (2018) Santo Fortunato, Carl T. Bergstrom, Katy Börner, James A. Evans, Dirk Helbing, Staša Milojević, Alexander M. Petersen, Filippo Radicchi, Roberta Sinatra, Brian Uzzi, Alessandro Vespignani, Ludo Waltman, Dashun Wang, and Albert-László Barabási. 2018.  Science of science.  _Science_ , 359(6379):eaao0185.
  * Habibi and Popescu-Belis (2015) Maryam Habibi and Andrei Popescu-Belis. 2015.  Keyword extraction and clustering for document recommendation in conversations.  _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , 23(4):746–759.
  * Hoberg and Phillips (2010) Gerard Hoberg and Gordon Phillips. 2010.  Product market synergies and competition in mergers and acquisitions: A text-based analysis.  _Review of Financial Studies_ , 23(10):3773–3811.
  * Hoberg and Phillips (2016) Gerard Hoberg and Gordon Phillips. 2016.  Text-based network industries and endogenous product differentiation.  _Journal of Political Economy_ , 124(5):1423–1465.
  * Hu et al. (2021) Edward J Hu, Yelong Shen, Phillip Wallis, Zeyuan Allen-Zhu, Yuanzhi Li, Shean Wang, Lu Wang, and Weizhu Chen. 2021.  LoRA: Low-Rank adaptation of large language models.
  * Jiang et al. (2023) Albert Q. Jiang, Alexandre Sablayrolles, Arthur Mensch, Chris Bamford, Devendra Singh Chaplot, Diego de las Casas, Florian Bressand, Gianna Lengyel, Guillaume Lample, Lucile Saulnier, Lélio Renard Lavaud, Marie-Anne Lachaux, Pierre Stock, Teven Le Scao, Thibaut Lavril, Thomas Wang, Timothée Lacroix, and William El Sayed. 2023.  Mistral 7b.
  * Lee and Dernoncourt (2016) Ji Young Lee and Franck Dernoncourt. 2016.  Sequential short-text classification with recurrent and convolutional neural networks.  In _Proceedings of the 2016 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies_ , pages 515–520, San Diego, California. Association for Computational Linguistics.
  * Liu et al. (2017) Jingzhou Liu, Wei-Cheng Chang, Yuexin Wu, and Yiming Yang. 2017.  Deep learning for extreme multi-label text classification.  In _Proceedings of the 40th International ACM SIGIR Conference on Research and Development in Information Retrieval_ , SIGIR ’17, page 115–124, New York, NY, USA. Association for Computing Machinery.
  * Liu et al. (2016) Pengfei Liu, Xipeng Qiu, and Xuanjing Huang. 2016.  Recurrent neural network for text classification with multi-task learning.
  * Liu et al. (2019) Yinhan Liu, Myle Ott, Naman Goyal, Jingfei Du, Mandar Joshi, Danqi Chen, Omer Levy, Mike Lewis, Luke Zettlemoyer, and Veselin Stoyanov. 2019.  Roberta: A robustly optimized bert pretraining approach.
  * Ostendorff et al. (2022) Malte Ostendorff, Nils Rethmeier, Isabelle Augenstein, Bela Gipp, and Georg Rehm. 2022.  Neighborhood contrastive learning for scientific document representations with citation embeddings.
  * Reid Smith and Hammond (2021) Tanya Serry Reid Smith, Pamela Snow and Lorraine Hammond. 2021.  The role of background knowledge in reading comprehension: A critical review.  _Reading Psychology_ , 42(3):214–240.
  * Reimers and Gurevych (2019) Nils Reimers and Iryna Gurevych. 2019.  Sentence-bert: Sentence embeddings using siamese bert-networks.
  * Rios and Kavuluru (2015) Anthony Rios and Ramakanth Kavuluru. 2015.  Convolutional neural networks for biomedical text classification: application in indexing biomedical articles.  In _Proceedings of the 6th ACM Conference on Bioinformatics, Computational Biology and Health Informatics_ , BCB ’15, page 258–267, New York, NY, USA. Association for Computing Machinery.
  * Shen et al. (2021) Jiaming Shen, Wenda Qiu, Yu Meng, Jingbo Shang, Xiang Ren, and Jiawei Han. 2021.  TaxoClass: Hierarchical multi-label text classification using only class names.  In _Proceedings of the 2021 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies_ , pages 4239–4249, Online. Association for Computational Linguistics.
  * Shweta and Belz (2021) FNU Shweta and Andrea Belz. 2021.  Computational linguistic analysis of submitted sec information (classi).  In _2021 IEEE 33rd International Conference on Tools with Artificial Intelligence (ICTAI)_. IEEE.
  * Song et al. (2022) Dezhao Song, Andrew Vold, Kanika Madan, and Frank Schilder. 2022.  Multi-label legal document classification: A deep learning-based approach with label-attention and domain-specific pre-training.  _Information Systems_ , 106:101718.
  * Sun et al. (2019) Xingping Sun, Yibing Li, Hongwei Kang, and Yong Shen. 2019.  Automatic document classification using convolutional neural network.  _Journal of Physics: Conference Series_ , 1176(3):032029.
  * Tang et al. (2016) B. Tang, H. He, P. M. Baggenstoss, and S. Kay. 2016.  A bayesian classification approach using class-specific features for text categorization.  _IEEE Transactions on Knowledge & Data Engineering_, 28(06):1602–1606.
  * Taylor et al. (2022) Ross Taylor, Marcin Kardas, Guillem Cucurull, Thomas Scialom, Anthony Hartshorn, Elvis Saravia, Andrew Poulton, Viktor Kerkez, and Robert Stojnic. 2022.  Galactica: A large language model for science.
  * Xiao et al. (2019) Lin Xiao, Xin Huang, Boli Chen, and Liping Jing. 2019.  Label-specific document representation for multi-label text classification.  In _Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP)_ , pages 466–475, Hong Kong, China. Association for Computational Linguistics.
  * Xu et al. (2023) Ran Xu, Yue Yu, Joyce Ho, and Carl Yang. 2023.  Weakly-supervised scientific document classification via retrieval-augmented multi-stage training.  In _Proceedings of the 46th International ACM SIGIR Conference on Research and Development in Information Retrieval_. ACM.
  * Xun et al. (2019) Guangxu Xun, Kishlay Jha, Ye Yuan, Yaqing Wang, and Aidong Zhang. 2019.  MeSHProbeNet: a self-attentive probe net for MeSH indexing.  _Bioinformatics_ , 35(19):3794–3802.
  * Yin et al. (2019) Wenpeng Yin, Jamaal Hay, and Dan Roth. 2019.  Benchmarking zero-shot text classification: Datasets, evaluation and entailment approach.  _CoRR_ , abs/1909.00161.
  * Yu et al. (2022) Hsiang-Fu Yu, Kai Zhong, Jiong Zhang, Wei-Cheng Chang, and Inderjit S Dhillon. 2022.  Pecos: Prediction for enormous and correlated output spaces.  _Journal of Machine Learning Research_ , 23(98):1–32.
  * Zhang et al. (2023) Yu Zhang, Bowen Jin, Xiusi Chen, Yanzhen Shen, Yunyi Zhang, Yu Meng, and Jiawei Han. 2023.  Weakly supervised multi-label classification of full-text scientific papers.
  * Zhang et al. (2022a) Yu Zhang, Zhihong Shen, Chieh-Han Wu, Boya Xie, Junheng Hao, Ye-Yi Wang, Kuansan Wang, and Jiawei Han. 2022a.  Metadata-induced contrastive learning for zero-shot multi-label text classification.  In _Proceedings of the ACM Web Conference 2022_ , WWW ’22, page 3162–3173, New York, NY, USA. Association for Computing Machinery.
  * Zhang et al. (2022b) Yu Zhang, Zhihong Shen, Chieh-Han Wu, Boya Xie, Junheng Hao, Ye-Yi Wang, Kuansan Wang, and Jiawei Han. 2022b.  Metadata-induced contrastive learning for zero-shot multi-label text classification.
  * Zong et al. (2015) Wei Zong, Feng Wu, Lap-Keung Chu, and Domenic Sculli. 2015.  A discriminative and semantic feature selection method for text categorization.  _International Journal of Production Economics_ , 165:215–222.